Proposed structure for coupled solvers

[rrsettgast]

The expanding number of coupled solvers warrants a pause to consider how we would like to handle the relationship between coupled solvers and single physics solvers. In case it is not clear, there are quite a few coupled solvers now with varying levels of coupling between different components (e.g. FlowProppantTransportSolver, HydrofractureSolver, LagrangianContactSolver, MultiphasePoromechanicsSolver, PhaseFieldFractureSolver, SinglePhasePoromechanicsEFEMKernel, SinglePhasePoromechanicsSolver, etc.)

Some of these use composition, some use inheritance, some use a combination of composition and inheritance. In short, it is kind of a free-for-all. Now that we have several good use cases to drive design, time has come to enforce a structured design. We propose the following:

1. Propose composition/inheritance rules for physics solvers:

   • single physics solvers may use inheritance or composition
   • coupled solver must use composition (except that they inherit from SolverBase). No inheritance...more on this

2. Generic Templated Coupled solvers (GTCS). We define a generic coupled solver structure

class CoupledSolver< FLOW_SOVLER , SOLID_SOLVER, CONTACT_SOLVER, THERMAL_SOLVER, ...>
{
...
private:
   FLOW_SOVLER * m_flow;
   SOLID_SOLVER * m_mech;
   CONTACT_SOLVER * m_contact;
   THERMAL_SOLVER * m_thermal;
}

The goal is to provide a truly generic coupled solver strategy for tying different single physics components together...or even multiple coupled physics components. For instance a tight coupling may warrant a tight coupling and custom kernels be created, and reused by other coupled solvers that are composed of the tight coupled solver.

There are implications for clean/strict interfaces for solver classes in general. We should keep a tighter grip on public/protected/private interfaces.

3. Introduction of a NullSolver to avoid GTCS definitions for different numbers of component types. NullSolver will just be a no-op for all interface functions.

4. Inheritance specialization mechanism for GTCS. There will be cases where the generic class will work, but may not take advantages of fusing redundant kernel operations into a single algorithm. In cases like this we may specialize the GTCS

class SinglePhaseSolidMechanicsLagrangeContact : public CoupledSolver< SinglePhaseFlow, 
                                                                       SolidMechanics, 
                                                                       LagrangeContact, 
                                                                       NullSolver >
{ ... }

In the specialized solver, we can call custom kernels, or make whatever optimizations that we would like. This retains the generic capability while giving the opportunity for specialized improvements.

The good news is we can try something like this out and not break anything.

@CusiniM @klevzoff @corbett5 @castelletto1 @TotoGaz @joshua-white

[TotoGaz]

Is it scientifically possible to deal with solvers agnostically like

private:
   std::list< Solver * > m_solvers;

instead of

private:
   FLOW_SOLVER * m_flow;
   SOLID_SOLVER * m_mech;
   CONTACT_SOLVER * m_contact;
   THERMAL_SOLVER * m_thermal;

Or does that make no sense?

Another flavor: is it possible to rely on virtual classes here instead of template parameters?

private:
   FlowSolver * m_flow;
   SolidSolver * m_mech;
   ContactSolver * m_contact;
   ThermalSolver * m_thermal;

If so, that would surely help understanding what we want for a these Solvers.
There will be some dispatch when it comes to the kernel, so I do not know...

[corbett5]

If you know that there's no optimizations going on for a combination you could always instantiate it as CoupledSolver< SolverBase , SolverBase, SolverBase, SolverBase, ...> and that would help eliminate some template explosion. But if every combo has some optimizations then I think you do need to specify each type.

[CusiniM]

The good thing is that each specialization can be in a separate compilation unit so it should not be that bad for the compilation time, right?

[rrsettgast]

The good thing is that can each specialization ca be in a separate compilation unit so it should not be that bad for the compilation time, right?

we can also just stick a specific instantiation into a single source file. I think this is what you mean.

[rrsettgast]

@TotoGaz We did discuss the non-templated virtual interface. And it is certainly possible to have inherited specialization on the non-templated virtual class....but I would push for more control at the cost of dealing with the strong typing. I considered the container for some collection of SolverBase*, but I think there might be some confusion about what solver is what. i.e. we have to keep track of which entry in the list is the SolidMechanics and which is the FlowSolver...etc. Not that this would be hard, but it is just harder to use imo.

[TotoGaz]

IIUC, at some time, you will inevitably need to convert the virtual classes into more "static" ones.
I would think that we should wait for the last moment to do this.
What is this CoupledSolver supposed to do? Create a kernel and launch it? Something else? Would it make sense to make the virtual dispatch in one function of this class? Do we need to do it before the CoupledSolver is created?

