Adding sparse grid example code and doxygen documentation for inhomogeneous...

Adding sparse grid example code and doxygen documentation for inhomogeneous diffusion in the fluid phase of a CaCO3 particle packed bed.

example/SparseGrid/9_inhomogeneous_diffusion_porous_catalyst_CaCO3/config.cfg

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

example/SparseGrid/9_inhomogeneous_diffusion_porous_catalyst_CaCO3/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/9_inhomogeneous_diffusion_porous_catalyst_CaCO3/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/9_inhomogeneous_diffusion_porous_catalyst_CaCO3/main.cu

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