Crate corries

Corries - CORrosive RIEman Solver

A library/framework (honestly what’s the difference nowadays anyways) to run 1D-hydrodynamics simulations solved with Riemann solvers. As such, corries is grid-based (in opposition to doing smooth particle hydrodynamics).

Corries should be ready to use for simple shocktube-like simulations at this point.

Disclaimer about nomenclature

Just a couple of notes about how things are named:

• corries uses generalised coordinates, and names them xi, eta, and Phi
• the coordinate systems in corries are always right-handed
• corries may be specialised for 1D simulations, but the other dimensions are still important, even though they are only 1 cell wide
• xi is the primary coordinate, the one that has more than one cell. That’s why the coordinate vectors in Mesh are called xi_*, for example.

Dependencies

Corries has one dependency that you need to use in your apps, which is color_eyre. Thus, in order to use corries, you need add to lines to the dependencies in your Cargo.toml:

[dependencies]
color-eyre = { version = "0.6", default-features = false }
corries = { version = "0.4" }

corries also has one default feature: validation. This feature performs run-time checks on the state of your simulation, like checking for finite values, positive mass densities and pressures, and the like.

As useful as they are, they do slow down the simulation. Such you can disable it by turning off the the default-features by modifying the corries entry like this:

corries = { versoin = "0.4", default-features = false }

Usage

The main building blocks of preparing and running a simulation with corries are:

• Choosing a couple of type parameters, namely for the traits: Physics, NumFlux, and TimeSolver, as well as the Mesh size (also through a compile time constant)
• Building a CorriesConfig, which allows for fine-grained control over the simulation
• Constructing and applying initial conditions to the necessary objects through the CorriesConfig::init_corries
• Running the simulation with the CorriesComponents’s Runner::run_corries method

You can see examples of this in the docs for CorriesConfig::init_corries and Runner::run_corries, as well as when looking through the source code for the integration tests in the tests/ folder. In those examples I usually use a couple of default constructors for configuration structs, but in examples/template.rs you can also see how a full-blown CorriesConfig can be set up.

Building a Sod test

Let’s walk through building a Sod test.

The first thing we should do is setting some compile constants and type definitions. The type defs are not necessary, but I highly recommend adding them to make your code more flexible when you decide to switch out some of these type defs.

use color_eyre::Result;
use corries::prelude::*;

// The number of grid cells in the Mesh
const S: usize = 100;

// Expands to:
// type P = Euler1DAdiabatic<S>;
// const E: usize = P::NUM_EQ; // which in turn equals 3
set_Physics_and_E!(Euler1DAdiabatic);

// The type for the numerical flux scheme
type N = Kt<E, S>;

// The type for the time discretisation scheme
type T = RungeKuttaFehlberg<P, E, S>;

To break this down, first we have a couple of use statements. corries features a prelude module that imports everything you need, and then you also need to use color_eyre::Result, because that is the return type for a lot of functions in corries.

Then we need to set the size of your mesh, by definining const S: usize. In this case, we set the mesh to have 100 grid cells (including ghost cells).

Next, we set up our Physics type. This determines the differential equations that are being solved in the simulation. You can check out the current options in the Implementors section of Physics. This macro also sets up const E: usize to be the number of equations for the Physics implementor you chose.

Then we set N to the NumFlux implementor we want, and T to the TimeSolver implementor. Again, check the docs for those traits to see the options.

Next up, we need to set up our CorriesConfig. This is probably the most complicated part, because there are so many options. In case this piece of the docs ever gets outdated, I would implore you to the CorriesConfig docs to check out the options you have, as well as the default constructors for that struct as well as some fields in it.

Either way, a hand-rolled CorriesConfig for this type of simulation would look like this:

