src/BoundaryConditions/tests/libForTesting/HelpFunctionsForTestingNBCs.hpp

+//
+// Created by jstark on 07.06.21.
+//
+
+#ifndef OPENFPM_NUMERICS_HELPFUNCTIONSFORTESTINGNBCS_HPP
+#define OPENFPM_NUMERICS_HELPFUNCTIONSFORTESTINGNBCS_HPP
+
+#include <iostream>
+#include "level_set/redistancing_Sussman/NarrowBand.hpp"
+
+template <size_t Phi_SDF_grid, size_t Phi_SDF, size_t dPhi, size_t dPhi_magn, typename grid_in_type, typename
+vd_in_type>
+void get_diffusion_domain(grid_in_type & g_dist, vd_in_type & vd, const double b_low, const double b_up)
+{
+	const double EPSILON = std::numeric_limits<double>::epsilon();
+	const double _b_low = b_low + EPSILON;
+	const double _b_up  = b_up - EPSILON;
+	NarrowBand<grid_in_type> narrowBand(g_dist, _b_low, _b_up); // Instantiation of NarrowBand class
+	narrowBand.template get_narrow_band<Phi_SDF_grid, Phi_SDF, dPhi, dPhi_magn>(g_dist, vd);
+}
+
+/**@brief Fill pid_source with IDs of particles whose SDF is >= b_low and <= b_up.
+ *
+ * @param vd Input particle vector_dist of type vd_type.
+ * @param b_low Lower bound for the SDF value a particle must have in order to be a source particle.
+ * @param b_up Upper bound for the SDF value a particle must have in order to be a source particle.
+ */
+template <size_t Phi_SDF, typename vd_type>
+void get_source_particle_ids(vd_type & vd, openfpm::vector<vect_dist_key_dx> & keys_source, const double b_low, const
+double
+b_up)
+{
+	const double EPSILON = std::numeric_limits<double>::epsilon();
+	auto dom = vd.getDomainIterator();
+	while(dom.isNext())
+	{
+		auto key = dom.get();
+		if (vd.template getProp<Phi_SDF>(key) >= b_low + EPSILON && vd.template getProp<Phi_SDF>(key) <= b_up + EPSILON)
+		{
+			keys_source.add(key);
+		}
+		++dom;
+	}
+}
+
+
+template <size_t R, size_t U, typename vd_type>
+void get_IC_Eq28(vd_type &vd, const double a)
+{
+	auto dom = vd.getDomainIterator();
+	while(dom.isNext())
+	{
+		auto key = dom.get();
+		double r = vd.template getProp<R>(key);
+		if (abs(r) < a) vd.template getProp<U>(key) = cos(M_PI * r / (2*a)) * cos(M_PI * r / (2*a));
+		else vd.template getProp<U>(key) = 0;
+		++dom;
+	}
+}
+
+template <size_t R, size_t U, typename vd_type>
+void get_IC_Eq28(vd_type &vd, const double a, const openfpm::vector<aggregate<int>> & pids)
+{
+	for (int i = 0; i < pids.size(); i++)
+	{
+		auto key = pids.get<0>(i);
+		double r = vd.template getProp<R>(key);
+		if (abs(r) < a) vd.template getProp<U>(key) = cos(M_PI * r / (2*a)) * cos(M_PI * r / (2*a));
+		else vd.template getProp<U>(key) = 0;
+	}
+}
+
+template <size_t U, typename vd_type>
+void get_FS_Eq29(vd_type &vd, const double a, const double R_disk)
+{
+	const double R = R_disk;
+	double A = 0;
+	if (a < R) A = M_PI * a*a/2 - 2*R*R/M_PI;
+	if (a >= R) A = M_PI * R*R/2 + a*a*( cos(M_PI * R / a) -1 ) / M_PI + a * R * sin(M_PI * R / a);
+	
+	auto u = getV<U>(vd);
+	u = A / (M_PI * R*R);
+	
+	std::cout << "Analytical steady state concentration = " << A / (M_PI * R*R) << std::endl;
+}
+
+
+#endif //OPENFPM_NUMERICS_HELPFUNCTIONSFORTESTINGNBCS_HPP

src/BoundaryConditions/tests/method_of_images_cylinder_unit_test.cpp

