diff --git a/code_students/apps/CMakeLists.txt b/code_students/apps/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f88fdd9dca93574070ae021f60e3a2679352047b
--- /dev/null
+++ b/code_students/apps/CMakeLists.txt
@@ -0,0 +1,41 @@
+cmake_minimum_required(VERSION 3.10)
+
+
+add_executable (run_full_code
+ main_full_code.cpp
+ )
+
+target_include_directories(run_full_code
+ PUBLIC ${CODE_INCLUDE_DIR}
+)
+
+target_link_libraries(run_full_code
+ PUBLIC utilities
+ PUBLIC sim_setup
+ PUBLIC solver
+)
+
+install(TARGETS run_full_code
+ RUNTIME DESTINATION ${PROJECT_SOURCE_DIR}/bin
+)
+
+
+
+add_executable (run_singlestep
+ main_test_singlesteps.cpp
+ )
+
+target_include_directories(run_singlestep
+ PUBLIC ${CODE_INCLUDE_DIR}
+)
+
+target_link_libraries(run_singlestep
+ PUBLIC utilities
+ PUBLIC sim_setup
+ PUBLIC solver
+ PUBLIC IO
+)
+
+install(TARGETS run_singlestep
+ RUNTIME DESTINATION ${PROJECT_SOURCE_DIR}/bin
+)
\ No newline at end of file
diff --git a/code_students/apps/main_full_code.cpp b/code_students/apps/main_full_code.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5370f9745b19ae61102bb099e51fecb63d67b522
--- /dev/null
+++ b/code_students/apps/main_full_code.cpp
@@ -0,0 +1,88 @@
+#include "core/config.hpp"
+#include "setup/fluid.hpp"
+#include "setup/grid.hpp"
+#include "solver/finite_volume_solver.hpp"
+#include "solver/time_integrator.hpp"
+#include "util/matrix.hpp"
+#include "util/utility_functions.hpp"
+
+#include <cmath>
+#include <functional>
+#include <iostream>
+
+double Sedov_volume;
+
+void init_Sedov(fluid_cell &fluid, double x_position, double y_position, double z_position) {
+ double radius = sqrt(sim_util::square(x_position) + sim_util::square(y_position) + sim_util::square(z_position));
+ double radius_init = 0.05;
+ double volume_init = 4.0 / 3.0 * M_PI * radius_init * radius_init * radius_init;
+ volume_init = Sedov_volume;
+ double E_init = 1.0;
+ double e_dens_init = E_init / volume_init;
+ if (radius < 0.1) {
+ fluid.fluid_data[fluid.get_index_energy()] = e_dens_init;
+ fluid.fluid_data[fluid.get_index_tracer()] = 1.0;
+ } else {
+ fluid.fluid_data[fluid.get_index_energy()] = 1.e-5 / 0.4;
+ fluid.fluid_data[fluid.get_index_tracer()] = 0.0;
+ }
+ fluid.fluid_data[fluid.get_index_density()] = 1.0;
+ fluid.fluid_data[fluid.get_index_v_x()] = 0.0;
+ fluid.fluid_data[fluid.get_index_v_y()] = 0.0;
+ fluid.fluid_data[fluid.get_index_v_z()] = 0.0;
+}
+
+int main() {
+ std::vector<double> bound_low(3), bound_up(3);
+ bound_low[0] = -0.5;
+ bound_low[1] = -0.5;
+ bound_low[2] = -0.5;
+
+ bound_up[0] = 0.5;
+ bound_up[1] = 0.5;
+ bound_up[2] = 0.5;
+
+ std::vector<int> num_cells(3);
+ num_cells[0] = 128;
+ num_cells[1] = 128;
+ num_cells[2] = 128;
+
+ grid_3D my_grid(bound_low, bound_up, num_cells, 2);
+
+ // Get number of Sedov cells
+ Sedov_volume = 0.0;
+ int num_Sedov_cells = 0;
+ double volume_cell = my_grid.x_grid.get_dx() * my_grid.y_grid.get_dx() * my_grid.z_grid.get_dx();
+
+ for (int ix = 0; ix < my_grid.get_num_cells(0); ++ix) {
+ double x_position = my_grid.x_grid.get_center(ix);
+ for (int iy = 0; iy < my_grid.get_num_cells(1); ++iy) {
+ double y_position = my_grid.y_grid.get_center(iy);
+ for (int iz = 0; iz < my_grid.get_num_cells(2); ++iz) {
+ double z_position = my_grid.z_grid.get_center(iz);
+ double dist = sqrt(sim_util::square(x_position) + sim_util::square(y_position) + sim_util::square(z_position));
+ if (dist < 0.1) {
+ Sedov_volume += volume_cell;
+ num_Sedov_cells++;
+ }
+ }
+ }
+ }
+ std::cout << " Volume of Sedov region: " << Sedov_volume << " in " << num_Sedov_cells << " cells\n";
+
+ // Now, I will create a HD fluid
+ fluid hd_fluid(parallelisation::FluidType::adiabatic);
+ hd_fluid.setup(my_grid);
+
+ std::function<void(fluid_cell &, double, double, double)> function_init = init_Sedov;
+
+ finite_volume_solver solver(hd_fluid);
+ solver.set_init_function(function_init);
+
+ double t_final = 0.1;
+ double dt_out = 0.005;
+
+ solver.run(my_grid, hd_fluid, t_final, dt_out);
+
+ return 0;
+}
\ No newline at end of file
diff --git a/code_students/apps/main_test_singlesteps.cpp b/code_students/apps/main_test_singlesteps.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..06d94ce2b83418b43dc8528cbee78522205ff2be
--- /dev/null
+++ b/code_students/apps/main_test_singlesteps.cpp
@@ -0,0 +1,64 @@
+#include "IO/data_storage.hpp"
+#include "core/config.hpp"
+#include "setup/fluid.hpp"
+#include "setup/grid.hpp"
+#include "solver/finite_volume_solver.hpp"
+#include "solver/time_integrator.hpp"
+#include "util/matrix.hpp"
+
+#include <iostream>
+
+int main() {
+ std::cout << " Trying to create a matrix \n";
+ size_t dims_2D[2] = {10, 10};
+ size_t dims_1D[1] = {10};
+ matrix<double, 1> my_matr_1D(dims_1D);
+ matrix<double, 2> my_matr_2D(dims_2D);
+
+ // my_matr_2D(2,3) = 22.2;
+ // std::cout << " Values " << my_matr_2D(1,1) << " ";
+ // std::cout << my_matr_2D(2,3) << "\n";
+
+ std::vector<double> bound_low(3), bound_up(3);
+ bound_low[0] = 0.0;
+ bound_low[1] = 0.0;
+ bound_low[2] = 0.0;
+
+ bound_up[0] = 1.0;
+ bound_up[1] = 1.0;
+ bound_up[2] = 1.0;
+
+ std::vector<int> num_cells(3);
+ num_cells[0] = 40;
+ num_cells[1] = 30;
+ num_cells[2] = 20;
+
+ grid_3D my_grid(bound_low, bound_up, num_cells, 2);
+
+ // Now, I will create a HD fluid
+ fluid hd_fluid(parallelisation::FluidType::adiabatic);
+ hd_fluid.setup(my_grid);
+
+ fluid hd_changes(hd_fluid.get_fluid_type());
+ hd_changes.setup(my_grid);
+
+ finite_volume_solver solver(hd_fluid);
+
+ RungeKutta2 time_steper(my_grid, hd_fluid);
+
+ double delta_t = 0.1;
+
+ // Do two RK steps:
+ solver.singlestep(my_grid, hd_fluid, hd_changes);
+ time_steper.do_sub_step(my_grid, hd_changes, hd_fluid, delta_t, 0);
+ solver.singlestep(my_grid, hd_fluid, hd_changes);
+ time_steper.do_sub_step(my_grid, hd_changes, hd_fluid, delta_t, 1);
+
+ // Testing IO
+ data_storage my_storage("test.h5");
+ std::vector<double> raw_data = hd_fluid.fluid_data[0].get_raw_data();
+ std::vector<size_t> extent = hd_fluid.fluid_data[0].get_dims();
+ my_storage.write_dataset(raw_data, extent, "Density");
+
+ return 0;
+}
\ No newline at end of file
diff --git a/code_students/include/IO/data_storage.hpp b/code_students/include/IO/data_storage.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e2744a30a6622171e0adc9d112ebd89487c8cdc9
--- /dev/null
+++ b/code_students/include/IO/data_storage.hpp
@@ -0,0 +1,37 @@
+#ifndef DATA_STORAGE_HPP
+#define DATA_STORAGE_HPP
+
+#include <hdf5.h>
+#include <string>
+#include <type_traits>
+#include <vector>
+
+// Some definitions for convenience
+template <typename T> using Invoke = typename T::type;
+template <typename Condition> using negation = std::integral_constant<bool, !bool(Condition::value)>;
+template <typename Condition> using EnableIf = Invoke<std::enable_if<Condition::value>>;
+template <typename Condition> using DisableIf = EnableIf<negation<Condition>>;
+
+// Do not deduce from templates for pointers
+#define template_noPointer template <typename T, typename = DisableIf<std::is_pointer<T>>>
+
+template <typename T> static hid_t get_hdf5_data_type();
+
+template <> hid_t get_hdf5_data_type<int>() { return H5Tcopy(H5T_NATIVE_INT); }
+template <> hid_t get_hdf5_data_type<float>() { return H5Tcopy(H5T_NATIVE_FLOAT); }
+template <> hid_t get_hdf5_data_type<double>() { return H5Tcopy(H5T_NATIVE_DOUBLE); }
+
+class data_storage {
+public:
+ data_storage(std::string file_name);
+ void write_dataset(const std::vector<double> &raw_data, const std::vector<size_t> &extents, std::string dataset_name);
+ template_noPointer void AddGlobalAttr(const std::string &AttrName, T AttrData);
+
+private:
+ void open_file();
+ void close_file();
+ hid_t hdf5file, hdf5group;
+ std::string file_name;
+};
+
+#endif
\ No newline at end of file
diff --git a/code_students/include/core/config.hpp b/code_students/include/core/config.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..7f48d2c1005427c57708c2661e939ecc2f858732
--- /dev/null
+++ b/code_students/include/core/config.hpp
@@ -0,0 +1,14 @@
+#ifndef CONFIG_HPP
+#define CONFIG_HPP
+
+namespace parallelisation {
+
+enum class direction { x, y, z };
+
+enum class FluidType { isothermal, adiabatic, none };
+
+enum class BoundaryType { constant_extrapolation };
+
+} // namespace parallelisation
+
+#endif
\ No newline at end of file
diff --git a/code_students/include/setup/fluid.hpp b/code_students/include/setup/fluid.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..7d8e4d23f6294ea702fb8fb61ff9bfaeedd1cf2f
--- /dev/null
+++ b/code_students/include/setup/fluid.hpp
@@ -0,0 +1,90 @@
+#ifndef FLUID_HPP
+#define FLUID_HPP
+
+#include "core/config.hpp"
+#include "setup/grid.hpp"
+#include "util/matrix.hpp"
+
+#include <map>
+#include <memory>
+#include <vector>
+
+class base_fluid {
+
+ enum class FieldName { density, v_x, v_y, v_z, energy, tracer, none };
+
+public:
+ base_fluid();
+ base_fluid(const parallelisation::FluidType &type);
+ void setup();
+ int get_field_index(const FieldName &field);
+ int get_index_density() const;
+ int get_index_v_x() const;
+ int get_index_v_y() const;
+ int get_index_v_z() const;
+ int get_index_energy() const;
+ int get_index_tracer() const;
+ int get_num_fields() const;
+ bool is_adiabatic() const;
+ void set_characteristic();
+ void set_conservative();
+ bool is_conservative() const;
+ parallelisation::FluidType get_fluid_type() const;
+
+private:
+ parallelisation::FluidType type;
+ std::map<FieldName, int> map_fields;
+ int index_density, index_v_x, index_v_y, index_v_z;
+ int index_energy, index_tracer;
+ bool uses_characteristic_quantities;
+};
+
+class fluid_cell : public base_fluid {
+public:
+ fluid_cell(const parallelisation::FluidType &type);
+
+ void setup();
+
+ /**
+ * Array holding fluid data
+ */
+ std::vector<double> fluid_data;
+
+private:
+};
+
+class fluxes_cell : base_fluid {
+public:
+ fluxes_cell(const parallelisation::FluidType &type);
+
+ void setup();
+
+ /**
+ * Array holding fluxes
+ */
+ std::vector<double> flux_data;
+};
+
+class fluid : public base_fluid {
+public:
+ fluid(const parallelisation::FluidType &type);
+
+ void setup(grid_3D &spatial_grid);
+ /**
+ * Array of matrices holding the actual data
+ */
+ std::vector<matrix<double, 3>> fluid_data;
+
+ /**
+ * Get number of fields stored in fluid_data
+ */
+ size_t get_number_fields() const;
+ void get_fluid_cell(fluid_cell &local_fluid, int index_x, int index_y, int index_z) const;
+ void set_cell_values(const fluid_cell &local_fluid, int index_x, int index_y, int index_z);
+ void get_fluid_cell_raw(fluid_cell &local_fluid, int index_1D) const;
+ void set_fluid_cell_raw(fluid_cell &local_fluid, int index_1D);
+
+private:
+};
+
+#endif
\ No newline at end of file
diff --git a/code_students/include/setup/grid.hpp b/code_students/include/setup/grid.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f87df4097962f067f0e0a50a94683f888f37f58a
--- /dev/null
+++ b/code_students/include/setup/grid.hpp
@@ -0,0 +1,23 @@
+#ifndef GRID_HPP
+#define GRID_HPP
+
+#include "setup/grid1D.hpp"
+#include "util/matrix.hpp"
+
+class grid_3D {
+public:
+ grid_3D(std::vector<double> &lower_boundary, std::vector<double> &upper_boundary, std::vector<int> &num_cells, int num_ghostcells);
+ int get_num_cells(int i_dir) const;
+ grid_1D x_grid, y_grid, z_grid;
+
+private:
+ /**
+ * Compute positions of gridpoints in all directions
+ */
+ void make_axes(std::vector<double> &lower_boundary, std::vector<double> &upper_boundary);
+ matrix<double, 1> xCen, yCen, zCen;
+ std::vector<double> delta;
+ std::vector<int> num_cells;
+ int dim, num_ghostcells;
+};
+#endif
\ No newline at end of file
diff --git a/code_students/include/setup/grid1D.hpp b/code_students/include/setup/grid1D.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..12ca01f5e87303b9d5c1db403dd6b9957c08d2fa
--- /dev/null
+++ b/code_students/include/setup/grid1D.hpp
@@ -0,0 +1,86 @@
+#ifndef GRID1D_HPP
+#define GRID1D_HPP
+
+#include "util/matrix.hpp"
+
+class grid_1D {
+public:
+ /**
+ * List of available grid types
+ */
+ enum class GridType { linear, none };
+
+ grid_1D();
+
+ /**
+ * Constructor which directly produces a grid
+ * @param type type of grid selected from GridType enum
+ * @param min left boundary of physical domain
+ * @param max right boundary of physical domain
+ * @param num_cells number of cells to be distributed between min and max
+ * @param num_ghost_cells number of boundary cells outside physical domain
+ */
+ grid_1D(const GridType &type, double min, double max, size_t num_cells, int num_ghost_cells = 0);
+
+ /**
+ * make a grid arrays
+ * @param type type of grid selected from GridType enum
+ * @param min left boundary of physical domain
+ * @param max right boundary of physical domain
+ * @param num_cells number of cells to be distributed between min and max
+ * @param num_ghost_cells number of boundary cells outside physical domain
+ */
+ void make_grid(const GridType &type, double min, double max, size_t num_cells, int num_ghost_cells = 0);
+
+ /**
+ * get position of cell center / left edge
+ * @param i_cell index of cell
+ */
+ double get_center(int i_cell) const;
+ double get_left(int i_cell) const;
+
+ /**
+ * get first / last index of grid cells
+ */
+ int get_index_lowest() const;
+ int get_index_highest() const;
+
+ /**
+ * Return size of grid cell for homogeneous grid
+ */
+ double get_dx() const;
+
+ /**
+ * Return inverse of size of grid cell for homogeneous grid
+ */
+ double get_inv_dx() const;
+
+ /**
+ * Return size of grid cell for arbitrary grid
+ */
+ double get_dx(int i_cell) const;
+
+ /**
+ * Return inverse of size of grid cell for arbitrary grid
+ */
+ double get_inv_dx(int i_cell) const;
+
+protected:
+ int num_ghostcells;
+
+private:
+ void build_lin_axis(double min, double max, size_t num_cells, int num_ghost_cells = 0);
+ /**
+ * get grid centers from left-handed grid boundareis
+ */
+ void get_centers();
+
+ matrix<double, 1> grid_centers, grid_left;
+
+ double delta_x, inv_delta_x;
+
+ GridType my_type;
+ bool is_set;
+};
+
+#endif
\ No newline at end of file
diff --git a/code_students/include/setup/physics.hpp b/code_students/include/setup/physics.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..07d70209bee62e90859ab1ffb529a99a6463d5cf
--- /dev/null
+++ b/code_students/include/setup/physics.hpp
@@ -0,0 +1,24 @@
+#ifndef PHYSICS_HPP
+#define PHYSICS_HPP
+
+#include "core/config.hpp"
+#include "setup/fluid.hpp"
+
+class physics {
+public:
+ physics();
+ void get_physical_fluxes(const fluid_cell &fluid, fluxes_cell &fluxes, const parallelisation::direction &local_direction);
+ void get_lambda_min_max(double &lambda_min, double &lambda_max, const fluid_cell &fluid_left, const fluid_cell &fluid_right,
+ const parallelisation::direction &local_direction);
+ double get_sound_speed(double density, double pressure);
+ double get_lambda_abs_max(const fluid_cell &fluid);
+ double get_pressure(const fluid_cell &fluid);
+ void transform_characteristic_to_conservative(fluid_cell &fluid);
+ void transform_conservative_to_characteristic(fluid_cell &fluid);
+ double get_e_total(const fluid_cell &fluid);
+
+private:
+ double adiabatic_index;
+};
+
+#endif
\ No newline at end of file
diff --git a/code_students/include/solver/Riemann_solvers.hpp b/code_students/include/solver/Riemann_solvers.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4813be68c6493d0d7a533417d972307ce809812c
--- /dev/null
+++ b/code_students/include/solver/Riemann_solvers.hpp
@@ -0,0 +1,29 @@
+#ifndef RIEMANN_SOLVERS_HPP
+#define RIEMANN_SOLVERS_HPP
+
+#include "setup/fluid.hpp"
+
+#include <vector>
+
+class Riemann_solver {
+public:
+ Riemann_solver(std::size_t num_fields_);
+ virtual ~Riemann_solver();
+ virtual void get_num_flux(fluid_cell &fluid_left_cell, fluid_cell &fluid_right_cell, const std::vector<double> &phys_flux_left_cell,
+ const std::vector<double> &phys_flux_right_cell, std::vector<double> &num_flux, double v_char_slowest, double v_char_fastest) = 0;
+
+protected:
+ std::size_t num_fields;
+};
+
+class HLL_solver : public Riemann_solver {
+public:
+ HLL_solver(std::size_t num_fields_);
+ void get_num_flux(fluid_cell &fluid_left_cell, fluid_cell &fluid_right_cell, const std::vector<double> &phys_flux_left_cell,
+ const std::vector<double> &phys_flux_right_cell, std::vector<double> &num_flux, double v_char_slowest, double v_char_fastest);
+
+private:
+ double epsilon;
+};
+
+#endif
\ No newline at end of file
diff --git a/code_students/include/solver/finite_volume_solver.hpp b/code_students/include/solver/finite_volume_solver.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..986d4f23c8f1adffe6d9c04e245f9a4cffaabc14
--- /dev/null
+++ b/code_students/include/solver/finite_volume_solver.hpp
@@ -0,0 +1,52 @@
+#ifndef FINITE_VOLUME_SOLVER_HPP
+#define FINITE_VOLUME_SOLVER_HPP
+
+#include "core/config.hpp"
+#include "setup/fluid.hpp"
+#include "setup/grid.hpp"
+#include "setup/physics.hpp"
+#include "solver/Riemann_solvers.hpp"
+#include "solver/reconstruction.hpp"
+
+#include <functional>
+#include <memory>
+
+class finite_volume_solver {
+public:
+ finite_volume_solver(fluid ¤t_fluid);
+ double singlestep(grid_3D &spatial_grid, fluid ¤t_fluid, fluid ¤t_changes);
+ int run(grid_3D &spatial_grid, fluid ¤t_fluid, double time_final, double delta_t_output);
+ void set_init_function(std::function<void(fluid_cell &, double, double, double)>);
+ matrix<double, 3> get_data_computational_volume(grid_3D &spatial_grid, fluid ¤t_fluid, int index_data);
+ void set_verbosity(int verbosity);
+
+private:
+ double get_CFL(grid_3D &spatial_grid, fluid ¤t_fluid);
+ void set_initial_conditions(grid_3D &spatial_grid, fluid ¤t_fluid);
+ void transform_fluid_to_conservative(fluid ¤t_fluid);
+ void transform_fluid_to_characteristic(fluid ¤t_fluid);
+ void apply_boundary_conditions(grid_3D &spatial_grid, fluid ¤t_fluid);
+ void store_timestep(grid_3D &spatial_grid, fluid ¤t_fluid);
+ double compute_delta_t_next(grid_3D &spatial_grid, fluid ¤t_fluid);
+ physics fluid_physics;
+ std::unique_ptr<HLL_solver> Riemann;
+ reconsctruction_second_order reconst;
+ // std::unique_ptr<Riemann_solver> Riemann;
+ fluxes_cell num_flux_left_x, num_flux_right_x, num_flux_left_y, num_flux_right_y;
+ fluxes_cell num_flux_left_z, num_flux_right_z;
+ fluxes_cell fluxes_local, fluxes_previous;
+ fluid_cell quantities_local, quantities_previous;
+ fluid_cell quantities_temp_left, quantities_temp_right;
+ std::function<void(fluid_cell &, double, double, double)> init_function;
+ parallelisation::BoundaryType boundary_type;
+ /**
+ * Maximum allowd CFL number for scheme
+ */
+ double CFL_max;
+ double time;
+ double time_output_next, delta_t_output;
+ int num_time_steps, verbosity;
+ bool init_set, write_next_step;
+};
+
+#endif
\ No newline at end of file
diff --git a/code_students/include/solver/limiter.hpp b/code_students/include/solver/limiter.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..924d390b742eb583c96ccc0c140ad3f123c8ab2f
--- /dev/null
+++ b/code_students/include/solver/limiter.hpp
@@ -0,0 +1,21 @@
+#ifndef LIMITER_HPP
+#define LIMITER_HPP
+
+class limiter_base {
+public:
+ limiter_base();
+ // TBD by students -> use sensible parameter list
+ virtual double compute(double first, double second, double third) = 0;
+};
+
+class limiter_minmod : public limiter_base {
+public:
+ limiter_minmod(double theta);
+ // TBD by students -> use sensible parameter list
+ double compute(double first, double second, double third);
+
+private:
+ double theta;
+};
+
+#endif
\ No newline at end of file
diff --git a/code_students/include/solver/reconstruction.hpp b/code_students/include/solver/reconstruction.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..2582afc8b3b1c648cc70ec8ec4c4218ac22f43e6
--- /dev/null
+++ b/code_students/include/solver/reconstruction.hpp
@@ -0,0 +1,29 @@
+#ifndef RECONSTRUCTION_HPP
+#define RECONSTRUCTION_HPP
+
+#include "core/config.hpp"
+#include "setup/fluid.hpp"
+#include "solver/limiter.hpp"
+
+class reconstruction {
+public:
+ reconstruction();
+ virtual void compute_point_values(const fluid &fluid3D, fluid_cell &values_left, fluid_cell &values_right, const parallelisation::direction &local_direction,
+ int index_x, int index_y, int index_z) = 0;
+ void get_derivatives(const fluid &fluid, fluid_cell &derivatives, const parallelisation::direction &local_direction, int index_x, int index_y, int index_z);
+
+protected:
+ limiter_minmod limiter;
+};
+
+class reconsctruction_second_order : public reconstruction {
+public:
+ reconsctruction_second_order(const fluid &fluid3D);
+ void compute_point_values(const fluid &fluid3D, fluid_cell &values_left, fluid_cell &values_right, const parallelisation::direction &local_direction,
+ int index_x, int index_y, int index_z);
+
+private:
+ fluid_cell derivatives;
+};
+
+#endif
diff --git a/code_students/include/solver/time_integrator.hpp b/code_students/include/solver/time_integrator.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..3db04ce363095e2d7d639b709306511d1efb7e4c
--- /dev/null
+++ b/code_students/include/solver/time_integrator.hpp
@@ -0,0 +1,19 @@
+#ifndef TIME_INTEGRATOR_HPP
+#define TIME_INTEGRATOR_HPP
+
+#include "setup/fluid.hpp"
+#include "setup/grid.hpp"
+
+class RungeKutta2 {
+public:
+ RungeKutta2(const grid_3D &grid, const fluid &fluid3D);
+ void save_data(const fluid &fluid3D);
+ void do_sub_step(const grid_3D &grid, const fluid &changes, fluid &fluid3D, double delta_t, int sub_step);
+ int get_number_substeps() const;
+
+private:
+ grid_3D simulation_grid;
+ fluid fluid_storage;
+};
+
+#endif
\ No newline at end of file
diff --git a/code_students/include/util/matrix.hpp b/code_students/include/util/matrix.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..6608a0a357024bad03b676c07757e24473f12913
--- /dev/null
+++ b/code_students/include/util/matrix.hpp
@@ -0,0 +1,71 @@
+#ifndef MATRIX_HPP
+#define MATRIX_HPP
+
+#include <array>
+#include <memory>
+#include <stddef.h>
+#include <vector>
+
+template <class T, int rank> class matrix;
+
+template <class T, int dim> class matrix {
+public:
+ matrix();
+ matrix(size_t *size_dim);
+ matrix(int index_lo[dim], int index_hi[dim]);
+
+ void resize(size_t *size_dim);
+ void resize(int index_lo[dim], int index_hi[dim]);
+
+ /** index operator, writing - 1D version */
+ T &operator()(int i);
+ /** index operator, reading - 1D version */
+ T operator()(int i) const;
+ /** index operator, writing - 3D version */
+ T &operator()(int i, int j, int k);
+ /** index operator, reading - 3D version */
+ T operator()(int i, int j, int k) const;
+
+ /** Setter that accesses the global 1D array */
+ void set_raw_data(size_t index_1D, T value);
+
+ /** Getter that accesses the global 1D array */
+ T get_raw_data(size_t index_1D) const;
+
+ /**
+ * Read access to raw data
+ */
+ // const T* raw_data() const;
+ std::vector<T> get_raw_data() const;
+
+ /**
+ * get lowest / highest index in direction i_dir
+ * @param i_dir direction
+ */
+ int get_lowest(size_t i_dir) const;
+ int get_highest(size_t i_dir) const;
+
+ /**
+ * get total number of entries
+ */
+ size_t get_size() const;
+
+ /**
+ * get extent of all dimensions
+ */
+ std::vector<size_t> get_dims() const;
+
+ /**
+ * set all entries to zero
+ */
+ void clear();
+
+private:
+ void reset_size(int index_lo[dim], int index_hi[dim]);
+ std::vector<T> data;
+ // std::unique_ptr<T> data;
+ std::vector<int> index_lo, index_hi, extent;
+ size_t size;
+};
+
+#endif
\ No newline at end of file
diff --git a/code_students/include/util/utility_functions.hpp b/code_students/include/util/utility_functions.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..66504c859b91bec2096f9f38e2ad12ef75e824c7
--- /dev/null
+++ b/code_students/include/util/utility_functions.hpp
@@ -0,0 +1,8 @@
+#ifndef UTILITY_FUNCTIONS_HPP
+#define UTILITY_FUNCTIONS_HPP
+
+namespace sim_util {
+inline double square(double value) { return value * value; }
+} // namespace sim_util
+
+#endif
\ No newline at end of file
diff --git a/code_students/src/CMakeLists.txt b/code_students/src/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..0a52267591d34f4f829c226d619e08ce435a35f1
--- /dev/null
+++ b/code_students/src/CMakeLists.txt
@@ -0,0 +1,6 @@
+cmake_minimum_required(VERSION 3.10)
+
+add_subdirectory(util)
+add_subdirectory(setup)
+add_subdirectory(solver)
+add_subdirectory(IO)
diff --git a/code_students/src/IO/CMakeLists.txt b/code_students/src/IO/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..dcd288aec8992ed7a9fe6b7726451e530f7ec3bf
--- /dev/null
+++ b/code_students/src/IO/CMakeLists.txt
@@ -0,0 +1,12 @@
+add_library(IO
+ data_storage.cpp
+)
+
+target_include_directories(IO
+ PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
+ PUBLIC ${PROJECT_SOURCE_DIR}/include
+)
+
+target_link_libraries(IO
+ PUBLIC SerialHdf5
+)
\ No newline at end of file
diff --git a/code_students/src/IO/data_storage.cpp b/code_students/src/IO/data_storage.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1913d45defd511df5701bc9b8594b81f5878f8d7
--- /dev/null
+++ b/code_students/src/IO/data_storage.cpp
@@ -0,0 +1,84 @@
+#include "IO/data_storage.hpp"
+
+data_storage::data_storage(std::string file_name) {
+ this->file_name = file_name;
+ hdf5file = H5Fcreate(file_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+ hdf5group = H5Gcreate2(hdf5file, "/Data", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ H5Gclose(hdf5group);
+ H5Fclose(hdf5file);
+}
+
+void data_storage::open_file() {
+ hdf5file = H5Fopen(file_name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
+ hdf5group = H5Gopen2(hdf5file, "/Data", H5P_DEFAULT);
+}
+
+void data_storage::close_file() {
+ H5Gclose(hdf5group);
+ H5Fclose(hdf5file);
+}
+
+template <typename T, typename> void data_storage::AddGlobalAttr(const std::string &AttrName, T AttrData) {
+ // Open the file
+ open_file();
+
+ // little endian of size 1
+ hid_t datatype = get_hdf5_data_type<T>();
+ if (!std::is_same<T, std::string>::value) {
+ // This is not allowed for strings
+ H5Tset_order(datatype, H5T_ORDER_LE);
+ }
+
+ // Create dataspace
+ hid_t AttrSpace = H5Screate(H5S_SCALAR);
+
+ // Create Attribute
+ hid_t info = H5Acreate2(hdf5group, AttrName.c_str(), datatype, AttrSpace, H5P_DEFAULT, H5P_DEFAULT);
+
+ if (std::is_same<T, std::string>::value) {
+ std::string *data = (std::string *)(&AttrData); // Workaround under c++11, since we dont have if constexpr
+ const char *temp = data->c_str();
+ H5Awrite(info, datatype, &temp);
+ } else {
+ H5Awrite(info, datatype, &AttrData);
+ }
+
+ // Close Attribute
+ H5Aclose(info);
+ // Close Dataspace
+ H5Sclose(AttrSpace);
+
+ close_file();
+}
+
+template void data_storage::AddGlobalAttr(const std::string &AttrName, int AttrData);
+template void data_storage::AddGlobalAttr(const std::string &AttrName, float AttrData);
+template void data_storage::AddGlobalAttr(const std::string &AttrName, double AttrData);
+
+void data_storage::write_dataset(const std::vector<double> &raw_data, const std::vector<size_t> &extents, std::string dataset_name) {
+ // Open the file
+ open_file();
+
+ int dataset_dimension = extents.size();
+ std::vector<hsize_t> DimsData;
+ // hsize_t DimsData[dataset_dimension];
+ for (int i_dir = 0; i_dir < dataset_dimension; ++i_dir) {
+ DimsData.push_back(extents[i_dir]);
+ }
+
+ // Make dataspace:
+ hid_t dataspace = H5Screate_simple(dataset_dimension, &DimsData[0], NULL);
+
+ // Set datatype to float
+ hid_t datatype = H5Tcopy(H5T_NATIVE_DOUBLE);
+
+ // Create dataset
+ hid_t dataset = H5Dcreate2(hdf5group, dataset_name.c_str(), datatype, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+
+ H5Dwrite(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &raw_data[0]);
+
+ H5Dclose(dataset);
+ H5Sclose(dataspace);
+
+ close_file();
+}
diff --git a/code_students/src/setup/CMakeLists.txt b/code_students/src/setup/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..85043de56e5d4ca405ab46a2705dfb43cdbf671b
--- /dev/null
+++ b/code_students/src/setup/CMakeLists.txt
@@ -0,0 +1,11 @@
+add_library(sim_setup
+ grid1D.cpp
+ grid.cpp
+ fluid.cpp
+ physics.cpp
+)
+
+target_include_directories(sim_setup
+ PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
+ PUBLIC ${PROJECT_SOURCE_DIR}/include
+)
\ No newline at end of file
diff --git a/code_students/src/setup/fluid.cpp b/code_students/src/setup/fluid.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..39aafa7937f7804d63c24bf758694bd3cdb21e83
--- /dev/null
+++ b/code_students/src/setup/fluid.cpp
@@ -0,0 +1,177 @@
+#include "setup/fluid.hpp"
+#include "util/matrix.hpp"
+
+#include <cassert>
+#include <iostream>
+
+base_fluid::base_fluid() {
+ this->type = parallelisation::FluidType::none;
+ set_characteristic();
+}
+
+base_fluid::base_fluid(const parallelisation::FluidType &type) {
+ this->type = type;
+ this->setup();
+ set_characteristic();
+}
+
+void base_fluid::setup() {
+ assert(this->type != parallelisation::FluidType::none);
+
+ // Set up isothermal or adiabatic fluid:
+
+ map_fields.clear(); // Make sure that map is empty
+
+ map_fields.insert(std::pair<FieldName, int>(FieldName::density, 0));
+ map_fields.insert(std::pair<FieldName, int>(FieldName::v_x, 1));
+ map_fields.insert(std::pair<FieldName, int>(FieldName::v_y, 2));
+ map_fields.insert(std::pair<FieldName, int>(FieldName::v_z, 3));
+
+ if (type == parallelisation::FluidType::adiabatic) {
+ map_fields.insert(std::pair<FieldName, int>(FieldName::energy, 4));
+ }
+
+ map_fields.insert(std::pair<FieldName, int>(FieldName::tracer, 5));
+
+ // Set individual indices:
+ index_density = get_field_index(FieldName::density);
+ index_v_x = get_field_index(FieldName::v_x);
+ index_v_y = get_field_index(FieldName::v_y);
+ index_v_z = get_field_index(FieldName::v_z);
+
+ if (type == parallelisation::FluidType::adiabatic) {
+ index_energy = get_field_index(FieldName::energy);
+ }
+
+ index_tracer = get_field_index(FieldName::tracer);
+}
+
+int base_fluid::get_field_index(const FieldName &field) {
+ assert(field != FieldName::none);
+
+ if (map_fields.find(field) != map_fields.end()) {
+ int map_index = map_fields.at(field);
+ return map_index;
+ } else {
+ return -1;
+ }
+}
+
+parallelisation::FluidType base_fluid::get_fluid_type() const { return type; }
+
+void base_fluid::set_characteristic() { uses_characteristic_quantities = true; }
+
+void base_fluid::set_conservative() { uses_characteristic_quantities = false; }
+
+bool base_fluid::is_conservative() const { return !uses_characteristic_quantities; }
+
+int base_fluid::get_index_density() const { return index_density; }
+
+int base_fluid::get_index_v_x() const { return index_v_x; }
+
+int base_fluid::get_index_v_y() const { return index_v_y; }
+
+int base_fluid::get_index_v_z() const { return index_v_z; }
+
+int base_fluid::get_index_energy() const {
+ assert(type == parallelisation::FluidType::adiabatic);
+ return index_energy;
+}
+
+int base_fluid::get_index_tracer() const { return index_tracer; }
+
+int base_fluid::get_num_fields() const {
+ // std::cerr << " This still needs to be tested\n" << map_fields.size() << "\n";
+ return map_fields.size();
+}
+
+bool base_fluid::is_adiabatic() const {
+ if (type == parallelisation::FluidType::adiabatic) {
+ return true;
+ } else {
+ return false;
+ }
+}
+
+fluid::fluid(const parallelisation::FluidType &type) : base_fluid(type) {}
+
+void fluid::setup(grid_3D &spatial_grid) {
+
+ // Get number of fields:
+ size_t num_fields = get_num_fields();
+
+ // Make fluid arrays
+ // Append spatial grids
+ for (size_t i_field = 0; i_field < num_fields; ++i_field) {
+ fluid_data.push_back(matrix<double, 3>());
+ }
+
+ // fluid_data.push_back(matrix<double, 3>()); // density
+ // fluid_data.push_back(matrix<double, 3>()); // velocity (x)
+ // fluid_data.push_back(matrix<double, 3>()); // velocity (y)
+ // fluid_data.push_back(matrix<double, 3>()); // velocity (z)
+ // fluid_data.push_back(matrix<double, 3>()); // energy density
+ // fluid_data.push_back(matrix<double, 3>()); // Tracer
+
+ // Get extent of grid
+ int index_low[3], index_hi[3];
+ index_low[0] = spatial_grid.x_grid.get_index_lowest();
+ index_hi[0] = spatial_grid.x_grid.get_index_highest();
+
+ index_low[1] = spatial_grid.y_grid.get_index_lowest();
+ index_hi[1] = spatial_grid.y_grid.get_index_highest();
+
+ index_low[2] = spatial_grid.z_grid.get_index_lowest();
+ index_hi[2] = spatial_grid.z_grid.get_index_highest();
+
+ // Resize grids
+ for (size_t i_field = 0; i_field < fluid_data.size(); ++i_field) {
+ std::cout << " resizing field " << i_field << ":\n";
+ fluid_data[i_field].resize(index_low, index_hi);
+ std::cout << " -> " << fluid_data[i_field].get_size() << "\n";
+ }
+}
+
+void fluid::get_fluid_cell(fluid_cell &local_fluid, int index_x, int index_y, int index_z) const {
+ for (size_t i_field = 0; i_field < fluid_data.size(); ++i_field) {
+ local_fluid.fluid_data[i_field] = fluid_data[i_field](index_x, index_y, index_z);
+ }
+ if (is_conservative()) {
+ local_fluid.set_conservative();
+ } else {
+ local_fluid.set_characteristic();
+ }
+}
+
+void fluid::set_cell_values(const fluid_cell &local_fluid, int index_x, int index_y, int index_z) {
+ for (size_t i_field = 0; i_field < fluid_data.size(); ++i_field) {
+ fluid_data[i_field](index_x, index_y, index_z) = local_fluid.fluid_data[i_field];
+ }
+}
+
+void fluid::get_fluid_cell_raw(fluid_cell &local_fluid, int index_1D) const {
+ for (size_t i_field = 0; i_field < fluid_data.size(); ++i_field) {
+ local_fluid.fluid_data[i_field] = fluid_data[i_field].get_raw_data(index_1D);
+ }
+ if (is_conservative()) {
+ local_fluid.set_conservative();
+ } else {
+ local_fluid.set_characteristic();
+ }
+}
+
+void fluid::set_fluid_cell_raw(fluid_cell &local_fluid, int index_1D) {
+ for (size_t i_field = 0; i_field < fluid_data.size(); ++i_field) {
+ fluid_data[i_field].set_raw_data(index_1D, local_fluid.fluid_data[i_field]);
+ }
+}
+
+size_t fluid::get_number_fields() const { return fluid_data.size(); }
+
+fluid_cell::fluid_cell(const parallelisation::FluidType &type) : base_fluid(type) { setup(); }
+
+void fluid_cell::setup() { fluid_data.resize(get_num_fields()); }
+
+fluxes_cell::fluxes_cell(const parallelisation::FluidType &type) : base_fluid(type) { setup(); }
+
+void fluxes_cell::setup() { flux_data.resize(get_num_fields()); }
diff --git a/code_students/src/setup/grid.cpp b/code_students/src/setup/grid.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1f93d0989f31990aa40209739e246b9ce02e36ed
--- /dev/null
+++ b/code_students/src/setup/grid.cpp
@@ -0,0 +1,41 @@
+#include "setup/grid.hpp"
+
+#include <cassert>
+#include <iostream>
+
+grid_3D::grid_3D(std::vector<double> &lower_boundary, std::vector<double> &upper_boundary, std::vector<int> &num_cells, int num_ghostcells) {
+
+ dim = 3;
+ delta.resize(3);
+
+ std::cout << " Constructor\n";
+
+ this->num_cells = num_cells;
+ this->num_ghostcells = num_ghostcells;
+
+ std::cout << " axes\n";
+
+ make_axes(lower_boundary, upper_boundary);
+}
+
+int grid_3D::get_num_cells(int i_dir) const { return num_cells[i_dir]; }
+
+void grid_3D::make_axes(std::vector<double> &lower_boundary, std::vector<double> &upper_boundary) {
+ std::cout << "wtf\n";
+
+ const grid_1D::GridType type_x_grid = grid_1D::GridType::linear;
+ x_grid.make_grid(type_x_grid, lower_boundary[0], upper_boundary[0], num_cells[0], num_ghostcells);
+
+ const grid_1D::GridType type_y_grid = grid_1D::GridType::linear;
+ y_grid.make_grid(type_y_grid, lower_boundary[1], upper_boundary[1], num_cells[1], num_ghostcells);
+
+ const grid_1D::GridType type_z_grid = grid_1D::GridType::linear;
+ z_grid.make_grid(type_z_grid, lower_boundary[2], upper_boundary[2], num_cells[2], num_ghostcells);
+
+ // // Test x-axis
+ // std::cout << " X-Grid: \n";
+ // for(int i_x=-num_ghostcells; i_x<num_cells[0]+num_ghostcells; ++i_x) {
+ // std::cout << " Writing point " << i_x << "\n";
+ // std::cout << i_x << " " << x_grid.get_center(i_x) << "\n";
+ // }
+}
\ No newline at end of file
diff --git a/code_students/src/setup/grid1D.cpp b/code_students/src/setup/grid1D.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2d4699551dc060c81373ac692349377538166d7b
--- /dev/null
+++ b/code_students/src/setup/grid1D.cpp
@@ -0,0 +1,104 @@
+#include "setup/grid1D.hpp"
+
+#include <cassert>
+#include <iostream>
+
+grid_1D::grid_1D() {
+ is_set = false;
+ my_type = GridType::none;
+}
+
+grid_1D::grid_1D(const GridType &type, double min, double max, size_t num_cells, int num_ghost_cells) { make_grid(type, min, max, num_cells, num_ghost_cells); }
+
+void grid_1D::make_grid(const GridType &type, double min, double max, size_t num_cells, int num_ghost_cells) {
+
+ my_type = type;
+ if (type == GridType::linear) {
+ build_lin_axis(min, max, num_cells, num_ghost_cells);
+ delta_x = grid_left(1) - grid_left(0);
+ inv_delta_x = 1. / delta_x;
+ }
+
+ is_set = true;
+}
+
+double grid_1D::get_center(int i_cell) const {
+ assert(is_set);
+ return grid_centers(i_cell);
+}
+
+double grid_1D::get_left(int i_cell) const {
+ assert(is_set);
+ return grid_left(i_cell);
+}
+
+int grid_1D::get_index_lowest() const {
+ assert(is_set);
+ return grid_centers.get_lowest(0);
+}
+
+int grid_1D::get_index_highest() const {
+ assert(is_set);
+ return grid_centers.get_highest(0);
+}
+
+double grid_1D::get_dx() const {
+ assert(is_set);
+ assert(my_type == GridType::linear);
+ // return grid_left(1) - grid_left(0);
+ return delta_x;
+}
+
+double grid_1D::get_inv_dx() const {
+ assert(is_set);
+ assert(my_type == GridType::linear);
+ return inv_delta_x;
+}
+
+double grid_1D::get_dx(int i_cell) const {
+ assert(is_set);
+ return grid_left(i_cell + 1) - grid_left(i_cell);
+}
+
+double grid_1D::get_inv_dx(int i_cell) const {
+ assert(is_set);
+ return 1. / (grid_left(i_cell + 1) - grid_left(i_cell));
+}
+
+void grid_1D::build_lin_axis(double min, double max, size_t num_cells, int num_ghost_cells) {
+
+ assert(min < max); // Min must be smaller than max
+ assert(num_cells > 1); // size must be larger than a single cell
+ assert(num_ghost_cells >= 0); // number of ghoscells must not be smaller than zero
+
+ this->num_ghostcells = num_ghost_cells;
+
+ int grid_beg = -num_ghostcells;
+ int grid_end = num_cells + num_ghostcells;
+
+ std::cout << " Making linear grid from ";
+ std::cout << min << " " << max << " with " << num_cells << " cells\n";
+ grid_left.resize(&grid_beg, &grid_end);
+
+ const auto dx = (max - min) / static_cast<double>(num_cells);
+ for (int i = grid_beg; i <= grid_end; ++i) {
+ const auto value = min + dx * i;
+ grid_left(i) = value;
+ }
+
+ // Compute positions of cell centers
+ get_centers();
+}
+
+void grid_1D::get_centers() {
+ int cen_grid_beg = grid_left.get_lowest(0);
+ int cen_grid_end = grid_left.get_highest(0) - 1;
+ grid_centers.resize(&cen_grid_beg, &cen_grid_end);
+
+ for (int i = cen_grid_beg; i <= cen_grid_end; ++i) {
+ const auto pos_left = grid_left(i);
+ const auto pos_right = grid_left(i + 1);
+ const auto pos_center = 0.5 * (pos_left + pos_right);
+ grid_centers(i) = pos_center;
+ }
+}
\ No newline at end of file
diff --git a/code_students/src/setup/physics.cpp b/code_students/src/setup/physics.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f734c40b727a68156d49654ad4fbd33efc78fd0c
--- /dev/null
+++ b/code_students/src/setup/physics.cpp
@@ -0,0 +1,143 @@
+#include "setup/physics.hpp"
+#include "setup/fluid.hpp"
+#include "util/utility_functions.hpp"
+
+#include <cassert>
+#include <cmath>
+#include <iostream>
+
+using namespace sim_util;
+
+physics::physics() { adiabatic_index = 1.4; }
+
+double physics::get_pressure(const fluid_cell &fluid) { return (adiabatic_index - 1.0) * fluid.fluid_data[fluid.get_index_energy()]; }
+
+void physics::transform_characteristic_to_conservative(fluid_cell &fluid) {
+ assert(!fluid.is_conservative());
+ double density = fluid.fluid_data[fluid.get_index_density()];
+ double v_x = fluid.fluid_data[fluid.get_index_v_x()];
+ double v_y = fluid.fluid_data[fluid.get_index_v_y()];
+ double v_z = fluid.fluid_data[fluid.get_index_v_z()];
+ double e_thermal = fluid.fluid_data[fluid.get_index_energy()];
+
+ double e_kin = 0.5 * density * (square(v_x) + square(v_y) + square(v_z));
+ double momentum_x = density * v_x;
+ double momentum_y = density * v_y;
+ double momentum_z = density * v_z;
+ double e_total = e_thermal + e_kin;
+
+ fluid.fluid_data[fluid.get_index_v_x()] = momentum_x;
+ fluid.fluid_data[fluid.get_index_v_y()] = momentum_y;
+ fluid.fluid_data[fluid.get_index_v_z()] = momentum_z;
+ fluid.fluid_data[fluid.get_index_energy()] = e_total;
+ fluid.set_conservative();
+}
+
+void physics::transform_conservative_to_characteristic(fluid_cell &fluid) {
+ assert(fluid.is_conservative());
+ double density = fluid.fluid_data[fluid.get_index_density()];
+ double momentum_x = fluid.fluid_data[fluid.get_index_v_x()];
+ double momentum_y = fluid.fluid_data[fluid.get_index_v_y()];
+ double momentum_z = fluid.fluid_data[fluid.get_index_v_z()];
+ double e_total = fluid.fluid_data[fluid.get_index_energy()];
+
+ double v_x = momentum_x / density;
+ double v_y = momentum_y / density;
+ double v_z = momentum_z / density;
+ double e_kin = 0.5 * density * (square(v_x) + square(v_y) + square(v_z));
+ double e_thermal = e_total - e_kin;
+
+ fluid.fluid_data[fluid.get_index_v_x()] = v_x;
+ fluid.fluid_data[fluid.get_index_v_y()] = v_y;
+ fluid.fluid_data[fluid.get_index_v_z()] = v_z;
+ fluid.fluid_data[fluid.get_index_energy()] = e_thermal;
+ fluid.set_conservative();
+}
+
+double physics::get_e_total(const fluid_cell &fluid) {
+ double density = fluid.fluid_data[fluid.get_index_density()];
+ double v_x = fluid.fluid_data[fluid.get_index_v_x()];
+ double v_y = fluid.fluid_data[fluid.get_index_v_y()];
+ double v_z = fluid.fluid_data[fluid.get_index_v_z()];
+ double e_therm = fluid.fluid_data[fluid.get_index_energy()];
+
+ double e_total = 0.5 * density * (sim_util::square(v_x) + sim_util::square(v_y) + sim_util::square(v_z)) + e_therm;
+ return e_total;
+}
+
+void physics::get_physical_fluxes(const fluid_cell &fluid, fluxes_cell &fluxes, const parallelisation::direction &local_direction) {
+ // Set physical fluxes for hydrodynamics
+ assert(!fluid.is_conservative());
+ double density = fluid.fluid_data[fluid.get_index_density()];
+ double v_x = fluid.fluid_data[fluid.get_index_v_x()];
+ double v_y = fluid.fluid_data[fluid.get_index_v_y()];
+ double v_z = fluid.fluid_data[fluid.get_index_v_z()];
+ double tracer = fluid.fluid_data[fluid.get_index_tracer()];
+
+ double pressure = get_pressure(fluid);
+ double e_total = get_e_total(fluid);
+
+ // fluxes.flux_data[fluid.get_index_density()] = fluid_cell.fluid_data[]
+ if (local_direction == parallelisation::direction::x) {
+ fluxes.flux_data[fluid.get_index_density()] = density * v_x;
+ fluxes.flux_data[fluid.get_index_v_x()] = 0.0; // TBD by students
+ fluxes.flux_data[fluid.get_index_v_y()] = 0.0; // TBD by students
+ fluxes.flux_data[fluid.get_index_v_z()] = 0.0; // TBD by students
+ fluxes.flux_data[fluid.get_index_energy()] = 0.0; // TBD by students
+ fluxes.flux_data[fluid.get_index_tracer()] = tracer * v_x;
+ } else if (local_direction == parallelisation::direction::y) {
+ // TBD by students
+ } else {
+ // TBD by students
+ }
+}
+
+double physics::get_sound_speed(double density, double pressure) {
+ // TBD by students
+ return 42.0;
+}
+
+double physics::get_lambda_abs_max(const fluid_cell &fluid) {
+ assert(!fluid.is_conservative());
+ double density = fluid.fluid_data[fluid.get_index_density()];
+ double v_x = fluid.fluid_data[fluid.get_index_v_x()];
+ double v_y = fluid.fluid_data[fluid.get_index_v_y()];
+ double v_z = fluid.fluid_data[fluid.get_index_v_z()];
+ double pressure = get_pressure(fluid);
+
+ double sound_speed = get_sound_speed(density, pressure);
+ double velo_abs = sqrt(sim_util::square(v_x) + sim_util::square(v_y) + sim_util::square(v_z));
+ double lambda_abs_max = velo_abs + sound_speed;
+ // std::cout << " lambda " << density << "\n";
+ return lambda_abs_max;
+}
+
+void physics::get_lambda_min_max(double &lambda_min, double &lambda_max, const fluid_cell &fluid_left_cell, const fluid_cell &fluid_right_cell,
+ const parallelisation::direction &local_direction) {
+
+ assert(!fluid_left_cell.is_conservative());
+ assert(!fluid_right_cell.is_conservative());
+
+ int index_density = fluid_left_cell.get_index_density();
+ int index_velocity_parallel = fluid_left_cell.get_index_v_x();
+ if (local_direction == parallelisation::direction::y) {
+ // TBD by students
+ } else if (local_direction == parallelisation::direction::z) {
+ // TBD by students
+ }
+
+ double density_left = 42.0; // TBD by students
+ double density_right = 42.0; // TBD by students
+
+ double v_parallel_left = fluid_left_cell.fluid_data[index_velocity_parallel];
+ double v_parallel_right = 42.0; // TBD by students
+
+ double pressure_left = get_pressure(fluid_left_cell);
+ double pressure_right = 42.0; // TBD by students
+
+ double sound_speed_left = get_sound_speed(density_left, pressure_left);
+ double sound_speed_right = 42.0; // TBD by students
+
+ lambda_max = std::max(v_parallel_left + sound_speed_left, v_parallel_right + sound_speed_right);
+ lambda_min = 42.0; // TBD by students
+}
diff --git a/code_students/src/solver/CMakeLists.txt b/code_students/src/solver/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..22b31cf64b91451065c49d9f4d7904d514bc4dd5
--- /dev/null
+++ b/code_students/src/solver/CMakeLists.txt
@@ -0,0 +1,21 @@
+add_library(solver
+ singlestep.cpp
+ Riemann_solvers.cpp
+ Riemann_solver_HLL.cpp
+ finite_volume_solver.cpp
+ reconstruction.cpp
+ limiter.cpp
+ time_integrator.cpp
+)
+
+target_include_directories(solver
+ PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
+ PUBLIC ${PROJECT_SOURCE_DIR}/include
+)
+
+
+target_link_libraries(solver
+ PUBLIC utilities
+ PUBLIC sim_setup
+ PUBLIC IO
+)
diff --git a/code_students/src/solver/Riemann_solver_HLL.cpp b/code_students/src/solver/Riemann_solver_HLL.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..efb2ea813c9dc49cf2eefdec276ed366b65f7d62
--- /dev/null
+++ b/code_students/src/solver/Riemann_solver_HLL.cpp
@@ -0,0 +1,23 @@
+#include "solver/Riemann_solvers.hpp"
+
+HLL_solver::HLL_solver(std::size_t num_fields_) : Riemann_solver(num_fields_) { epsilon = 1.e-50; }
+
+void HLL_solver::get_num_flux(fluid_cell &fluid_left_cell, fluid_cell &fluid_right_cell, const std::vector<double> &phys_flux_left_cell,
+ const std::vector<double> &phys_flux_right_cell, std::vector<double> &num_flux, double v_char_slowest, double v_char_fastest) {
+
+ // Apply HLL fluxes
+ if (v_char_slowest > 0.0) {
+
+ for (std::size_t i_field = 0; i_field < num_fields; i_field++) {
+ num_flux[i_field] = phys_flux_left_cell[i_field];
+ }
+
+ } else if (v_char_fastest < 0.0) {
+
+ // TBD by students
+
+ } else {
+
+ // TBD by students
+ }
+}
diff --git a/code_students/src/solver/Riemann_solvers.cpp b/code_students/src/solver/Riemann_solvers.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b7db6ff5a13442402a51db7308a5aa730d224d90
--- /dev/null
+++ b/code_students/src/solver/Riemann_solvers.cpp
@@ -0,0 +1,5 @@
+#include "solver/Riemann_solvers.hpp"
+
+Riemann_solver::Riemann_solver(std::size_t num_fields_) { this->num_fields = num_fields_; }
+
+Riemann_solver::~Riemann_solver() {}
\ No newline at end of file
diff --git a/code_students/src/solver/finite_volume_solver.cpp b/code_students/src/solver/finite_volume_solver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..25a27ac4fd32c0734293df428975e2b3c626b474
--- /dev/null
+++ b/code_students/src/solver/finite_volume_solver.cpp
@@ -0,0 +1,329 @@
+#include "solver/finite_volume_solver.hpp"
+#include "IO/data_storage.hpp"
+#include "core/config.hpp"
+#include "setup/fluid.hpp"
+#include "setup/physics.hpp"
+#include "solver/time_integrator.hpp"
+
+#include <cassert>
+#include <iostream>
+#include <sstream>
+
+finite_volume_solver::finite_volume_solver(fluid ¤t_fluid)
+: reconst(current_fluid.get_fluid_type()), num_flux_left_x(current_fluid.get_fluid_type()), num_flux_right_x(current_fluid.get_fluid_type()),
+ num_flux_left_y(current_fluid.get_fluid_type()), num_flux_right_y(current_fluid.get_fluid_type()), num_flux_left_z(current_fluid.get_fluid_type()),
+ num_flux_right_z(current_fluid.get_fluid_type()), fluxes_local(current_fluid.get_fluid_type()), fluxes_previous(current_fluid.get_fluid_type()),
+ quantities_local(current_fluid.get_fluid_type()), quantities_previous(current_fluid.get_fluid_type()), quantities_temp_left(current_fluid.get_fluid_type()),
+ quantities_temp_right(current_fluid.get_fluid_type()) {
+
+ CFL_max = 0.4;
+
+ init_set = false;
+ write_next_step = false;
+ verbosity = 0;
+
+ boundary_type = parallelisation::BoundaryType::constant_extrapolation;
+
+ // Get number of fields:
+ size_t num_fields = quantities_local.get_num_fields();
+
+ // Choose a specific Riemann solver (here: HLL)
+ Riemann = std::make_unique<HLL_solver>(num_fields);
+}
+
+void finite_volume_solver::set_init_function(std::function<void(fluid_cell &, double, double, double)> init_function_user) {
+ init_function = init_function_user;
+ init_set = true;
+}
+
+void finite_volume_solver::set_verbosity(int verbosity) { this->verbosity = verbosity; }
+
+double finite_volume_solver::compute_delta_t_next(grid_3D &spatial_grid, fluid ¤t_fluid) {
+ double CFL_number = get_CFL(spatial_grid, current_fluid);
+ double delta_t = CFL_max / CFL_number;
+
+ write_next_step = false;
+
+ double delta_t_output_next = time_output_next - time;
+ if (delta_t_output_next < delta_t) {
+ delta_t = delta_t_output_next;
+ time_output_next += delta_t_output;
+ write_next_step = true;
+ }
+
+ return delta_t;
+}
+
+int finite_volume_solver::run(grid_3D &spatial_grid, fluid ¤t_fluid, double time_final, double delta_t_output) {
+ assert(init_set);
+
+ this->delta_t_output = delta_t_output;
+
+ // Set initial conditions
+ set_initial_conditions(spatial_grid, current_fluid);
+
+ // prepare
+ fluid fluid_changes(current_fluid.get_fluid_type());
+ fluid_changes.setup(spatial_grid);
+ RungeKutta2 time_stepper(spatial_grid, current_fluid);
+
+ time = 0.0;
+ num_time_steps = 0;
+
+ store_timestep(spatial_grid, current_fluid);
+
+ time_output_next = delta_t_output;
+
+ // Store initial time step
+ if (num_time_steps % 10 == 0) {
+ store_timestep(spatial_grid, current_fluid);
+ }
+
+ while (time < time_final) {
+
+ double delta_t = compute_delta_t_next(spatial_grid, current_fluid);
+
+ std::cout << " Integrating time step " << num_time_steps << " at t =" << time << "\n";
+ std::cout << " \t\t step size " << delta_t << "\n";
+
+ // Do individual Runge-Kutta steps
+ int n_RK_steps = time_stepper.get_number_substeps();
+
+ for (int i_RK_step = 0; i_RK_step < n_RK_steps; ++i_RK_step) {
+
+ singlestep(spatial_grid, current_fluid, fluid_changes);
+
+ // Transform to conservative variables
+ transform_fluid_to_conservative(current_fluid);
+ time_stepper.do_sub_step(spatial_grid, fluid_changes, current_fluid, delta_t, i_RK_step);
+ // Transform back to characteristic variables
+ transform_fluid_to_characteristic(current_fluid);
+
+ apply_boundary_conditions(spatial_grid, current_fluid);
+ }
+
+ num_time_steps += 1;
+ time += delta_t;
+
+ // Write some test outputs
+ // if(num_time_steps%10 == 0) {
+ if (write_next_step) {
+ store_timestep(spatial_grid, current_fluid);
+ }
+ }
+
+ return num_time_steps;
+}
+
+void finite_volume_solver::set_initial_conditions(grid_3D &spatial_grid, fluid ¤t_fluid) {
+ assert(init_set);
+
+ for (int ix = 0; ix < spatial_grid.get_num_cells(0); ++ix) {
+ double x_position = spatial_grid.x_grid.get_center(ix);
+ for (int iy = 0; iy < spatial_grid.get_num_cells(1); ++iy) {
+ double y_position = spatial_grid.y_grid.get_center(iy);
+ for (int iz = 0; iz < spatial_grid.get_num_cells(2); ++iz) {
+ double z_position = spatial_grid.z_grid.get_center(iz);
+
+ init_function(quantities_local, x_position, y_position, z_position);
+ current_fluid.set_cell_values(quantities_local, ix, iy, iz);
+ }
+ }
+ }
+
+ apply_boundary_conditions(spatial_grid, current_fluid);
+}
+
+void finite_volume_solver::apply_boundary_conditions(grid_3D &spatial_grid, fluid ¤t_fluid) {
+
+ if (boundary_type == parallelisation::BoundaryType::constant_extrapolation) {
+
+ for (size_t i_field = 0; i_field < current_fluid.fluid_data.size(); ++i_field) {
+
+ // Lower x boundary
+ for (int ix = -1; ix >= -2; --ix) {
+ for (int iy = 0; iy < spatial_grid.get_num_cells(1); ++iy) {
+ for (int iz = 0; iz < spatial_grid.get_num_cells(2); ++iz) {
+ current_fluid.fluid_data[i_field](ix, iy, iz) = current_fluid.fluid_data[i_field](ix + 1, iy, iz);
+ }
+ }
+ }
+ // Upper x boundary
+ for (int ix = spatial_grid.get_num_cells(0); ix < spatial_grid.get_num_cells(0) + 2; ++ix) {
+ for (int iy = 0; iy < spatial_grid.get_num_cells(1); ++iy) {
+ for (int iz = 0; iz < spatial_grid.get_num_cells(2); ++iz) {
+ current_fluid.fluid_data[i_field](ix, iy, iz) = current_fluid.fluid_data[i_field](ix - 1, iy, iz);
+ }
+ }
+ }
+
+ // Lower y boundary
+ for (int ix = 0; ix < spatial_grid.get_num_cells(0); ++ix) {
+ for (int iy = -1; iy >= -2; --iy) {
+ for (int iz = 0; iz < spatial_grid.get_num_cells(2); ++iz) {
+ current_fluid.fluid_data[i_field](ix, iy, iz) = current_fluid.fluid_data[i_field](ix, iy + 1, iz);
+ }
+ }
+ }
+ // Upper y boundary
+ for (int ix = 0; ix < spatial_grid.get_num_cells(0); ++ix) {
+ for (int iy = spatial_grid.get_num_cells(1); iy < spatial_grid.get_num_cells(1) + 2; ++iy) {
+ for (int iz = 0; iz < spatial_grid.get_num_cells(2); ++iz) {
+ current_fluid.fluid_data[i_field](ix, iy, iz) = current_fluid.fluid_data[i_field](ix, iy - 1, iz);
+ }
+ }
+ }
+
+ // Lower z boundary
+ for (int ix = 0; ix < spatial_grid.get_num_cells(0); ++ix) {
+ for (int iy = 0; iy < spatial_grid.get_num_cells(1); ++iy) {
+ for (int iz = -1; iz >= -2; --iz) {
+ current_fluid.fluid_data[i_field](ix, iy, iz) = current_fluid.fluid_data[i_field](ix, iy, iz + 1);
+ }
+ }
+ }
+ // Upper z boundary
+ for (int ix = 0; ix < spatial_grid.get_num_cells(0); ++ix) {
+ for (int iy = 0; iy < spatial_grid.get_num_cells(1); ++iy) {
+ for (int iz = spatial_grid.get_num_cells(2); iz < spatial_grid.get_num_cells(2) + 2; ++iz) {
+ current_fluid.fluid_data[i_field](ix, iy, iz) = current_fluid.fluid_data[i_field](ix, iy, iz - 1);
+ }
+ }
+ }
+ }
+ }
+}
+
+void finite_volume_solver::transform_fluid_to_conservative(fluid ¤t_fluid) {
+ size_t num_cells = current_fluid.fluid_data[0].get_size();
+
+ for (size_t i_cell = 0; i_cell < num_cells; ++i_cell) {
+ current_fluid.get_fluid_cell_raw(quantities_local, i_cell);
+ fluid_physics.transform_characteristic_to_conservative(quantities_local);
+ current_fluid.set_fluid_cell_raw(quantities_local, i_cell);
+ }
+ current_fluid.set_conservative();
+}
+
+void finite_volume_solver::transform_fluid_to_characteristic(fluid ¤t_fluid) {
+ size_t num_cells = current_fluid.fluid_data[0].get_size();
+
+ for (size_t i_cell = 0; i_cell < num_cells; ++i_cell) {
+ current_fluid.get_fluid_cell_raw(quantities_local, i_cell);
+ fluid_physics.transform_conservative_to_characteristic(quantities_local);
+ current_fluid.set_fluid_cell_raw(quantities_local, i_cell);
+ }
+ current_fluid.set_characteristic();
+}
+
+double finite_volume_solver::get_CFL(grid_3D &spatial_grid, fluid ¤t_fluid) {
+ double CLF_number = 0.0;
+ double delta_min = std::min(spatial_grid.x_grid.get_dx(), std::min(spatial_grid.y_grid.get_dx(), spatial_grid.z_grid.get_dx()));
+ std::cout << " delta: " << delta_min << "\n";
+ // Loop over full grid and compute local CFL value
+ for (int ix = 0; ix < spatial_grid.get_num_cells(0); ++ix) {
+ for (int iy = 0; iy < spatial_grid.get_num_cells(1); ++iy) {
+ for (int iz = 0; iz < spatial_grid.get_num_cells(2); ++iz) {
+ current_fluid.get_fluid_cell(quantities_local, ix, iy, iz);
+ double v_max = fluid_physics.get_lambda_abs_max(quantities_local);
+ CLF_number = std::max(CLF_number, v_max / delta_min);
+ }
+ }
+ }
+ return CLF_number;
+}
+
+void finite_volume_solver::store_timestep(grid_3D &spatial_grid, fluid ¤t_fluid) {
+ std::ostringstream oss_file_name;
+ oss_file_name << "output_step" << num_time_steps << ".h5";
+ std::string file_name = oss_file_name.str();
+
+ data_storage storage(file_name);
+ storage.AddGlobalAttr<int>("Output step", num_time_steps);
+ storage.AddGlobalAttr<int>("Output time", time);
+
+ // Store the computational grid
+ std::vector<double> grid_positions_center, grid_positions_left;
+ std::vector<size_t> extent_grid;
+
+ // x direction
+ for (int index_x = 0; index_x < spatial_grid.get_num_cells(0); ++index_x) {
+ grid_positions_center.push_back(spatial_grid.x_grid.get_center(index_x));
+ grid_positions_left.push_back(spatial_grid.x_grid.get_left(index_x));
+ }
+ grid_positions_left.push_back(spatial_grid.x_grid.get_left(spatial_grid.get_num_cells(0)));
+ extent_grid.push_back(spatial_grid.get_num_cells(0));
+ storage.write_dataset(grid_positions_center, extent_grid, "x_grid_centers");
+ extent_grid[0] += 1;
+ storage.write_dataset(grid_positions_left, extent_grid, "x_grid_left");
+
+ // y direction
+ grid_positions_center.clear();
+ grid_positions_left.clear();
+ for (int index_y = 0; index_y < spatial_grid.get_num_cells(1); ++index_y) {
+ grid_positions_center.push_back(spatial_grid.y_grid.get_center(index_y));
+ grid_positions_left.push_back(spatial_grid.y_grid.get_left(index_y));
+ }
+ grid_positions_left.push_back(spatial_grid.y_grid.get_left(spatial_grid.get_num_cells(2)));
+ extent_grid[0] = spatial_grid.get_num_cells(1);
+ storage.write_dataset(grid_positions_center, extent_grid, "y_grid_centers");
+ extent_grid[0] += 1;
+ storage.write_dataset(grid_positions_left, extent_grid, "y_grid_left");
+
+ // z direction
+ grid_positions_center.clear();
+ grid_positions_left.clear();
+ for (int index_z = 0; index_z < spatial_grid.get_num_cells(2); ++index_z) {
+ grid_positions_center.push_back(spatial_grid.z_grid.get_center(index_z));
+ grid_positions_left.push_back(spatial_grid.z_grid.get_left(index_z));
+ }
+ grid_positions_left.push_back(spatial_grid.z_grid.get_left(spatial_grid.get_num_cells(2)));
+ extent_grid[0] = spatial_grid.get_num_cells(2);
+ storage.write_dataset(grid_positions_center, extent_grid, "z_grid_centers");
+ extent_grid[0] += 1;
+ storage.write_dataset(grid_positions_left, extent_grid, "z_grid_left");
+
+ std::vector<size_t> extents;
+ // extents.push_back(current_fluid.fluid_data[0].get_highest(0)-current_fluid.fluid_data[0].get_lowest(0)+1);
+ // extents.push_back(current_fluid.fluid_data[0].get_highest(1)-current_fluid.fluid_data[0].get_lowest(1)+1);
+ // extents.push_back(current_fluid.fluid_data[0].get_highest(2)-current_fluid.fluid_data[0].get_lowest(2)+1);
+
+ extents.push_back(spatial_grid.get_num_cells(0));
+ extents.push_back(spatial_grid.get_num_cells(1));
+ extents.push_back(spatial_grid.get_num_cells(2));
+
+ // Store the actual data
+ matrix<double, 3> storage_data = get_data_computational_volume(spatial_grid, current_fluid, current_fluid.get_index_density());
+ storage.write_dataset(storage_data.get_raw_data(), extents, "Density");
+ storage_data = get_data_computational_volume(spatial_grid, current_fluid, current_fluid.get_index_v_x());
+ storage.write_dataset(storage_data.get_raw_data(), extents, "v_x");
+ storage_data = get_data_computational_volume(spatial_grid, current_fluid, current_fluid.get_index_v_y());
+ storage.write_dataset(storage_data.get_raw_data(), extents, "v_y");
+ storage_data = get_data_computational_volume(spatial_grid, current_fluid, current_fluid.get_index_v_z());
+ storage.write_dataset(storage_data.get_raw_data(), extents, "v_z");
+ storage_data = get_data_computational_volume(spatial_grid, current_fluid, current_fluid.get_index_energy());
+ storage.write_dataset(storage_data.get_raw_data(), extents, "Energy");
+ storage_data = get_data_computational_volume(spatial_grid, current_fluid, current_fluid.get_index_tracer());
+ storage.write_dataset(storage_data.get_raw_data(), extents, "Tracer");
+}
+
+matrix<double, 3> finite_volume_solver::get_data_computational_volume(grid_3D &spatial_grid, fluid ¤t_fluid, int index_data) {
+ matrix<double, 3> return_data;
+ int index_low[3], index_hi[3];
+ index_low[0] = 0;
+ index_low[1] = 0;
+ index_low[2] = 0;
+ index_hi[0] = spatial_grid.get_num_cells(0) - 1;
+ index_hi[1] = spatial_grid.get_num_cells(1) - 1;
+ index_hi[2] = spatial_grid.get_num_cells(2) - 1;
+ return_data.resize(index_low, index_hi);
+
+ for (int ix = 0; ix < spatial_grid.get_num_cells(0); ++ix) {
+ for (int iy = 0; iy < spatial_grid.get_num_cells(1); ++iy) {
+ for (int iz = 0; iz < spatial_grid.get_num_cells(2); ++iz) {
+ return_data(ix, iy, iz) = current_fluid.fluid_data[index_data](ix, iy, iz);
+ }
+ }
+ }
+ return return_data;
+}
\ No newline at end of file
diff --git a/code_students/src/solver/limiter.cpp b/code_students/src/solver/limiter.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1644b64e3eb707898d924fab262610aa0fa33171
--- /dev/null
+++ b/code_students/src/solver/limiter.cpp
@@ -0,0 +1,14 @@
+#include "solver/limiter.hpp"
+
+#include <algorithm>
+#include <iostream>
+
+limiter_base::limiter_base() {}
+
+limiter_minmod::limiter_minmod(double theta) { this->theta = theta; }
+
+double limiter_minmod::compute(double first, double second, double third) {
+
+ // TBD by students
+ return 42.0;
+}
\ No newline at end of file
diff --git a/code_students/src/solver/reconstruction.cpp b/code_students/src/solver/reconstruction.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9dcbe27eb73361b85406a691b41e3b44e6a54f08
--- /dev/null
+++ b/code_students/src/solver/reconstruction.cpp
@@ -0,0 +1,48 @@
+#include "solver/reconstruction.hpp"
+#include "core/config.hpp"
+#include "solver/limiter.hpp"
+
+reconstruction::reconstruction() : limiter(1.0) {}
+
+void reconstruction::get_derivatives(const fluid &fluid3D, fluid_cell &derivatives, const parallelisation::direction &local_direction, int index_x, int index_y,
+ int index_z) {
+
+ size_t num_fields = fluid3D.get_number_fields();
+
+ for (size_t index_field = 0; index_field < num_fields; ++index_field) {
+ double delta_left(0.0), delta_right(0.0);
+ if (local_direction == parallelisation::direction::x) {
+ delta_left = fluid3D.fluid_data[index_field](index_x, index_y, index_z) - fluid3D.fluid_data[index_field](index_x - 1, index_y, index_z);
+ delta_right = fluid3D.fluid_data[index_field](index_x + 1, index_y, index_z) - fluid3D.fluid_data[index_field](index_x, index_y, index_z);
+ } else if (local_direction == parallelisation::direction::y) {
+ delta_left = fluid3D.fluid_data[index_field](index_x, index_y, index_z) - fluid3D.fluid_data[index_field](index_x, index_y - 1, index_z);
+ delta_right = fluid3D.fluid_data[index_field](index_x, index_y + 1, index_z) - fluid3D.fluid_data[index_field](index_x, index_y, index_z);
+ } else {
+ delta_left = fluid3D.fluid_data[index_field](index_x, index_y, index_z) - fluid3D.fluid_data[index_field](index_x, index_y, index_z - 1);
+ delta_right = fluid3D.fluid_data[index_field](index_x, index_y, index_z + 1) - fluid3D.fluid_data[index_field](index_x, index_y, index_z);
+ }
+ double delta_central = delta_left + delta_right;
+
+ // Now, apply the limiter
+ // TBD by students: here you have to pass sensible parameters
+ double derivative = limiter.compute(42.0, 42.0, 42.0);
+ derivatives.fluid_data[index_field] = derivative;
+ }
+}
+
+reconsctruction_second_order::reconsctruction_second_order(const fluid &fluid3D) : derivatives(fluid3D.get_fluid_type()) {}
+
+void reconsctruction_second_order::compute_point_values(const fluid &fluid3D, fluid_cell &values_left, fluid_cell &values_right,
+ const parallelisation::direction &local_direction, int index_x, int index_y, int index_z) {
+
+ // Start by computing the derivatives:
+ get_derivatives(fluid3D, derivatives, local_direction, index_x, index_y, index_z);
+
+ // Now, compute point values
+ for (size_t index_field = 0; index_field < fluid3D.get_number_fields(); ++index_field) {
+ double value_local = fluid3D.fluid_data[index_field](index_x, index_y, index_z);
+
+ values_left.fluid_data[index_field] = value_local - 0.5 * derivatives.fluid_data[index_field];
+ values_right.fluid_data[index_field] = value_local + 0.5 * derivatives.fluid_data[index_field];
+ }
+}
\ No newline at end of file
diff --git a/code_students/src/solver/singlestep.cpp b/code_students/src/solver/singlestep.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4390feddf9d84f4f33604065a1c389223edda10c
--- /dev/null
+++ b/code_students/src/solver/singlestep.cpp
@@ -0,0 +1,166 @@
+#include "core/config.hpp"
+#include "setup/physics.hpp"
+#include "solver/finite_volume_solver.hpp"
+
+#include <iostream>
+
+double finite_volume_solver::singlestep(grid_3D &spatial_grid, fluid ¤t_fluid, fluid ¤t_changes) {
+
+ size_t num_fields = current_fluid.get_number_fields();
+
+ for (size_t i_field = 0; i_field < num_fields; ++i_field) {
+ current_changes.fluid_data[i_field].clear();
+ }
+
+ for (int ix = 0; ix <= spatial_grid.get_num_cells(0); ++ix) {
+ for (int iy = 0; iy <= spatial_grid.get_num_cells(1); ++iy) {
+ for (int iz = 0; iz <= spatial_grid.get_num_cells(2); ++iz) {
+
+ // x-dimension
+
+ // Get local values using a reconstruction
+ // previous cell
+ reconst.compute_point_values(current_fluid, quantities_temp_left, quantities_temp_right, parallelisation::direction::x, ix - 1, iy, iz);
+ quantities_previous = quantities_temp_right;
+ // local cell
+ reconst.compute_point_values(current_fluid, quantities_temp_left, quantities_temp_right, parallelisation::direction::x, ix, iy, iz);
+ quantities_local = quantities_temp_left;
+
+ // compute physical fluxes
+ fluid_physics.get_physical_fluxes(quantities_previous, fluxes_previous, parallelisation::direction::x);
+ fluid_physics.get_physical_fluxes(quantities_local, fluxes_local, parallelisation::direction::x);
+
+ // Compute characteristic velocities from point values
+ double lambda_min_x(0.0), lambda_max_x(0.0);
+ fluid_physics.get_lambda_min_max(lambda_min_x, lambda_max_x, quantities_previous, quantities_local, parallelisation::direction::x);
+
+ // Need conservative variables for numerical fluxes
+ fluid_physics.transform_characteristic_to_conservative(quantities_previous);
+ fluid_physics.transform_characteristic_to_conservative(quantities_local);
+
+ // Compute numerical fluxes from point values and physical fluxes
+ // Compute left-handed flux in x-direction
+ Riemann->get_num_flux(quantities_previous, quantities_local, fluxes_previous.flux_data, fluxes_local.flux_data, num_flux_left_x.flux_data, lambda_min_x,
+ lambda_max_x);
+
+ // Compute changes of cell - left handed flux applies to previous and local cell
+ for (size_t i_field = 0; i_field < num_fields; ++i_field) {
+ if (ix > 0) {
+ current_changes.fluid_data[i_field](ix - 1, iy, iz) -= num_flux_left_x.flux_data[i_field] * spatial_grid.x_grid.get_inv_dx();
+ }
+ if (ix < spatial_grid.get_num_cells(0)) {
+ current_changes.fluid_data[i_field](ix, iy, iz) += num_flux_left_x.flux_data[i_field] * spatial_grid.x_grid.get_inv_dx();
+ }
+ }
+
+ // y-dimension
+
+ // Get local values using a reconstruction
+ // previous cell
+ reconst.compute_point_values(current_fluid, quantities_temp_left, quantities_temp_right, parallelisation::direction::y, ix, iy - 1, iz);
+ quantities_previous = quantities_temp_right;
+ // local cell
+ reconst.compute_point_values(current_fluid, quantities_temp_left, quantities_temp_right, parallelisation::direction::y, ix, iy, iz);
+ quantities_local = quantities_temp_left;
+
+ // compute physical fluxes
+ fluid_physics.get_physical_fluxes(quantities_previous, fluxes_previous, parallelisation::direction::y);
+ fluid_physics.get_physical_fluxes(quantities_local, fluxes_local, parallelisation::direction::y);
+
+ // Compute characteristic velocities from point values
+ double lambda_min_y(0.0), lambda_max_y(0.0);
+ fluid_physics.get_lambda_min_max(lambda_min_y, lambda_max_y, quantities_previous, quantities_local, parallelisation::direction::y);
+
+ // Need conservative variables for numerical fluxes
+ fluid_physics.transform_characteristic_to_conservative(quantities_previous);
+ fluid_physics.transform_characteristic_to_conservative(quantities_local);
+
+ // Compute numerical fluxes from point values and physical fluxes
+ // Compute left-handed flux in x-direction
+ Riemann->get_num_flux(quantities_previous, quantities_local, fluxes_previous.flux_data, fluxes_local.flux_data, num_flux_left_y.flux_data, lambda_min_y,
+ lambda_max_y);
+
+ // Compute changes of cell
+ for (size_t i_field = 0; i_field < num_fields; ++i_field) {
+ if (iy > 0) {
+ current_changes.fluid_data[i_field](ix, iy - 1, iz) -= num_flux_left_y.flux_data[i_field] * spatial_grid.y_grid.get_inv_dx();
+ }
+ if (iy < spatial_grid.get_num_cells(1)) {
+ current_changes.fluid_data[i_field](ix, iy, iz) += num_flux_left_y.flux_data[i_field] * spatial_grid.y_grid.get_inv_dx();
+ }
+ }
+
+ // z-dimension
+
+ // Get local values using a reconstruction
+ // previous cell
+ reconst.compute_point_values(current_fluid, quantities_temp_left, quantities_temp_right, parallelisation::direction::z, ix, iy, iz - 1);
+ quantities_previous = quantities_temp_right;
+ // local cell
+ reconst.compute_point_values(current_fluid, quantities_temp_left, quantities_temp_right, parallelisation::direction::z, ix, iy, iz);
+ quantities_local = quantities_temp_left;
+
+ // compute physical fluxes
+ fluid_physics.get_physical_fluxes(quantities_previous, fluxes_previous, parallelisation::direction::z);
+ fluid_physics.get_physical_fluxes(quantities_local, fluxes_local, parallelisation::direction::z);
+
+ // Compute characteristic velocities from point values
+ double lambda_min_z(0.0), lambda_max_z(0.0);
+ fluid_physics.get_lambda_min_max(lambda_min_z, lambda_max_z, quantities_previous, quantities_local, parallelisation::direction::z);
+
+ // Need conservative variables for numerical fluxes
+ fluid_physics.transform_characteristic_to_conservative(quantities_previous);
+ fluid_physics.transform_characteristic_to_conservative(quantities_local);
+
+ // Compute numerical fluxes from point values and physical fluxes
+ // Compute left-handed flux in z-direction
+ Riemann->get_num_flux(quantities_previous, quantities_local, fluxes_previous.flux_data, fluxes_local.flux_data, num_flux_left_z.flux_data, lambda_min_z,
+ lambda_max_z);
+
+ // Compute changes of cell
+ for (size_t i_field = 0; i_field < num_fields; ++i_field) {
+ if (iz > 0) {
+ current_changes.fluid_data[i_field](ix, iy, iz - 1) -= num_flux_left_z.flux_data[i_field] * spatial_grid.z_grid.get_inv_dx();
+ }
+ if (iz < spatial_grid.get_num_cells(2)) {
+ current_changes.fluid_data[i_field](ix, iy, iz) += num_flux_left_z.flux_data[i_field] * spatial_grid.z_grid.get_inv_dx();
+ }
+ }
+ }
+ }
+ }
+
+ if (verbosity > 10) {
+ std::cout << " Values: ";
+ std::cout << current_fluid.fluid_data[0](14, 16, 16) << " ";
+ std::cout << current_fluid.fluid_data[4](13, 16, 16) << " ";
+ std::cout << current_fluid.fluid_data[4](14, 16, 16) << " ";
+ std::cout << current_fluid.fluid_data[4](15, 16, 16) << " ";
+ std::cout << "\n";
+ std::cout << " Changes ";
+ std::cout << current_changes.fluid_data[0](13, 16, 16) << " ";
+ std::cout << current_changes.fluid_data[0](14, 16, 16) << " ";
+ std::cout << current_changes.fluid_data[0](15, 16, 16) << " ";
+ std::cout << "\n";
+ std::cout << current_changes.fluid_data[1](12, 16, 16) << " ";
+ std::cout << current_changes.fluid_data[1](13, 16, 16) << " ";
+ std::cout << current_changes.fluid_data[1](14, 16, 16) << " ";
+ std::cout << current_changes.fluid_data[1](15, 16, 16) << " ";
+ std::cout << "\n";
+ std::cout << current_changes.fluid_data[2](13, 16, 16) << " ";
+ std::cout << current_changes.fluid_data[2](14, 16, 16) << " ";
+ std::cout << current_changes.fluid_data[2](15, 16, 16) << " ";
+ std::cout << "\n";
+ std::cout << current_changes.fluid_data[3](13, 16, 16) << " ";
+ std::cout << current_changes.fluid_data[3](14, 16, 16) << " ";
+ std::cout << current_changes.fluid_data[3](15, 16, 16) << " ";
+ std::cout << "\n";
+ std::cout << current_changes.fluid_data[4](12, 16, 16) << " ";
+ std::cout << current_changes.fluid_data[4](13, 16, 16) << " ";
+ std::cout << current_changes.fluid_data[4](14, 16, 16) << " ";
+ std::cout << current_changes.fluid_data[4](15, 16, 16) << " ";
+ std::cout << "\n";
+ }
+
+ return 0.0;
+}
\ No newline at end of file
diff --git a/code_students/src/solver/time_integrator.cpp b/code_students/src/solver/time_integrator.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fa8f7e02b57ff129c5e4e82366c7049cbeb040ab
--- /dev/null
+++ b/code_students/src/solver/time_integrator.cpp
@@ -0,0 +1,41 @@
+#include "solver/time_integrator.hpp"
+#include "setup/grid.hpp"
+
+#include <iostream>
+
+RungeKutta2::RungeKutta2(const grid_3D &grid, const fluid &fluid3D) : simulation_grid(grid), fluid_storage(fluid3D.get_fluid_type()) {
+ fluid_storage.setup(simulation_grid);
+}
+
+void RungeKutta2::save_data(const fluid &fluid3D) {
+ for (int i_field = 0; i_field < fluid3D.get_num_fields(); ++i_field) {
+ fluid_storage.fluid_data[i_field] = fluid3D.fluid_data[i_field];
+ }
+}
+
+void RungeKutta2::do_sub_step(const grid_3D &grid, const fluid &changes, fluid &fluid3D, double delta_t, int sub_step) {
+
+ int num_fields = fluid3D.get_num_fields();
+
+ std::cout << " Doing Runge-Kutta step " << sub_step << "\n";
+ if (sub_step == 0) {
+ save_data(fluid3D);
+ }
+
+ for (int i_field = 0; i_field < num_fields; ++i_field) {
+ if (sub_step == 0) {
+
+ // TBD by students -> do first Runge-Kutta step
+
+ } else if (sub_step == 1) {
+
+ // TBD by students -> do second Runge-Kutta step
+
+ } else {
+ std::cerr << " Error no such substep: " << sub_step << "\n";
+ exit(3);
+ }
+ }
+}
+
+int RungeKutta2::get_number_substeps() const { return 2; }
\ No newline at end of file
diff --git a/code_students/src/util/CMakeLists.txt b/code_students/src/util/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a0e85e5e266cee7201d73bafd123cd2d57bc9d0a
--- /dev/null
+++ b/code_students/src/util/CMakeLists.txt
@@ -0,0 +1,8 @@
+add_library(utilities
+ matrix.cpp
+)
+
+target_include_directories(utilities
+ PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
+ PUBLIC ${PROJECT_SOURCE_DIR}/include
+)
\ No newline at end of file
diff --git a/code_students/src/util/matrix.cpp b/code_students/src/util/matrix.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9ad73df44f00fbcc110260e78dd7ab83537c6fe0
--- /dev/null
+++ b/code_students/src/util/matrix.cpp
@@ -0,0 +1,139 @@
+#include "util/matrix.hpp"
+
+#include <cassert>
+#include <iostream>
+
+template <class T, int dim> matrix<T, dim>::matrix() {
+ data.resize(0);
+ extent.resize(dim);
+ index_lo.resize(dim);
+ index_hi.resize(dim);
+ size = 0;
+ for (size_t i_dir = 0; i_dir < dim; ++i_dir) {
+ extent[i_dir] = 0;
+ index_lo[i_dir] = 0;
+ index_hi[i_dir] = 0;
+ }
+}
+
+template <class T, int dim> matrix<T, dim>::matrix(size_t *size_dim) {
+
+ int ind_lo[dim] = {0};
+ int ind_hi[dim];
+
+ for (size_t i_dir = 0; i_dir < dim; ++i_dir) {
+ ind_hi[i_dir] = size_dim[i_dir];
+ }
+ reset_size(ind_lo, ind_hi);
+}
+
+template <class T, int dim> matrix<T, dim>::matrix(int index_lo[dim], int index_hi[dim]) { reset_size(index_lo, index_hi); }
+
+template <class T, int dim> void matrix<T, dim>::resize(size_t *size_dim) {
+ int ind_lo[dim] = {0};
+ int ind_hi[dim];
+
+ for (size_t i_dir = 0; i_dir < dim; ++i_dir) {
+ ind_hi[i_dir] = size_dim[i_dir];
+ }
+ reset_size(ind_lo, ind_hi);
+}
+
+template <class T, int dim> void matrix<T, dim>::resize(int _index_lo[dim], int index_hi[dim]) { reset_size(_index_lo, index_hi); }
+
+template <class T, int dim> void matrix<T, dim>::reset_size(int index_lo[dim], int index_hi[dim]) {
+ size = 1;
+ extent.resize(dim);
+ this->index_lo.resize(dim);
+ this->index_hi.resize(dim);
+ for (size_t i_dir = 0; i_dir < dim; ++i_dir) {
+ this->index_lo[i_dir] = index_lo[i_dir];
+ this->index_hi[i_dir] = index_hi[i_dir];
+ extent[i_dir] = index_hi[i_dir] - index_lo[i_dir] + 1;
+ size *= extent[i_dir];
+ }
+
+ data.resize(size);
+}
+
+template <class T, int dim> T &matrix<T, dim>::operator()(int i) {
+ assert(dim == 1);
+ assert(i >= index_lo[0] && i <= index_hi[0]);
+ return data[i - index_lo[0]];
+}
+
+template <class T, int dim> T matrix<T, dim>::operator()(int i) const {
+ assert(dim == 1);
+ assert(i >= index_lo[0] && i <= index_hi[0]);
+ return data[i - index_lo[0]];
+}
+
+template <class T, int dim> T &matrix<T, dim>::operator()(int i, int j, int k) {
+ assert(dim == 3);
+ assert(i >= index_lo[0] && i <= index_hi[0]);
+ assert(j >= index_lo[1] && j <= index_hi[1]);
+ assert(k >= index_lo[2] && k <= index_hi[2]);
+ // TBD by students - check if this is the best implementation
+ return data[((k - index_lo[2]) * extent[1] + (j - index_lo[1])) * extent[0] + (i - index_lo[0])];
+}
+
+template <class T, int dim> T matrix<T, dim>::operator()(int i, int j, int k) const {
+ assert(dim == 3);
+ assert(i >= index_lo[0] && i <= index_hi[0]);
+ assert(j >= index_lo[1] && j <= index_hi[1]);
+ assert(k >= index_lo[2] && k <= index_hi[2]);
+ // TBD by students - check if this is the best implementation
+ return data[((k - index_lo[2]) * extent[1] + (j - index_lo[1])) * extent[0] + (i - index_lo[0])];
+}
+
+template <class T, int dim> void matrix<T, dim>::set_raw_data(size_t index_1D, T value) {
+ assert(index_1D < size);
+ data[index_1D] = value;
+}
+
+template <class T, int dim> T matrix<T, dim>::get_raw_data(size_t index_1D) const {
+ assert(index_1D < size);
+ return data[index_1D];
+}
+
+template <class T, int dim> int matrix<T, dim>::get_lowest(size_t i_dir) const {
+ assert(i_dir < dim); // index of direction must be within bounds of dim
+ return index_lo[i_dir];
+}
+
+template <class T, int dim> int matrix<T, dim>::get_highest(size_t i_dir) const {
+ assert(i_dir < dim); // index of direction must be within bounds of dim
+ return index_hi[i_dir];
+}
+
+template <class T, int dim> size_t matrix<T, dim>::get_size() const { return size; }
+
+template <class T, int dim> std::vector<size_t> matrix<T, dim>::get_dims() const {
+ std::vector<size_t> dim_extents;
+ for (size_t i_dir = 0; i_dir < dim; ++i_dir) {
+ dim_extents.push_back(extent[i_dir]);
+ }
+ return dim_extents;
+}
+
+template <class T, int dim> void matrix<T, dim>::clear() {
+ for (size_t i_entry = 0; i_entry < size; ++i_entry) {
+ data[i_entry] = 0;
+ }
+}
+
+// template<class T, int dim>
+// const T* matrix<T,dim>::raw_data() const {
+// return &data[0];
+// }
+
+template <class T, int dim> std::vector<T> matrix<T, dim>::get_raw_data() const {
+ std::vector<T> raw_data;
+ raw_data = data;
+ return raw_data;
+}
+
+// // Make instances of matrix:
+template class matrix<double, 1>;
+template class matrix<double, 2>;
+template class matrix<double, 3>;
\ No newline at end of file
diff --git a/code_students/tests/CMakeLists.txt b/code_students/tests/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..7ee6cb15dd85771ab4aeac1ca0c0fa5a4cbd0292
--- /dev/null
+++ b/code_students/tests/CMakeLists.txt
@@ -0,0 +1,22 @@
+set(EXECUTABLE_OUTPUT_PATH ${PROJECT_SOURCE_DIR}/bin/tests )
+MESSAGE(STATUS "Installing tests in ${EXECUTABLE_OUTPUT_PATH}")
+
+add_executable(test_limiter test_limiter.cpp)
+target_include_directories(test_limiter PRIVATE ${PROJECT_SOURCE_DIR}/include)
+target_link_libraries(test_limiter PRIVATE solver PRIVATE gtest PRIVATE gtest_main PRIVATE pthread)
+add_test(NAME test_limiter COMMAND test_limiter WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH})
+
+add_executable(test_physics test_physics.cpp)
+target_include_directories(test_physics PRIVATE ${PROJECT_SOURCE_DIR}/include)
+target_link_libraries(test_physics PRIVATE sim_setup PRIVATE solver PRIVATE utilities PRIVATE gtest PRIVATE gtest_main PRIVATE pthread)
+add_test(NAME test_physics COMMAND test_physics WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH})
+
+add_executable(test_Riemann_solvers tests_Riemann_solvers.cpp)
+target_include_directories(test_Riemann_solvers PRIVATE ${PROJECT_SOURCE_DIR}/include)
+target_link_libraries(test_Riemann_solvers PRIVATE sim_setup PRIVATE solver PRIVATE utilities PRIVATE gtest PRIVATE gtest_main PRIVATE pthread)
+add_test(NAME test_Riemann_solvers COMMAND test_Riemann_solvers WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH})
+
+add_executable(test_time_integrator test_time_integrator.cpp)
+target_include_directories(test_time_integrator PRIVATE ${PROJECT_SOURCE_DIR}/include)
+target_link_libraries(test_time_integrator PRIVATE sim_setup PRIVATE solver PRIVATE utilities PRIVATE gtest PRIVATE gtest_main PRIVATE pthread)
+add_test(NAME test_time_integrator COMMAND test_time_integrator WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH})
diff --git a/code_students/tests/test_limiter.cpp b/code_students/tests/test_limiter.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2cf481cd820b87aed1b6238be3da22b4c6121eb5
--- /dev/null
+++ b/code_students/tests/test_limiter.cpp
@@ -0,0 +1,54 @@
+#include "solver/limiter.hpp"
+
+#include <gtest/gtest.h>
+#include <iostream>
+#include <memory>
+
+// Choice of accuracy requirement for passing a test
+const double eps = 1.e-8;
+
+class TestLimiter : public ::testing::Test {
+protected:
+ virtual void SetUp() { ptr_limiter = std::make_unique<limiter_minmod>(1.0); }
+ virtual void TearDown() {}
+
+ std::unique_ptr<limiter_minmod> ptr_limiter;
+};
+
+TEST_F(TestLimiter, same_sign) {
+ // delta_right > 0.0
+ double derivative = ptr_limiter->compute(1.5, 3.0, 2.0);
+ EXPECT_NEAR(derivative, 1.5, eps);
+
+ derivative = ptr_limiter->compute(2.0, 3.0, 1.5);
+ EXPECT_NEAR(derivative, 1.5, eps);
+
+ derivative = ptr_limiter->compute(3.0, 2.0, 1.5);
+ EXPECT_NEAR(derivative, 1.5, eps);
+
+ // delta_right < 0.0
+ derivative = ptr_limiter->compute(-1.5, -3.0, -2.0);
+ EXPECT_NEAR(derivative, -1.5, eps);
+
+ ptr_limiter->compute(-2.0, -3.0, -1.5);
+ EXPECT_NEAR(derivative, -1.5, eps);
+
+ derivative = ptr_limiter->compute(-3.0, -2.0, -1.5);
+ EXPECT_NEAR(derivative, -1.5, eps);
+}
+
+TEST_F(TestLimiter, different_sign) {
+ double derivative = ptr_limiter->compute(1.5, 3.0, -2.0);
+ EXPECT_NEAR(derivative, 0.0, eps);
+
+ derivative = ptr_limiter->compute(-2.0, 3.0, 1.5);
+ EXPECT_NEAR(derivative, 0.0, eps);
+
+ derivative = ptr_limiter->compute(-3.0, -2.0, 1.5);
+ EXPECT_NEAR(derivative, 0.0, eps);
+}
+
+int main(int argc, char **argv) {
+ ::testing::InitGoogleTest(&argc, argv);
+ return RUN_ALL_TESTS();
+}
diff --git a/code_students/tests/test_physics.cpp b/code_students/tests/test_physics.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e0986a9c5fed9c03ee9aebec0d7f3f9cadd97ad5
--- /dev/null
+++ b/code_students/tests/test_physics.cpp
@@ -0,0 +1,169 @@
+#include "core/config.hpp"
+#include "setup/fluid.hpp"
+#include "setup/physics.hpp"
+
+#include <cmath>
+#include <gtest/gtest.h>
+#include <iostream>
+#include <memory>
+
+// Choice of accuracy requirement for passing a test
+const double eps = 1.e-8;
+
+class TestLimiter : public ::testing::Test {
+protected:
+ virtual void SetUp() { ptr_physics = std::make_unique<physics>(); }
+ virtual void TearDown() {}
+
+ std::unique_ptr<physics> ptr_physics;
+};
+
+TEST_F(TestLimiter, variable_transforms) {
+ std::cout << " Testing computation of thermal pressure\n";
+
+ fluid_cell quantities_local(parallelisation::FluidType::adiabatic);
+ quantities_local.fluid_data[quantities_local.get_index_density()] = 2.0;
+ quantities_local.fluid_data[quantities_local.get_index_v_x()] = 2.5;
+ quantities_local.fluid_data[quantities_local.get_index_v_y()] = 4.5;
+ quantities_local.fluid_data[quantities_local.get_index_v_z()] = 8.5;
+ quantities_local.fluid_data[quantities_local.get_index_energy()] = 0.2;
+ quantities_local.fluid_data[quantities_local.get_index_tracer()] = -3.2;
+
+ std::cout << " Testing computation of pressure\n";
+ double pressure = ptr_physics->get_pressure(quantities_local);
+
+ EXPECT_NEAR(pressure, 0.08, eps);
+
+ std::cout << " Testing computation of total energy\n";
+ double e_total = ptr_physics->get_e_total(quantities_local);
+
+ EXPECT_NEAR(e_total, 98.95, eps);
+
+ std::cout << " Testing transform to conservative variables\n";
+ ptr_physics->transform_characteristic_to_conservative(quantities_local);
+ EXPECT_NEAR(quantities_local.fluid_data[quantities_local.get_index_density()], 2.0, eps);
+ EXPECT_NEAR(quantities_local.fluid_data[quantities_local.get_index_v_x()], 5.0, eps);
+ EXPECT_NEAR(quantities_local.fluid_data[quantities_local.get_index_v_y()], 9.0, eps);
+ EXPECT_NEAR(quantities_local.fluid_data[quantities_local.get_index_v_z()], 17.0, eps);
+ EXPECT_NEAR(quantities_local.fluid_data[quantities_local.get_index_energy()], 98.95, eps);
+ EXPECT_NEAR(quantities_local.fluid_data[quantities_local.get_index_tracer()], -3.2, eps);
+
+ std::cout << " Testing transform to characteristic variables\n";
+ ptr_physics->transform_conservative_to_characteristic(quantities_local);
+ EXPECT_NEAR(quantities_local.fluid_data[quantities_local.get_index_density()], 2.0, eps);
+ EXPECT_NEAR(quantities_local.fluid_data[quantities_local.get_index_v_x()], 2.5, eps);
+ EXPECT_NEAR(quantities_local.fluid_data[quantities_local.get_index_v_y()], 4.5, eps);
+ EXPECT_NEAR(quantities_local.fluid_data[quantities_local.get_index_v_z()], 8.5, eps);
+ EXPECT_NEAR(quantities_local.fluid_data[quantities_local.get_index_energy()], 0.2, eps);
+ EXPECT_NEAR(quantities_local.fluid_data[quantities_local.get_index_tracer()], -3.2, eps);
+}
+
+TEST_F(TestLimiter, physical_fluxes) {
+ std::cout << " Testing physical fluxes\n";
+ fluid_cell quantities_local(parallelisation::FluidType::adiabatic);
+ quantities_local.fluid_data[quantities_local.get_index_density()] = 1.0;
+ quantities_local.fluid_data[quantities_local.get_index_v_x()] = 2.5;
+ quantities_local.fluid_data[quantities_local.get_index_v_y()] = 4.5;
+ quantities_local.fluid_data[quantities_local.get_index_v_z()] = 8.5;
+ quantities_local.fluid_data[quantities_local.get_index_energy()] = 0.2;
+ quantities_local.fluid_data[quantities_local.get_index_tracer()] = -3.2;
+
+ fluxes_cell fluxes_local(quantities_local.get_fluid_type());
+
+ ptr_physics->get_physical_fluxes(quantities_local, fluxes_local, parallelisation::direction::x);
+
+ EXPECT_NEAR(fluxes_local.flux_data[quantities_local.get_index_density()], 2.5, eps);
+ EXPECT_NEAR(fluxes_local.flux_data[quantities_local.get_index_v_x()], 6.33, eps);
+ EXPECT_NEAR(fluxes_local.flux_data[quantities_local.get_index_v_y()], 11.25, eps);
+ EXPECT_NEAR(fluxes_local.flux_data[quantities_local.get_index_v_z()], 21.25, eps);
+ EXPECT_NEAR(fluxes_local.flux_data[quantities_local.get_index_energy()], 124.1375, eps);
+ EXPECT_NEAR(fluxes_local.flux_data[quantities_local.get_index_tracer()], -8.0, eps);
+
+ ptr_physics->get_physical_fluxes(quantities_local, fluxes_local, parallelisation::direction::y);
+
+ EXPECT_NEAR(fluxes_local.flux_data[quantities_local.get_index_density()], 4.5, eps);
+ EXPECT_NEAR(fluxes_local.flux_data[quantities_local.get_index_v_x()], 11.25, eps);
+ EXPECT_NEAR(fluxes_local.flux_data[quantities_local.get_index_v_y()], 20.33, eps);
+ EXPECT_NEAR(fluxes_local.flux_data[quantities_local.get_index_v_z()], 38.25, eps);
+ EXPECT_NEAR(fluxes_local.flux_data[quantities_local.get_index_energy()], 223.4475, eps);
+ EXPECT_NEAR(fluxes_local.flux_data[quantities_local.get_index_tracer()], -14.4, eps);
+
+ ptr_physics->get_physical_fluxes(quantities_local, fluxes_local, parallelisation::direction::z);
+
+ EXPECT_NEAR(fluxes_local.flux_data[quantities_local.get_index_density()], 8.5, eps);
+ EXPECT_NEAR(fluxes_local.flux_data[quantities_local.get_index_v_x()], 21.25, eps);
+ EXPECT_NEAR(fluxes_local.flux_data[quantities_local.get_index_v_y()], 38.25, eps);
+ EXPECT_NEAR(fluxes_local.flux_data[quantities_local.get_index_v_z()], 72.33, eps);
+ EXPECT_NEAR(fluxes_local.flux_data[quantities_local.get_index_energy()], 422.0675, eps);
+ EXPECT_NEAR(fluxes_local.flux_data[quantities_local.get_index_tracer()], -27.2, eps);
+}
+
+TEST_F(TestLimiter, characteristic_speeds) {
+ std::cout << " Testing characteristic speeds\n";
+ fluid_cell quantities_local(parallelisation::FluidType::adiabatic);
+ quantities_local.fluid_data[quantities_local.get_index_density()] = 2.0;
+ quantities_local.fluid_data[quantities_local.get_index_v_x()] = 2.5;
+ quantities_local.fluid_data[quantities_local.get_index_v_y()] = 4.5;
+ quantities_local.fluid_data[quantities_local.get_index_v_z()] = 8.5;
+ quantities_local.fluid_data[quantities_local.get_index_energy()] = 0.2;
+ quantities_local.fluid_data[quantities_local.get_index_tracer()] = -3.2;
+
+ double pressure = ptr_physics->get_pressure(quantities_local);
+ double c_sound = ptr_physics->get_sound_speed(quantities_local.fluid_data[quantities_local.get_index_density()], pressure);
+
+ EXPECT_NEAR(c_sound * c_sound, 0.056, eps);
+
+ double lambda_abs_max = ptr_physics->get_lambda_abs_max(quantities_local);
+
+ double v_abs = sqrt(2.5 * 2.5 + 4.5 * 4.5 + 8.5 * 8.5);
+ EXPECT_NEAR(lambda_abs_max, v_abs + sqrt(0.056), eps);
+
+ fluid_cell quantities_previous(parallelisation::FluidType::adiabatic);
+ quantities_previous.fluid_data[quantities_previous.get_index_density()] = 1.0;
+ quantities_previous.fluid_data[quantities_previous.get_index_v_x()] = 1.5;
+ quantities_previous.fluid_data[quantities_previous.get_index_v_y()] = 2.5;
+ quantities_previous.fluid_data[quantities_previous.get_index_v_z()] = 4.5;
+ quantities_previous.fluid_data[quantities_previous.get_index_energy()] = 0.1;
+ quantities_previous.fluid_data[quantities_previous.get_index_tracer()] = -2.2;
+
+ double lambda_min, lambda_max;
+ ptr_physics->get_lambda_min_max(lambda_min, lambda_max, quantities_previous, quantities_local, parallelisation::direction::x);
+ EXPECT_NEAR(lambda_max, 2.5 + sqrt(0.056), eps);
+ // EXPECT_NEAR( lambda_min, 0.0, eps);
+ EXPECT_NEAR(lambda_min, 1.5 - sqrt(0.056), eps);
+
+ ptr_physics->get_lambda_min_max(lambda_min, lambda_max, quantities_previous, quantities_local, parallelisation::direction::y);
+ EXPECT_NEAR(lambda_max, 4.5 + sqrt(0.056), eps);
+ // EXPECT_NEAR( lambda_min, 0.0, eps);
+ EXPECT_NEAR(lambda_min, 2.5 - sqrt(0.056), eps);
+
+ ptr_physics->get_lambda_min_max(lambda_min, lambda_max, quantities_previous, quantities_local, parallelisation::direction::z);
+ EXPECT_NEAR(lambda_max, 8.5 + sqrt(0.056), eps);
+ // EXPECT_NEAR( lambda_min, 0.0, eps);
+ EXPECT_NEAR(lambda_min, 4.5 - sqrt(0.056), eps);
+
+ quantities_previous.fluid_data[quantities_local.get_index_density()] = 1.0;
+ quantities_previous.fluid_data[quantities_local.get_index_v_x()] = 1.5;
+ quantities_previous.fluid_data[quantities_local.get_index_v_y()] = 2.5;
+ quantities_previous.fluid_data[quantities_local.get_index_v_z()] = 4.5;
+ quantities_previous.fluid_data[quantities_local.get_index_energy()] = 22.0;
+ quantities_previous.fluid_data[quantities_local.get_index_tracer()] = -2.2;
+
+ // p=0.4*22=8.8 c_s**2 = 8.8/1.0*1.4=12.32
+ ptr_physics->get_lambda_min_max(lambda_min, lambda_max, quantities_previous, quantities_local, parallelisation::direction::x);
+ EXPECT_NEAR(lambda_max, 1.5 + sqrt(12.32), eps);
+ EXPECT_NEAR(lambda_min, 1.5 - sqrt(12.32), eps);
+
+ ptr_physics->get_lambda_min_max(lambda_min, lambda_max, quantities_previous, quantities_local, parallelisation::direction::y);
+ EXPECT_NEAR(lambda_max, 2.5 + sqrt(12.32), eps);
+ EXPECT_NEAR(lambda_min, 2.5 - sqrt(12.32), eps);
+
+ ptr_physics->get_lambda_min_max(lambda_min, lambda_max, quantities_previous, quantities_local, parallelisation::direction::z);
+ EXPECT_NEAR(lambda_max, 8.5 + sqrt(0.056), eps); // left-handed value still higher
+ EXPECT_NEAR(lambda_min, 4.5 - sqrt(12.32), eps);
+}
+
+int main(int argc, char **argv) {
+ ::testing::InitGoogleTest(&argc, argv);
+ return RUN_ALL_TESTS();
+}
diff --git a/code_students/tests/test_time_integrator.cpp b/code_students/tests/test_time_integrator.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0197f79d1b207b80cf17cc889c50f28efba89f83
--- /dev/null
+++ b/code_students/tests/test_time_integrator.cpp
@@ -0,0 +1,71 @@
+#include "setup/grid.hpp"
+#include "solver/time_integrator.hpp"
+
+#include <gtest/gtest.h>
+#include <iostream>
+#include <memory>
+
+// Choice of accuracy requirement for passing a test
+const double eps = 1.e-8;
+
+class TestTimeIntegrator : public ::testing::Test {
+protected:
+ virtual void SetUp() {
+ std::vector<double> bound_low(3), bound_up(3);
+ bound_low[0] = 0.0;
+ bound_low[1] = 0.0;
+ bound_low[2] = 0.0;
+
+ bound_up[0] = 1.0;
+ bound_up[1] = 1.0;
+ bound_up[2] = 1.0;
+
+ std::vector<int> num_cells(3);
+ num_cells[0] = 2;
+ num_cells[1] = 2;
+ num_cells[2] = 2;
+
+ grid = std::make_unique<grid_3D>(bound_low, bound_up, num_cells, 0);
+
+ std::cout << " Making fluid\n";
+ hd_fluid = std::make_unique<fluid>(parallelisation::FluidType::adiabatic);
+ std::cout << " setup of fluid\n";
+ hd_fluid->setup(*grid);
+
+ std::cout << " Making time stepper\n";
+ time_stepper = std::make_unique<RungeKutta2>(*grid, *hd_fluid);
+ std::cout << " Finished making time stepper\n";
+ }
+ virtual void TearDown() {}
+ std::unique_ptr<grid_3D> grid;
+ std::unique_ptr<fluid> hd_fluid;
+ std::unique_ptr<RungeKutta2> time_stepper;
+};
+
+TEST_F(TestTimeIntegrator, exponential) {
+ std::cout << " Testing time integration of an exponential\n";
+ fluid fluid_changes(hd_fluid->get_fluid_type());
+ fluid_changes.setup(*grid);
+ hd_fluid->fluid_data[0](0, 0, 0) = 1.0;
+ double lambda = -2.5;
+ double delta_t = 0.2;
+
+ // RK step 0:
+ fluid_changes.fluid_data[0](0, 0, 0) = lambda * hd_fluid->fluid_data[0](0, 0, 0);
+
+ time_stepper->do_sub_step(*grid, fluid_changes, *hd_fluid, delta_t, 0);
+ EXPECT_NEAR(hd_fluid->fluid_data[0](0, 0, 0), 1.0 + delta_t * lambda, eps);
+ std::cout << " Result(RK step 1): " << hd_fluid->fluid_data[0](0, 0, 0) << "\n";
+
+ // RK step 1:
+ fluid_changes.fluid_data[0](0, 0, 0) = lambda * hd_fluid->fluid_data[0](0, 0, 0);
+
+ time_stepper->do_sub_step(*grid, fluid_changes, *hd_fluid, delta_t, 1);
+ EXPECT_NEAR(hd_fluid->fluid_data[0](0, 0, 0), 1.0 + delta_t * lambda + 0.5 * delta_t * delta_t * lambda * lambda, eps);
+ std::cout << " Result(RK step 2): " << hd_fluid->fluid_data[0](0, 0, 0) << "\n";
+}
+
+int main(int argc, char **argv) {
+ ::testing::InitGoogleTest(&argc, argv);
+ return RUN_ALL_TESTS();
+}
diff --git a/code_students/tests/tests_Riemann_solvers.cpp b/code_students/tests/tests_Riemann_solvers.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d39903b2a17a7838830a0fe9fdc8349035f349ec
--- /dev/null
+++ b/code_students/tests/tests_Riemann_solvers.cpp
@@ -0,0 +1,102 @@
+#include "core/config.hpp"
+#include "setup/fluid.hpp"
+#include "setup/physics.hpp"
+#include "solver/Riemann_solvers.hpp"
+
+#include <cmath>
+#include <gtest/gtest.h>
+#include <iostream>
+#include <memory>
+
+// Choice of accuracy requirement for passing a test
+const double eps = 1.e-8;
+
+class TestRiemannSolvers : public ::testing::Test {
+protected:
+ virtual void SetUp() {
+ fluid_cell hd_fluid(parallelisation::FluidType::adiabatic);
+ Riemann_HLL = std::make_unique<HLL_solver>(hd_fluid.get_num_fields());
+ ptr_physics = std::make_unique<physics>();
+ fluid_type = parallelisation::FluidType::adiabatic;
+ }
+ virtual void TearDown() {}
+ std::unique_ptr<HLL_solver> Riemann_HLL;
+ std::unique_ptr<physics> ptr_physics;
+ parallelisation::FluidType fluid_type;
+};
+
+TEST_F(TestRiemannSolvers, hll) {
+ std::cout << " Testing the Riemann solver\n";
+ fluid_cell quantities_previous(fluid_type);
+ quantities_previous.fluid_data[quantities_previous.get_index_density()] = 1.0;
+ quantities_previous.fluid_data[quantities_previous.get_index_v_x()] = 1.5;
+ quantities_previous.fluid_data[quantities_previous.get_index_v_y()] = 2.5;
+ quantities_previous.fluid_data[quantities_previous.get_index_v_z()] = 4.5;
+ quantities_previous.fluid_data[quantities_previous.get_index_energy()] = 0.1;
+ quantities_previous.fluid_data[quantities_previous.get_index_tracer()] = -2.2;
+
+ fluid_cell quantities_local(fluid_type);
+ quantities_local.fluid_data[quantities_local.get_index_density()] = 2.0;
+ quantities_local.fluid_data[quantities_local.get_index_v_x()] = 2.5;
+ quantities_local.fluid_data[quantities_local.get_index_v_y()] = 4.5;
+ quantities_local.fluid_data[quantities_local.get_index_v_z()] = 8.5;
+ quantities_local.fluid_data[quantities_local.get_index_energy()] = 0.2;
+ quantities_local.fluid_data[quantities_local.get_index_tracer()] = -3.2;
+
+ fluxes_cell fluxes_previous(fluid_type), fluxes_local(fluid_type);
+ fluxes_previous.flux_data[quantities_local.get_index_density()] = 1.5;
+ fluxes_previous.flux_data[quantities_local.get_index_v_x()] = 2.5;
+ fluxes_previous.flux_data[quantities_local.get_index_v_y()] = 3.5;
+ fluxes_previous.flux_data[quantities_local.get_index_v_z()] = 4.5;
+ fluxes_previous.flux_data[quantities_local.get_index_energy()] = 5.5;
+ fluxes_previous.flux_data[quantities_local.get_index_tracer()] = 6.5;
+
+ fluxes_local.flux_data[quantities_local.get_index_density()] = 6.5;
+ fluxes_local.flux_data[quantities_local.get_index_v_x()] = 7.5;
+ fluxes_local.flux_data[quantities_local.get_index_v_y()] = 8.5;
+ fluxes_local.flux_data[quantities_local.get_index_v_z()] = 9.5;
+ fluxes_local.flux_data[quantities_local.get_index_energy()] = 10.5;
+ fluxes_local.flux_data[quantities_local.get_index_tracer()] = 11.5;
+
+ // Test Riemann fan tipped to the right
+ double lambda_min_x = 2.0;
+ double lambda_max_x = 3.0;
+
+ fluxes_cell num_flux(quantities_local.get_fluid_type());
+ Riemann_HLL->get_num_flux(quantities_previous, quantities_local, fluxes_previous.flux_data, fluxes_local.flux_data, num_flux.flux_data, lambda_min_x,
+ lambda_max_x);
+
+ EXPECT_NEAR(num_flux.flux_data[0], 1.5, eps);
+ EXPECT_NEAR(num_flux.flux_data[3], 4.5, eps);
+
+ std::cout << " Left-handed flux " << num_flux.flux_data[0] << "\n";
+
+ // Test Riemann fan tipped to the left
+ lambda_min_x = -3.0;
+ lambda_max_x = -2.0;
+
+ Riemann_HLL->get_num_flux(quantities_previous, quantities_local, fluxes_previous.flux_data, fluxes_local.flux_data, num_flux.flux_data, lambda_min_x,
+ lambda_max_x);
+
+ EXPECT_NEAR(num_flux.flux_data[1], 7.5, eps);
+ EXPECT_NEAR(num_flux.flux_data[4], 10.5, eps);
+
+ std::cout << " Right-handed flux " << num_flux.flux_data[1] << "\n";
+
+ // Test local Riemann fan
+ lambda_min_x = -2.0;
+ lambda_max_x = 3.0;
+
+ Riemann_HLL->get_num_flux(quantities_previous, quantities_local, fluxes_previous.flux_data, fluxes_local.flux_data, num_flux.flux_data, lambda_min_x,
+ lambda_max_x);
+
+ EXPECT_NEAR(num_flux.flux_data[2], 3.1, eps);
+ EXPECT_NEAR(num_flux.flux_data[5], 9.7, eps);
+
+ std::cout << " Riemann-fan-internal flux " << num_flux.flux_data[2] << "\n";
+}
+
+int main(int argc, char **argv) {
+ ::testing::InitGoogleTest(&argc, argv);
+ return RUN_ALL_TESTS();
+}