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?
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.
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.