Adding a new Physics Solver

In this tutorial, you will learn how to construct a new GEOS Physics Solver class. We will use LaplaceFEM solver, computing the solution of the Laplace problem in a specified material, as a starting point.

D^* \Delta X = f \quad \mbox{in\;} \Omega \\
X = X^g \quad \mbox{on\;} \Gamma_g

It is advised to read XML Input preliminary to this tutorial. The goal of this document is to explain how to develop a new solver that solves Laplace’s equation with a constant diffusion coefficient that is specified by users in the XML input.

For readability, member functions in the text will be referenced by their names but their arguments will be omitted.

LaplaceFEM overview

The LaplaceFEM solver can be found in ./src/coreComponents/physicsSolvers/simplePDE/. Let us inspect declarations in LaplaceFEM.hpp and implementations in LaplaceFEM.cpp before diving into specifying a new solver class that meets our needs.

Declaration file (reference)

The included header is physicsSolvers/simplePDE/LaplaceBaseH1.hpp which declares the base class LaplaceBaseH1, shared by all Laplace solvers. Moreover, physicsSolver/simplePDE/LaplaceBaseH1.hpp includes the following headers:

  • common/EnumStrings.hpp which includes facilities for enum-string conversion (useful for reading enum values from input);

  • physicsSolver/SolverBase.hpp which declares the abstraction class shared by all physics solvers.

  • managers/FieldSpecification/FieldSpecificationManager.hpp which declares a manager used to access and to set field on the discretized domain.

Let us jump forward to the class enum and variable as they contain the data used specifically in the implementation of LaplaceFEM.

class enums and variables (reference)

The class exhibits two member variables:

  • m_fieldName which stores the name of the diffused variable (e.g. the temperature) as a string;

  • m_timeIntegrationOption an enum value allowing to dispatch with respect to the transient treatment.

TimeIntegrationOption is an enum specifying the transient treatment which can be chosen respectively between SteadyState and ImplicitTransient depending on whether we are interested in the transient state.

  enum class TimeIntegrationOption : integer

In order to register an enumeration type with the Data Repository and have its value read from input, we must define stream insertion/extraction operators. This is a common task, so GEOS provides a facility for automating it. Upon including common/EnumStrings.hpp, we can call the following macro at the namespace scope (in this case, right after the LaplaceBaseH1 class definition is complete):

ENUM_STRINGS( LaplaceBaseH1::TimeIntegrationOption,
              "ImplicitTransient" );

Once explained the main variables and enum, let us start reading through the different member functions:

class LaplaceFEM : public LaplaceBaseH1
  /// The default nullary constructor is disabled to avoid compiler auto-generation:
  LaplaceFEM() = delete;

  /// The constructor needs a user-defined "name" and a parent Group (to place this instance in the
  /// tree structure of classes)
  LaplaceFEM( const string & name,
              Group * const parent );

  /// Destructor
  virtual ~LaplaceFEM() override;

  /// "CatalogName()" return the string used as XML tag in the input file.  It ties the XML tag with
  /// this C++ classes. This is important.
  static string catalogName() { return "LaplaceFEM"; }
   * @copydoc SolverBase::getCatalogName()
  string getCatalogName() const override { return catalogName(); }

Start looking at the class LaplaceFEM constructor and destructor declarations shows the usual string name and Group* pointer to parent that are required to build the global file-system like structure of GEOS (see Group : the base class of GEOS for details). It can also be noted that the nullary constructor is deleted on purpose to avoid compiler automatic generation and user misuse.

The next method catalogName() is static and returns the key to be added to the Catalog for this type of solver (see A few words about the ObjectCatalog for details). It has to be paired with the following macro in the implementation file.

REGISTER_CATALOG_ENTRY( SolverBase, LaplaceFEM, string const &, Group * const )

Finally, the member function registerDataOnMesh() is declared in the LaplaceBaseH1 class as

  /// This method ties properties with their supporting mesh
  virtual void registerDataOnMesh( Group & meshBodies ) override final;

It is used to assign fields onto the discretized mesh object and will be further discussed in the Implementation File (reference) section.

