Adding a new Physics Solver¶

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

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

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 four included headers are:

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;

linearAlgebra/interfaces/InterfaceTypes.hpp which declares an interface to linear solvers and linear algebra libraries.

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 three 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, ImplicitTransient and ExplicitTransient depending on whether we are interested in the transient state and whether we want it to be discretized as a backward or a forward Euler scheme.

  enum class TimeIntegrationOption : integer
  {
    SteadyState,
    ImplicitTransient,
    ExplicitTransient
  };

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 GEOSX 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 LaplaceFEM class definition is complete):

ENUM_STRINGS( LaplaceFEM::TimeIntegrationOption, "SteadyState", "ImplicitTransient", "ExplicitTransient" )

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

class LaplaceFEM : public SolverBase
{
public:
  // 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 std::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"; }

  // This method ties properties with their supporting mesh
  virtual void RegisterDataOnMesh( Group * const MeshBodies ) override final;

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 GEOSX (see Group : the base class of GEOSX 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 return 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, std::string const &, Group * const )

The following member function RegisterDataOnMesh() 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 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
  SetupSystem( DomainPartition & domain,
               DofManager & dofManager,
               CRSMatrix< real64, globalIndex > & localMatrix,
               array1d< real64 > & localRhs,
               array1d< real64 > & localSolution,
               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;

  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
  SolveSystem( DofManager const & dofManager,
               ParallelMatrix & matrix,
               ParallelVector & rhs,
               ParallelVector & solution ) override;

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

  virtual void
    ResetStateToBeginningOfStep( DomainPartition & GEOSX_UNUSED_PARAM( domain ) ) override;

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

Eventually, ApplyDirichletBC_implicit() 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 BaseSolver 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
  {
    dataRepository::ViewKey timeIntegrationOption = { "timeIntegrationOption" };
    dataRepository::ViewKey fieldVarName = { "fieldName" };

  } laplaceFEMViewKeys;

We can check that in the LaplaceFEM companion integratedTest

    <LaplaceFEM
      name="laplace"
      discretization="FE1"
      timeIntegrationOption="SteadyState"
      fieldName="Temperature"
      logLevel="0"
      targetRegions="{ Region1 }">
      <LinearSolverParameters
        solverType="direct"
        directParallel="0"
        logLevel="0"/>
    </LaplaceFEM>

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.

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

  registerWrapper( laplaceFEMViewKeys.fieldVarName.Key(), &m_fieldName )->
    setInputFlag( InputFlags::REQUIRED )->
    setDescription( "Name of field variable" );
}

Checking out the constructor, we can see that the use of a registerWrapper<T>(...) allow 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 LaplaceFEM::RegisterDataOnMesh( Group * const MeshBodies )
{
  for( auto & mesh : MeshBodies->GetSubGroups() )
  {
    NodeManager * const nodes = mesh.second->group_cast< MeshBody * >()->getMeshLevel( 0 )->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 GEOSX_UNUSED_PARAM( time_n ),
                                 real64 const GEOSX_UNUSED_PARAM( dt ),
                                 DomainPartition & domain,
                                 DofManager const & dofManager,
                                 CRSMatrixView< real64, globalIndex const > const & localMatrix,
                                 arrayView1d< real64 > const & localRhs )
{
  MeshLevel * const mesh = domain.getMeshBodies()->GetGroup< MeshBody >( 0 )->getMeshLevel( 0 );

  NodeManager & nodeManager = *(mesh->getNodeManager());

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


  finiteElement::
    regionBasedKernelApplication< parallelDevicePolicy< 32 >,
                                  constitutive::NullModel,
                                  CellElementSubRegion,
                                  LaplaceFEMKernel >( *mesh,
                                                      targetRegionNames(),
                                                      this->getDiscretizationName(),
                                                      arrayView1d< string const >(),
                                                      dofIndex,
                                                      dofManager.rankOffset(),
                                                      localMatrix,
                                                      localRhs,
                                                      m_fieldName );

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.

Note

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 geosx
{

class LaplaceDiffFEM : public LaplaceFEM
{
public:

  LaplaceDiffFEM() = delete;

  LaplaceDiffFEM( const std::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;

  protected:
  virtual void PostProcessInput() override final;

private:
  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 std::string& name,
                                Group * const parent ):
LaplaceFEM( name, parent ), m_diffusion(0.0)
{
  registerWrapper<string>(laplaceDiffFEMViewKeys.diffusionCoeff.Key())->
    setInputFlag(InputFlags::REQUIRED)->
    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()
{
  LaplaceFEM::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.