let config: CorriesConfig = CorriesConfig {
    // prints a welcome banner upon running [init_corries]
    print_banner: true,

    // Sets up a cartesian mesh, with the edges of the computational area (not counting ghost
    // cells) are at the coordinates `xi = 1.0` and `xi = 2.0`. Remember that `xi` is the name
    // of the primary coordinate in corries.
    mesh_config: MeshConfig {
        mode: MeshMode::Cartesian,
        xi_in: 1.0,
        xi_out: 2.0,
    },

    // Since we want to run an adiabatic simulaion, we need to set the adiabatic index.
    // For the Sod test, let's set it to `1.4`.
    // We can also set the type of units we want, but there is not much reason for it in
    // corries right now (yes, I should just remove it for now...). Just set it to
    // [UnitsMode::SI].
    physics_config: PhysicsConfig {
        adiabatic_index: 1.4,
        units_mode: UnitsMode::SI,
    },

    // Sets the boundary conditions for the west (inner) edge of the computational area.
    // [BoundaryMode::Custom] allows you to set individual boundary conditions for
    // every equation by putting pairs of the equation index the the custom boundary
    // mode in the embedded vector. Check out the documentation for the [Physics]
    // implementor you chose to see which indexes correspond to which equations, and check out
    // the documentation for [CustomBoundaryMode] to see what the available conditions are.
    //
    // Alternatively, this field could have been set to:
    // `boundary_condition_west: BoundaryMode::NoGradients`
    // This would set the `NoGradients` condition for all equations, without setting up this
    // vector. See `boundary_condition_east` to see how that would look like.
    boundary_condition_west: BoundaryMode::Custom(vec![
        (0, CustomBoundaryMode::NoGradients),
        (1, CustomBoundaryMode::NoGradients),
        (2, CustomBoundaryMode::NoGradients),
    ]),

    // Sets the boundary conditions for the east (outer) edge of the computational area.
    boundary_condition_east: BoundaryMode::NoGradients,

    // Sets up everything regarding numerics.
    numerics_config: NumericsConfig {

        // Sets up configuration for the numerical flux scheme.
        // [NumFluxConfig] is an enum that needs to correspond to the [NumFlux] implementor you
        // chose up top, otherwise the constructor for the [NumFlux] object will panic.
        //
        // `Kt` should be your go-to. It's fast and accurate, without running into Godunov's
        // order barrier since we only use linear limiting functions.
        numflux_config: NumFluxConfig::Kt {
            // When choosing `Kt`, you also need to set the limiter function here. Your go-to
            // should be the `VanLeer` limiter.
            limiter_mode: LimiterMode::VanLeer,
        },

        // Sets up the time integration scheme.
        // Currently, corries only supports Runge-Kutta-Fehlberg schemes, which are set here.
        time_integration_config: TimeIntegrationConfig::Rkf(RkfConfig {
            // The only important bit about this config is the exact scheme you want to use.
            // SSPRK5 and RKF4 are the go-to choices, though I would recommend the first
            // usually.
            rkf_mode: RKFMode::SSPRK5,

            // Whether or not to use automatic step control. Not really needed for this test,
            // so we set this to false.
            asc: false,

            // The relative tolerance when calculating the error between high and low order
            // solutions, used by the automatic step control.
            asc_relative_tolerance: 0.001,

            // The absolute tolerance when calculating the error between high and low order
            // solutions, used by the automatic step control.
            asc_absolute_tolerance: 0.1,

            // The timestep friction factor. This dampens new time steps widths calculated by
            // the automatic step control.
            asc_timestep_friction: 0.08,
        }),

        // How many full loop iterations corries is allowed to run before aborting the
        // simulation, if it has not ended due to other break conditions yet.
        iter_max: usize::MAX - 2,

        // The time coordinate at which the simulation starts.
        t0: 0.0,

        // The time coordinate at which the simulation ends. We set this to 0.25, so that the
        // the simulation stays contained within our spatial bounds.
        t_end: 0.25,

        // If the calculated time step width falls below this value, the simulation ends with
        // an error. This is to prevent simulations that would run on too long because
        // something is driving the time step width into the ground.
        dt_min: 1.0e-12,

        // The maximum value for the time step width. If the calculated time step width exceeds
        // this value, former will simply be capped to this.
        dt_max: f64::MAX,

        /// The CFL (Courant-Friedrichs-Lewy) condition parameter. Basically a safety factor to
        // dampen the time step widths calculated by the CFL condition.
        dt_cfl_param: 0.4,
    },

    // How many times should [Writer] write outputs during the simulation (not counting the
    // initial output).
    // These are evenly distributed throughout the simulation time, so for example, if you set
    // `numerics_config.t0 = 0.0` and `numerics_config.t_end = 10.0`, then setting
    // `output_counter_max = 10` would mean that corries writes output at
    // `t = [1.0, 2.0, ..., 10.0]` in addition to the extra output at `t = 0.0`.
    output_counter_max: 10,

    // Sets up how outputs are generated.
    // This is a vector, where every entry is a instance of [OutputConfig], and each of these
    // instances sets up one output stream. In this example, we set up to streams: one for
    // stdout (the first entry), and one for files (the second entry). You are free to mix and
    // match different output streams as you like, but this is probably a good go-to.
    writer_config: vec![

        // Sets up a stdout output stream, with TSV formatting (tab separated values) for
        // scalar values, i.e. values that do not differ for different parts of the
        // computational area.
        //
        // In most cases, setting this stream up like this is probably too verbose, so
        // OutputConfig provides constructor that you only need to pass the data_names field
        // to.
        // This same setup could have been achieved by replacing this initialiser with:
        // OutputConfig::default_stdout_with_names(
        //     vec![DataName::Iter, DataName::T, DataName::Dt, DataName::DtKind]
        // ),
        //
        // You can also replicate this setup with the default_stout constructor:
        //
        // OutputConfig::default_stdout(),
        //
        // that sets up data_names exactly like this too.
        OutputConfig {
            // Sets up that this stream writes to stdout
            stream_mode: StreamMode::Stdout,

            // Formats the output vales such that the values are separated by tabs
            formatting_mode: FormattingMode::TSV,

            // Tells the stream to only print "global" data, i.e. values that are the same
            // everywhere in the computational area.
            string_conversion_mode: ToStringConversionMode::Scalar,

            // Ignored for stdout output; for file output this field sets the folder the
            // simulation data should be written to.
            folder_name: "".to_string(),

            // Ignored for stdout output; Whether to clear the contents of that folder during
            // initialisation.
            // THIS CAN LEAD TO LOSS OF DATA SO BE VERY CAREFUL WHAT FOLDER_NAME POINTS TO WHEN
            // SETTING THIS VALUE TO true!
            should_clear_out_folder: true,

            // Irrelevant for stdout output; for file output this field sets the base file name
            // for the simulation data.
            file_name: "".to_string(),

            // Floating point precision when printing floating point numbers
            precision: 3,

            // Ignored for scalar output; for vector output, this sets whether to include the
            // ghost cells in the output.
            should_print_ghostcells: false,

            // Whether to print the configuration metadata. For file based output, this
            // generates a <file_name>__metadata.json next to your simulation data, which is the
            // deserialised [CorriesConfig] for this simulation.
            should_print_metadata: true,

            // Defines which values should be printed to the output.
            // These values will be printed from left to right in the output. Check out the
            // docs for [DataName] for the options you have.
            data_names: vec![DataName::Iter, DataName::T, DataName::Dt, DataName::DtKind],
        },

        // Sets up a file output stream.
        //
        // There is also a handy constructor that sets up a default file output.
        // You could also define this exact same setup by constructing the OutputConfig like
        // this:
        // OutputConfig::default_file_with_names(
        //     "results/sod",
        //     "sod",
        //     vec![DataName::T, DataName::XiCent, DataName::Prim(0), DataName::Prim(1), DataName::Prim(2)],
        // ),
        //
        // There is also OutputConfig::default_file() that sets the data_names field like this
        // automatically. You would call it like this:
        //
        // OutputConfig::default_file("results/sod", "sod", E),
        //
        // where E is the number of equations (we set this through the set_Physics_and_E macro
        // up top. This function would also replicate this exact setup.
        OutputConfig {
            // Configures this object to write to a file
            stream_mode: StreamMode::File,

            // Formats the data into csv files
            formatting_mode: FormattingMode::CSV,

            // Defines that the data is written as vectors. This allows the output to write
            // data like DataName::XiCent or DataName::Prim(n), which are vector-like values.
            string_conversion_mode: ToStringConversionMode::Vector,

            // Which folder to write the data to, relative to the cwd when calling the
            // executable.
            folder_name: "results/sod".to_string(),

            // Whether to clear out the folder defined in `folder_name` during initialisation.
            // THIS CAN LEAD TO LOSS OF DATA SO BE VERY CAREFUL WHAT FOLDER_NAME POINTS TO WHEN
            // SETTING THIS VALUE TO true!
            should_clear_out_folder: true,

            // The base file name for output files.
            // In this setup, your output files would be named:
            // `results/sod/sod_{00,01,02,..,10}.csv`
            file_name: "sod".to_string(),

            // Floating point precision for floating point numbers
            precision: 7,

            // Whether to include ghostcells when writing vector-like values.
            should_print_ghostcells: true,

            // Whether to write the metadata dump. In this setup, this file would be called:
            // `results/sod/sod__metadata.json`
            should_print_metadata: false,

            // Which values to include in this output, read from left to right.
            data_names: vec![
                DataName::T, DataName::XiCent, DataName::Prim(0), DataName::Prim(1), DataName::Prim(2)
            ],
        },
    ],
};

