Refactoring Legacy Code Using Object-Oriented Numeric

Introduction

The migration project of the numerical library to modern C++20 is completed, and I am working on the validation of the physics algorithms. I started the rafctoring of legacy code, and I want to share a design or a prototype.   The original version of the code was written with the C-function approach; the new version of our library must be backward compatible.

More details about Object-Oriented Numeric see  Object-Oriented Numeric: A New Approach

First, I present a quick outlook of the class I have designed for the numerical model that we are solving. Briefly, the application is simulating a Dam Break wave propagation using the St-Venant equations and the finite volume method. For more details about the mathematical and numerical method follow the link below 

Mathematical and Numerical Model

To be backward compatibility with legacy code, I decided to create a class that encapsulates all those C-functions and hides them from the user. Users should see a method and not bother with details. 

Since we are solving the explicit scheme, RHS(Right-Hand_Side) terms are concerned with spatial terms, the ones that contain derivatives and other mathematical matter. 

C-functions 

Static class “BaseNumericalTreatment” (rhs discretization of each term), I put all functions in this class and made it static, because the sole purpose of the function is to calculate and return result in a numerical array. 

Convective flux algorithms (based on physics considerations) are free functions. These are particularly complex (the convective term in the equation can give rise to discontinuity in the solution and require special treatment such as shock-capturing method). 

I didn’t wrap this function in a class like others because their implementation can be very different and the signature of these functions is as well. From my point of view physics algorithms should be implemented as functions (template functions). Modern C++ added a lot of new features for generic programming (variadic template, traits, concept, perfect forwarding …).

 

—– DamBreakSim  simulator class (original code)

Below is a snippet of the original code (legacy) C-functions that compute rhs(right-hand-side) terms. Solving an explicit numerical scheme in a conservative way. It first call the convective flux algorithm, then the pressure term if desired, finally compute slope energy and then bed slope.

Then an integrator computes the solution at the next step in time.   

...

CalculFF (FF1, FF2, U1, U2, dU1, dU2, NbSections, B);
TraitementTermeP (PF2, P2, U1, NbSections, B);
TraitementTermeSource2 (S, U2, U1, H, n, dx, NbSections, B);

//	n+1/2 time step

U1p[0] = U1[0];
U2p[0] = U2[0];

U1p[NbSections-1] = U1[NbSections-1];
U2p[NbSections-1] = U2[NbSections-1];
		
for (j = 1; j < NbSections-1; j++)
{
	U1p[j] = U1[j] - dt/dx * (FF1[j] - FF1[j-1]);
	U2p[j] = U2[j] - dt/dx * (FF2[j] - FF2[j-1])
		           - dt/dx * (PF2[j] - PF2[j-1])
				   + dt * S[j];
}

...

—– Ptr2FuncFluxAlgorithm class

This class is responsible for wrapping flux algorithm pointer-2-function or C-function and calling it. Constructor takes argument on the pointer. All other methods return information about the algorithm such as reconstruction procedure of state variables and so on.  