The next block consists in solver interface functions. These member functions set up and specialize every time step from the system matrix assembly to the solver stage.

  virtual void
  setupSystem( DomainPartition & domain,
               DofManager & dofManager,
               CRSMatrix< real64, globalIndex > & localMatrix,
               ParallelVector & rhs,
               ParallelVector & solution,
               bool const setSparsity = false ) override;

  virtual void
  assembleSystem( real64 const time,
                  real64 const dt,
                  DomainPartition & domain,
                  DofManager const & dofManager,
                  CRSMatrixView< real64, globalIndex const > const & localMatrix,
                  arrayView1d< real64 > const & localRhs ) override;

Furthermore, the following functions are inherited from the base class.

  virtual real64 solverStep( real64 const & time_n,
                             real64 const & dt,
                             integer const cycleNumber,
                             DomainPartition & domain ) override;

  virtual void
  implicitStepSetup( real64 const & time_n,
                     real64 const & dt,
                     DomainPartition & domain ) override;

  virtual void
  setupDofs( DomainPartition const & domain,
             DofManager & dofManager ) const override;

  virtual void
  applyBoundaryConditions( real64 const time,
                           real64 const dt,
                           DomainPartition & domain,
                           DofManager const & dofManager,
                           CRSMatrixView< real64, globalIndex const > const & localMatrix,
                           arrayView1d< real64 > const & localRhs ) override;

  virtual void
  applySystemSolution( DofManager const & dofManager,
                       arrayView1d< real64 const > const & localSolution,
                       real64 const scalingFactor,
                       real64 const dt,
                       DomainPartition & domain ) override;

  virtual void updateState( DomainPartition & domain ) override final;

  virtual void
    resetStateToBeginningOfStep( DomainPartition & GEOS_UNUSED_PARAM( domain ) ) override;

  virtual void
  implicitStepComplete( real64 const & time,
                        real64 const & dt,
                        DomainPartition & domain ) override;

  /// This method is specific to this Laplace solver.
  /// It is used to apply Dirichlet boundary condition
  /// and called when the base class applyBoundaryConditions() is called.
  virtual void applyDirichletBCImplicit( real64 const time,
                                         DofManager const & dofManager,
                                         DomainPartition & domain,
                                         CRSMatrixView< real64, globalIndex const > const & localMatrix,
                                         arrayView1d< real64 > const & localRhs );

Eventually, applyDirichletBCImplicit() is the working specialized member functions called when applyBoundaryConditions() is called in this particular class override.

Browsing the base class SolverBase, it can be noted that most of the solver interface functions are called during either SolverBase::linearImplicitStep() or SolverBase::nonlinearImplicitStep() depending on the solver strategy chosen.

Switching to protected members, postProcessInput() is a central member function and will be called by Group object after input is read from XML entry file. It will set and dispatch solver variables from the base class SolverBase to the most derived class. For LaplaceFEM, it will allow us to set the right time integration scheme based on the XML value as will be further explored in the next Implementation File (reference) section.

Let us focus on a struct that plays an important role: the viewKeyStruct structure.

viewKeyStruct structure (reference)

This embedded instantiated structure is a common pattern shared by all solvers. It stores dataRepository::ViewKey type objects that are used as binding data between the input XML file and the source code.

  struct viewKeyStruct : public SolverBase::viewKeyStruct
    static constexpr char const * timeIntegrationOption() { return "timeIntegrationOption"; }
    static constexpr char const * fieldVarName() { return "fieldName"; }

We can check that in the LaplaceFEM companion integratedTest

      targetRegions="{ Domain }">

In the following section, we will see where this binding takes place.

Implementation File (reference)

Switching to implementation, we will focus on few implementations, leaving details to other tutorials. The LaplaceFEM constructor is implemented as follows.

LaplaceFEM::LaplaceFEM( const string & name,
                        Group * const parent ):
  LaplaceBaseH1( name, parent )

As we see, it calls the LaplaceBaseH1 constructor, that is implemented as follows.

