Programming Physics Algorithm Using DamBreak++

DamBreak++ Simulation Framework

Programming physics algorithm can be incredibly complex due to the mathematical model complexity, precision, and performance demands involved. DamBreak++ is a programming environment designed to accelerate the prototyping phase in the validation of physics-based algorithms on the so-called Dam Break problem. Efforts are focused on reusability and extensibility of software components. The One-Dimensional dam break problem is a very useful benchmark to assess the numerical stability of the model in extreme conditions. See technical note below

Scientific Framework (ODE Solvers)

A key element in numerical simulation is the selection of numerical algorithm for the problem at hand. This is particularly true in industrial projects; it is often necessary to change a model or algorithm for the sake of correction or improvements during the development or execution of a numerical simulation. The numerical treatment of the mathematical terms of the governing equations is of paramount importance to get accurate results. We want to combine different numerical treatment for each of these terms as well as numerical methods. 

Implementing such a programming environment requires an architecture that specifies flexible data sharing and method invocation relationships between components while maintaining performance. However, despite the differences, most algorithmic developments share many commons aspects.  

The framework I have developed is essentially a collection of C++ classes organized in libraries. This framework implemented contains basic building blocks for the numerical solution, corresponding to the explicit conservative finite difference (based on finite volume discretization). The code uses inheritance to group similar data structures and procedures into families that share some elements and that are, in appropriate contexts, interchangeable.

Physical simulations typically involve several different types of physics. Each type of physics is best handled with a different numerical scheme (discretization).  For example, derivative can be approximated by finite difference approximation, numerical flux approximated by complex algorithms (such as shock-capturing).

I want to show a physical algorithm prototype I have implemented using the framework.

Physical Algorithm Components Based to Build Numerical Method

Finite volume method the computational domain is represented by a collection of simple domains, called cells or elements. The differential system is integrated over each control volume to produce the discrete equivalent of the conservation law. Main task is to evaluate the numerical fluxes at cell face by some algorithm, but also to evaluate source terms which contains derivative.

The main purpose of the PhysicalAlgorithm abstract class is to provide the guts of the simulation (using components to program complex algorithms quickly). The framework provides a rich set of basic types from which the user can use or inherit from. Details are encapsulated in classes and can be interchangeable. As required by the contract of PhysicalAlgorithm, the derived class provides an implementation for the calculate method as shown below:   

3 arguments: PhysicalSystem, FiniteVolumeDiscretization and Simulation bean

Physical algorithm based on HLL Riemann solver for convective flux evaluation and Nujic source term discretization (Nujic 1995). This algorithm was proposed by E. McNeil (Hydro-Quebec).