+//
+// Created by Abhinav Singh & Justina Stark on 10.06.21.
+//
+#define BOOST_TEST_DYN_LINK
+#include <iostream>
+#include "config.h"
+#include "util/common.hpp"
+#include "util/util_debug.hpp"
+#include <boost/test/unit_test.hpp>
+
+#include "Vector/vector_dist_subset.hpp"
+#include "Operators/Vector/vector_dist_operators.hpp"
+
+// For the Neumann BCs (Method of images)
+#include "BoundaryConditions/MethodOfImages.hpp"
+
+
+constexpr int x = 0;
+constexpr int y = 1;
+constexpr int z = 2;
+constexpr size_t CONCENTRATION = 0;
+constexpr size_t NORMAL        = 1;
+constexpr size_t IS_SOURCE     = 2;
+constexpr size_t subset_id_real   = 0;
+constexpr size_t subset_id_mirror = 1;
+
+
+
+BOOST_AUTO_TEST_SUITE(MethodOfImagesTestSuite)
+	BOOST_AUTO_TEST_CASE(mirror_cylinder_base_test) {
+		const size_t sz[3] = {18, 18, 5};
+		double boxsize = 20.0;
+		double spacing_p = 10.0 / (sz[x]-1);
+		double spacing_np = 10.0 / (sz[x]-1);
+		Box<3, double> box({-2.0, -2.0, 0}, {boxsize, boxsize, (sz[z])*spacing_p});
+		size_t bc[3] = {NON_PERIODIC, NON_PERIODIC, PERIODIC};
+//		Ghost<3, double> ghost(2.9*spacing_np+spacing_np/8.0);
+		const double ghost_layer_thickness = 2.9*spacing_np+spacing_np/8.0;
+//		const double ghost_layer_thickness = 5.0 - 3*spacing_np/4;
+		Ghost<3, double> ghost(ghost_layer_thickness);
+		
+		typedef vector_dist_ws<3, double, aggregate<double, VectorS<3,double>, int>> vd_type;
+		vd_type Particles(0,box,bc,ghost);
+		Particles.setPropNames({"Concentration", "Normal", "is_source"});
+		
+		//Grid Like Particles
+		auto it = Particles.getGridIterator(sz);
+		while (it.isNext())
+		{
+			auto key = it.get();
+			double x_coord = key.get(x) * spacing_np;
+			double y_coord = key.get(y) * spacing_np;
+			double z_coord = key.get(z) * spacing_p;
+			
+			if(sqrt((5.0-x_coord)*(5.0-x_coord)+(5.0-y_coord)*(5.0-y_coord))<5.0-7*spacing_np/6.0)
+			{
+				Particles.add();
+				if (key.get(x)==sz[x]-1)
+				{
+					Particles.getLastPos()[x] = boxsize-1e-6;
+				}
+				else
+				{
+					Particles.getLastPos()[x] = x_coord;
+				}
+				
+				if (key.get(y)==sz[y]-1)
+				{
+					Particles.getLastPos()[y] = boxsize-1e-6;
+				}
+				else
+				{
+					Particles.getLastPos()[y] = y_coord;
+				}
+				Particles.getLastPos()[z] = z_coord;
+				Particles.getLastProp<NORMAL>()[x] = 0;
+				Particles.getLastProp<NORMAL>()[y] = 0;
+				Particles.getLastProp<NORMAL>()[z] = 0;
+				Particles.getLastProp<CONCENTRATION>() = x_coord+y_coord+z_coord;
+				
+				Particles.getLastProp<IS_SOURCE>() = 0;
+			}
+			++it;
+		}
+		
+		//Adding a layer to represent geometry better. (Mirroring should only have this layer as it is slightly inside radius and the gridlike particles have been set to 0 normal)
+		int n_b=int(sz[0])*5;
+		double radius = 5.0 - 3*spacing_np/4;
+		//double Golden_angle=M_PI * (3.0 - sqrt(5.0));
+		double Golden_angle=2.0*M_PI/double(n_b);
+		int number_of_border_particles = 0;
+		for (int j=0;j<int(sz[z]);j++)
+		{
+			for(int i=1;i<=n_b;i++)
+			{
+				double Golden_theta = Golden_angle * i;
+				double x_coord = 5.0+cos(Golden_theta) * radius;
+				double y_coord = 5.0+sin(Golden_theta) * radius;
+				double z_coord = j*spacing_p;
+				Particles.add();
+				Particles.getLastPos()[x] = x_coord;
+				Particles.getLastPos()[y] = y_coord;
+				Particles.getLastPos()[z] = z_coord;
+
+				Particles.getLastProp<NORMAL>()[x] = (x_coord-5.0)/sqrt((x_coord-5.0)*(x_coord-5.0)+(y_coord-5.0)*(y_coord-5.0));
+				Particles.getLastProp<NORMAL>()[y] = (y_coord-5.0)/sqrt((x_coord-5.0)*(x_coord-5.0)+(y_coord-5.0)*(y_coord-5.0));
+				Particles.getLastProp<NORMAL>()[z] = 0.0;
+				Particles.getLastProp<CONCENTRATION>() = x_coord+y_coord+z_coord;
+				
+				Particles.getLastProp<IS_SOURCE>() = 1;
+
+				number_of_border_particles++;
+			}
+		}
+		auto &v_cl = create_vcluster();
+		v_cl.sum(number_of_border_particles);
+		v_cl.execute();
+		
+		if (v_cl.rank() == 0) std::cout << "Number of particles with surface normal = " << number_of_border_particles
+		<< std::endl;
+		Particles.map();
+		Particles.ghost_get<NORMAL,IS_SOURCE>();
+		//We write the particles to check if the initialization is correct.
+//		Particles.write("Init");
+		
+		//Here Mirror Particle and do method of Images and check if it matches  property 0 mirroring (x+y+z of the mirror).
+		
+		openfpm::vector<vect_dist_key_dx> keys_source;
+		
+		auto dom = Particles.getDomainIterator();
+		while(dom.isNext())
+		{
+			auto key = dom.get();
+			if (Particles.template getProp<IS_SOURCE>(key) == 1)
+			{
+				keys_source.add(key);
+			}
+			++dom;
+		}
+		size_t number_of_source_particles = keys_source.size();
+		v_cl.sum(number_of_source_particles);
+		v_cl.execute();
+		size_t number_of_real_particle_no_ghost = Particles.size_local();
+		size_t number_of_real_particle_with_ghost = Particles.size_local_with_ghost();
+		
+		/*
+		if (v_cl.rank() == 0)
+		{
+			std::cout << "number_of_source_particles = " << number_of_source_particles << std::endl;
+			std::cout << "number_of_real_particle_no_ghost = " << number_of_real_particle_no_ghost << std::endl;
+			std::cout << "number_of_real_particle_with_ghost before mirroring = " << number_of_real_particle_with_ghost << std::endl;
+		}
+		*/
+		
+		
+		// Apply Method of images to impose noflux Neumann Boundary Conditions
+		MethodOfImages<NORMAL, vd_type> NBCs(Particles, keys_source, subset_id_real, subset_id_mirror);
+		NBCs.get_mirror_particles(Particles);
+		NBCs.apply_noflux<CONCENTRATION>(Particles);
+		
+		size_t number_of_mirror_particles = Particles.size_local() - number_of_real_particle_no_ghost;
+		v_cl.sum(number_of_mirror_particles);
+		v_cl.execute();
+		
+		if (v_cl.rank() == 0) std::cout << "Number of mirror particles = " << number_of_mirror_particles << std::endl;
+		
+		/*
+		if (v_cl.rank() == 0)
+		{
+			std::cout << "number_of_real_particle_with_ghost + mirror particles = " << Particles.size_local_with_ghost() <<
+					std::endl;
+			std::cout << "Total number of particles expected after mirroring = " << number_of_real_particle_with_ghost +
+					keys_source.size() << std::endl;
+		}
+		*/
+		
+		Particles.write("Cylinder_with_mirror_particles");
+		
+		BOOST_CHECK(number_of_source_particles == number_of_border_particles);
+		BOOST_CHECK(number_of_mirror_particles == number_of_source_particles);
+		
+		for (int i = 0; i < NBCs.key_map_source_mirror.size(); ++i)
+		{
+			openfpm::vector<size_t> row = NBCs.key_map_source_mirror.get(i);
+			vect_dist_key_dx key_source, key_mirror;
+			key_source.setKey(row.get(0));
+			key_mirror.setKey(row.get(1));
+			BOOST_CHECK(Particles.template getProp<CONCENTRATION>(key_mirror) == Particles.template getProp<CONCENTRATION>(key_source));
+		}
+		
+	}
+BOOST_AUTO_TEST_SUITE_END()
+

