Add main application template

code_students/apps/CMakeLists.txt

+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

code_students/apps/main_full_code.cpp

+#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

code_students/apps/main_test_singlesteps.cpp

+#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

code_students/include/IO/data_storage.hpp

+#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

code_students/include/core/config.hpp

+#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

code_students/include/setup/fluid.hpp

+#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

code_students/include/setup/grid.hpp

+#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

code_students/include/setup/grid1D.hpp

+#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

code_students/include/setup/physics.hpp

+#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

code_students/include/solver/Riemann_solvers.hpp

+#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

code_students/include/solver/finite_volume_solver.hpp

+#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 &current_fluid);
+	double singlestep(grid_3D &spatial_grid, fluid &current_fluid, fluid &current_changes);
+	int run(grid_3D &spatial_grid, fluid &current_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 &current_fluid, int index_data);
+	void set_verbosity(int verbosity);
+
+private:
+	double get_CFL(grid_3D &spatial_grid, fluid &current_fluid);
+	void set_initial_conditions(grid_3D &spatial_grid, fluid &current_fluid);
+	void transform_fluid_to_conservative(fluid &current_fluid);
+	void transform_fluid_to_characteristic(fluid &current_fluid);
+	void apply_boundary_conditions(grid_3D &spatial_grid, fluid &current_fluid);
+	void store_timestep(grid_3D &spatial_grid, fluid &current_fluid);
+	double compute_delta_t_next(grid_3D &spatial_grid, fluid &current_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

code_students/include/solver/limiter.hpp

+#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

code_students/include/solver/reconstruction.hpp

+#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

code_students/include/solver/time_integrator.hpp

+#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

code_students/include/util/matrix.hpp

+#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

code_students/include/util/utility_functions.hpp

+#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

code_students/src/CMakeLists.txt

+cmake_minimum_required(VERSION 3.10)
+
+add_subdirectory(util)
+add_subdirectory(setup)
+add_subdirectory(solver)
+add_subdirectory(IO)

code_students/src/IO/CMakeLists.txt

+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

code_students/src/IO/data_storage.cpp

+#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();
+}

code_students/src/setup/CMakeLists.txt

+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