From dbdc366042c7ecdc10f04dc026b9630fb93b6a64 Mon Sep 17 00:00:00 2001
From: jstark <jstark@mpi-cbg.de>
Date: Mon, 3 Jan 2022 14:32:02 +0100
Subject: [PATCH] Adding generic types.
---
src/FiniteDifference/Upwind_gradient.hpp | 12 ++-
.../HelpFunctionsForGrid.hpp | 18 ++++
.../RedistancingSussman.hpp | 32 ++++---
.../tests/convergence_test.cpp | 90 ++++++++++---------
.../tests/redistancingSussman_unit_test.cpp | 2 +-
5 files changed, 91 insertions(+), 63 deletions(-)
diff --git a/src/FiniteDifference/Upwind_gradient.hpp b/src/FiniteDifference/Upwind_gradient.hpp
index 8f8f29f3..5d7047d7 100644
--- a/src/FiniteDifference/Upwind_gradient.hpp
+++ b/src/FiniteDifference/Upwind_gradient.hpp
@@ -35,9 +35,10 @@
* @param sign: Sign of the velocity with which the wave front is moving.
* @return Scalar double upwind gradient approximation in the dimension given.
*/
-static double upwinding(double dplus, double dminus, int sign)
+template <typename field_type>
+static field_type upwinding(field_type dplus, field_type dminus, int sign)
{
- double grad_upwind = 0;
+ field_type grad_upwind = 0.0;
if (dplus * sign < 0
&& (dminus + dplus) * sign < 0) grad_upwind = dplus;
@@ -66,9 +67,12 @@ static double upwinding(double dplus, double dminus, int sign)
* @return Upwind finite difference in one dimension of the property under index Field on the current node with index key.
*/
template <size_t Field, size_t Sign, typename gridtype, typename keytype>
-double FD_upwind(gridtype &grid, keytype &key, size_t d, size_t order)
+auto FD_upwind(gridtype &grid, keytype &key, size_t d, size_t order)
{
- double dplus = 0, dminus = 0;
+// std::cout << boost::typeindex::type_id_with_cvr<std::decay_t<decltype(first_node)>>() << std::endl;
+
+ typedef typename std::decay_t<decltype(grid.template get<Field>(key))> field_type;
+ field_type dplus = 0.0, dminus = 0.0;
int sign_phi0 = grid.template get<Sign> (key);
switch(order)
diff --git a/src/level_set/redistancing_Sussman/HelpFunctionsForGrid.hpp b/src/level_set/redistancing_Sussman/HelpFunctionsForGrid.hpp
index ee47be5f..86e599f8 100644
--- a/src/level_set/redistancing_Sussman/HelpFunctionsForGrid.hpp
+++ b/src/level_set/redistancing_Sussman/HelpFunctionsForGrid.hpp
@@ -295,4 +295,22 @@ void get_vector_magnitude(gridtype & grid)
}
}
+/**@brief Computes the magnitude of the gradient (L2-norm of gradient vector).
+ *
+ * @tparam Phi_grad_in Size_t index of property that contains the gradient.
+ * @tparam gridtype Type of input grid.
+ * @param grid Grid, on which the magnitude of gradient should be computed.
+ */
+template <size_t Vector_in, typename magnitude_type, typename key_type, typename gridtype>
+magnitude_type get_vector_magnitude(gridtype & grid, key_type & key)
+{
+ magnitude_type sum = 0;
+ for(size_t d = 0; d < gridtype::dims; d++)
+ {
+ sum += grid.template get<Vector_in> (key)[d] * grid.template get<Vector_in> (key)[d];
+ }
+ return sqrt(sum);
+}
+
+
#endif //REDISTANCING_SUSSMAN_HELPFUNCTIONSFORGRID_HPP
diff --git a/src/level_set/redistancing_Sussman/RedistancingSussman.hpp b/src/level_set/redistancing_Sussman/RedistancingSussman.hpp
index a95d0643..899b6ef7 100644
--- a/src/level_set/redistancing_Sussman/RedistancingSussman.hpp
+++ b/src/level_set/redistancing_Sussman/RedistancingSussman.hpp
@@ -136,15 +136,16 @@ struct Redist_options
*
* @see RedistancingSussman::get_residual_and_change_NB()
*/
+template <typename phi_type>
struct DistFromSol
{
- auto change; ///< Variable that contains the absolute value of the change of \a φ between the
+ phi_type change; ///< Variable that contains the absolute value of the change of \a φ between the
///< current
///< iteration \a i and the previous iteration \a i-1: @f[ abs( \phi_i - \phi_{i-1} ) @f] It is
///< computed for all grid
///< points that lie within the narrow band and normalized to the number of grid points that
///< lie in that narrow band.
- auto residual; ///< Variable that contains the absolute value of how far the gradient magnitude of
+ phi_type residual; ///< Variable that contains the absolute value of how far the gradient magnitude of
///< the current \a φ of iteration number \a i is away from being equal to 1: @f[ abs(|\nabla\phi_i| - 1 ) @f]
///< It is computed for all grid points that lie within the narrow band and
///< normalized to the number of grid points that lie in that narrow band.
@@ -169,7 +170,7 @@ public:
* @param redistOptions User defined options for the Sussman redistancing process
*
*/
- RedistancingSussman(grid_in_type &grid_in, Redist_options &redistOptions) : redistOptions(redistOptions),
+ RedistancingSussman(grid_in_type &grid_in, Redist_options<phi_type> &redistOptions) : redistOptions(redistOptions),
r_grid_in(grid_in),
g_temp(grid_in.getDecomposition(),
grid_in.getGridInfoVoid().getSize(),
@@ -189,7 +190,7 @@ public:
* gradient of Phi_{n+1},
* sign of the original input Phi_0 (for the upwinding).
*/
- typedef aggregate<phi_type, Point<grid_in_type::dims, phi_type>, int>
+ typedef aggregate<phi_type, phi_type[grid_in_type::dims], int>
props_temp;
/** @brief Type definition for the temporary grid.
*/
@@ -211,11 +212,12 @@ public:
* having more properties. Computes the gradients. Runs the redistancing on the internal temporary grid. Copies
* resulting signed distance function to the Phi_SDF_out property of the input grid.
*/
- template<size_t Phi_0_in, size_t Phi_SDF_out> void run_redistancing()
+ template <size_t Phi_0_in, size_t Phi_SDF_out>
+ void run_redistancing()
{
init_temp_grid<Phi_0_in>();
- init_sign_prop<Phi_n_temp, Phi_0_sign_temp>(
- g_temp); // initialize Phi_0_sign_temp with the sign of the initial (pre-redistancing) Phi_0
+ init_sign_prop<Phi_n_temp, Phi_0_sign_temp>(g_temp); // Initialize Phi_0_sign_temp with the sign of the
+ // initial (pre-redistancing) Phi_0
// Get initial gradients
get_upwind_gradient<Phi_n_temp, Phi_0_sign_temp, Phi_grad_temp>(g_temp, order_upwind_gradient, true);
iterative_redistancing(g_temp); // Do the redistancing on the temporary grid
@@ -271,10 +273,11 @@ private:
static constexpr size_t Phi_0_sign_temp = 2; ///< Property index of sign of initial (input) Phi_0 (temp. grid).
// Member variables
- Redist_options redistOptions; ///< Instantiate redistancing options.
+ Redist_options<phi_type> redistOptions; ///< Instantiate redistancing options.
grid_in_type &r_grid_in; ///< Define reference to input grid.
- DistFromSol distFromSol; ///< Instantiate distance from solution in terms of change, residual, numb. point in NB.
+ DistFromSol<phi_type> distFromSol; ///< Instantiate distance from solution in terms of change, residual, numb.
+ // point in NB.
int final_iter = 0; ///< Will be set to the final iteration when redistancing ends.
/// Transform the half-bandwidth in no_of_grid_points into physical half-bandwidth kappa.
@@ -282,7 +285,7 @@ private:
/**@brief Artificial timestep for the redistancing iterations.
* @see get_time_step_CFL(g_temp_type &grid), get_time_step(), set_user_time_step()
*/
- grid_type::stype time_step;
+ typename grid_in_type::stype time_step;
int order_upwind_gradient;
// Member functions
#ifdef SE_CLASS1
@@ -306,7 +309,7 @@ private:
template<size_t Phi_0_in>
void init_temp_grid()
{
- phi_type min_value = get_min_val<Phi_0_in>(r_grid_in); // get minimum Phi_0 value on the input grid
+ phi_type min_value = get_min_val<Phi_0_in, phi_type>(r_grid_in); // get minimum Phi_0 value on the input grid
init_grid_and_ghost<Phi_n_temp>(g_temp, min_value); // init. Phi_n_temp (incl. ghost) with min. Phi_0
copy_gridTogrid<Phi_0_in, Phi_n_temp>(r_grid_in, g_temp); // Copy Phi_0 from the input grid to Phi_n_temp
}
@@ -321,7 +324,8 @@ private:
* @return Phi_n+1 which is the Phi of the next time step on current node.
*
*/
- phi_type get_phi_nplus1(phi_type phi_n, phi_type phi_n_magnOfGrad, grid_in_type::stype dt, phi_type sgn_phi_n)
+ phi_type get_phi_nplus1(phi_type phi_n, phi_type phi_n_magnOfGrad, typename grid_in_type::stype dt, phi_type
+ sgn_phi_n)
{
return phi_n + dt * sgn_phi_n * (1 - phi_n_magnOfGrad);
}
@@ -339,7 +343,7 @@ private:
{
auto key = dom.get();
const phi_type phi_n = grid.template get<Phi_n_temp>(key);
- const phi_type phi_n_magnOfGrad = grid.template get<Phi_grad_temp>(key).norm();
+ const phi_type phi_n_magnOfGrad = get_vector_magnitude<Phi_grad_temp, phi_type>(grid, key);
phi_type epsilon = phi_n_magnOfGrad * grid.getSpacing()[0];
grid.template get<Phi_n_temp>(key) = get_phi_nplus1(phi_n, phi_n_magnOfGrad, time_step,
smooth_S(phi_n, epsilon));
@@ -375,7 +379,7 @@ private:
if (lays_inside_NB(grid.template get<Phi_n_temp>(key)))
{
count++;
- phi_type phi_n_magnOfGrad = grid.template get<Phi_grad_temp>(key).norm();
+ phi_type phi_n_magnOfGrad = get_vector_magnitude<Phi_grad_temp, phi_type>(grid, key);
phi_type epsilon = phi_n_magnOfGrad * grid.getSpacing()[0];
phi_type phi_nplus1 = get_phi_nplus1(grid.template get<Phi_n_temp>(key), phi_n_magnOfGrad, time_step,
smooth_S(grid.template get<Phi_n_temp>(key), epsilon));
diff --git a/src/level_set/redistancing_Sussman/tests/convergence_test.cpp b/src/level_set/redistancing_Sussman/tests/convergence_test.cpp
index 1b18701c..ccd7f25b 100644
--- a/src/level_set/redistancing_Sussman/tests/convergence_test.cpp
+++ b/src/level_set/redistancing_Sussman/tests/convergence_test.cpp
@@ -14,24 +14,27 @@
#include "l_norms/LNorms.hpp"
#include "analytical_SDF/AnalyticalSDF.hpp"
-double get_min_required_time_step(size_t Nmax, double b_low, double b_up, int dims)
-{
- double dx = (b_up - b_low) / (double) Nmax;
- double sum = 0;
- for (int d = 0; d < dims; d++)
- {
- sum += 1 / (dx * dx);
- }
- return 0.5 / sum;
-}
+
BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
+ typedef double phi_type;
+ typedef double space_type;
+ const phi_type EPSILON = std::numeric_limits<phi_type>::epsilon();
+
+ space_type get_min_required_time_step(size_t Nmax, space_type b_low, space_type b_up, int dims)
+ {
+ space_type dx = (b_up - b_low) / (space_type) Nmax;
+ space_type sum = 0;
+ for (int d = 0; d < dims; d++)
+ {
+ sum += 1 / (dx * dx);
+ }
+ return 0.5 / sum;
+ }
BOOST_AUTO_TEST_CASE(RedistancingSussman_convergence_test_1D)
{
// CFL dt in 1D for N=64 and c=1 is 0.000503905
- const double EPSILON = std::numeric_limits<double>::epsilon();
- std::cout << "epsilon = " << EPSILON << std::endl;
const size_t grid_dim = 1;
// some indices
@@ -42,17 +45,17 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
const size_t SDF_exact_grid = 2;
const size_t Error_grid = 3;
- const double box_lower = -1.0;
- const double box_upper = 1.0;
- Box<grid_dim, double> box({box_lower}, {box_upper});
+ const space_type box_lower = -1.0;
+ const space_type box_upper = 1.0;
+ Box<grid_dim, space_type> box({box_lower}, {box_upper});
Ghost<grid_dim, long int> ghost(0);
- typedef aggregate<double, double, double, double> props;
- typedef grid_dist_id<grid_dim, double, props> grid_in_type;
+ typedef aggregate<phi_type, phi_type, phi_type, phi_type> props;
+ typedef grid_dist_id<grid_dim, space_type, props> grid_in_type;
-// double dt = 0.000503905;
+// space_type dt = 0.000503905;
size_t Nmin = 32;
size_t Nmax = 4096;
- double dt = get_min_required_time_step(Nmax, box_lower, box_upper, grid_dim);
+ space_type dt = get_min_required_time_step(Nmax, box_lower, box_upper, grid_dim);
for (size_t N = Nmin; N <= Nmax; N*=2)
{
const size_t sz[grid_dim] = {N};
@@ -64,14 +67,14 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
while(dom.isNext())
{
auto key = dom.get();
- Point<grid_dim, double> coords = g_dist.getPos(key);
+ Point<grid_dim, space_type> coords = g_dist.getPos(key);
g_dist.template get<Phi_0_grid>(key) = sgn(coords.get(x));
g_dist.template get<SDF_exact_grid>(key) = coords.get(x);
++dom;
}
{
- Redist_options redist_options;
+ Redist_options<phi_type> redist_options;
redist_options.min_iter = 1e3;
redist_options.max_iter = 1e16;
@@ -89,7 +92,7 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
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> redist_obj(g_dist,
+ RedistancingSussman<grid_in_type, phi_type> redist_obj(g_dist,
redist_options); // Instantiation of Sussman-redistancing class
std::cout << "CFL dt = " << redist_obj.get_time_step() << std::endl;
redist_obj.set_user_time_step(dt);
@@ -99,14 +102,14 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
// Compute the absolute error between analytical and numerical solution at each grid point
get_absolute_error<SDF_sussman_grid, SDF_exact_grid, Error_grid>(g_dist);
- g_dist.write("grid_postRedistancing", FORMAT_BINARY);
+// g_dist.write("grid_postRedistancing", FORMAT_BINARY);
/////////////////////////////////////////////////////////////////////////////////////////////
// Get narrow band: Place particles on interface (narrow band width e.g. 4 grid points on each side of the
// interface)
size_t bc[grid_dim] = {NON_PERIODIC};
- typedef aggregate<double> props_nb;
- typedef vector_dist<grid_dim, double, props_nb> vd_type;
- Ghost<grid_dim, double> ghost_vd(0);
+ typedef aggregate<phi_type> props_nb;
+ typedef vector_dist<grid_dim, space_type, props_nb> vd_type;
+ Ghost<grid_dim, space_type> ghost_vd(0);
vd_type vd_narrow_band(0, box, bc, ghost_vd);
vd_narrow_band.setPropNames({"error"});
const size_t Error_vd = 0;
@@ -115,7 +118,7 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
NarrowBand<grid_in_type> narrowBand_points(g_dist, narrow_band_width); // 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_error_N" + std::to_string(N), FORMAT_BINARY);
+// vd_narrow_band.write("vd_error_N" + std::to_string(N), 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
L_norms lNorms_vd;
@@ -128,8 +131,6 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
BOOST_AUTO_TEST_CASE(RedistancingSussman_convergence_test_3D)
{
- const double EPSILON = std::numeric_limits<double>::epsilon();
-
const size_t grid_dim = 3;
// some indices
const size_t x = 0;
@@ -143,24 +144,24 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
for (size_t N=32; N <=128; N*=2)
{
- const double dt = 0.000165334;
+ const space_type dt = 0.000165334;
const size_t sz[grid_dim] = {N, N, N};
- const double radius = 1.0;
- const double box_lower = 0.0;
- const double box_upper = 4.0 * radius;
- Box<grid_dim, double> box({box_lower, box_lower, box_lower}, {box_upper, box_upper, box_upper});
+ const space_type radius = 1.0;
+ const space_type box_lower = 0.0;
+ const space_type box_upper = 4.0 * radius;
+ Box<grid_dim, space_type> box({box_lower, box_lower, box_lower}, {box_upper, box_upper, box_upper});
Ghost<grid_dim, long int> ghost(0);
- typedef aggregate<double, double, double, double> props;
- typedef grid_dist_id<grid_dim, double, props > grid_in_type;
+ typedef aggregate<phi_type, phi_type, phi_type, phi_type> props;
+ typedef grid_dist_id<grid_dim, space_type, props > grid_in_type;
grid_in_type g_dist(sz, box, ghost);
g_dist.setPropNames({"Phi_0", "SDF_sussman", "SDF_exact", "Relative error"});
- const double center[grid_dim] = {0.5*(box_upper+box_lower), 0.5*(box_upper+box_lower), 0.5*(box_upper+box_lower)};
+ const space_type center[grid_dim] = {0.5*(box_upper+box_lower), 0.5*(box_upper+box_lower), 0.5*(box_upper+box_lower)};
init_grid_with_sphere<Phi_0_grid>(g_dist, radius, center[x], center[y], center[z]); // Initialize sphere onto grid
int order = 1;
{
- Redist_options redist_options;
+ Redist_options<phi_type> redist_options;
redist_options.min_iter = 1e4;
redist_options.max_iter = 1e4;
@@ -174,7 +175,8 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
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> redist_obj(g_dist, redist_options); // Instantiation of Sussman-redistancing class
+ 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;
// 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.
@@ -191,9 +193,9 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
// Get narrow band: Place particles on interface (narrow band width e.g. 4 grid points on each side of the
// interface)
size_t bc[grid_dim] = {PERIODIC, PERIODIC, PERIODIC};
- typedef aggregate<double> props_nb;
- typedef vector_dist<grid_dim, double, props_nb> vd_type;
- Ghost<grid_dim, double> ghost_vd(0);
+ typedef aggregate<phi_type> props_nb;
+ typedef vector_dist<grid_dim, space_type, props_nb> vd_type;
+ Ghost<grid_dim, space_type> ghost_vd(0);
vd_type vd_narrow_band(0, box, bc, ghost_vd);
vd_narrow_band.setPropNames({"error"});
const size_t Error_vd = 0;
@@ -202,8 +204,8 @@ BOOST_AUTO_TEST_SUITE(ConvergenceTestSuite)
NarrowBand<grid_in_type> narrowBand_points(g_dist, narrow_band_width); // 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_nb8p_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
L_norms lNorms_vd;
diff --git a/src/level_set/redistancing_Sussman/tests/redistancingSussman_unit_test.cpp b/src/level_set/redistancing_Sussman/tests/redistancingSussman_unit_test.cpp
index c3d9516e..34950d59 100644
--- a/src/level_set/redistancing_Sussman/tests/redistancingSussman_unit_test.cpp
+++ b/src/level_set/redistancing_Sussman/tests/redistancingSussman_unit_test.cpp
@@ -51,7 +51,7 @@ BOOST_AUTO_TEST_SUITE(RedistancingSussmanTestSuite)
init_grid_with_sphere<Phi_0_grid>(g_dist, radius, center[x], center[y], center[z]); // Initialize sphere onto grid
- Redist_options redist_options;
+ Redist_options<phi_type> redist_options;
redist_options.min_iter = 1e4;
redist_options.max_iter = 1e4;
--
GitLab