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, physicsSolver/simplePDE/LaplaceBaseH1.hpp includes the following headers:
common/EnumStrings.hppwhich includes facilities for enum-string conversion (useful for reading enum values from input);
physicsSolver/PhysicsSolverBase.hppwhich declares the abstraction class shared by all physics solvers.
managers/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 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/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 );
/// 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 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
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 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.
Switching to protected members, postInputInitialization() 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 PhysicsSolverBase 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 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,
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";
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
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 postInputInitialization() 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
{
public:
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;
protected:
virtual void postInputInitialization() 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 , postInputInitialization() 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)
{
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 postInputInitialization().
void LaplaceDiffFEM::postInputInitialization()
{
LaplaceFEM::postInputInitialization();
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 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 integratedTests for the solver class.