Interfacing with 3rd-party PDE libraries

Summary:

Allow 3rd-party PDE(partial differential equation) libraries to be used to perform inference that involve PDEs.

Description:

The design involves cmdstan, stan, and math. This issue ticket solicits the discussion regarding the related design in all the repos. See https://github.com/stan-dev/math/issues/931 for corresponding issue for Math library.

Design goals

• Stan interface with any 3rd-party PDE library, so that QoI(quantity of interest, PDE jargon) can be connected to input parameters.
• Use external library in Stan for PDE computing. The external library is provided by the user through user_header and Stan's external C++ interface.
• Hide stan internals from PDE developers.

This design allows the user to use his own PDE solver in Stan, by using solve_pde function and make with STANCFLAGS=--allow_undefined. Similar to ODE solver in Stan, the PDE is provided through a functor(Stan's user-defined function). The actual action of the functor is provided through the external PDE solver in user_header.hpp(see example below). In this design the user is asked to do the following to use his external PDE solver:

• write a stan user-function, in the following example, this function is laplace_model.
• declare the external functions used in the above stan user-function.
• in Stan call solve_pde, with the above stan user-function as the first argument.
• provide Makefile for the external PDE solver
• build the Stan model with

make WITH_EXTERN_PDE=1 EXTERN_PDE_MAKEFILE=/path/to/makefile stan_model

An additional design goal is to allow stan to coordinate MPI computation with the external lib. This belongs to next stage, will be not be discussed in the rest of the discussion.

Also note that the design works on any 3rd-party library that can do sensitivity analysis, but my interest is mostly in PDE.

The rest subsections introduce the proposed design.

forward_pde solve_pde function

Implementation here
https://github.com/yizhang-cae/math/tree/forward_pde

forward_pde solve_pde in math maps input PDE model functor and parameters to QoI:

  template<typename F_pde, typename T>
  inline std::vector<T> solve_pde(const F_pde& pde,
                                         const std::vector<T>& theta,
                                         const std::vector<double>& x_r,
                                         const std::vector<int>& x_i,
                                         std::ostream* msgs = nullptr);

This is similar to ODE solver, except that all model-related information is provide by functor pde, which as the following signature

  inline std::vector<std::vector<double> >
  operator()(const std::vector<double>& theta,
             const bool need_sens,
             const std::vector<double>& x_r,
             const std::vector<int>& x_i,
             std::ostream* msgs = nullptr);

The rationale is to relieve PDE lib user from working on vars. He needs to make decision on

• Whether to compute sensitivity(need_sens?),
• How the sensitivity is computed from theta,
  and doesn't need to know how the result is incorporated into math.

Linking forward_pde to Stan language

Implementation is at
https://github.com/yizhang-cae/stan/tree/forward_pde
Nothing new here, just exposing a high-order function.

Using forward_pde solve_pde in Stan would be something like the following.

functions {
  real[,] solve_with_sensitivity(real[] theta);
  real[,] solve(real[] theta);
  real[,] laplace_model(real[] theta, int need_sens, real[] x_r, int[] x_i){
    if(need_sens)
      return solve_with_sensitivity(theta);
    else
      return solve(theta);
    }
}
/* .... */
parameters {
  real<lower = 0> k[2];
} 

transformed parameters{
  real QoI[1];
  QoI = solve_pde(laplace_model, k, x_r, x_i);
}
/* .... */

make in cmdstan

In order to make the above model, incmdstan two switches are added to Makefile

• WITH_EXTERN_PDE: =1 indicates user needs to supply user_header.hpp.
• EXTERN_PDE_MAKEFILE: with WITH_EXTERN_PDE=1 user needs to provide make info regarding his library, so that the above Laplace model can be built with

make WITH_EXTERN_PDE=1 EXTERN_PDE_MAKEFILE=./examples/forward_pde_laplace/Makefile ./examples/forward_pde_laplace/forward_pde_laplace

The example and the Makefile change can be found at
https://github.com/yizhang-cae/cmdstan/tree/forward_pde

Reproducible Steps:

Final model in ./examples/forward_pde_lapace
https://github.com/yizhang-cae/cmdstan/tree/forward_pde
Unit tests in the above math branch.

Current Output:

n/a

Expected Output:

Three external libraries has been tested and included in unit tests. To verify you need to install these libraries and provide make info indicated in the Makefile of the unit test directories in math repo.

• http://libmesh.github.io/
• http://mfem.org/
• http://opensees.berkeley.edu/

Pressure contour of the Darcy's flow from the MFEM unit test

Parameter(normalized porosity k) from Stan, based on simulated data.

Additional Information:

Provide any additional information here.

Current Version:

v2.17.1

Review the proposed solve_pde implementation in the math solve_pde branch, the Stan forward_pde branch, and the cmdstan forward_pde example. Start with the referenced Makefile changes, user_header.hpp interface, and math unit tests. Done means the cross-repository design is agreed and the external PDE example and tests work with the documented build flags.