Programming Physics Algorithm Using DamBreak++

DamBreak++ Simulation Framework

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. 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.

Scientific Framework (ODE Solvers)

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

Physical Algorithm Components Based to Build Numerical Method

See the folowing document for details about Mathematical And Numerical Method

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).We want to combine different numerical treatment for each of these terms as well as numerical methods. 

The main purpose of the PhysicalAlgorithm abstract class is to provide the guts of the simulation (using components to program complex algorithms quickly).   

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();

...
 }