diff --git a/example/SparseGrid/9_inhomogeneous_diffusion_porous_catalyst_CaCO3/config.cfg b/example/SparseGrid/9_inhomogeneous_diffusion_porous_catalyst_CaCO3/config.cfg
new file mode 100755
index 0000000000000000000000000000000000000000..067b942df5494321fd7d5440df1b7dff32b4f11a
--- /dev/null
+++ b/example/SparseGrid/9_inhomogeneous_diffusion_porous_catalyst_CaCO3/config.cfg
@@ -0,0 +1,2 @@
+[pack]
+files = main.cu Makefile
\ No newline at end of file
diff --git a/example/SparseGrid/9_inhomogeneous_diffusion_porous_catalyst_CaCO3/include/DiffusionSpace_sparseGrid.hpp b/example/SparseGrid/9_inhomogeneous_diffusion_porous_catalyst_CaCO3/include/DiffusionSpace_sparseGrid.hpp
new file mode 100755
index 0000000000000000000000000000000000000000..d613e0a43d3d205f0cf966f20052d1dd8d3b10d9
--- /dev/null
+++ b/example/SparseGrid/9_inhomogeneous_diffusion_porous_catalyst_CaCO3/include/DiffusionSpace_sparseGrid.hpp
@@ -0,0 +1,295 @@
+//
+// Created by jstark on 28.09.21.
+//
+
+#ifndef SPARSEGRID_DIFFUSIONSPACE_SPARSEGRID_HPP
+#define SPARSEGRID_DIFFUSIONSPACE_SPARSEGRID_HPP
+
+#include "math.h"
+
+template <typename T>
+static bool is_diffusionSpace(const T & phi_sdf, const T & lower_bound, const T & upper_bound)
+{
+ const T EPSILON = std::numeric_limits<T>::epsilon();
+ const T _b_low = lower_bound + EPSILON;
+ const T _b_up = upper_bound - EPSILON;
+ return (phi_sdf > _b_low && phi_sdf < _b_up);
+}
+
+template <size_t PHI_SDF_ECS_FULL, size_t PHI_SDF_ECS_SPARSE, typename grid_type, typename sparse_in_type, typename T>
+void get_diffusion_domain_sparse_grid(grid_type & grid,
+ sparse_in_type & sparse_grid,
+ const T b_low,
+ const T b_up)
+{
+ auto dom = grid.getDomainIterator();
+ while(dom.isNext())
+ {
+ auto key = dom.get();
+ if(is_diffusionSpace(grid.template get<PHI_SDF_ECS_FULL>(key), b_low, b_up))
+ {
+ sparse_grid.template insertFlush<PHI_SDF_ECS_SPARSE>(key) = grid.template get<PHI_SDF_ECS_FULL>(key);
+ }
+ ++dom;
+ }
+ // Mapping?
+ // Load Balancing?
+}
+
+#if 0
+template <size_t PHI_SDF_ECS_FULL, size_t PHI_SDF_ECS_SPARSE, typename grid_type, typename grid_type_upscale, typename sparse_in_type, typename T>
+void get_diffusion_domain_sparse_grid_upscale(grid_type & grid,
+ grid_type_upscale & grid_upscale,
+ sparse_in_type & sparse_grid,
+ const int upscale_factor, // total upscale factor that is product all dimensions, e.g. 2^3 = 8
+ const T b_low,
+ const T b_up)
+{
+ auto dom = grid.getDomainIterator();
+ auto dom_upscale = grid_upscale.getDomainIterator();
+ while(dom.isNext())
+ {
+ auto key = dom.get();
+ auto key_upscale = dom_upscale.get();
+
+ for(int i = 0, i < upscale_factor, ++i) // Place upscale number of points for each key
+ {
+ if(is_diffusionSpace(grid.template get<PHI_SDF_ECS_FULL>(key), b_low, b_up))
+ {
+ sparse_grid.template insertFlush<PHI_SDF_ECS_SPARSE>(key_upscale) = grid.template get<PHI_SDF_ECS_FULL>(key);
+ }
+ ++dom_upscale
+ }
+
+
+ ++dom;
+ }
+
+}
+#endif
+
+template <
+size_t PHI_SDF_ECS_FULL,
+size_t PHI_SDF_SHELL_FULL,
+size_t PHI_SDF_ECS_SPARSE,
+size_t PHI_SDF_SHELL_SPARSE,
+typename grid_type,
+typename sparse_in_type,
+typename T>
+void get_diffusion_domain_sparse_grid_with_shell(grid_type & grid_ecs,
+ grid_type & grid_shell,
+ sparse_in_type & sparse_grid,
+ const T b_low,
+ const T b_up)
+{
+ auto dom = grid_ecs.getDomainIterator();
+ while(dom.isNext())
+ {
+ auto key = dom.get();
+ if(is_diffusionSpace(grid_ecs.template get<PHI_SDF_ECS_FULL>(key), b_low, b_up))
+ {
+ sparse_grid.template insertFlush<PHI_SDF_ECS_SPARSE>(key) = grid_ecs.template get<PHI_SDF_ECS_FULL>(key);
+ sparse_grid.template insertFlush<PHI_SDF_SHELL_SPARSE>(key) = grid_shell.template get<PHI_SDF_SHELL_FULL>(key);
+ }
+ ++dom;
+ }
+}
+
+
+#if 0
+// Init c0
+ auto dom = g_sparse.getDomainIterator();
+ while(dom.isNext())
+ {
+ auto key = dom.get();
+ if (g_sparse.template getProp<PHI_SDF>(key) < 4 * g_sparse.spacing(x))
+ {
+ g_sparse.template getProp<CONC_N>(key) =
+ (1 - 0.25 * g_sparse.template getProp<PHI_SDF>(key)) / g_sparse.spacing(x);
+ }
+ ++dom;
+ }
+ if(v_cl.rank() == 0) std::cout << "Initialized grid with smoothed c0 at the boundary." << std::endl;
+ g_sparse.write(path_output + "g_init", FORMAT_BINARY);
+#endif
+
+
+
+
+template <typename point_type, typename y_margin_type>
+auto distance_from_margin(point_type & coord, y_margin_type y_margin)
+{
+ return(coord.get(1) - y_margin);
+}
+
+template <typename point_type, typename y_margin_type, typename width_type>
+bool is_source(point_type & coord, y_margin_type y_margin, width_type width_source)
+{
+ return (coord.get(1) >= y_margin
+ && distance_from_margin(coord, y_margin) <= width_source);
+}
+
+template <typename T>
+static bool is_inner_surface(const T & phi_sdf, const T & lower_bound)
+{
+ const T EPSILON = std::numeric_limits<T>::epsilon();
+ const T _b_low = lower_bound + EPSILON;
+ return (phi_sdf > _b_low);
+}
+
+template <size_t PHI_SDF, size_t K_SOURCE, size_t K_SINK, typename grid_type, typename y_margin_type, typename width_type, typename k_type>
+void init_reactionTerms(grid_type & grid,
+ width_type width_source,
+ y_margin_type y_margin,
+ k_type k_source,
+ k_type k_sink,
+ int no_membrane_points=4)
+{
+ /*
+ * Setting the reaction terms according to their distance from the margin of the blastoderm along the dorsal-ventral axis
+ *
+ * */
+
+ auto dom = grid.getDomainIterator();
+ while(dom.isNext())
+ {
+ auto key = dom.get();
+ Point<grid_type::dims, y_margin_type> coord = grid.getPos(key);
+
+ if(grid.template get<PHI_SDF>(key) < no_membrane_points * grid.spacing(0)) // If point lies close to the cell surface
+ {
+ if(is_source(coord, y_margin, width_source)) // If point lies within width_source distance from the margin, set k_source
+ {
+ grid.template insertFlush<K_SOURCE>(key) = k_source;
+ }
+ else
+ {
+ grid.template insertFlush<K_SOURCE>(key) = 0.0; // If membrane point not in source, emission is zero
+ }
+ grid.template insertFlush<K_SINK>(key) = k_sink; // Have sink at all membranes
+
+ }
+
+ else // For the nodes that are further away from the membrane, set both reaction terms to zero
+ {
+ grid.template insertFlush<K_SOURCE>(key) = 0.0;
+ grid.template insertFlush<K_SINK>(key) = 0.0;
+ }
+ ++dom;
+ }
+}
+
+template <
+size_t PHI_SDF_ECS_SPARSE,
+size_t PHI_SDF_SHELL_SPARSE,
+size_t K_SOURCE,
+size_t K_SINK,
+typename grid_type,
+typename y_margin_type,
+typename width_type,
+typename k_type,
+typename boundary_type>
+void init_reactionTerms_with_shell(grid_type & grid,
+ width_type width_source,
+ y_margin_type y_margin,
+ k_type k_source,
+ k_type k_sink,
+ boundary_type b_low_shell,
+ int no_membrane_points=4)
+{
+ /*
+ * Setting the reaction terms according to their distance from the margin of the blastoderm along the dorsal-ventral axis
+ *
+ * */
+
+ auto dom = grid.getDomainIterator();
+ while(dom.isNext())
+ {
+ auto key = dom.get();
+ Point<grid_type::dims, y_margin_type> coord = grid.getPos(key);
+
+ if(
+ grid.template get<PHI_SDF_ECS_SPARSE>(key) < no_membrane_points * grid.spacing(0) // If point is a surface point
+ && is_inner_surface(grid.template get<PHI_SDF_SHELL_SPARSE>(key), b_low_shell) // and this surface is not towards the yolk or EVL
+ )
+ {
+ if(is_source(coord, y_margin, width_source)) // If point lies within width_source distance from the margin, set k_source
+ {
+ grid.template insertFlush<K_SOURCE>(key) = k_source;
+ }
+ else
+ {
+ grid.template insertFlush<K_SOURCE>(key) = 0.0; // If membrane point not in source, emission is zero
+ }
+ grid.template insertFlush<K_SINK>(key) = k_sink; // Have sink at all membranes
+
+ }
+
+ else // For the nodes that are further away from the membrane or belong to the outer surface, set both reaction terms to zero
+ {
+ grid.template insertFlush<K_SOURCE>(key) = 0.0;
+ grid.template insertFlush<K_SINK>(key) = 0.0;
+ }
+ ++dom;
+ }
+}
+
+
+template <size_t PHI_SDF, size_t K_SOURCE, size_t K_SINK, typename grid_type, typename coord_type, typename k_type>
+void init_reactionTerms_smoothed(grid_type & grid,
+ coord_type y_start,
+ coord_type y_peak,
+ coord_type y_stop,
+ k_type k_source,
+ k_type k_sink,
+ coord_type no_membrane_points=4,
+ coord_type smoothing=0.25)
+{
+ /*
+ * Setting the reaction terms according to the position along the animal-vegetal axis = y-axis
+ *
+ * */
+ const int x = 0;
+ const int y = 1;
+ const int z = 2;
+
+ auto dom = grid.getDomainIterator();
+ while(dom.isNext())
+ {
+ auto key = dom.get();
+ Point<grid_type::dims, coord_type> coord = grid.getPos(key);
+// k_type sdf_based_smoothing = (1 - smoothing * grid.template getProp<PHI_SDF>(key)) / grid.spacing(x);
+ k_type sdf_based_smoothing = 1.0;
+ if(grid.template get<PHI_SDF>(key) < no_membrane_points * grid.spacing(0)) // If grid node lies close to the
+ // membrane
+ {
+ // If membrane on animal side of the peak, emission strength increases with y
+ if (coord[y] > y_start
+ && coord[y] < y_peak + std::numeric_limits<coord_type>::epsilon())
+ {
+ grid.template insertFlush<K_SOURCE>(key) = k_source * (coord[y] - y_start) * sdf_based_smoothing;
+ }
+ // Else if membrane on vegetal side of the peak, emission strength decreases with y as moving towards the
+ // boundary to the yolk
+ else if (coord[y] >= y_peak - std::numeric_limits<coord_type>::epsilon()
+ && coord[y] < y_stop - std::numeric_limits<coord_type>::epsilon())
+ {
+ grid.template insertFlush<K_SOURCE>(key) = k_source * (y_stop - coord[y]) * sdf_based_smoothing;
+ }
+ // Else source is zero
+ else grid.template insertFlush<K_SOURCE>(key) = 0.0;
+ grid.template insertFlush<K_SINK>(key) = k_sink; // Have sink at all membranes
+ }
+
+ else // For the nodes that are further away from the membrane, set the reaction terms to zero
+ {
+ grid.template insertFlush<K_SOURCE>(key) = 0.0;
+ grid.template insertFlush<K_SINK>(key) = 0.0;
+ }
+ ++dom;
+ }
+}
+
+
+
+#endif //SPARSEGRID_DIFFUSIONSPACE_SPARSEGRID_HPP
diff --git a/example/SparseGrid/9_inhomogeneous_diffusion_porous_catalyst_CaCO3/include/HelpFunctions_diffusion.hpp b/example/SparseGrid/9_inhomogeneous_diffusion_porous_catalyst_CaCO3/include/HelpFunctions_diffusion.hpp
new file mode 100755
index 0000000000000000000000000000000000000000..7e3069c8549f142ae2f3b31a730533d196530365
--- /dev/null
+++ b/example/SparseGrid/9_inhomogeneous_diffusion_porous_catalyst_CaCO3/include/HelpFunctions_diffusion.hpp
@@ -0,0 +1,195 @@
+//
+// Created by jstark on 08.10.21
+//
+
+#ifndef DIFFUSION_HELPFUNCTIONSDIFFUSION_HPP
+#define DIFFUSION_HELPFUNCTIONSDIFFUSION_HPP
+#include "util/PathsAndFiles.hpp"
+
+#include "level_set/redistancing_Sussman/HelpFunctionsForGrid.hpp"
+
+/**@brief Get timestep that fulfills stability criterion nD-diffusion using a FTCS scheme
+ *
+ * @tparam grid_type Template type of input g_sparse.
+ * @tparam T Template type of diffusion coefficient and timestep.
+ * @param grid Input grid of type grid_type.
+ * @param k_max Max. diffusion coefficient in the domain.
+ * @return T type timestep.
+ */
+template <typename grid_type, typename T>
+T diffusion_time_step(grid_type & grid, T k_max)
+{
+ T sum = 0;
+ for(int d = 0; d < grid_type::dims; ++d)
+ {
+ sum += (T)1.0 / (T)(grid.spacing(d) * grid.spacing(d));
+ }
+
+ return 1.0 / ((T)4.0 * k_max * sum);
+}
+
+/**@brief Sigmoidal which can be shifted and stretched in order to get desired smooth inhomog. diffusion-coeff scalar
+ * field
+ *
+ * @tparam T Template type of coordinates.
+ * @param x T type abscissa axis of sigmoidal function.
+ * @param x_shift T type shift of sigmoidal along abscissa.
+ * @param x_stretch T type stretching of sigmoidal along acscissa - the smaller, the steeper!
+ * @param y_min T type asymptotic lowest y-value.
+ * @param y_max T type asymptotic maximal y-value.
+ * @return Function value of sigmoidal for input value x.
+ */
+template <typename T>
+T get_smooth_sigmoidal(T x, T x_shift, T x_stretch, T y_min, T y_max)
+{
+ return y_min + y_max / (1.0 + exp(- (x_shift + x_stretch * x)));
+}
+
+/**@brief Sums up property values over all particles in grid.
+ *
+ * @tparam Prop Index of property whose values will be summed up over all particles.
+ * @tparam grid_type Template type of input particle vector.
+ * @param grid Input particle vector. Can be a subset.
+ * @return Sum of values stored in Prop.
+ */
+template <size_t Prop, typename grid_type>
+auto sum_prop_over_grid(grid_type & grid)
+{
+ auto sum = 0.0;
+ auto dom = grid.getDomainIterator();
+ while(dom.isNext())
+ {
+ auto key = dom.get();
+ sum += grid.template getProp<Prop>(key);
+ ++dom;
+ }
+ auto &v_cl = create_vcluster();
+ v_cl.sum(sum);
+ v_cl.execute();
+ return sum;
+}
+
+/**@brief Writing out the total mass to a csv file.
+ *
+ * @tparam Conc Index of property containing the concentration.
+ * @tparam T_mass Template type of initial mass.
+ * @param grid Input particle vector_dist of type grid_type.
+ * @param m_initial Initial total mass.
+ * @param t Current time of the diffusion.
+ * @param i Current iteration of the diffusion.
+ * @param path_output std::string Path where csv file will be saved.
+ * @param filename std::string Name of csv file. Default: "/total_mass.csv"
+ * @Outfile 4 columns: Time, Total mass, Mass difference, Max mass
+ */
+template <size_t Conc, typename grid_type, typename T, typename T_mass>
+void monitor_total_mass(grid_type & grid, const T_mass m_initial, const T p_volume, const T t, const size_t i,
+ const std::string & path_output, const std::string & filename="total_mass.csv")
+{
+ auto &v_cl = create_vcluster();
+
+// if (m_initial == 0)
+// {
+// if (v_cl.rank() == 0) std::cout << "m_initial is zero! Normalizing the total mass with the initial mass will "
+// "not work. Aborting..." << std::endl;
+// abort();
+// }
+
+ T m_total = sum_prop_over_grid<Conc>(grid) * p_volume;
+ T m_diff = m_total - m_initial;
+ T m_max = get_max_val<Conc>(grid) * p_volume;
+
+ if (v_cl.rank() == 0)
+ {
+ std::string outpath = path_output + "/" + filename;
+ create_file_if_not_exist(outpath);
+ std::ofstream outfile;
+ outfile.open(outpath, std::ios_base::app); // append instead of overwrite
+
+ if (i == 0)
+ {
+ outfile << "Time, Total mass, Mass difference, Max mass" << std::endl;
+// std::cout << "Time, Total mass, Total mass in-/decrease, Max mass" << std::endl;
+ }
+
+ outfile
+ << to_string_with_precision(t, 6) << ","
+ << to_string_with_precision(m_total, 6) << ","
+ << to_string_with_precision(m_diff, 6) << ","
+ << to_string_with_precision(m_max, 6)
+ << std::endl;
+ outfile.close();
+#if 0
+ std::cout
+ << to_string_with_precision(t, 6) << ","
+ << to_string_with_precision(m_total, 6) << ","
+ << to_string_with_precision(m_diff, 6) << ","
+ << to_string_with_precision(m_max, 6)
+ << std::endl;
+#endif
+// if (m_domain_normalized > 100 || m_domain_normalized < 0)
+// {
+// std::cout << "Mass increases or is negative, something went wrong. Aborting..." << std::endl;
+// abort();
+// }
+ }
+}
+
+/**@brief Writing out the total concentration over time to a csv file.
+ *
+ * @tparam Conc Index of property containing the concentration.
+ * @param grid Input particle vector_dist of type grid_type.
+ * @param t Current time of the diffusion.
+ * @param i Current iteration of the diffusion.
+ * @param path_output std::string Path where csv file will be saved.
+ * @param filename std::string Name of csv file. Default: "/normalized_total_mass.csv"
+ * @Outfile 3 columns: timestep of writeout i, total concentration
+ */
+template <size_t Conc, typename grid_type, typename T>
+void monitor_total_concentration(grid_type & grid, const T t, const size_t i,
+ const std::string & path_output, const std::string & filename="total_conc.csv")
+{
+ auto &v_cl = create_vcluster();
+ T c_total = sum_prop_over_grid<Conc>(grid);
+
+ if (v_cl.rank() == 0)
+ {
+ std::string outpath = path_output + "/" + filename;
+ create_file_if_not_exist(outpath);
+ std::ofstream outfile;
+ outfile.open(outpath, std::ios_base::app); // append instead of overwrite
+
+ if (i == 0)
+ {
+ outfile << "Time, Total concentration" << std::endl;
+ }
+
+ outfile
+ << to_string_with_precision(t, 6) << ","
+ << to_string_with_precision(c_total, 6)
+ << std::endl;
+ outfile.close();
+ }
+}
+
+/**@brief Sets the emission term to 0 if concentration surpasses a threshold. This is to avoid accumulation of mass
+ * at unconnected islands.
+ *
+ * @tparam Conc Index of property containing the concentration.
+ * @tparam Emission Index of property containing the emission constant.
+ * @tparam grid_type Template type of input particle vector.
+ * @tparam key_vector_type Template type of particle keys.
+ * @param grid Input particle vector. Can be a subset.
+ * @param keys Keys for which necessity for adaptation is checked.
+ * @param threshold Maximal concentration that is tolerated for keeping emission.
+ */
+template <size_t Conc, size_t Emission, typename grid_type, typename key_vector_type, typename T>
+void adapt_emission(grid_type & grid, key_vector_type & keys, T threshold)
+{
+ for (int i = 0; i < keys.size(); ++i)
+ {
+ auto key = keys.get(i);
+ if (grid.template getProp<Conc>(key) > threshold) grid.template getProp<Emission>(key) = 0.0;
+ }
+}
+
+#endif //DIFFUSION_HELPFUNCTIONSDIFFUSION_HPP
diff --git a/example/SparseGrid/9_inhomogeneous_diffusion_porous_catalyst_CaCO3/main.cu b/example/SparseGrid/9_inhomogeneous_diffusion_porous_catalyst_CaCO3/main.cu
new file mode 100755
index 0000000000000000000000000000000000000000..b607f52d5cbd5e0f7022882eead120535a5746bc
--- /dev/null
+++ b/example/SparseGrid/9_inhomogeneous_diffusion_porous_catalyst_CaCO3/main.cu
@@ -0,0 +1,498 @@
+//
+// Created by jstark on 2023-04-24.
+//
+/**
+ * @file 9_inhomogeneous_diffusion_porous_catalyst/main.cu
+ * @page 9_inhomogeneous_diffusion_porous_catalyst diffusion_porous_media 3D
+ *
+ * # Simulate inhomogeneous diffusion in a CaCO\f$_3\f$ particle packed bed #
+ *
+ * In this example, we simulate inhomogeneous diffusion in the fluid phase of a CaCO\f$_3\f$ particle packed bed. For the complete image-based simulation pipeline and performance of this simulation, we refer to
+ * <a href="https://arxiv.org/abs/2304.11165">J. Stark, I. F. Sbalzarini "An open-source pipeline for solving continuous reaction-diffusion models in image-based geometries of porous media." arXiv preprint arXiv:2304.11165 (2023).</a>
+
+ * The geometry of the fluid phase is reconstructed based on $\upmu$CT images provided by kindly provided by Prof. Jörg Petrasch (Michigan State University, College of Engineering) (see: <a href="https://asmedigitalcollection.asme.org/heattransfer/article/131/7/072701/444766/Tomographic-Characterization-of-a-Semitransparent">S. Haussener et al. “Tomographic characterization of a semitransparent-particle packed bed and determination of its thermal radiative propertiesâ€. In: Journal of Heat Transfer 131.7 (2009)</a>).
+
+ * For image based reconstruction and redistancing see @ref example_sussman_images_3D.
+ *
+ */
+
+/**
+ * @page 9_inhomogeneous_diffusion_porous_catalyst diffusion_porous_media 3D
+ *
+ * ## Include ## {#e_CaC03_include}
+ *
+ * We start with inluding header files, setting some paths and indices:
+ *
+ * \snippet SparseGrid/9_inhomogeneous_diffusion_porous_catalyst/main.cu Include
+ *
+ */
+//! \cond [Include] \endcond
+
+
+#include <iostream>
+#include <typeinfo>
+
+#include "CSVReader/CSVReader.hpp"
+
+// Include redistancing files
+#include "util/PathsAndFiles.hpp"
+#include "level_set/redistancing_Sussman/RedistancingSussman.hpp"
+#include "RawReader/InitGridWithPixel.hpp"
+#include "RemoveLines.hpp" // For removing thin (diagonal or straight) lines
+
+#include "level_set/redistancing_Sussman/HelpFunctionsForGrid.hpp"
+#include "level_set/redistancing_Sussman/AnalyticalSDF.hpp"
+#include "FiniteDifference/FD_simple.hpp"
+#include "Decomposition/Distribution/BoxDistribution.hpp"
+#include "timer.hpp"
+
+// Help functions for simulating diffusion in the image-based geometry
+#include "include/DiffusionSpace_sparseGrid.hpp"
+#include "include/HelpFunctions_diffusion.hpp"
+
+
+
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+// input
+const std::string path_to_redistancing_result =
+"/INPUT_PATH/benchmarks/CaCO3/sussman_redistancing/build/output_sussman_maxIter6e3_CaCO3_fluidPhase_531x531x531/";
+
+
+const std::string redistancing_filename = "grid_CaCO3_post_redistancing.hdf5";
+const std::string path_to_size = path_to_redistancing_result;
+
+// output
+const std::string output_name = "output_inhomogDiffusion_CaCO3_fluid_phase";
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Property indices full grid
+const size_t PHI_FULL = 0;
+
+typedef aggregate<float> props_full;
+const openfpm::vector<std::string> prop_names_full = {"Phi_Sussman_Out"};
+
+
+// Property indices sparse grid
+const size_t PHI_PHASE = 0;
+const size_t CONC_N = 1;
+const size_t CONC_NPLUS1 = 2;
+const size_t DIFFUSION_COEFFICIENT = 3;
+
+
+typedef aggregate<float, float, float, float> props_sparse;
+const openfpm::vector<std::string> prop_names_sparse = {"PHI_PHASE",
+ "CONC_N",
+ "CONC_NPLUS1",
+ "DIFFUSION_COEFFICIENT"};
+
+// Space indices
+constexpr size_t x = 0, y = 1, z = 2;
+
+// Parameters for grid size
+const size_t dims = 3;
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+//! \cond [Include] \endcond
+
+/**
+ * @page 9_inhomogeneous_diffusion_porous_catalyst diffusion_porous_media 3D
+ *
+ * ## Initialization and output folder ## {#e_CaC03_init}
+ *
+ * We start with
+ * * Initializing OpenFPM
+ * * Setting the output path and creating an output folder
+ *
+ * \snippet SparseGrid/9_inhomogeneous_diffusion_porous_catalyst/main.cu Initialization and output folder
+ *
+ */
+//! \cond [Initialization and output folder] \endcond
+int main(int argc, char* argv[])
+{
+ // Initialize library.
+ openfpm_init(&argc, &argv);
+ auto & v_cl = create_vcluster();
+
+ timer t_total;
+ timer t_iterative_diffusion_total;
+ timer t_iteration_wct;
+ timer t_GPU;
+
+ t_total.start();
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ // Set current working directory, define output paths and create folders where output will be saved
+ std::string cwd = get_cwd();
+ const std::string path_output = cwd + "/" + output_name + "/";
+ create_directory_if_not_exist(path_output);
+
+ if(v_cl.rank()==0) std::cout << "Redistancing result will be loaded from " << path_to_redistancing_result << redistancing_filename << std::endl;
+
+ if(v_cl.rank()==0) std::cout << "Outputs will be saved to " << path_output << std::endl;
+
+ //! \cond [Initialization and output folder] \endcond
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ /**
+ * @page 9_inhomogeneous_diffusion_porous_catalyst diffusion_porous_media 3D
+ *
+ * ## Create a dense grid and load redistancing result ## {#e_CaC03_grid}
+ *
+ * \snippet SparseGrid/9_inhomogeneous_diffusion_porous_catalyst/main.cu Grid creation
+ *
+ */
+ //! \cond [Grid creation] \endcond
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ // Create full grid and load redistancing result
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ // Read size of sussman redistancing output grid from csv files
+ openfpm::vector<size_t> v_sz;
+ openfpm::vector<float> v_L;
+
+ size_t m, n;
+ read_csv_to_vector(path_to_size + "/N.csv", v_sz, m, n);
+ read_csv_to_vector(path_to_size + "/L.csv", v_L, m, n);
+
+ const float Lx_low = 0.0;
+ const float Lx_up = v_L.get(x);
+ const float Ly_low = 0.0;
+ const float Ly_up = v_L.get(y);
+ const float Lz_low = 0.0;
+ const float Lz_up = v_L.get(z);
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ size_t sz[dims] = {v_sz.get(x), v_sz.get(y), v_sz.get(z)};
+ Box<dims, float> box({Lx_low, Ly_low, Lz_low}, {Lx_up, Ly_up, Lz_up});
+ Ghost<dims, long int> ghost(1);
+ typedef grid_dist_id<dims, float, props_full> grid_in_type;
+ grid_in_type g_dist(sz, box, ghost);
+
+ g_dist.load(path_to_redistancing_result + "/" + redistancing_filename); // Load SDF
+ //! \cond [Grid creation] \endcond
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ /**
+ * @page 9_inhomogeneous_diffusion_porous_catalyst diffusion_porous_media 3D
+ *
+ * ## Create a sparse grid of fluid phase ## {#e_CaC03_sparse_grid}
+ *
+ * Obtain sparse grid by placing grid points inside the fluid phase only using the signed distance function
+ *
+ * \snippet SparseGrid/9_inhomogeneous_diffusion_porous_catalyst/main.cu Sparse grid creation
+ *
+ */
+ //! \cond [Sparse grid creation] \endcond
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ // Get sparse grid of diffusion domain
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ // Create sparse grid
+ std::cout << "Rank " << v_cl.rank() << " starts creating sparse grid." << std::endl;
+ typedef CartDecomposition<dims,float, CudaMemory, memory_traits_inte, BoxDistribution<dims,float> > Dec;
+ typedef sgrid_dist_id_gpu<dims, float, props_sparse, CudaMemory, Dec> sparse_grid_type;
+ sparse_grid_type g_sparse(sz, box, ghost);
+ g_sparse.setPropNames(prop_names_sparse);
+
+ // Lower and upper bound for SDF indicating fluid phase
+ const float d_low = 0.0, d_up = Lx_up;
+
+ // Alternatively: load sparse grid of diffusion domain from HDF5 file
+// g_sparse.load(path_to_redistancing_result + "/" + redistancing_filename);
+
+ // Obtain sparse grid
+ get_diffusion_domain_sparse_grid<PHI_FULL, PHI_PHASE>(g_dist, g_sparse, d_low, d_up);
+ std::cout << "Rank " << v_cl.rank() << " finished creating sparse grid." << std::endl;
+
+ //! \cond [Sparse grid creation] \endcond
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ /**
+ * @page 9_inhomogeneous_diffusion_porous_catalyst diffusion_porous_media 3D
+ *
+ * ## Defining the inhomogeneous diffusion coefficient ## {#e_CaC03_diff_coeff}
+ *
+ * The inhomogeneous diffusion coefficient is smoothed around the interface by a sigmoidal function that depends on
+ * the signed distance function
+ *
+ * \snippet SparseGrid/9_inhomogeneous_diffusion_porous_catalyst/main.cu Diffusion coefficient
+ *
+ */
+ //! \cond [Diffusion coefficient] \endcond
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ // Initialize properties on sparse grid
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ // Diffusion coefficient
+ const float min_porosity = 0.0;
+ const float max_porosity = 1.0;
+
+ float D_molecular = 1.0;
+ float D_min = D_molecular * min_porosity;
+ float D_max = D_molecular * max_porosity;
+
+ float min_sdf_calcite = get_min_val<PHI_PHASE>(g_sparse);
+
+ float x_stretch = 4.0 * g_sparse.spacing(x);
+ float x_shift = min_sdf_calcite * x_stretch;
+ float y_min = D_min;
+ float y_max = D_max - D_min;
+
+ const float c0 = 10.0;
+
+ auto dom = g_sparse.getDomainIterator();
+ while(dom.isNext())
+ {
+ auto key = dom.get();
+
+ Point<dims, float> coords = g_sparse.getPos(key);
+ if(coords.get(x) <= 0.05 * Lx_up)
+ {
+ g_sparse.template insertFlush<CONC_N>(key) = c0;
+ }
+ else g_sparse.template insertFlush<CONC_N>(key) = 0.0;
+
+ // Compute diffusion coefficient and store on grid
+ g_sparse.template insertFlush<DIFFUSION_COEFFICIENT>(key) =
+ get_smooth_sigmoidal(
+ g_sparse.template getProp<PHI_PHASE>(key) - min_sdf_calcite,
+ x_shift,
+ x_stretch,
+ y_min,
+ y_max);
+ ++dom;
+ }
+
+ //! \cond [Diffusion coefficient] \endcond
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ /**
+ * @page 9_inhomogeneous_diffusion_porous_catalyst diffusion_porous_media 3D
+ *
+ * ## Defining the functor for the inhomogeneous diffusion ## {#e_CaC03_functor}
+ *
+ * The functor contains the forward time central space discretization of the inhomogeneous diffusion equation.
+ * The stencil size is 1.
+ * It will be passed to the GPU and convolved with the sparse grid.
+ *
+ * \snippet SparseGrid/9_inhomogeneous_diffusion_porous_catalyst/main.cu Functor
+ *
+ */
+ //! \cond [Functor] \endcond
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ // Defining the functor that will be passed to the GPU to simulate reaction-diffusion
+ // GetCpBlockType<typename SGridGpu, unsigned int prp, unsigned int stencil_size>
+ typedef typename GetCpBlockType<decltype(g_sparse),0,1>::type CpBlockType;
+
+ const float dx = g_sparse.spacing(x), dy = g_sparse.spacing(y), dz = g_sparse.spacing(z);
+ // Stability criterion for FTCS
+ const float dt = diffusion_time_step(g_sparse, D_max);
+
+ std::cout << "dx = " << dx << "dy = " << dy << "dz = " << dz << std::endl;
+ std::cout << "dt = " << dt << std::endl;
+
+
+ auto func_inhomogDiffusion = [dx, dy, dz, dt, d_low] __device__ (
+ float & u_out, // concentration output
+ float & D_out, // diffusion coefficient output (dummy)
+ float & phi_out, // sdf of domain output (dummy)
+ CpBlockType & u, // concentration input
+ CpBlockType & D, // diffusion coefficient input (to get also the half step diffusion coefficients)
+ CpBlockType & phi, // sdf of domain (for boundary conditions)
+ auto & block,
+ int offset,
+ int i,
+ int j,
+ int k
+ )
+ {
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ // Stencil
+ // Concentration
+ float u_c = u(i, j, k);
+ float u_px = u(i+1, j, k);
+ float u_mx = u(i-1, j, k);
+ float u_py = u(i, j+1, k);
+ float u_my = u(i, j-1, k);
+ float u_pz = u(i, j, k+1);
+ float u_mz = u(i, j, k-1);
+
+ // Signed distance function
+ float phi_c = phi(i, j, k);
+ float phi_px = phi(i+1, j, k);
+ float phi_mx = phi(i-1, j, k);
+ float phi_py = phi(i, j+1, k);
+ float phi_my = phi(i, j-1, k);
+ float phi_pz = phi(i, j, k+1);
+ float phi_mz = phi(i, j, k-1);
+
+ // Diffusion coefficient
+ float D_c = D(i, j, k);
+ float D_px = D(i+1, j, k);
+ float D_mx = D(i-1, j, k);
+ float D_py = D(i, j+1, k);
+ float D_my = D(i, j-1, k);
+ float D_pz = D(i, j, k+1);
+ float D_mz = D(i, j, k-1);
+
+ // Impose no-flux boundary conditions
+ if (phi_px <= d_low + std::numeric_limits<float>::epsilon()) {u_px = u_c; D_px = D_c;}
+ if (phi_mx <= d_low + std::numeric_limits<float>::epsilon()) {u_mx = u_c; D_mx = D_c;}
+ if (phi_py <= d_low + std::numeric_limits<float>::epsilon()) {u_py = u_c; D_py = D_c;}
+ if (phi_my <= d_low + std::numeric_limits<float>::epsilon()) {u_my = u_c; D_my = D_c;}
+ if (phi_pz <= d_low + std::numeric_limits<float>::epsilon()) {u_pz = u_c; D_pz = D_c;}
+ if (phi_mz <= d_low + std::numeric_limits<float>::epsilon()) {u_mz = u_c; D_mz = D_c;}
+
+ // Interpolate diffusion constant of half-step neighboring points
+ float D_half_px = (D_c + D_px)/2.0;
+ float D_half_mx = (D_c + D_mx)/2.0;
+ float D_half_py = (D_c + D_py)/2.0;
+ float D_half_my = (D_c + D_my)/2.0;
+ float D_half_pz = (D_c + D_pz)/2.0;
+ float D_half_mz = (D_c + D_mz)/2.0;
+
+ // Compute concentration of next time point
+ u_out = u_c + dt *
+ (1/(dx*dx) * (D_half_px * (u_px - u_c) - D_half_mx * (u_c - u_mx)) +
+ 1/(dy*dy) * (D_half_py * (u_py - u_c) - D_half_my * (u_c - u_my)) +
+ 1/(dz*dz) * (D_half_pz * (u_pz - u_c) - D_half_mz * (u_c - u_mz)));
+
+ // dummy outs
+ D_out = D_c;
+ phi_out = phi_c;
+ };
+
+
+ //! \cond [Functor] \endcond
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ /**
+ * @page 9_inhomogeneous_diffusion_porous_catalyst diffusion_porous_media 3D
+ *
+ * ## Run inhomogeneous diffusion ## {#e_CaC03_run}
+ *
+ *
+ * \snippet SparseGrid/9_inhomogeneous_diffusion_porous_catalyst/main.cu Run
+ *
+ */
+ //! \cond [Run] \endcond
+ // Create text file to which wall clock time for iteration incl. ghost communication will be saved
+ std::string path_time_filewct = path_output + "/time_wct_incl_ghostCommunication" + std::to_string(v_cl.rank()) + ".txt";
+ std::ofstream out_wct_file(path_time_filewct, std::ios_base::app);
+
+ // Create text file to which GPU time will be saved
+ std::string path_time_filegpu = path_output + "/time_GPU" + std::to_string(v_cl.rank()) + ".txt";
+ std::ofstream out_GPUtime_file(path_time_filegpu, std::ios_base::app);
+
+ // Iterative diffusion
+ const size_t iterations = 103;
+ const size_t number_write_outs = 10;
+ const size_t interval_write = iterations / number_write_outs; // Find interval for writing to reach
+ size_t iter = 0;
+ float t = 0;
+
+ // Copy from host to GPU for simulation
+ g_sparse.template hostToDevice<CONC_N, CONC_NPLUS1, DIFFUSION_COEFFICIENT, PHI_PHASE>();
+ t_iterative_diffusion_total.start();
+ while(iter <= iterations)
+ {
+ if (iter % 2 == 0)
+ {
+ t_iteration_wct.start();
+ g_sparse.template ghost_get<PHI_PHASE, CONC_N, DIFFUSION_COEFFICIENT>(RUN_ON_DEVICE);
+ t_GPU.startGPU();
+ g_sparse.template conv3_b<
+ CONC_N,
+ DIFFUSION_COEFFICIENT,
+ PHI_PHASE,
+ CONC_NPLUS1,
+ DIFFUSION_COEFFICIENT,
+ PHI_PHASE, 1>
+ ({0, 0, 0},
+ {(long int) sz[x]-1, (long int) sz[y]-1, (long int) sz[z]-1},
+ func_inhomogDiffusion);
+ t_GPU.stopGPU();
+ t_iteration_wct.stop();
+
+ // Write out time to text-file
+ out_wct_file << t_iteration_wct.getwct() << ",";
+ out_GPUtime_file << t_GPU.getwctGPU() << ",";
+ }
+ else
+ {
+ t_iteration_wct.start();
+ g_sparse.template ghost_get<PHI_PHASE, CONC_NPLUS1, DIFFUSION_COEFFICIENT>(RUN_ON_DEVICE);
+ t_GPU.startGPU();
+ g_sparse.template conv3_b<
+ CONC_NPLUS1,
+ DIFFUSION_COEFFICIENT,
+ PHI_PHASE,
+ CONC_N,
+ DIFFUSION_COEFFICIENT,
+ PHI_PHASE, 1>
+ ({0, 0, 0},
+ {(long int) sz[x]-1, (long int) sz[y]-1, (long int) sz[z]-1},
+ func_inhomogDiffusion);
+ t_GPU.stopGPU();
+ t_iteration_wct.stop();
+ // Write out time to text-file
+ out_wct_file << t_iteration_wct.getwct() << ",";
+ out_GPUtime_file << t_GPU.getwctGPU() << ",";
+ }
+
+
+ if (iter % interval_write == 0)
+ {
+ // Copy from GPU to host for writing
+ g_sparse.template deviceToHost<CONC_N, CONC_NPLUS1, DIFFUSION_COEFFICIENT, PHI_PHASE>();
+
+ // Write g_sparse to vtk
+ g_sparse.write_frame(path_output + "g_sparse_diffuse", iter, FORMAT_BINARY);
+
+ // Write g_sparse to hdf5
+ g_sparse.save(path_output + "g_sparse_diffuse_" + std::to_string(iter) + ".hdf5");
+
+ if (iter % 2 == 0)
+ {
+ monitor_total_concentration<CONC_N>(g_sparse, t, iter, path_output,"total_conc.csv");
+ }
+ else
+ {
+ monitor_total_concentration<CONC_NPLUS1>(g_sparse, t, iter, path_output,"total_conc.csv");
+ }
+ }
+
+
+ t += dt;
+ iter += 1;
+ }
+ t_iterative_diffusion_total.stop();
+ out_GPUtime_file.close();
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ t_total.stop();
+ std::cout << "Rank " << v_cl.rank() << " total time for the whole programm : " << t_total.getwct() << std::endl;
+
+ //! \cond [Run] \endcond
+
+ /**
+ * @page 9_inhomogeneous_diffusion_porous_catalyst diffusion_porous_media 3D
+ *
+ * ## Terminate ## {#e_CaC03_terminate}
+ *
+ * We end with terminating OpenFPM.
+ *
+ * \snippet SparseGrid/9_inhomogeneous_diffusion_porous_catalyst/main.cu Terminate
+ */
+ //! \cond [Terminate] \endcond
+ openfpm_finalize(); // Finalize openFPM library
+ return 0;
+}
+//! \cond [Terminate] \endcond
+
+/**
+ * @page 9_inhomogeneous_diffusion_porous_catalyst diffusion_porous_media 3D
+ *
+ * ## Full code ## {#e_CaC03_full}
+ *
+ * @include SparseGrid/9_inhomogeneous_diffusion_porous_catalyst/main.cu
+ */