Now with our configuration done, we can initialise the objects we need. The CorriesConfig::init_corries() method does exactly what you need here, and all you need to pass it are the generic parameters we set up top, as well as a function for your initial conditions.

That function has the signature (barring generic parameters):

fn(&mut State, &mut Solver, &Mesh) -> Result<()>

As was already stated, the purpose of it is to apply your initial conditions to the State object.

The Sod test for this set of equations is set up like this:

• net zero velocities
• a density and pressure jump in the horizontal centre of the shocktube:
  • To the left of that shock the mass density and pressure are set to 1.0
  • To the right of that shock the mass density is set to 0.125 and the pressure is set to 1.0

When setting up initial conditions we generally only have to set up the primitive variables for the cell centric variable set. The CorriesConfig::init_corries() method handles applying the boundary conditions and making sure that State::west and State::east are set up too.

config.init_corries::<P, N, T, E, S>(|u, _, _| {
    // this function does need the [Solver] and [Mesh] objects, hence the two `_` in the
    // arguments.
    let breakpoint_index = (S as f64 * 0.5) as usize;
    for i in 0..breakpoint_index {
        // set mass density (rho) and pressure of the cell central values left of the shock to
        // 1.0
        u.cent.prim[[P::JRHO, i]] = 1.0;
        u.cent.prim[[P::JPRESSURE, i]] = 1.0;

        // set the xi velocity to 0.0
        u.cent.prim[[P::JXI, i]] = 0.0;
    }
    for i in breakpoint_index..S {
        // set the mass density to the right of the shock to 0.125
        u.cent.prim[[P::JRHO, i]] = 0.125;

        // set the pressure to the right of the shock to 0.1
        u.cent.prim[[P::JPRESSURE, i]] = 0.1;

        // set the xi velocity to 0.0
        u.cent.prim[[P::JXI, i]] = 0.0;
    }
    Ok(())
}).unwrap();

