Generic elliptic solver

From Einstein Toolkit Documentation
Jump to: navigation, search

Notes

Requirements:

  • Steve: FunWave, need parallelism and mesh refinement;
  • Ian: 6-variable, linear elliptic equation, may need mesh refinement and definitely parallelism;
  • Eloisa: generic solver, easy to use and to experiment with more important than efficiency, no restriction on topology, mesh refinement would be good but doesn't need it for everything;
  • Do these requirements limit the type of solver? Can we use multigrid? Spectral?
  • Boundary conditions play now a more important role than for hyperbolic problems. Are the current mechanisms for specifying them enough? Do we need, say, inner boundaries?

Existing tools:

  • Scott's multigrid AMR elliptic solver, second order with extension to fourth order coming soon; integrate with Cactus via own data conversions, but currently making it talk to Carpet. Available immediately via SVN;
  • Eloisa has experimented with OpenFOAM: nice and flexible, not great for accuracy, need to import data to Cactus afterwards (not complicated, but unfeasible to do at each timestep).
  • Current CactusElliptic: EllBase gives interface to register elliptic solvers, currently not much implemented (SOR), equation type is a little restrictive (linear), compatibility with AMR unknown;
  • CarpetPETSc (can PETSc work with AMR)? LSUPETSc? Note: PETSc doesn't work with OpenMP;
  • Lorene/KADATH;
  • Existing implementations:
  • There's currently nothing that allows for AMR!
  • Chombo multigrid AMR?
  • HYPRE?
  • PAMR?

Possible directions:

  • EllBase road, make type of equation more generic:
    • Add more terms and coefficients;
    • Ian: specify the equation via a callback for the residual, a registration method for the variables to solve for and possibly the derivatives/jacobians, much like in MoL (these could be generated via Kranc); this would also really only work with relaxation methods, not for direct inversion methods. Many solvers, all working along these guidelines? EllBase would provide the interface;
  • Petsc
    • Use petsc directly
    • Write an interface between petsc and cactus

Discussion and conclusions

  • There is a distinction between inversion methods and relaxation methods; the former use matrix inversion to solve a linear system (possibly coming from the linearization of a non-linear problem), whereas the latter reduce the residuals to zero by implementing the corresponding parabolic equation and let the solution relax to a static one. PETSc has some infrastructure for relaxation methods (namely, multigrid solvers), but has much more in the way of inversion methods (along with linearization methods). Additionally, AMR is easy to cope with in the inversion methods, where the grid (in real space or in frequency space) is mapped to a matrix and its "connectivity" is irrelevant, but is harder to have around in relaxation methods, where the jacobians of all the involved operators (including prolongation, ...) have to appear;
  • Erik has some experience implementing solvers based on the PETSc library, and has confidence that this method would yield robust tools as long as one is willing to code up all the equation jacobians; these will be variably complicated depending on the type of mesh refinement, i.e. relatively simple for unigrid or old-school multipatch, more complicated for new multipatch and AMR. This means one could complete this project step-by-step, progressively adding AMR capabilities by coding up the appropriate jacobians;
  • TexMEx, a multigrid solver with AMR capabilities, is currently available in Cactus. The integration with Cactus (and Carpet) is currently minimal: TexMEx is used as a black box that produces the result (with its own data structures, etc.) and only at the end communicates it to Cactus. This prevents the framework from seeing inside its algorithm, which is not ideal for debugging, and additionally makes it hard to extend the package, e.g. by letting the user specify a callback function to calculate the residual. As of now, TexMEx implements the hamiltonian constraint, which doesn't help who (Steve, Ian) is interested in solving other types of equations.
  • Peter's NoExcision thorn uses [Conjugate Gradient] methods with finite differencing. This seems like a good method to use and should be easy to implement. A generic infrastructure could be written for this, where a user provides just the RHS function as in MoL. This could also be done in Kranc. Not sure what the drawback would be - maybe it would not be efficient enough, or maybe there are subtleties we are not aware of due to lack of knowledge. Assuming it would work with mesh-refinement, but have no knowledge of any possible problems.

Plan

  • Take LSUPETSc as a basis, and try to add the missing pieces to make it work with the grids we're interested in;
  • Check with the TexMEx developers what their roadmap is for integration with Carpet, and see if we can make this tool solve generic equations.
  • Implement a conjugate-gradient solver in Kranc and then try to generalise to a generic elliptic interface.