Fixed cfl condition for timestep during redistancing.

src/level_set/redistancing_Sussman/HelpFunctionsForGrid.hpp

@@ -19,31 +19,14 @@
 #include "VCluster/VCluster.hpp"
 #include "Grid/grid_dist_id.hpp"
 
-/**@brief Computes the time step for the iterative re-distancing fulfilling CFL condition.
- *
- * @tparam grid_type Template type of the input grid.
- * @param grid Input OpenFPM grid.
- * @return Time step.
- */
-template <typename grid_type>
-typename grid_type::stype get_time_step_CFL(grid_type &grid)
-{
-	typename grid_type::stype sum = 0.0;
-	for (size_t d = 0; d < grid_type::dims; d++)
-	{
-		sum += 1.0 / (grid.spacing(d) * grid.spacing(d));
-	}
-	return 0.5 / sum;
-}
 
-#if 0
 /**@brief Computes the time step size fulfilling CFL condition according to https://www.cfd-online
  * .com/Wiki/Courant–Friedrichs–Lewy_condition for arbitrary dimensionality.
  *
  * @tparam grid_type Template type of the input grid.
  * @param grid Input OpenFPM grid.
  * @param u Array of size grid_type::dims containing the velocity in each dimension.
- * @param Cmax Courant number.
+ * @param C Courant number.
  * @return Time step.
  */
 template <typename grid_type>
@@ -62,7 +45,7 @@ typename grid_type::stype get_time_step_CFL(grid_type & grid, typename grid_type
  * @tparam grid_type Template type of the input grid.
  * @param grid Input OpenFPM grid.
  * @param u Velocity of propagating wave if isotropic for each direction.
- * @param Cmax Courant number.
+ * @param C Courant number.
  * @return Time step.
  */
 template <typename grid_type>
@@ -75,7 +58,7 @@ typename grid_type::stype get_time_step_CFL(grid_type & grid, typename grid_type
 	}
 	return C / sum;
 }
