Class Guides

Most classes in Yggdrasil are purely virtual. This allows for the interaction between class species (like physics classes and integrators) to be abstracted and predictable.

Physics Classes

The primary class that you’ll interact with in Yggdrasil is the Physics class. This class contains the prescription for the State vector and all of the necessary logic for computing derivatives of the State vector. The general concept for how physics classes work in Yggdrasil is that the NodeList holds immutable quantities that are valid for the beginning of an integration step, while the State vector (usually a subset copy of some of the Fields on the NodeList) may update throughout a step. Yggdrasil’s integrators expect assigned (derived) physics classes to override the physics base class method EvaluateDerivatives at least, and some may also override PrestepInitialize and FinalizeStep as well.

Let’s use constantGravity.cc in the src/Physics directory as an example of a derived class that overrides only the derivatives.

template <int dim>
class ConstantGravity : public Physics<dim> {
protected:
    Lin::Vector<dim> gravityVector;
    double dtmin;
public:
    using Vector = Lin::Vector<dim>;
    using VectorField = Field<Vector>;
    using ScalarField = Field<double>;

    ConstantGravity(NodeList* nodeList, PhysicalConstants& constants, Vector& gravityVector) :
        Physics<dim>(nodeList,constants),
        gravityVector(gravityVector) {

        int numNodes = nodeList->size();
        this->template EnrollFields<Vector>({"acceleration", "velocity", "position"});
        
        for (int i=0; i<numNodes; ++i)
            nodeList->getField<Vector>("acceleration")->setValue(i,gravityVector);

        this->template EnrollStateFields<Vector>({"velocity", "position"});
    }

    ~ConstantGravity() {}

In the constructor for ConstantGravity, we pass in the NodeList pointer, a PhysicalConstants object, and a Vector that prescribes the direction and magnitude of the acceleration. Then we enroll the acceleration, velocity, and position Vector Fields using the type templated EnrollFields method from the Physics base class. This ensures that these fields exist on the NodeList, and if they don’t, this method will create them.

Next we enroll the position and velocity Fields to the State. In the case of the constant gravity package, the only derivatives are the position derivative and velocity derivative, and as a standard, we keep Fields whose derivatives we want to iterate throughout an integration step on the State object. For that reason, we do not assign acceleration to the State object for this physics class since it never changes.

