diff --git a/example/SparseGrid/10_heat_conduction_reticulate_porous_ceramics/config.cfg b/example/SparseGrid/10_heat_conduction_reticulate_porous_ceramics/config.cfg
new file mode 100755
index 0000000000000000000000000000000000000000..067b942df5494321fd7d5440df1b7dff32b4f11a
--- /dev/null
+++ b/example/SparseGrid/10_heat_conduction_reticulate_porous_ceramics/config.cfg
@@ -0,0 +1,2 @@
+[pack]
+files = main.cu Makefile
\ No newline at end of file
diff --git a/example/SparseGrid/10_heat_conduction_reticulate_porous_ceramics/include/DiffusionSpace_sparseGrid.hpp b/example/SparseGrid/10_heat_conduction_reticulate_porous_ceramics/include/DiffusionSpace_sparseGrid.hpp
new file mode 100755
index 0000000000000000000000000000000000000000..d613e0a43d3d205f0cf966f20052d1dd8d3b10d9
--- /dev/null
+++ b/example/SparseGrid/10_heat_conduction_reticulate_porous_ceramics/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/10_heat_conduction_reticulate_porous_ceramics/include/HelpFunctions_diffusion.hpp b/example/SparseGrid/10_heat_conduction_reticulate_porous_ceramics/include/HelpFunctions_diffusion.hpp
new file mode 100755
index 0000000000000000000000000000000000000000..7e3069c8549f142ae2f3b31a730533d196530365
--- /dev/null
+++ b/example/SparseGrid/10_heat_conduction_reticulate_porous_ceramics/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/10_heat_conduction_reticulate_porous_ceramics/main.cu b/example/SparseGrid/10_heat_conduction_reticulate_porous_ceramics/main.cu
new file mode 100644
index 0000000000000000000000000000000000000000..96bc781fd46b8c42ce93b6a0a56be8659ca837f3
--- /dev/null
+++ b/example/SparseGrid/10_heat_conduction_reticulate_porous_ceramics/main.cu
@@ -0,0 +1,386 @@
+//
+// Created by jstark on 2023-05-05.
+//
+/*!
+ *
+ * \page 10_heat_conduction_RPC Heat conduction in reticulate porous ceramics using sparse grids on GPU
+ *
+ * [TOC]
+ *
+ * # Solving heat conduction in the image-based geometry of reticulate porous ceramics # {#e10_heat_conduction_RPC_gpu}
+ *
+ * In this example, we simulate heat conduction in the solid phase of reticulate porous ceramics with heat dissipation at the surface. 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 solid 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://doi.org/10.1111/j.1551-2916.2008.02308.x">J. Petrasch et al., "Tomography-based multiscale analyses of the 3D geometrical morphology of reticulated porous ceramics", Journal of the American Ceramic Society (2008)</a>
+ * and <a href="https://doi.org/10.1115/1.4000226">S. Haussener et al., "Tomography-based heat and mass transfer characterization of reticu- late porous ceramics for high-temperature processing", Journal of Heat Transfer (2010)</a>).
+ *
+ *
+ * For image based reconstruction and redistancing see @ref example_sussman_images_3D.
+ *
+ *
+ * First, we initialize and generate sparse grid of solid phase. In this example, an indicator funtion (color field)
+ * would also be sufficient for defining the diffusion domain as we do not scale any parameter by their distance to the
+ * surface here, as opposed to the inhomogeneous diffusion in the CaCO\f$_3\f$ fluid phase,
+ * see @ref 9_inhomogeneous_diffusion_porous_catalyst.
+ * We anyways load the SDF as this gives us more flexibility when we want to extend the code towards applications
+ * for which we need the distance to the surface.
+ *
+ * \snippet SparseGrid/10_heat_conduction_reticulate_porous_ceramics/main.cu initialize
+ *
+ * Define the functor for heat conduction
+ * \snippet SparseGrid/10_heat_conduction_reticulate_porous_ceramics/main.cu functor
+ *
+ * Iterative heat conduction
+ * \snippet SparseGrid/10_heat_conduction_reticulate_porous_ceramics/main.cu iteration
+ *
+ */
+
+//! \cond [initialize] \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 "level_set/redistancing_Sussman/HelpFunctionsForGrid.hpp" // For the ghost initialization
+#include "RemoveLines.hpp" // For removing thin (diagonal or straight) lines
+
+#include "FiniteDifference/FD_simple.hpp"
+
+#include "DiffusionSpace_sparseGrid.hpp"
+#include "HelpFunctions_diffusion.hpp"
+
+#include "Decomposition/Distribution/BoxDistribution.hpp"
+
+const size_t dims = 3;
+
+// 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 U_N = 1;
+const size_t U_NPLUS1 = 2;
+const size_t K_SINK = 3;
+
+
+typedef aggregate<float, float, float, float> props_sparse;
+const openfpm::vector<std::string> prop_names_sparse = {"PHI_PHASE",
+ "U_N",
+ "U_NPLUS1",
+ "K_SINK"};
+
+
+// Space indices
+constexpr size_t x = 0, y = 1, z = 2;
+
+// input
+const std::string path_to_redistancing_result =
+ "/MY_PATH/porous_ceramics/sussman_with_cuda/build/output_sussman_sparse_grid_porousCeramics_1216x1016x941/";
+
+
+const std::string redistancing_filename = "sparseGrid_initial.hdf5";
+
+const std::string path_to_size = path_to_redistancing_result;
+
+const std::string output_name = "output_heat_conduction_porous_ceramics_1216x1016x941_sparse_grid";
+
+
+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;
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////
+ // 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);
+
+ std::cout << "Rank " << v_cl.rank() << " starts loading redistancing result..." << std::endl;
+ g_dist.load(path_to_redistancing_result + "/" + redistancing_filename); // Load SDF
+ std::cout << "Rank " << v_cl.rank() << " finished loading redistancing result." << std::endl;
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ // Get sparse grid of heat conduction 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> sparse_grid_type;
+ sparse_grid_type g_sparse(sz, box, ghost);
+ g_sparse.setPropNames(prop_names_sparse);
+
+
+ const float d_low = 0.0, d_up = Lx_up;
+ get_diffusion_domain_sparse_grid<PHI_FULL, PHI_PHASE>(g_dist, g_sparse, d_low, d_up);
+ // g_sparse.load(path_to_redistancing_result + "/" + redistancing_filename);
+ std::cout << "Rank " << v_cl.rank() << " finished creating sparse grid." << std::endl;
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ // Define parameters for heat conduction
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ const float u0 = 1.0; // Initial heat from one side
+ const float D = 0.1; // Heat diffusivity
+ const float sink_rate = 0.1; // Rate of heat loss at solid-air interface
+ std::cout << "Diffusion coefficient in s/(mm^2): " << D << std::endl;
+ std::cout << "Sink rate in 1/s: " << sink_rate << std::endl;
+
+ // Set initial condition: initial heat source is a sphere of radius R at the domain center
+ const float center[dims] = {(Lx_up+Lx_low)/(float)2,
+ (Ly_up+Ly_low)/(float)2,
+ (Lz_up+Lz_low)/(float)2};
+ const float R = Lz_up / 4.0;
+
+ auto dom = g_sparse.getDomainIterator();
+ while(dom.isNext())
+ {
+ auto key = dom.get();
+
+ // Initial condition: heat comes from sphere at center
+ Point<dims, float> coords = g_sparse.getPos(key);
+
+ if(((coords.get(x) - center[x]) * (coords.get(x) - center[x])
+ + (coords.get(y) - center[y]) * (coords.get(y) - center[y])
+ + (coords.get(z) - center[z]) * (coords.get(z) - center[z])) <= R * R)
+ {
+ g_sparse.template insertFlush<U_N>(key) = u0;
+ }
+ else g_sparse.template insertFlush<U_N>(key) = 0.0;
+
+ ++dom;
+ }
+
+ //! \cond [initialize] \endcond
+
+ //! \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 : dt = 1/4 * 1/(1/dx^2 + 1/dy^2)
+ const float dt = diffusion_time_step(g_sparse, D);
+ std::cout << "dx, dy, dz = " << g_sparse.spacing(x) << "," << g_sparse.spacing(y) << "," << g_sparse.spacing(z) << std::endl;
+ std::cout << "dt = 1/(4D) * 1/(1/dx^2 + 1/dy^2 + 1/dz^2) = " << dt << std::endl;
+
+ auto func_heatDiffusion_withSink = [dx, dy, dz, dt, d_low, D, sink_rate] __device__ (
+ float & u_out, // concentration output
+ float & phi_out, // sdf of domain output (dummy)
+ CpBlockType & u, // concentration input
+ 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);
+
+ // Get sink term
+ float k_sink = 0; // Sink equal to zero in bulk and equal to sink_rate only at the surfaces
+
+ // Impose no-flux boundary conditions and sink term at the surface (i.e., when neighbor is outside)
+ if (phi_px <= d_low + std::numeric_limits<float>::epsilon()) {u_px = u_c; k_sink = sink_rate;}
+ if (phi_mx <= d_low + std::numeric_limits<float>::epsilon()) {u_mx = u_c; k_sink = sink_rate;}
+ if (phi_py <= d_low + std::numeric_limits<float>::epsilon()) {u_py = u_c; k_sink = sink_rate;}
+ if (phi_my <= d_low + std::numeric_limits<float>::epsilon()) {u_my = u_c; k_sink = sink_rate;}
+ if (phi_pz <= d_low + std::numeric_limits<float>::epsilon()) {u_pz = u_c; k_sink = sink_rate;}
+ if (phi_mz <= d_low + std::numeric_limits<float>::epsilon()) {u_mz = u_c; k_sink = sink_rate;}
+
+ block.template get<K_SINK>()[offset] = k_sink; // Just for checking in paraview
+
+ // Compute concentration of next time point
+ u_out = u_c + dt * ( D * ((u_px + u_mx - 2 * u_c)/(dx * dx)
+ + (u_py + u_my - 2 * u_c)/(dy * dy)
+ + (u_pz + u_mz - 2 * u_c)/(dz * dz))
+ - k_sink * u_c);
+
+
+ // dummy outs
+ phi_out = phi_c;
+ };
+
+ //! \cond [functor] \endcond
+
+
+ //! \cond [iteration] \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 = 1e6;
+ const size_t iterations = 103;
+ const size_t number_write_outs = 1;
+ 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<U_N, U_NPLUS1, K_SINK, 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, U_N, K_SINK>(RUN_ON_DEVICE);
+ t_GPU.startGPU();
+ g_sparse.template conv2_b<
+ U_N,
+ PHI_PHASE,
+ U_NPLUS1,
+ PHI_PHASE, 1>
+ ({0, 0, 0},
+ {(long int) sz[x]-1, (long int) sz[y]-1, (long int) sz[z]-1},
+ func_heatDiffusion_withSink);
+ 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() << ",";
+
+ t_iteration_wct.reset();
+ }
+ else
+ {
+ t_iteration_wct.start();
+ g_sparse.template ghost_get<PHI_PHASE, U_NPLUS1, K_SINK>(RUN_ON_DEVICE);
+ t_GPU.startGPU();
+ g_sparse.template conv2_b<
+ U_NPLUS1,
+ PHI_PHASE,
+ U_N,
+ PHI_PHASE, 1>
+ ({0, 0, 0},
+ {(long int) sz[x]-1, (long int) sz[y]-1, (long int) sz[z]-1},
+ func_heatDiffusion_withSink);
+ 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() << ",";
+
+ t_iteration_wct.reset();
+ }
+#if 0
+
+ if (iter % interval_write == 0)
+ {
+ // Copy from GPU to host for writing
+ g_sparse.template deviceToHost<U_N, U_NPLUS1, K_SINK, PHI_PHASE>();
+
+ // Write g_sparse to vtk
+ g_sparse.write_frame(path_output + "g_sparse_diffuse", iter, FORMAT_BINARY);
+
+ if (iter % 2 == 0)
+ {
+ monitor_total_concentration<U_N>(g_sparse, t, iter, path_output,"total_conc.csv");
+ }
+ else
+ {
+ monitor_total_concentration<U_NPLUS1>(g_sparse, t, iter, path_output,"total_conc.csv");
+ }
+ }
+#endif
+ t += dt;
+ iter += 1;
+ }
+ t_iterative_diffusion_total.stop();
+ std::cout << "Rank " << v_cl.rank() << " total time for iterative diffusion incl. write outs: " << t_iterative_diffusion_total.getwct() << std::endl;
+ std::cout << "Rank " << v_cl.rank() << " finished diffusion." << std::endl;
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ t_total.stop();
+ std::cout << "Rank " << v_cl.rank() << " total time for the whole programm : " << t_total.getwct() << std::endl;
+
+
+ openfpm_finalize();
+ return 0;
+}
+
+//! \cond [iteration] \endcond
diff --git a/example/SparseGrid/9_inhomogeneous_diffusion_porous_catalyst_CaCO3/main.cu b/example/SparseGrid/9_inhomogeneous_diffusion_porous_catalyst_CaCO3/main.cu
index b607f52d5cbd5e0f7022882eead120535a5746bc..c0d1136e580ae228069a7ea7e28d87e49853b2e4 100755
--- a/example/SparseGrid/9_inhomogeneous_diffusion_porous_catalyst_CaCO3/main.cu
+++ b/example/SparseGrid/9_inhomogeneous_diffusion_porous_catalyst_CaCO3/main.cu
@@ -4,15 +4,15 @@
/**
* @file 9_inhomogeneous_diffusion_porous_catalyst/main.cu
* @page 9_inhomogeneous_diffusion_porous_catalyst diffusion_porous_media 3D
- *
+ * [TOC]
* # 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>).
+ * 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.
+ * For image based reconstruction and redistancing, see @ref example_sussman_images_3D.
*
*/