LaplaceBaseH1::LaplaceBaseH1( const string & name,
                              Group * const parent ):
  SolverBase( name, parent ),
  m_fieldName( "primaryField" ),
  m_timeIntegrationOption( TimeIntegrationOption::ImplicitTransient )
  this->registerWrapper( viewKeyStruct::timeIntegrationOption(), &m_timeIntegrationOption ).
    setInputFlag( InputFlags::REQUIRED ).
    setDescription( "Time integration method. Options are:\n* " + EnumStrings< TimeIntegrationOption >::concat( "\n* " ) );

  this->registerWrapper( viewKeyStruct::fieldVarName(), &m_fieldName ).
    setRTTypeName( rtTypes::CustomTypes::groupNameRef ).
    setInputFlag( InputFlags::REQUIRED ).
    setDescription( "Name of field variable" );


Checking out the constructor, we can see that the use of a registerWrapper<T>(...) allows us to register the key value from the enum viewKeyStruct defining them as:

  • InputFlags::OPTIONAL if they are optional and can be provided;

  • InputFlags::REQUIRED if they are required and will throw error if not;

and their associated descriptions for auto-generated docs.

void LaplaceBaseH1::registerDataOnMesh( Group & meshBodies )
  meshBodies.forSubGroups< MeshBody >( [&] ( MeshBody & meshBody )
    NodeManager & nodes = meshBody.getBaseDiscretization().getNodeManager();

    nodes.registerWrapper< real64_array >( m_fieldName ).
      setApplyDefaultValue( 0.0 ).
      setPlotLevel( PlotLevel::LEVEL_0 ).
      setDescription( "Primary field variable" );
  } );

registerDataOnMesh() is browsing all subgroups in the mesh Group object and for all nodes in the sub group:

  • register the observed field under the chosen m_fieldName key;

  • apply a default value;

  • set the output verbosity level (here PlotLevel::LEVEL_0);

  • set the field associated description for auto generated docs.

void LaplaceFEM::assembleSystem( real64 const GEOS_UNUSED_PARAM( time_n ),
                                 real64 const dt,
                                 DomainPartition & domain,
                                 DofManager const & dofManager,
                                 CRSMatrixView< real64, globalIndex const > const & localMatrix,
                                 arrayView1d< real64 > const & localRhs )
  forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
                                                                MeshLevel & mesh,
                                                                arrayView1d< string const > const & regionNames )
    NodeManager & nodeManager = mesh.getNodeManager();
    string const dofKey = dofManager.getKey( m_fieldName );
    arrayView1d< globalIndex const > const &
    dofIndex =  nodeManager.getReference< array1d< globalIndex > >( dofKey );

    LaplaceFEMKernelFactory kernelFactory( dofIndex, dofManager.rankOffset(), localMatrix, localRhs, dt, m_fieldName );

    string const dummyString = "dummy";
      regionBasedKernelApplication< parallelDevicePolicy< >,
                                    CellElementSubRegion >( mesh,
                                                            kernelFactory );

  } );


assembleSystem() will be our core focus as we want to change the diffusion coefficient from its hard coded value to a XML read user-defined value. One can see that this method is in charge of constructing in a parallel fashion the FEM system matrix. Bringing nodeManager and ElementRegionManager from domain local MeshLevel object together with FiniteElementDiscretizationManager from the NumericalMethodManager, it uses nodes embedded loops on degrees of freedom in a local index embedded loops to fill a matrix and a rhs container.

As we spotted the place to change in a code to get a user-defined diffusion coefficient into the game, let us jump to writing our new LaplaceDiffFEM solver.


We might want to remove final keyword from postProcessInput() as it will prevent you from overriding it.

Start doing your own Physic solver

As we will extend LaplaceFEM capabilities, we will derive publicly from it.

Declaration File

As there is only few places where we have to change, the whole declaration file is reported below and commented afterwards.

#include "physicsSolvers/simplePDE/LaplaceFEM.hpp"

namespace geos