/**
  * @brief Numerical Face Flux Algorithm (Legacy code support Pointer-2-Function) 
 */
  class Ptr2FLegacyFluxAlgorithm : public FluxAlgorithm
  {
  private:
    std::string m_ptr2fName; /**< function name*/
    int32 m_varOrder;        /**< variable order*/
    CalculFF m_ptr2func;     /**< pointer function*/
  protected:
    /** @brief Legacy code support of pointer-to-function flux algorithm.
     *  Reconstr procedure of state variables use MUSCL type at second-order (support legacy code)
    */
    ff12pair calculFF( const std::valarray<float64>& aU1, const std::valarray<float64>& aU2) const; 
  public:
    /**
     * @brief support legacy code (flux algorithm scheme pointer-to-function)
     * @param aPtr2Func pointer-to-function 
     *  Use default value variable order (support 2)
    */
    Ptr2FLegacyFluxAlgorithm(std::string aFluxAlgoName, CalculFF aPtr2Func);
    /**
     * @brief numerical face flux calculation 
     * @param aFaceVariables cell face variables
     * @return flux tensor
    */
    FluxTensor calculFF( const Sfx::cellFaceVariables& aCellFaceVar) override final;
    /**
     * @brief numerical face flux calculation
     * @param U1 first state variable
     * @param U2 second state variable
     * @param aDomain computational domain
     * @param aPhysbc physical boundary
     * @return flux tensor map
    */
    FluxTensorMap
      calculFF( const Sfx::scalarField1D& U1, const Sfx::scalarField1D& U2, //computational dommain 
        const Omega& aDomain, const PhysicalBoundaryCnd& aPhysbc) override final;
...

std::pair<valarray64_t, valarray64_t>
    Ptr2FLegacyFluxAlgorithm::calculFF(const Sfx::scalarField1D& U1, 
    const Sfx::scalarField1D& U2, const PhysicalBoundaryCnd& aPhysbc)
  {
//adapt function-to-pointer signature (takes numerical array)
    std::valarray<float64> w_U1(Sfx::DIM::value + 1);  // allocate space for ghost node
    std::valarray<float64> w_U2(Sfx::DIM::value + 1); // ditto

    // scalar field satisfy concept of a range
    if( std::ranges::range<Sfx::scalarField1D>)
    {
      if( U1.values().size() == Sfx::DIM::value)
      {
        // copy file valuestoarray
        std::ranges::copy(U1.values(), std::ranges::begin(w_U1));
        std::ranges::copy(U2.values(), std::ranges::begin(w_U2));

        // apply physical b.c. cond. at both end
        float64 w_Aupstrm{};
        float64 w_Qupstrm{};
        std::tie(w_Aupstrm, w_Qupstrm, std::ignore) = aPhysbc.getLeftEnd().Values();
        w_U1[0] = w_Aupstrm;
        w_U2[0] = w_Qupstrm;
        float64 w_Adwnstrm{};
        float64 w_Qdwnstrm{};
        std::tie(w_Adwnstrm, w_Qdwnstrm, std::ignore) = aPhysbc.getRightEnd().Values();
        w_U1[U1.values().size()] = w_Adwnstrm; // ghost node
        w_U2[U2.values().size()] = w_Qdwnstrm; // ditto
      }// if
    }

    return calculFF(w_U1, w_U1);
  }

ff12pair
    Ptr2FLegacyFluxAlgorithm::calculFF( const std::valarray<float64>& aU1,
      const std::valarray<float64>& aU2) const
  {
    // sanity check (physical boundary cond applied include ghost node)
    assert( Sfx::DIM::value+1==aU1.size());
    assert( aU2.size() == aU1.size());

    // numerical flux at cell face (computational domain is 100 nodes)
    std::valarray<float64> w_FF1(Sfx::DIM::value); //zero initialization
    std::valarray<float64> w_FF2(Sfx::DIM::value); // ditto

    // Numerical Flux basic algorithm (support legacy code pointer-to-function)
    assert( nullptr != m_ptr2func);

    // sanity check (legacy code pointer-2-function are not declare no except)
    static_assert(!std::is_nothrow_invocable_v<CalculFF, decltype(w_FF1), decltype(w_FF2), decltype(aU1), decltype(aU2)>);
    // check if can be call with the args
    if( std::is_invocable_v<CalculFF, decltype(w_FF1), decltype(w_FF2), decltype(aU1), decltype(aU2) >)
    {
      std::invoke(m_ptr2func, w_FF1, w_FF2, aU1, aU2);
    }

    return std::make_pair(std::move(w_FF1),std::move(w_FF2));
  }

—– RhsSrcFlux class

Class responsible to create the disretization of each terms (spatial derivatives, convective flux and energy slope). This class inherit from an abstract concept call “SweRhsAlgorithm“. The idea behind this concept is to set the discretization of the RHS terms and call them accordingly (see DamBreak Simulator code snippet above). Now these calls are wrapped in the caculate method (see second code snippet below) and return to user the result in a type ‘SweRhsData’ (templated aggregate). The code snippet below show class layout.

 /**
   * @brief Wrappper class to support legacy code (C-functions).
  */
  class RhsPtr2FuncFlux final : public SweRhsAlgorithm 
  {
  public:
    // should pass the physical system, retrieve physical boundary condition
    // and list of physical objects under study, section flow (next version) 
    RhsPtr2FuncFlux( std::string aFFname, CalculFF aPtr2Func, const ListSectionsFlow& aSectionList, bool aUsePressure=false);

    // shared data?? is that make sense??
    RhsPtr2FuncFlux( std::string aFFname, CalculFF aPtr2Func, SrcNumericalTreatment* aSrcTreatment,
    bool aUsePressure);

...

SweRhsData<> calculate( const Sfx::StateVectorField& aU) override final; 
/**
 * @brief compute RHS ... 
 * @param aU state vector field
 * @param aList 
 * @return RHS data (...)
*/
SweRhsData<> calculate( const Sfx::StateVectorField& aU, const ListSectionsFlow& aList);  

// Source term: pressure is not considered (i am not sure what we are doing exactly)
// incomplete flux
SweRhsData<> calculate(const Sfx::StateVectorField& aU, 
  const std::shared_ptr<FiniteVolumeDiscretization>& aGblDiscr) override final;

...

SweRhsData<> RhsPtr2FuncFlux::calculate( const Sfx::StateVectorField& aU, const ListSectionsFlow& aListSect)
{
   // Numerical flux (convective part)
    auto w_FF12 = m_ptr2FuncFF.calculFF( aU.getAfield(), aU.getQfield(),
      PhysicalBoundaryCnd{ m_leftBcNode,m_rightBcNode });

    // pressure term discretization
    [[maybe_unused]] std::valarray<float64> w_PF2(aU.Avalues().size()+1); // first-order derivative
    if( m_usePressure) // for now not active
    {
      std::valarray<float64> w_U1(aU.Avalues().size() + 1);  // allocate space for ghost node
      std::ranges::copy(aU.getAfield(), std::ranges::begin(w_U1));
      w_U1[0] = std::get<0>(m_leftBcNode.Values()); //A state variable
      w_U1[aU.Avalues().size()] = std::get<0>(m_rightBcNode.Values()); // ghost node
      std::valarray<float64> w_P2(w_U1.size());  // hydrostatic pressure
      // legacy code call pressure numerical algorithm
      BaseNumericalTreatmemt::TraitementTermeP(w_PF2, w_P2, w_U1, static_cast<int>(w_U1.size()));
    }//if-pressure

    // Source treatment (legacy code wrapper) energy slope and bottom slope 
    SbSfTermsEvaluation w_SfSb;
    auto w_sfsbArray = w_SfSb.TraitementTermeSource2(aU, aListSect,
      PhysicalBoundaryCnd{ m_leftBcNode, m_rightBcNode });

    using enum FluxTensorMap::eFFcomp;
    // RVO (Return Value Optimization)
    return SweRhsData<> { std::move(w_FF12.first), std::move(w_FF12.second), std::move(w_sfsbArray) };
}

—– Refactoring Code Simulator

Below a code snippet of the end product. It shows the calls of the different components of the API to solve the numerical model.

The simulation Manager is responsible for managing simulation configuration read from the config file. Users need to fill in the different fields of the config file before launching the simulation, such as simulation start/stop time, pointer-to-function name and so on…

Once this is completed, create an RHS discretization and call calculate(compute right hand side terms). This is equivalent to the 3 calls of the original code (see first code snippet above) and it’s all done by ‘calculate’ method. On ce this done, an explicit numerical integrator is use to advance solution in time step.

Main advantage of this is you switch between different easily. Let’s say you want to call another convective flux algorithm, all you must do is to go in the config file and set a new function name and restart the simulation.

...

// sanity check
auto w_FFname = Sfx::SimulationMgr::getSingleton().getPtrFluxAlgorithmName();
if( auto cmpName = w_FFname <=> mMethodName; cmpName != 0)
{
  Sfx::Logger::instance()->OutputError("Ptr-2-Func are not the same"s.data());
}

m_calculFF = Sfx::SimulationMgr::getSingleton().getPtr2FAlgorithm(w_FFname);
// Create RHS discretization
RhsPtr2FuncFlux w_rhsAlgo{ Sfx::SimulationMgr::getSingleton().getPtrFluxAlgorithmName(),
m_calculFF, nullptr, false }; // no source treatment/pressure, included in pointer-to-function
// set physical boundary cnd
w_rhsAlgo.setPhysicalBoundaryCnd(w_phyBCond);

// compute RHS (use std::span to initialize scalar field, cheap to move)
auto retAlgop = w_rhsAlgo.calculate( // StateVectorField
	{ std::make_shared<Sfx::scalarField1D>(grid1D,"A-field"s, 
	  std::span{rng::cbegin(mU1),Sfx::DIM::value}),
	  std::make_shared<Sfx::scalarField1D>(grid1D,"Q-field"s, 
	  std::span{rng::cbegin(mU2),Sfx::DIM::value}) }, m_listDamBreak);

// time stepping n+1/2
using enum NujicIntegrator::eIntegratorStep;
auto w_leftBCValues = w_phyBCond.getLeftEnd();
const auto& [Abc,Qbc,Hbc] = w_leftBCValues.Values();
w_twoStepsIntegrator.setLeftNodeBCValues(std::make_pair(Abc,Qbc));
w_twoStepsIntegrator.setIntegratorStep(predictor);
w_twoStepsIntegrator.step(std::move(retAlgop), m_Tip.Delta());
auto&& [U1p,U2p] = w_twoStepsIntegrator.getSolutionAsVarray();

...

Conclusion what do we gain from the refactoring? 

• Less code in the main loop (couple of calls) details done by class
• Config done externally avoid modifying the core
• Easier to extend to new types or add new algorithms without modifying the core 
• more flexibility and enhances the object-oriented separation of concerns