    virtual void
    EvaluateDerivatives(const State<dim>* initialState, State<dim>& deriv, const double time, const double dt) override {
        NodeList* nodeList = this->nodeList;
        PhysicalConstants constants = this->constants;
        int numNodes = nodeList->size();

        VectorField* position       = initialState->template getField<Vector>("position");
        VectorField* acceleration   = nodeList->getField<Vector>("acceleration");
        // since the acceleration is literally constant, we keep it on the nodeList
        VectorField* velocity       = initialState->template getField<Vector>("velocity");

        VectorField* dxdt           = deriv.template getField<Vector>("position");
        VectorField* dvdt           = deriv.template getField<Vector>("velocity");

        dxdt->copyValues(velocity);
        dxdt->operator+(*acceleration*dt);
        dvdt->copyValues(acceleration);

        double local_dtmin = 1e30;

        #pragma omp parallel for reduction(min:local_dtmin)
        for (int i=0; i<numNodes ; ++i) {
            double amag = acceleration->getValue(i).mag2();
            double vmag = velocity->getValue(i).mag2();
            local_dtmin = std::min(local_dtmin,vmag/amag);
        }
        dtmin = local_dtmin;
        this->lastDt = dt;
    }

EvaluateDerivatives is used to compute the derivatives of the State object at each node. In this case, we are computing the change in velocity and position due to a constant acceleration, so we have a pair of ODEs:

\[\begin{split}\begin{aligned} \dot{x} &= v(t) + a \Delta t\\ \dot{v} &= a\\ \end{aligned}\end{split}\]

The integrator invokes this method with intial and derivatives State vectors and a time step, and the physics class stores the value of those derivatives back to the deriv State object.

Note

The integrator will call this method multiple times per time step if the particular integrator in question is of high temporal order. Thus, the value of dt here may be a partial step and the initial State object may change within a time step, but in the case of a simple, forward Euler integration (\(\mathcal{O}(\Delta t)\)), dt=0 means that we are computing the derivatives using State vector quantities evaluated at the current time.

The final bit of code in EvaluateDerivatives merely calculates the minimum timestep based on a velocity condition.

Generally, the base class method for FinalizeStep should suffice for most physics packages as the base class takes care of copying the final State after an integration step back to the NodeList on the derived physics class.

Note

If you don’t want to recapitulate what is already in the base class FinalizeStep method, but you’d still like to tie up some loose ends at the end of an integration step (e.g. you may want to set floor limits on certain quantities), you can do so by overriding the FinalChecks() method. By dafault, this method doesn’t do anything, but is always called after FinalizeStep.

Integrator Classes

Nihoggr’s integrators follow a basic pattern of initializing the state of the physics objects at each step, evaluating derivatives, advancing the state of each physics object from those derivatives, applying boundary conditions, and then finalizing each physics object’s state. An example of this pattern for simple forward Euler integration is shown below:

template <int dim>
Integrator<dim>::Integrator(std::vector<Physics<dim>*> packages, double dtmin, bool verbose)
    : packages(packages), dt(dtmin), dtmin(dtmin), cycle(0), time(0.0), verbose(verbose) {}

template <int dim>
Integrator<dim>::~Integrator() {}

template <int dim>
void Integrator<dim>::Step() {  
    if (cycle == 0) {
        for (Physics<dim>* physics : packages)
            physics->ZeroTimeInitialize();
    }

    for (Physics<dim>* physics : packages)
    {
        physics->UpdateState();
        physics->PreStepInitialize();
        
        State<dim> finalState = Integrate(physics);

        physics->ApplyBoundaries(&finalState);
        physics->FinalizeStep(&finalState);
    }

    time += dt;
    cycle += 1;

    VoteDt();
}

template <int dim>
State<dim> 
Integrator<dim>::Integrate(Physics<dim>* physics) {
    const State<dim>* state = physics->getState();
    State<dim> derivatives(state->size());
    derivatives.ghost(state);

    physics->EvaluateDerivatives(state, derivatives, time, 0);

    derivatives *= dt;
    State<dim> newState(state->size());
    newState = state->deepCopy();
    newState += derivatives;
    return newState;
}

template <int dim>
void Integrator<dim>::VoteDt() {
    double smallestDt = 1e30;
    for (Physics<dim>* physics : packages) {
        double newdt = physics->EstimateTimestep();
        if (newdt < smallestDt) {
            smallestDt = newdt;
            if (verbose)
                std::cout << physics->name() << " requested timestep of " << newdt << "\n";
        }
    }

    dt = (dt < smallestDt ?  dt + 0.2 * (smallestDt - dt) : smallestDt);

    this->dt = std::max(dt, this->dtmin) * this->dtMultiplier;
}

Note

While this class is a foward Euler integrator, it is also the base class for all integrators in Yggdrasil, and as such, its Step() method is rarely overridden by derived classes. Instead, the most common override will be to the Integrate() method.

Equations of State

Equations of state are possibly the simplest classes in Yggdrasil. They consume Field objects (or referenced doubles) and set the values of other Fields according to their respective closure equations. Most of the logic for how they work is self-described by the base class interface file equationOfState.hh.

// Base class for Equation of State
class EquationOfState {
public:
    EquationOfState(PhysicalConstants& constants) {}
    
    virtual ~EquationOfState() {}

    // Field-based methods

    virtual void 
    setPressure(Field<double>* pressure, Field<double>* density, Field<double>* internalEnergy) const = 0;

    virtual void 
    setInternalEnergy(Field<double>* internalEnergy, Field<double>* density, Field<double>* pressure) const = 0;

    virtual void 
    setSoundSpeed(Field<double>* soundSpeed, Field<double>* density, Field<double>* internalEnergy) const = 0;

    virtual void 
    setTemperature(Field<double>* temperature, Field<double>* density, Field<double>* internalEnergy) const = 0;

    virtual void 
    setInternalEnergyFromTemperature(Field<double>* internalEnergy, Field<double>* density, Field<double>* temperature) const = 0;

    // Scalar-based methods

    virtual void 
    setPressure(double* pressure, double* density, double* internalEnergy) const = 0;

    virtual void 
    setInternalEnergy(double* internalEnergy, double* density, double* pressure) const = 0;

    virtual void 
    setSoundSpeed(double* soundSpeed, double* density,double* internalEnergy) const = 0;

    virtual void 
    setTemperature(double* temperature, double* density, double* internalEnergy) const = 0;

    virtual void 
    setInternalEnergyFromTemperature(double* internalEnergy, double* density, double* temperature) const = 0;

    virtual std::string 

Note

Equations of state are not typically templated in Yggdrasil since they act only on scalar Fields.