-#endif
+
 /**@brief Initializes given property \p Prop of an OpenFPM grid including ghost layer with a given value from \p
  * init_value.
  *

src/level_set/redistancing_Sussman/RedistancingSussman.hpp

@@ -175,7 +175,8 @@ public:
 	                                                                                   grid_in.getGridInfoVoid().getSize(),
 	                                                                                   Ghost<grid_in_type::dims, long int>(3))
 	{
-		time_step = get_time_step_CFL(grid_in);
+		// Get timestep fulfilling CFL condition for velocity=1.0 and courant number=0.1
+		time_step = get_time_step_CFL(grid_in, 1.0, 0.1);
 		order_upwind_gradient = 1;
 #ifdef SE_CLASS1
 		assure_minimal_thickness_of_NB();

src/level_set/redistancing_Sussman/tests/convergence_test.cpp

@@ -114,8 +114,7 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
 				vd_narrow_band.setPropNames({"error"});
 				const size_t Error_vd = 0;
 				// Compute the L_2- and L_infinity-norm and save to file
-				size_t narrow_band_width = redist_options.width_NB_in_grid_points - 2;
-				NarrowBand<grid_in_type, phi_type> narrowBand_points(g_dist, narrow_band_width); // Instantiation of
+				NarrowBand<grid_in_type, phi_type> narrowBand_points(g_dist, redist_options.width_NB_in_grid_points ); // Instantiation of
 				// NarrowBand class
 				narrowBand_points.get_narrow_band_copy_specific_property<SDF_sussman_grid, Error_grid, Error_vd>(g_dist,
 				                                                                                                 vd_narrow_band);
@@ -163,8 +162,8 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
 			int order = 1;
 			{
 				Redist_options<phi_type> redist_options;
-				redist_options.min_iter                             = 1e4;
-				redist_options.max_iter                             = 1e4;
+				redist_options.min_iter                             = 1e3;
+				redist_options.max_iter                             = 1e3;
 				
 				// set both convergence criteria to false s.t. termination only when max_iterations reached
 				redist_options.convTolChange.check                  = false;    // define here which of the convergence criteria above should be used. If both are true, termination only occurs when both are fulfilled or when iter > max_iter
@@ -172,14 +171,15 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
 				
 				redist_options.interval_check_convergence           = 100;        // interval of #iterations at which
 				// convergence is checked (default: 100)
-				redist_options.width_NB_in_grid_points              = 8;        // width of narrow band in number of grid points. Must be at least 4, in order to have at least 2 grid points on each side of the interface. (default: 4)
+				redist_options.width_NB_in_grid_points              = 4;        // width of narrow band in number of
+				// grid points. Must be at least 4, in order to have at least 2 grid points on each side of the interface. (default: 4)
 				redist_options.print_current_iterChangeResidual     = true;     // if true, prints out every current iteration + corresponding change from the previous iteration + residual from SDF (default: false)
 				redist_options.print_steadyState_iter               = true;     // if true, prints out the final iteration number when steady state was reached + final change + residual (default: true)
 				
 				RedistancingSussman<grid_in_type, phi_type> redist_obj(g_dist, redist_options);   // Instantiation of
 				// Sussman-redistancing class
-				redist_obj.set_user_time_step(dt);
-				std::cout << "dt set to = " << dt << std::endl;
+				//redist_obj.set_user_time_step(dt);
+				//std::cout << "dt set to = " << dt << std::endl;
 				// Run the redistancing. in the <> brackets provide property-index where 1.) your initial Phi is stored and 2.) where the resulting SDF should be written to.
 				redist_obj.run_redistancing<Phi_0_grid, SDF_sussman_grid>();
 				
@@ -201,33 +201,36 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
 				vd_narrow_band.setPropNames({"error"});
 				const size_t Error_vd = 0;
 				// Compute the L_2- and L_infinity-norm and save to file
-				size_t narrow_band_width = 8;
-				NarrowBand<grid_in_type, phi_type> narrowBand_points(g_dist, narrow_band_width); // Instantiation of NarrowBand class
+				NarrowBand<grid_in_type, phi_type> narrowBand_points(g_dist, redist_options.width_NB_in_grid_points); // Instantiation of NarrowBand class
 				narrowBand_points.get_narrow_band_copy_specific_property<SDF_sussman_grid, Error_grid, Error_vd>(g_dist,
 				                                                                                                 vd_narrow_band);
-//				vd_narrow_band.write("vd_nb8p_error_N" + std::to_string(N) + "_order" +
-//				std::to_string(order), FORMAT_BINARY);
+				vd_narrow_band.write("vd_nb4p_error_N" + std::to_string(N) + "_order" +
+				std::to_string(order), FORMAT_BINARY);
 //				vd_narrow_band.save("test_data/output/vd_nb8p_error" + std::to_string(N) + ".bin");
 				// Compute the L_2- and L_infinity-norm and save to file
 				LNorms<phi_type> lNorms_vd;
 				lNorms_vd.get_l_norms_vector<Error_vd>(vd_narrow_band);
-				lNorms_vd.write_to_file(N, 6, "l_norms_vd_absError_8p_order" + std::to_string(order), "./");
+				lNorms_vd.write_to_file(N, 6, "l_norms_vd_absError_4p_order" + std::to_string(order), "./");
 				
-//				switch(order)
-//				{
-//					case 1:
-//						BOOST_CHECK(lNorms_vd.l2   < 0.03369 + EPSILON);
-//						BOOST_CHECK(lNorms_vd.linf < 0.06307 + EPSILON);
-//						break;
-//					case 3:
-//						BOOST_CHECK(lNorms_vd.l2   < 0.02794 + EPSILON);
-//						BOOST_CHECK(lNorms_vd.linf < 0.0586704 + EPSILON);
-//						break;
-//					case 5:
-//						BOOST_CHECK(lNorms_vd.l2   < 0.0187199 + EPSILON);
-//						BOOST_CHECK(lNorms_vd.linf < 0.0367638 + EPSILON);
-//						break;
-//				}
+//				32,0.044692,0.066994
+//				64,0.010542,0.018396
+//				128,0.003301,0.009491
+
+				switch(N)
+				{
+					case 32:
+						BOOST_CHECK(lNorms_vd.l2   < 0.044693);
+						BOOST_CHECK(lNorms_vd.linf < 0.066995);
+						break;
+					case 64:
+						BOOST_CHECK(lNorms_vd.l2   < 0.010543);
+						BOOST_CHECK(lNorms_vd.linf < 0.018397);
+						break;
+					case 128:
+						BOOST_CHECK(lNorms_vd.l2   < 0.003302);
+						BOOST_CHECK(lNorms_vd.linf < 0.009492);
+						break;
+				}
 			
 			}
 		}

src/level_set/redistancing_Sussman/tests/help_functions_unit_test.cpp

@@ -23,18 +23,28 @@ BOOST_AUTO_TEST_SUITE(HelpFunctionsTestSuite)
 		const size_t sz[grid_dim] = {32};
 		grid_type g_dist(sz, box, ghost);
 		
-		double i = 3.3;
+		// Each processor assigns smaller_value to the first grid node in its domain
+		double smaller_value = 0.1;
 		auto dom = g_dist.getDomainIterator();
+		auto key_0 = dom.get();
+		g_dist.template get<Field>(key_0) = smaller_value;
+		++dom;
+		
+		// Afterwards we loop over the other grid nodes and assign them another bigger number
+		double bigger_value = 0.5;
 		while(dom.isNext())
 		{
-			i -= 0.1;
 			auto key = dom.get();
-			g_dist.template get<Field>(key) = i;
+			g_dist.template get<Field>(key) = bigger_value;
 			++dom;
 		}
-		auto correct_value = i;
+		
+		// Now we check if get_min_value returns smaller_value
 		auto min_value = get_min_val<Field>(g_dist);
-//		BOOST_CHECK_MESSAGE(min_value == correct_value, "Checking if minimum value stored in grid is returned.");
+		double tolerance = 1e-12;
+//BOOST_CHECK_MESSAGE(isApproxEqual(min_value, smaller_value, tolerance), "Checking if smallest value stored "
+//		"in grid "
+//																		  "is returned.");
 	}
 	
 	BOOST_AUTO_TEST_CASE(get_max_val_test)
@@ -49,17 +59,27 @@ BOOST_AUTO_TEST_SUITE(HelpFunctionsTestSuite)
 		const size_t sz[grid_dim] = {32};
 		grid_type g_dist(sz, box, ghost);
 		
-		double i = 0.0;
+		// Each processor assigns smaller_value to the first grid node in its domain
+		double smaller_value = 0.1;
 		auto dom = g_dist.getDomainIterator();
+		auto key_0 = dom.get();
+		g_dist.template get<Field>(key_0) = smaller_value;
+		++dom;
+		
+		// Afterwards we loop over the other grid nodes and assign them another bigger number
+		double bigger_value = 0.5;
 		while(dom.isNext())
 		{
-			i += 0.1;
 			auto key = dom.get();
-			g_dist.template get<Field>(key) = i;
+			g_dist.template get<Field>(key) = bigger_value;
 			++dom;
 		}
-		auto correct_value = i;
+		
+		// Now we check if get_max_value returns bigger_value
 		auto max_value = get_max_val<Field>(g_dist);
-//		BOOST_CHECK_MESSAGE(max_value == correct_value, "Checking if maximum value stored in grid is returned.");
+		double tolerance = 1e-12;
+//		BOOST_CHECK_MESSAGE(isApproxEqual(max_value, bigger_value, tolerance), "Checking if smallest value stored in "
+//																			   "grid "
+//																		  "is returned.");
 	}
 BOOST_AUTO_TEST_SUITE_END()