Implementation of a problem: example

This page will explain briefly what is going on in CEPS when running a simulation. This will be illustrated with the pseudo wave front problem, which is a manufactured problem that mimicks the planar propagation of an action potential, but with an analytic expression of the solution. The solution is basically the product of two hyperbolic tangent functions.

Note that the code snippets below are manually extracted. The actual source files may have changed by the time you read this page. If the API didn't change, the class and method names are clickable, so that you can access their documentation.

The main program: instanciation of the problem

# applications/ceps.cpp:38
void
doSingleRun(InputParameters* p, Geometry* g)
{
  AbstractPdeProblem* pb = nullptr;
  CepsString s = ceps::toKey(p->getString("problem type"));
  // ...
  else if (s=="PSEUDOWAVEFRONT")
    pb = new PseudoWaveFrontProblem(g,p);
  // ...
  pb->run();
  ceps::destroyObject(pb);
  return;
}

Here, the input parameters have been read from the file passed in the command line, and the geometry is already initialized and partitioned. The type of problem to instantiate is deduced from parameters (what isn't?).

Let's see what happens in the constructor of the problem:

# pde/problem/PseudoWaveFrontProblem.cpp:32
 
PseudoWaveFrontProblem::PseudoWaveFrontProblem(Geometry* g, InputParameters *p) :
  HeatProblem(g,p),
  m_frontX0  (CepsMathVertex( 0,0,0)),
  m_frontX1  (CepsMathVertex(-1,0,0)),
  m_frontN   (CepsMathVertex( 1,0,0)),
  m_frontVel (1.),
  m_frontEps0(1.e-2),
  m_frontEps1(1.e-2)
{
  if (ceps::isValidPtr(p))
    PseudoWaveFrontProblem::setupWithParameters(p);
  m_functions->addEntry("dicoSrc",new PseudoWaveFrontSourceTerm(this));
  m_functions->addEntry("dicoSol",new PseudoWaveFrontSolution  (this));
}

which calls

# pde/problem/HeatProblem.cpp:33
 
HeatProblem::HeatProblem(Geometry* g, InputParameters* p, const CepsMathTensor& k) 
: AbstractTimedPdeProblem (g,p), m_k(k)
{
  if (ceps::isValidPtr(p))
    HeatProblem::setupWithParameters(p);
}

# pde/problem/AbstractPdeProblem.cpp:602
 
AbstractTimedPdeProblem::AbstractTimedPdeProblem(
  Geometry*        geom, 
  InputParameters* params
):
  AbstractPdeProblem(geom,params),
  m_pdeStartTime(0.),
  m_pdeEndTime(1.),
  m_pdeTimeStep(0.1),
  m_pdeSnapshotTime(0.1)
{
  if (m_parameters)
    AbstractTimedPdeProblem::setupWithParameters(m_parameters);
}

# pde/problem/AbstractPdeProblem.cpp:370
 
AbstractPdeProblem::AbstractPdeProblem(
  Geometry*        geom, 
  InputParameters* params
):
  m_geom              (geom),
  m_parameters        (params),
  m_discr             (nullptr),
  m_ownedDiscr        (false),
  m_analyticSolution  (nullptr),
  m_ownedRefSol       (false),
  m_refSolFiles       (""),
  m_refSolSnapDt      (1.),
  m_functions         (new FunctionDictionary()),
  m_boundaryConditions(new BoundaryConditionManager(m_functions)),
  m_sourceTerms       (new SourceTermManager       (m_functions)),
  m_ownedFunctions    (true),
  m_ownedBCs          (true),
  m_ownedSrcs         (true),
  m_outputFileBase    ("./result"),
  m_outputFormat      (CepsOutputFormat::PSERIES),
  m_writeGlobalIDs    (false),
  m_probePoints       ({})
{
  CEPS_ABORT_IF(ceps::isNullPtr(geom),
    "When instantiating pde problem : passed ptr to Geometry is null"
  );
  if (m_parameters)
    AbstractPdeProblem::setupWithParameters(m_parameters);
}

As you can see, the structures of the problem can be found at different levels of heritage as we tried to factorize in a reasonable way what could be. More importantly, each one of these constructors call a method named setupWithParameters which extracts the values required to initialize the problem from the input file.

Note that, at this stage, the spatial discretization is not created. m_discr is still a nullptr. Even if we have the geometry, we still do not know what defines the problems at first : its unknowns.

Defining functions for the problem

After the abstract problem, abstract timed problem and heat problem constructors have done their job and all parameters are read, you can see that we add two entries into the dictionary of functions:

# pde/problem/PseudoWaveFrontProblem.cpp:32
 
PseudoWaveFrontProblem::PseudoWaveFrontProblem(Geometry* g, InputParameters *p) :
  HeatProblem(g,p),
  m_frontX0  (CepsMathVertex( 0,0,0)),
  m_frontX1  (CepsMathVertex(-1,0,0)),
  m_frontN   (CepsMathVertex( 1,0,0)),
  m_frontVel (1.),
  m_frontEps0(1.e-2),
  m_frontEps1(1.e-2)
{
  if (ceps::isValidPtr(p))
    PseudoWaveFrontProblem::setupWithParameters(p);
-->  m_functions->addEntry("dicoSrc",new PseudoWaveFrontSourceTerm(this)); <--
-->  m_functions->addEntry("dicoSol",new PseudoWaveFrontSolution  (this)); <--
}

These are the functions that are the source term of the problem, and its analytic solution, respectively. They are declared within the PseudoWaveFrontProblem class:

# pde/problem/PseudoWaveFrontProblem.hpp
 
class PseudoWaveFrontProblem : public HeatProblem
{
  // ...
 
    class PseudoWaveFrontSourceTerm: public ScalarSAFunc
    {
 
      public:
 
        explicit PseudoWaveFrontSourceTerm(PseudoWaveFrontProblem* p);
 
        CepsReal
        eval(CepsStandardArgs args) override;
 
        CepsEnum
        getFlags() const override;
 
      private:
 
        PseudoWaveFrontProblem* m_p; 
 
    };
 
    class PseudoWaveFrontSolution: public ScalarSAFunc
    {
 
      public:
 
        explicit PseudoWaveFrontSolution(PseudoWaveFrontProblem* p);
 
        CepsReal
        eval(CepsStandardArgs args) override;
 
        CepsEnum
        getFlags() const override;
 
      private:
 
        PseudoWaveFrontProblem* m_p; 
 
    };
// ...
};

Both these classes are derived from ScalarSAFunc. SAFunc are functors that can return anything, but take only one argument for their evaluation : a variable of structured type CepsStandardArgs. You can extract the needed information from these arguments:

CepsReal
PseudoWaveFrontProblem::PseudoWaveFrontSolution::eval(CepsStandardArgs args)
{
  CepsReal dotProd0 = 0.e0;
  CepsReal dotProd1 = 0.e0;
  for (uint i = 0; i < 3; i++)
  {
    dotProd0 += (args.x[i]-m_p->m_frontX0(i))*m_p->m_frontN(i);
    dotProd1 += (args.x[i]-m_p->m_frontX1(i))*m_p->m_frontN(i);
  }
 
  return 0.5*(1+tanh((dotProd0-m_p->m_frontVel*args.t)/m_p->m_frontEps0))
         0.5*(1-tanh((dotProd1-m_p->m_frontVel*args.t)/m_p->m_frontEps1));
}

Now run !

Once the problem is instantiated, the main application calls its run method. Since for this problem, we do exactly the same steps as for the heat equation, you can find the implementation in HeatProblem. Here it is:

# pde/problem/HeatProblem.cpp
 
void
HeatProblem::run()
{
 
  if (ceps::isProfiling())
    getProfiler()->start("whole","solving the PDE (HeatProblem class)");
 
  initializeEquation();
 
  HeatSolver solver(this);
  
  solver.solve(); 
 
  if (hasAnalyticSolution())
    m_errors = solver.getErrors();
 
  if (ceps::isProfiling())
    getProfiler()->stop("whole");  
 
}

Side note: you can see that we use a profiler to track the time spent on a specific task !

The structure of the run is simple

1. Initialize the equation : we will see later what it does
2. Then a solver is created, and we directly call its solve method. It... solves the equation. But it also invokes a VtkWriter that will output the solution as the user asked. We won't go into the details of the numerical solver in this page, since there are too many possibilities and requires jumps back and forth between many files.
3. Finally, since this problem has an analytic solution, we fetch the errors that were computed by the solver all along the computation.

Defining the equation

Let's see the contents of initializeEquation

void
AbstractTimedPdeProblem::initializeEquation()
{
  AbstractPdeProblem::initializeEquation();
  this->defineInitialCondition();
}

I said, the real initializeEquation()

void
AbstractPdeProblem::initializeEquation()
{
  this->defineUnknowns();
  this->createSpatialDiscretization();
  this->defineSourceTerms();
  this->defineBoundaryConditions();
  this->defineAnalyticSolution();
}

Perfect. You can see that it is mostly a succession of defines. These methods are virtual and can be overloaded at will. By default, they do nothing. There is one intruder : you can see that right after the unknowns of the problem are defined, we can determine how many degrees of freedom there will be on each element of the mesh. This is why we create the spatial discretization right after that. It is quite a heavy routine ! For finite elements, it creates the nodal points from which the basis functions are built (they can be different from geom nodes!), then the finite elements, then the degrees of freedom, then the distributed factory from which we can build the linear system.

Let's see the defines for our pseudo wave front:

# pde/problem/heatProblem.cpp:58 Yes it is the same as Heat Problem
 
void
HeatProblem::defineUnknowns()
{
  // A single spatial unknown defined everywhere
  addUnknown("u");
}

Our problem has a single spatial unknown, and it is defined everywhere, so it is quite easy.

void
PseudoWaveFrontProblem::defineSourceTerms()
{
  m_sourceTerms->add("manufacturedSource dicoSrc UNKNOWN 0");
}
void
PseudoWaveFrontProblem::defineBoundaryConditions()
{
  m_boundaryConditions->add("dirichlet DIRICHLET dicoSol UNKNOWN 0");
}
 
void
PseudoWaveFrontProblem::defineInitialCondition()
{
  m_initialCondition = m_functions->getScalar("dicoSol");
}
 
void
PseudoWaveFrontProblem::defineAnalyticSolution()
{
  m_analyticSolution = m_functions->getScalar("dicoSol");
}

Here, note that the analytic solution and initial condition are directly picked from the functions dictionary of the problem. On the contrary, we must give additional information for the source term and the boundary condition (where do they apply, what is their type, which unknown)... The options that must be given to properly create boundary conditions and source terms are given here.

That being said, we covered most of the aspects you need to know to create a new problem in CEPS. But we barely scratched the surface, and there are of course many technical caveats. For that, you need to put on your spelunker hat and dig into the code.