Multi-physics simulation involves coupling multiple "physics", often with different space and/or time scales. Additionally, the single-physics codes are often written by different teams. The most commonly used coupling technique is first-order operator splitting, but this has poor accuracy and stability properties. How do I determine which algorithms will be effective for a problem of interest, and how should I structure my software to make these algorithms available?
2 Answers
Since this is a very active area of research, I hesitate to attempt an answer to this but I have some experience on what not to try.
Do Not:
- Take old application code A and old application B, then try to couple them together
- Use archaic (implying unusable in hindsight) code, instead of building a new application
- Require a huge framework (> 10 required dependencies) on new users beginning to contribute
- Assume the data layout (meshes, matrices, vectors, etc.) is easy to write yourself
Do:
- Use standard programming practices and, hopefully, good design patterns
- Read the literature on operator splitting to understand the limitations of accuracy and stability
I strongly advocate a fully coupled assembly since this can easily reproduce the operator split versions. Specifically, the routines which calculate the residual and Jacobian for different physics can be separate, but the framework should be able to combine them to form a unified residual for the entire system. This is how PETSc works.
Then, the operator split solutions can be used as a preconditioner for the fully coupled system, or as a solver in its own right, all from the command line. Moreover, some couplings can be preserved, while others can be disgarded. PETSc handles this through the PC FieldSplit interface. This allows replication of hybrid schemes, such as semi-implicit ICE for fluid dynamics.