Adding doxymentation for example for heat conduction in reticulate porous ceramics.

example/SparseGrid/10_heat_conduction_reticulate_porous_ceramics/config.cfg

+[pack]
+files = main.cu Makefile
\ No newline at end of file

example/SparseGrid/10_heat_conduction_reticulate_porous_ceramics/include/DiffusionSpace_sparseGrid.hpp

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

example/SparseGrid/10_heat_conduction_reticulate_porous_ceramics/include/HelpFunctions_diffusion.hpp

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

example/SparseGrid/10_heat_conduction_reticulate_porous_ceramics/main.cu

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

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.
  *
  */