class LaplaceDiffFEM : public LaplaceFEM

  LaplaceDiffFEM() = delete;

  LaplaceDiffFEM( const string& name,
                  Group * const parent );

  virtual ~LaplaceDiffFEM() override;

  static string catalogName() { return "LaplaceDiffFEM"; }

  virtual void
  assembleSystem( real64 const time,
                  real64 const dt,
                  DomainPartition * const domain,
                  DofManager const & dofManager,
                  ParallelMatrix & matrix,
                  ParallelVector & rhs ) override;

  struct viewKeyStruct : public LaplaceFEM::viewKeyStruct
    dataRepository::ViewKey diffusionCoeff = { "diffusionCoeff" };
  } laplaceDiffFEMViewKeys;

  virtual void postProcessInput() override final;

  real64 m_diffusion;


We intend to have a user-defined diffusion coefficient, we then need a real64 class variable m_diffusion to store it.

Consistently with LaplaceFEM, we will also delete the nullary constructor and declare a constructor with the same arguments for forwarding to Group master class. Another mandatory step is to override the static CatalogName() method to properly register any data from the new solver class.

Then as mentioned in Implementation File (reference), the diffusion coefficient is used when assembling the matrix coefficient. Hence we will have to override the assembleSystem() function as detailed below.

Moreover, if we want to introduce a new binding between the input XML and the code we will have to work on the three struct viewKeyStruct , postProcessInput() and the constructor.

Our new solver viewKeyStruct will have its own structure inheriting from the LaplaceFEM one to have the timeIntegrationOption and fieldName field. It will also create a diffusionCoeff field to be bound to the user defined homogeneous coefficient on one hand and to our m_diffusion class variable on the other.

Implementation File

As we have seen in Implementation File (reference), the first place where to implement a new register from XML input is in the constructor. The diffusionCoeff entry we have defined in the laplaceDiffFEMViewKeys will then be asked as a required input. If not provided, the error thrown will ask for it described asked an “input uniform diffusion coefficient for the Laplace equation”.

LaplaceDiffFEM::LaplaceDiffFEM( const string& name,
                                Group * const parent ):
LaplaceFEM( name, parent ), m_diffusion(0.0)
    setDescription("input uniform diffusion coeff for the laplace equation");

Another important spot for binding the value of the XML read parameter to our m_diffusion is in postProcessInput().

void LaplaceDiffFEM::postProcessInput()

  string sDiffCoeff = this->getReference<string>(laplaceDiffFEMViewKeys.diffusionCoeff);
  this->m_diffusion = std::stof(sDiffCoeff);

Now that we have required, read and bind the user-defined diffusion value to a variable, we can use it in the construction of our matrix into the overridden assembleSystem().

// begin element loop, skipping ghost elements
for( localIndex k=0 ; k<elementSubRegion->size() ; ++k )
  if(elemGhostRank[k] < 0)
    element_rhs = 0.0;
    element_matrix = 0.0;
    for( localIndex q=0 ; q<n_q_points ; ++q)
      for( localIndex a=0 ; a<numNodesPerElement ; ++a)
        elemDofIndex[a] = dofIndex[ elemNodes( k, a ) ];

        for( localIndex b=0 ; b<numNodesPerElement ; ++b)
          element_matrix(a,b) += detJ[k][q] *
                                 m_diffusion *
                               + Dot( dNdX[k][q][a], dNdX[k][q][b] );

    matrix.add( elemDofIndex, elemDofIndex, element_matrix );
    rhs.add( elemDofIndex, element_rhs );

This completes the implementation of our new solver LaplaceDiffFEM.

Nonetheless, the compiler should complain that m_fieldName is privately as inherited from LaplaceFEM. One should then either promote m_fieldName to protected or add a getter in LaplaceFEM class to correct the error. The getter option has been chosen and the fix in our solver is then:

array1d<globalIndex> const & dofIndex =
  nodeManager->getReference< array1d<globalIndex> >( dofManager.getKey( getFieldName() ) );

Note: For consistency do not forget to change LaplaceFEM to LaplaceDiffFEM in the guards comments

Last steps

After assembling both declarations and implementations for our new solver, the final steps go as:

  • add declarations to parent CMakeLists.txt (here add to physicsSolvers_headers );

  • add implementations to parent CMakeLists.txt (here add to physicsSolvers_sources);

  • check that Doxygen comments are properly set in our solver class;

  • uncrustify it to match the code style;

  • write unit tests for each new features in the solver class;

  • write an integratedTests for the solver class.