[klevzoff]

2. Generic Templated Coupled solvers (GTCS). We define a generic coupled solver structure

class CoupledSolver< FLOW_SOVLER , SOLID_SOLVER, CONTACT_SOLVER, THERMAL_SOLVER, ...>
{
...
private:
   FLOW_SOVLER * m_flow;
   SOLID_SOLVER * m_mech;
   CONTACT_SOLVER * m_contact;
   THERMAL_SOLVER * m_thermal;
}

The goal is to provide a truly generic coupled solver strategy for tying different single physics components together...or even multiple coupled physics components. For instance a tight coupling may warrant a tight coupling and custom kernels be created, and reused by other coupled solvers that are composed of the tight coupled solver.

There are implications for clean/strict interfaces for solver classes in general. We should keep a tighter grip on public/protected/private interfaces.

1. Introduction of a NullSolver to avoid GTCS definitions for different numbers of component types. NullSolver will just be a no-op for all interface functions.

(Q to everyone) What are the pros and cons vs an even more generic definition

template< typename ... SOLVERS >
class CoupledSolver
{
private:
  std::tuple< SOLVERS ... * > m_solvers;
}

The one major drawback I can think of is having to re-define the ordering in every derived coupled solver:

class SinglePhaseSolidMechanicsLagrangeContact : public CoupledSolver< SinglePhaseFlow, 
                                                                       SolidMechanics, 
                                                                       LagrangeContact >
{
  enum SolverTypes
  {
    Flow,
    Mechanics,
    Contact
  };
  ...
  void foo()
  {
    getSolver< Mechanics >()->foo();
    // instead of m_mech->foo();
  }
}

On the plus side, generic for-loop over tuple elements is a solved problem. This will come up a lot when calling operations on sub-solvers. We can still do it with the original definition:

template< typename LAMBDA >
void forEachSolver( LAMBDA && lambda ) const
{
  std::tuple< ... > solvers { m_flow, m_mech, m_contact, m_thermal }; // must be modified for each pointer added
  for_each_element( solvers, [lambda = std::forward< LAMBDA >( lambda )] ( auto * solver )
  {
    lambda( solver );
  } );
}

But now think about e.g. printing the nonlinear residuals for each solver. We don't want to see them printed for null solvers. So we'll have to work around it:

forEachSolver( []( auto * solver )
{
  real64 const norm = solver->calculateResidualNorm(); // can be outside the branch, since it's a no-op for NullSolver
  if( !std::is_same( TYPEOFPTR( solver ), NullSolver )
  {
    GEOSX_LOG_RANK0( GEOSX_FMT( "R[{}] = {}" ), solver->getName(), norm ) );
  }
} );

All in all, this just leaves a slightly sour taste (but I'll admit there's no perfect solution here in any case).

Moreover, NullSolver has to physically exist in the data repository tree. We already made users add NullModel to their input files, will NullSolver be required too, or do we always forcibly instantiate it?

[corbett5]

Sounds good to me.

[CusiniM]

This sounds like a good idea to me.

[rrsettgast]

The only con to a more generic definition is having to find what solver is what...which wouldn't really be a problem if we treated all solvers on equal ground. I don't know that we do this now...should we? As is pointed out, we can have both the specific and generic looping interface without drawback...and there is no reason that we couldn't have both members:

class CoupledSolver< FLOW_SOVLER , SOLID_SOLVER, CONTACT_SOLVER, THERMAL_SOLVER>
{

...
private:
   FLOW_SOVLER * m_flow;
   SOLID_SOLVER * m_mech;
   CONTACT_SOLVER * m_contact;
   THERMAL_SOLVER * m_thermal;
   std::tuple< FLOW_SOVLER * , SOLID_SOLVER * , CONTACT_SOLVER * , THERMAL_SOLVER *  > m_solvers;
...

// call this at the end of init after all the pointers are set.
void setSolverCollection()
{
  m_solvers = make_tuple( m_flow, m-mech, m_contact, m_thermal ) ;
}
...
};

At least this would avoid always having to create the tuple. Don't know if it matters between having to create the tuple every time, or having multiple ways to access the same data.

Regarding the residual output issue. I think this is actually a separate issue in that we need properly formatted output for coupled solvers, which probably collects the residual for each solver and outputs them in order on the same line to give component residuals, with an aggregate residual. In this respect seeing a 0 for the NullSolver component of the residual doesn't bother me too much. But this argument simplifies the problem since a component solver can be a coupled solver on its own.

I don't think we should make the user define a NullSolver or a NullModel. They can just be created by default. There is only ever one of each, so it isn't a problem...is it?