src/BoundaryConditions/tests/method_of_images_unit_test.cpp

+//
+// Created by jstark on 01.06.21.
+//
+
+#define BOOST_TEST_DYN_LINK
+
+#include <boost/test/unit_test.hpp>
+
+// Include header files for testing
+#include "BoundaryConditions/MethodOfImages.hpp"
+#include "BoundaryConditions/SurfaceNormal.hpp"
+#include "libForTesting/HelpFunctionsForTestingNBCs.hpp"
+
+BOOST_AUTO_TEST_SUITE(MethodOfImagesTestSuite)
+	
+	BOOST_AUTO_TEST_CASE(DiffusionWithNoFluxNBCs)
+	{
+		const size_t dims     = 2;
+		const double L_low    = 0.0;
+		const double L_up     = 1.0;
+		const double R_disk   = 0.25;
+		const double D = 0.005; // Diffusion constant.
+		const double a = R_disk;
+		const double t0 = 0; // Must be 0 in order to fit the IC from Eq. 28
+		const double tmax = 20;
+		constexpr size_t Phi_SDF = 0, dPhi = 1, dPhi_magn = 2, SurfaceNormal = 3, Conc_n = 4, Conc_nplus1 = 5, Lap_conc = 6,
+				Radius = 7, Conc_exact = 8, Error = 9;
+		constexpr size_t x = 0, y = 1;
+		
+//		for (size_t N = 32; N <= 512; N *= 2)
+		size_t N = 32;
+		{
+			auto &v_cl = create_vcluster();
+			
+			
+			////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+			// Some indices for the grid / grid properties
+			const size_t Phi_0_grid             = 0;
+			const size_t Phi_SDF_grid           = 1;
+			////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+			// Here we create a 2D grid that stores 2 properties:
+			// Prop1: store the initial Phi;
+			// Prop2: here the re-initialized Phi (signed distance function) will be written to in the re-distancing step
+			size_t sz[dims] = {N, N}; // Grid size in terms of number of grid points per dimension
+			Box<dims, double> box({L_low, L_low}, {L_up, L_up}); // 2D
+			Ghost<dims, long int> ghost(0);
+			typedef aggregate<double, double> props;
+			typedef grid_dist_id<dims, double, props > grid_in_type;
+			grid_in_type g_dist(sz, box, ghost);
+			g_dist.setPropNames({"Phi_0", "Phi_SDF"});
+			g_dist.load("/Users/jstark/Desktop/diffusion/simple_geometries/get_diffusion_space/disk_analytical_sdf"
+			            "/output_circle_" + std::to_string(N) + "/grid_disk_analyticalSdf.bin"); // Load SDF of a disk
+			// from hdf5 file.
+			double p_spacing = g_dist.spacing(x);
+			
+			////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+			// Place particles into diffusion space according to SDF stored in grid. Copy Phi_SDF, dPhi and dPhi_magn from the
+			// grid onto the particles.
+			////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+			size_t bc[dims] = {NON_PERIODIC, NON_PERIODIC};
+			Ghost<dims, double> ghost_vd(p_spacing * 4);
+			std::cout << "Ghost layer thickness: " << p_spacing * 10 << std::endl;
+			typedef aggregate<double, double[dims], double, Point<dims, double>, double, double, double, double, double, double>
+					props_diffSpace;
+			typedef vector_dist_ws<dims, double, props_diffSpace> vd_type;
+			vd_type vd(0, box, bc, ghost_vd);
+			vd.setPropNames({"Phi_SDF", "dPhi", "dPhi_magn", "SurfaceNormal", "Conc_n", "Conc_n+1", "Lap_conc", "Radius",
+			                 "Conc_exact", "Error"});
+			
+			
+			// Get diffusion domain
+			double dlow = 0, dup = 2*L_up;
+			get_diffusion_domain<Phi_SDF_grid, Phi_SDF, dPhi, dPhi_magn>(g_dist, vd, dlow, dup);
+//	vd.write(path_output + "/vd_diffusionSpace_real", FORMAT_BINARY); // VTK file
+//	vd.save(path_output + "/vd_diffusionSpace_real.bin"); // HDF5 file
+			
+			///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+			//	Initialize the mass
+			////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+			// Get polar radius from cartesian coordinates
+			double center_disk [] = {0.5, 0.5};
+			get_radius<Radius>(vd, center_disk);
+			
+			// Initialize with exact solution for diffusion on disk at t=0
+			// get_Eq27<Radius, Conc_n>(vd, a, D, t0);
+			get_IC_Eq28<Radius, Conc_n>(vd, a);
+			
+			vd.template ghost_get<Conc_n>(KEEP_PROPERTIES);
+			
+			///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+			// Place mirror particles for noflux boundary conditions.
+			////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+			size_t mirror_width = 3; // Mirror width in number of particles
+//		size_t mirror_width = std::round(0.1/p_spacing); // Mirror width in number of particles
+			
+			if (v_cl.rank() == 0) std::cout << "Width of mirror layer: " << mirror_width << " particles." << std::endl;
+			double b_low = 0, b_up = mirror_width * p_spacing;
+			openfpm::vector<vect_dist_key_dx> keys_source;
+			get_source_particle_ids<Phi_SDF>(vd, keys_source, b_low, b_up);
+			
+			get_surface_normal_sdf_subset<Phi_SDF, dPhi, SurfaceNormal>(vd, keys_source);
+			
+			MethodOfImages<SurfaceNormal, vd_type> NBCs(vd, keys_source, 0, 1);
+			NBCs.get_mirror_particles(vd);
+			
+		}
+		
+		
+	}
+BOOST_AUTO_TEST_SUITE_END()