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.
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, physicsSolvers/simplePDE/LaplaceBaseH1.hpp includes the following headers:
common/format/EnumStrings.hppwhich includes facilities for enum-string conversion (useful for reading enum values from input);
physicsSolvers/PhysicsSolverBase.hppwhich declares the abstraction class shared by all physics solvers.
fieldSpecification/FieldSpecificationManager.hppwhich 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 (protected) member variables:
m_fieldNamewhich stores the name of the diffused variable (e.g. the temperature) as a string;
m_timeIntegrationOptionan 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
{
SteadyState,
ImplicitTransient
};
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/format/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,
"SteadyState",
"ImplicitTransient" );
Once explained the main variables and enum, let us start reading through the different member functions:
class LaplaceFEM : public LaplaceBaseH1
{
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 string & name,
Group * const parent );
/// "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 PhysicsSolverBase::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( PhysicsSolverBase, 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 setSparsityPattern( DomainPartition & domain,
DofManager & dofManager,
CRSMatrix< real64, globalIndex > & localMatrix,
SparsityPattern< globalIndex > & pattern ) 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 PhysicsSolverBase, it can be noted that most of the solver interface functions are called during
either PhysicsSolverBase::linearImplicitStep() or PhysicsSolverBase::nonlinearImplicitStep() depending on the solver strategy chosen.
Another central member function is postInputInitialization(), inherited from Group and overridden in PhysicsSolverBase.
It is called once the input has been read from the XML file, and is the place where a solver can check or post-process its input values.
LaplaceFEM does not need to override it: as shown in the Implementation File (reference) section, its input values are registered
directly on its member variables, which are therefore filled automatically when the XML file is parsed.
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 PhysicsSolverBase::viewKeyStruct
{
static constexpr char const * timeIntegrationOption() { return "timeIntegrationOption"; }
static constexpr char const * fieldVarName() { return "fieldName"; }
};
We can check that in the LaplaceFEM companion integratedTest
<LaplaceFEM
name="laplace"
discretization="FE1"
timeIntegrationOption="SteadyState"
fieldName="Temperature"
targetRegions="{ Domain }">
<LinearSolverParameters
directParallel="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. 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 ):
PhysicsSolverBase( 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::OPTIONALif they are optional and can be provided;
InputFlags::REQUIREDif 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_fieldNamekey;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,
string_array 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";
finiteElement::
regionBasedKernelApplication< parallelDevicePolicy< >,
constitutive::NullModel,
CellElementSubRegion >( mesh,
regionNames,
this->getDiscretizationName(),
dummyString,
kernelFactory );
} );
}
assembleSystem() will be our core focus as we want to change the diffusion coefficient from its
implicit unit value to a user-defined value read from the XML file. This method is in charge of constructing
the FEM system matrix in parallel. For each mesh level targeted by the solver, it retrieves the degree-of-freedom numbers
of the nodes and builds a LaplaceFEMKernelFactory. The generic function finiteElement::regionBasedKernelApplication
then launches the kernel on all the target cell sub-regions, for the finite element space of the solver discretization.
The element-level computation is implemented in the kernel class LaplaceFEMKernel (file LaplaceFEMKernels.hpp), whose
quadraturePointKernel() adds the contribution of each quadrature point to the element matrix:
void quadraturePointKernel( localIndex const k,
localIndex const q,
StackVariables & stack ) const
{
GEOS_UNUSED_VAR( k );
real64 dNdX[ maxNumTestSupportPointsPerElem ][ 3 ];
real64 const detJ = FE_TYPE::calcGradN( q, stack.xLocal, stack.feStack, dNdX );
for( localIndex a = 0; a < stack.numRows; ++a )
{
for( localIndex b = 0; b < stack.numCols; ++b )
{
stack.localJacobian[ a ][ b ] += LvArray::tensorOps::AiBi< 3 >( dNdX[a], dNdX[b] ) * detJ;
}
}
}
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.
Start doing your own Physics solver
As we will extend LaplaceFEM capabilities, we will derive publicly from it, and write a kernel
deriving from LaplaceFEMKernel.
Note
The code snippets below illustrate the procedure; they are not compiled as part of GEOS.
Always refer to the current sources of LaplaceFEM for the exact signatures.
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
{
public:
LaplaceDiffFEM() = delete;
LaplaceDiffFEM( string const & name,
Group * const parent );
static string catalogName() { return "LaplaceDiffFEM"; }
string getCatalogName() const override { return catalogName(); }
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;
struct viewKeyStruct : public LaplaceFEM::viewKeyStruct
{
static constexpr char const * diffusionCoefficient() { return "diffusionCoefficient"; }
};
protected:
real64 m_diffusion;
};
} /* namespace geos */
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 the Group base class. Another mandatory step is to define the static catalogName() method (and the virtual
getCatalogName()), which gives the XML tag of the new solver and is used to register it in the catalog.
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
struct viewKeyStruct and on the constructor.
Our new solver viewKeyStruct will have its own structure inheriting from the LaplaceFEM one to have the timeIntegrationOption
and fieldName keys. It will also define a diffusionCoefficient key 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. Registering the wrapper with a pointer to m_diffusion binds the XML attribute to the member variable:
its value is read and stored automatically when the input file is parsed, so that there is no need to override postInputInitialization().
The diffusionCoefficient entry will be asked as a required input. If not provided, the error thrown will ask for it.
LaplaceDiffFEM::LaplaceDiffFEM( string const & name,
Group * const parent ):
LaplaceFEM( name, parent ),
m_diffusion( 0.0 )
{
registerWrapper( viewKeyStruct::diffusionCoefficient(), &m_diffusion ).
setInputFlag( dataRepository::InputFlags::REQUIRED ).
setDescription( "Uniform diffusion coefficient for the Laplace equation" );
}
REGISTER_CATALOG_ENTRY( PhysicsSolverBase, LaplaceDiffFEM, string const &, Group * const )
Now that we have required, read and bound the user-defined diffusion value to a variable, we can use it in the construction of our
matrix. The element-level computation takes place in the kernel, so we first derive a new kernel from LaplaceFEMKernel,
which receives the diffusion coefficient as an additional constructor argument and uses it in quadraturePointKernel().
Since the kernels are launched through templates on their exact type, defining quadraturePointKernel() in the derived kernel
is enough for it to be used instead of the base class one.
template< typename SUBREGION_TYPE,
typename CONSTITUTIVE_TYPE,
typename FE_TYPE >
class LaplaceDiffFEMKernel : public LaplaceFEMKernel< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >
{
public:
using Base = LaplaceFEMKernel< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE >;
using typename Base::StackVariables;
using Base::maxNumTestSupportPointsPerElem;
LaplaceDiffFEMKernel( NodeManager const & nodeManager,
EdgeManager const & edgeManager,
FaceManager const & faceManager,
localIndex const targetRegionIndex,
SUBREGION_TYPE const & elementSubRegion,
FE_TYPE const & finiteElementSpace,
CONSTITUTIVE_TYPE & inputConstitutiveType,
arrayView1d< globalIndex const > const inputDofNumber,
globalIndex const rankOffset,
CRSMatrixView< real64, globalIndex const > const inputMatrix,
arrayView1d< real64 > const inputRhs,
real64 const inputDt,
string const fieldName,
real64 const diffusion ):
Base( nodeManager, edgeManager, faceManager, targetRegionIndex, elementSubRegion,
finiteElementSpace, inputConstitutiveType, inputDofNumber, rankOffset,
inputMatrix, inputRhs, inputDt, fieldName ),
m_diffusion( diffusion )
{}
GEOS_HOST_DEVICE
inline
void quadraturePointKernel( localIndex const k,
localIndex const q,
StackVariables & stack ) const
{
GEOS_UNUSED_VAR( k );
real64 dNdX[ maxNumTestSupportPointsPerElem ][ 3 ];
real64 const detJ = FE_TYPE::calcGradN( q, stack.xLocal, stack.feStack, dNdX );
for( localIndex a = 0; a < stack.numRows; ++a )
{
for( localIndex b = 0; b < stack.numCols; ++b )
{
stack.localJacobian[ a ][ b ] += m_diffusion * LvArray::tensorOps::AiBi< 3 >( dNdX[a], dNdX[b] ) * detJ;
}
}
}
protected:
real64 const m_diffusion;
};
using LaplaceDiffFEMKernelFactory = finiteElement::KernelFactory< LaplaceDiffFEMKernel,
arrayView1d< globalIndex const > const,
globalIndex const,
CRSMatrixView< real64, globalIndex const > const,
arrayView1d< real64 > const,
real64 const,
string const,
real64 const >;
The overridden assembleSystem() is then identical to the one of LaplaceFEM, except that it builds the new kernel factory
and passes m_diffusion to it. Note that m_fieldName is a protected member of LaplaceBaseH1, so it is directly accessible
from the derived solver.
void LaplaceDiffFEM::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,
string_array const & regionNames )
{
NodeManager & nodeManager = mesh.getNodeManager();
string const dofKey = dofManager.getKey( m_fieldName );
arrayView1d< globalIndex const > const &
dofIndex = nodeManager.getReference< array1d< globalIndex > >( dofKey );
LaplaceDiffFEMKernelFactory kernelFactory( dofIndex, dofManager.rankOffset(), localMatrix, localRhs, dt,
m_fieldName, m_diffusion );
string const dummyString = "dummy";
finiteElement::
regionBasedKernelApplication< parallelDevicePolicy< >,
constitutive::NullModel,
CellElementSubRegion >( mesh,
regionNames,
this->getDiscretizationName(),
dummyString,
kernelFactory );
} );
}
This completes the implementation of our new solver LaplaceDiffFEM.
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 the new header files to the
simplePDESolvers_headerslist ofphysicsSolvers/simplePDE/CMakeLists.txt;add the new source file to the
simplePDESolvers_sourceslist of the same file;check that Doxygen comments are properly set in our solver class;
uncrustify it to match the code style by going to the build folder and running the command:
make uncrustify_style;write unit tests for each new features in the solver class;
write an input file and register it as an integrated test for the solver class (see Contributing Input Files and Integrated Tests).