Obviously, you should be using the ? operator instead of unwrap, but then this doc test does not run correctly… ¯_(ツ)_/¯

The object returned here is a CorriesComponents, which is simply a tuple around:

• u: the State object that holds the state of the physical variables
• solver: the Solver object that generates new solutions
• mesh: the Mesh object that models the grid the simulation runs on
• writer: the Writer object that handles writing output

Almost done. Now we can call Runner::run_corries to run the simulation!

config.init_corries::<P, N, T, E, S>(|u, _, _| {
    // this function does need the [Solver] and [Mesh] objects, hence the two `_` in the
    // arguments.
    let breakpoint_index = (S as f64 * 0.5) as usize;
    for i in 0..breakpoint_index {
        // set mass density (rho) and pressure of the cell central values left of the shock to
        // 1.0
        u.cent.prim[[P::JRHO, i]] = 1.0;
        u.cent.prim[[P::JPRESSURE, i]] = 1.0;

        // set the xi velocity to 0.0
        u.cent.prim[[P::JXI, i]] = 0.0;
    }
    for i in breakpoint_index..S {
        // set the mass density to the right of the shock to 0.125
        u.cent.prim[[P::JRHO, i]] = 0.125;

        // set the pressure to the right of the shock to 0.1
        u.cent.prim[[P::JPRESSURE, i]] = 0.1;

        // set the xi velocity to 0.0
        u.cent.prim[[P::JXI, i]] = 0.0;
    }
    Ok(())
}).unwrap().run_corries().unwrap();

Plans

Currently, corries only supports cartesian meshes, though non-cartesian meshes are coming very soon.

• Add cylindrical geometry + Source trait + GeometricSource + Sedov test
• Add spherical geometry + new Sedov case
• Add 2D adiabatic Euler and 2D isothermal Euler
• Add adiabatic and isothermal Navier Stokes physics
• Add logcylindrical geometry + Vortex test
• Add gravitational source + Bondi test
• Add viscosity source + Pringle test
• investigate the possibility of putting boundary conditions into State
  • not really possible currently, because I cannot split the borrows between the boundary conditions and the variables without using Cell for interior mutability.

Most importantly exports the module prelude.

Re-exports

pub use prelude::*;

Modules

components

Exports the CorriesComponents type alias, as well as Runner trait.

config

Exports the CorriesConfig structs and its nested structs for configuring Corries simulations.

directions

Exports the Direction enum

macros

Exports useful macros for public use

mesh

Exports the Mesh struct that provides coordinates and geometric information

prelude

Exports everything you need to run a corries simulation. This includes the following modules

rhs

Exports the Rhs struct that carries objects and methods for solving the right-hand side of a set of equations.

solver

Exports the Solver struct, responsible for generating new solutions for State objects

state

Exports the State struct

time

Exports the TimeSolver trait and the DtKind enum.

units

Exports the UnitsMode and Units structs for configuring and using unit systems.

writer

Exports the Writer struct that deals with writing data to multiple different output streams

Macros

check_finite_arrayd

Macro to check that each Array in a list is finite.

check_finite_double

Macro to check that each double in a list is finite.

check_positive_arrayd

Macro to check that each Array has only positive non-zero elements.

check_positive_double

Macro to check that each double in a list is finite.

set_Physics_and_E

Expands to a type alias P for the type of physics you are using, and a constant E that represents the number of equations in that system. This macro is used to make sure that E is always set correctly for your type of physics, while also giving you a useful type alias to make your code look more generic than it actually is.