void EMcNeilAlgorithm::calculate( dbpp::PhysicalSystem* aPhysys,
    const std::shared_ptr<dbpp::FiniteVolumeDiscretization>& aGlbdiscr, 
    Sfx::Simulation* aSim)
  {
    // Harten-Lax-van Leer flux algorithm
    dbpp::HLLFluxAlgorithm w_hllAlgo; 
    if( w_hllAlgo.useReconstr())
    {
      w_hllAlgo.setReconstrVarOrder(2);
      w_hllAlgo.setReconstrType(dbpp::eReconstrType::MUSCL);
    }
    w_hllAlgo.setApprRiemanSolver(nullptr); // default is HLL solver Riemann solver
    w_hllAlgo.setFluxAlgoPrms(); // flux parameters algorithm

    using enum dbpp::SrcNumericalTreatment::eDerivativeType;
    
    // Source terms evaluation Sf (energy slope) and Sb (bottom slope)
    dbpp::SrcNumericalTreatment w_treatmentSrc{};
    w_treatmentSrc.useManningFormula();
    w_treatmentSrc.setManningRoughnessCoeff(0.); //frictionless as default 
    w_treatmentSrc.setSfSbNumericalTreatment("Nujic(1995)"); // default
    w_treatmentSrc.setHderivativeType(centered_2nd); //ditto
    w_treatmentSrc.setSpatialDerivativeOrder(2); //bottom term d_xZ  
    w_treatmentSrc.setHAverage(true); // Nujic (1995) algorithm

    // Create a RHS discretization algorithm to compute RHS terms 
    dbpp::RhsHLLFluxSrc w_swerhsAlgo{ std::string{"RhsHLLFluxSrc"}, aPhysys };
    w_swerhsAlgo.setFluxAlgorithm(&w_hllAlgo); // convective flux treatement
    w_swerhsAlgo.setSourceTermDiscr(&w_treatmentSrc); // energy slope and bottom

...

This code above shows how the different components of the physical algorithm can be set to solve numerically the dam break problem.

Next is to select a numerical method to simulate or evaluate numerical values at global nodes. In this example I select “SemiDiscreteMethod”. The System of equations to be solved can be represented as ODE (Ordinary Differential Equations).

For more details see the following technical note: Numerical and Mathematical

Again an abstract class NumericalMethod from which SemiDiscreteMethod inherit from. FiniteVolumeDiscretization hold a reference on the method, strong dependence between these 2 components. 

The method use an ODESolver to solve numerically and advance in time the solution.

...

// retrieve numerical method (Numerical representation of math equations)
// use semi discrete method (time and space are discretized separetly) 
const auto& w_numethod = aGlbdiscr->method();
// numerical method in use for this algorithm
if( auto ptr2NumethodCast = 
  std::dynamic_pointer_cast<dbpp::SemiDiscreteMethod>(w_numethod); ptr2NumethodCast)
 {
   ptr2NumethodCast->setSweRhsAlgorithm(&w_swerhsAlgo);
  if( auto w_odeSolver1D = ptr2NumethodCast->getODESolver1D(); w_odeSolver1D)
  {
    w_odeSolver1D->registerPhysicalSystem(aPhysys);
    assert(w_odeSolver1D->getPhysicalSystem());
  }
}

...

Update global nodal values (numerical method) advance number in time by some method in this implementation a finite volume method is used to solve the math equations in their conservative form (numerical flux at cell face, mainly convective flux).

...

// Time parameters (start, stop and time step) read from simulation property file
Sfx::SfxTimePrm w_tip; // { 0., 1., 0.01);
using enum Sfx::SfxTimePrm::eTimeStepMode;
w_tip.setStartTime(Sfx::SimulationMgr::getSingleton().getStartTime());
w_tip.setStopTime(Sfx::SimulationMgr::getSingleton().getStopTime()); 
w_tip.setTimeStepMode(VARIABLE_TIME_STEP);
w_tip.initTimeLoop();
// Simulation time
aSim->setSimulationStartTime(Sfx::SimulationMgr::getSingleton().getStartTime());

...

std::cout << "Physical algorithm based on the following components:" << "\n";
std::cout << "Using HLL Riemann Solver (convective flux evaluation)" << "\n";
std::cout << "Using Nujic (1995) physics algorithm source terms" << "\n";
std::cout << "Using Time Two-steps Runge-Kutta solver family"<< "\n";
std::cout << "Using Semi-Discrete numerical method at second-order"<<"\n";
    
// update all gobal nodes
w_numethod->mainLoop(aGlbdiscr, w_tip);

...

Wave Simulator Example

Below a piece of code of the simulator showing basics steps to be coded by the user

void Wave1DSimulator::run()
  {  
    std::shared_ptr<dbpp::NumericalMethod> w_odeMethod{ new dbpp::SemiDiscreteMethod };

    // Create the global discretization for this simulation 
    auto w_gblDiscr = std::make_shared<dbpp::FiniteVolumeDiscretization>(w_odeMethod);

    // physical system that we simulate (solve dynamic equations for physical object)
    dbpp::PhysicalSystem* w_phySystem = new dbpp::PhysicalSystem; 
    w_phySystem->attach(this); // attach observer

    // physics configuration (set initial condition)                                                  
    Testvs19::DamBreakSysConfigure w_phyConfig; // configure profile and ...
    w_phyConfig.configure( w_phySystem, nullptr/*w_glbDiscr*/);

    // shared list sections (Visualizer wave profile)
    m_ListSectFlow = w_phySystem->getListSectionsPtr();
    
    // Physical algorithm
   w_phyAlgo = std::make_unique<Testvs19::EMcNeilAlgorithm>("HLLScheme"s);
   
   dbpp::PhysicalMeasurement w_phyMeasure; // physical quantities measured on the system
   dbpp::FileDataStore w_fileDatStore;     // save physical measure to the data store

   // time loop (scheduler tick) physical system state updated as we move along time line
   while( Sfx::SimulationMgr::getSingleton().getNbIterations() <
          Sfx::SimulationMgr::getSingleton().getNbIterationsMax())
   {
     ...
   update physical system ...
   w_phyAlgo->calculate(w_phySystem, w_gblDiscr, &Sfx::Simulation::getSingleton());

   // take physical measure on the physical system such as energy, velocity, ...
   w_phyMeasure.take(w_phySystem);

   // save measure to data store to eventually do some manip 
   // try-catch file open exception (take a look at Swe proto version)
   // create a final report (simulation result)
   if (!w_fileDatStore.isOpen())
   {
     w_fileDatStore.open();
   }
   w_fileDatStore.save(&w_phyMeasure, &Sfx::Simulation::getSingleton());
   w_fileDatStore.close();

...
 }