Solver is called every timestep. Derived from ISolver. Most SPH solvers further derive from GenericSolver.

Code currently implements following solvers:

1. ContinuitySolver

   Standard sets of hydrodynamical equations

2. SummationSolver

   Density is solved by direct summation rather than using ContinuityEquation

3. DensityIndependentSolver

   Experimental solver, attempts to resolver some issues of SPH, such as artificial pressure at contact discontinuity

4. GravitySolver

   Extension of ContinuitySolver also including self-gravity

Specifying equations to solve

Unless specified otherwise, the solver and the solved set of equations is determined from the parameters of RunSettings. If needed, however, solver can set manually and assigned to Run::solver. This way we can specify the set of equations explicitly.

Modifying the set of equations is fairly straighforward. Usually, we wish to add an equation to the existing set. If ContinuitySolver or SummationSolver are used, we can specify additional equations in constructor. Note that this does not affect the default set of equations, it only adds additional terms.

If we wish to specify the whole set of equations, we can utilize GenericSolver. This solver is a 'template' (not in C++ sense) that solves the set of equations given in the constructor. It does not solve any equations by default, except for helper 'dummy' equation NeighbourCountTerm. It therefore allows to create any combination of equations, leave out certain terms, etc.

Attention

Currently, it is necessary to specify an evolution equation for smoothing length. Unlike other quantities that remain constain in time if no evolution equation is specified, the smoothing length will change incorrectly. To keep the smoothing length constant, add ConstSmoothingLength into the list of equations. Usually, the better alternative is AdaptiveSmoothingLength that allows for variable spatial resolution of the simulation. This problem will be resolted in the future.

The following code adds an internal friction and non-intertial accelerations (centripetal force and Coriolis force) into the set of equations. Add it to your Run::setUp function, note that it overrides previous solver, of course.

#include "sph/solvers/ContinuitySolver.h" 
#include "sph/equations/Friction.h"
#include "sph/equations/Potentials.h"
 
void Run::setUp() {
 
    // Container holding the additional equations. Do not add equations already used in the solver!
    EquationHolder equations;
    
    // Add internal friction. The kinematic viscosity is taken from parameter BodySettingsId::KINEMATIC_VISCOSITY.
    equations += makeTerm<InternalFriction>();
 
    // Add non-intertial acceleration. It takes the angular frequency as a parameter.
    // Here we use rotation with period P = 10s. Accelerated rotation is currently not implemented.
    const Vector omega(0._f, 0._f, 2._f * PI / 10._f);
    equations += makeTerm<NoninertialForce>(omega);
 
    // Construct the solver, passing the additional forces as the second argument.
    solver = makeAuto<ContinuitySolver>(settings, equations);
 
    // Now we must construct all the necessary quantities used by she solver. We can either use class
    // InitialConditions for this, using its constructor that takes the solver, or we can do it manually as shown here.
 
    // The quantities must be created for each material in the run separately, as different materials
    // might behave differently. Here we assume there is only a single material in the storage.
    IMaterial& material = storage->getMaterial(0);
    solver->create(*storage, material);
 
    // It is also necessary to construct the material itself. This creates quantities needed for equation
    // of state (in case these quantities are not created by the solver), and also sets up flaws in the 
    // material if a fragmentation model is used.
    MaterialInitialContext context;
    // Select a random number generator. Needed to distribute the flaws, can be skipped if no fragmentation
    // is used.
    context.rng = makeAuto<RngWrapper<UniformRng>>(1234);
 
    // Create the material. Note that this is done by InitialCondition object if used.
    material.create(*storage, context);
 
    // Also set up initial conditions, output files, ...
    // See cli/main.cpp for example
};