Commit e6fd81f1 authored by jstark's avatar jstark
Browse files

Backing up changed in loops.

parent 486fa9e6
......@@ -62,7 +62,7 @@ struct Conv_tol_change
///< between the iterations is considered and the redistancing finishes when the Conv_tol_change.value (or the
///< max-iteration) is reached. If false, the change is not considered as steady-state
///< criterion. Default: true.
double value = 1e-6; ///< Double variable that stores the tolerance value for the change at which the iterative
double value = 1e-15; ///< Double variable that stores the tolerance value for the change at which the iterative
///< redistancing process is considered as steady-state. (see also DistFromSol::change)
};
......@@ -73,10 +73,10 @@ struct Conv_tol_change
*/
struct Conv_tol_residual
{
bool check = false; ///< If true, the residual of Phi (see DistFromSol::residual) is considered and the
bool check = true; ///< If true, the residual of Phi (see DistFromSol::residual) is considered and the
///< redistancing finishes when the Conv_tol_residual.value (or the
///< max-iteration) is reached. If false, the change is not considered as steady-state criterion. Default: false.
double value = 1e-1; ///< Double variable that stores the tolerance value for the residual at which the
///< max-iteration) is reached. If false, the change is not considered as steady-state criterion. Default: true.
double value = 1e-3; ///< Double variable that stores the tolerance value for the residual at which the
///< iterative redistancing process is considered as steady-state. (see also
///< DistFromSol::residual)
};
......@@ -86,35 +86,36 @@ struct Conv_tol_residual
* @details For the redistancing, we can choose some options. These options will then be passed bundled as a structure to
* the redistancing function. Setting these options is optional, since they all have a Default value as well.
*
* @param min_iter: Minimum number of iterations before steady state in narrow band will be checked (Default: 100).
* @param min_iter: Minimum number of iterations before steady state in narrow band will be checked (Default: 1e5).
* @param max_iter: Maximum number of iterations you want to run the redistancing, even if steady state might not yet
* have been reached (Default: 1e6).
* have been reached (Default: 1e12).
* @param order_space_op: Order of accuracy of the upwind gradient computation when solving the eikonal equation
* during the redistancing. Options are: {1, 3, 5}. Default: 5;
* during the redistancing. Options are: {1, 3, 5} (Default: 5).
* @param convTolChange.value: Convolution tolerance for the normalized total change of Phi in the narrow band between
* two consecutive iterations (Default: 1e-6).
* two consecutive iterations (Default: 1e-15).
* @param convTolChange.check: Set true, if you want to use the normalized total change between two iterations as
* measure of how close you are to the steady state solution. Redistancing will then stop
* if convTolChange.value is reached or if the current iteration is bigger than max_iter.
* @param convTolResidual.value: Convolution tolerance for the residual, that is abs(magnitude gradient of phi - 1) of
* Phi in the narrow band (Default 1e-1).
* Phi in the narrow band (Default 1e-3).
* @param convTolResidual.check: Set true, if you want to use the residual of the current iteration as measure of how
* close you are to the steady state solution. Redistancing will then stop if
* convTolResidual.value is reached or if the current iteration is bigger than max_iter.
* @param interval_check_convergence: Interval of number of iterations at which convergence to steady state is checked
* (Default: 100).
* @param width_NB_in_grid_points: 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. Is automatically set
* to 4, if a value smaller than 4 is chosen (Default: 4).
* @param width_NB_in_grid_points: Width of narrow band in number of grid points. Convergence is checked for this
* area around the interface only, so don't choose too small! (Default: 8).
* @param print_current_iterChangeResidual: If true, the number of the current iteration, the corresponding change
* w.r.t the previous iteration and the residual is printed (Default: false).
* @param print_steadyState_iter: If true, the number of the steady-state-iteration, the corresponding change
* w.r.t the previous iteration and the residual is printed (Default: false).
* @param save_temp_grid: If true, save the temporary grid as hdf5 that can be reloaded onto a grid
*/
struct Redist_options
{
size_t min_iter = 1000;
size_t max_iter = 1e6;
size_t min_iter = 1e5;
size_t max_iter = 1e12;
size_t order_space_op = 5;
......@@ -122,7 +123,7 @@ struct Redist_options
Conv_tol_residual convTolResidual;
size_t interval_check_convergence = 100;
size_t width_NB_in_grid_points = 4;
size_t width_NB_in_grid_points = 8;
bool print_current_iterChangeResidual = false;
bool print_steadyState_iter = true;
bool save_temp_grid = false;
......@@ -173,8 +174,9 @@ public:
Ghost<grid_in_type::dims, long int>(3))
{
time_step = get_time_step_CFL(g_temp);
// assure_minimal_thickness_of_NB(); // overwrites user-set NB thickness in which convergence and residual is
// checked
#ifdef SE_CLASS1
assure_minimal_thickness_of_NB();
#endif // SE_CLASS1
}
/**@brief Aggregated properties for the temporary grid.
......@@ -209,15 +211,13 @@ public:
template<size_t Phi_0_in, size_t Phi_SDF_out> void run_redistancing()
{
init_temp_grid<Phi_0_in>();
init_sign_prop<Phi_0_temp, Phi_0_sign_temp>(
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_0_temp, Phi_0_sign_temp, Phi_grad_temp>(g_temp, redistOptions.order_space_op, true);
get_upwind_gradient<Phi_n_temp, Phi_0_sign_temp, Phi_grad_temp>(g_temp, redistOptions.order_space_op, true);
iterative_redistancing(g_temp); // Do the redistancing on the temporary grid
copy_gridTogrid<Phi_nplus1_temp, Phi_SDF_out>(g_temp, r_grid_in); // Copy resulting SDF function to input grid
}
/** @brief Overwrite the time_step found via CFL condition with an individual time_step.
*
* @details If the user wants to overwrite the time_step found via CFL condition with an individual time_step.
......@@ -243,26 +243,38 @@ public:
private:
// Some indices for better readability
static constexpr size_t Phi_0_temp = 0; ///< Property index of Phi_0 on the temporary grid.
static constexpr size_t Phi_n_temp = 0; ///< Property index of Phi_0 on the temporary grid.
static constexpr size_t Phi_nplus1_temp = 1; ///< Property index of Phi_n+1 on the temporary grid.
static constexpr size_t Phi_grad_temp = 2; ///< Property index of gradient of Phi_n on the temporary grid.
static constexpr size_t Phi_0_sign_temp = 3; ///< Property index of sign of initial (input) Phi_0 (temp. grid).
// Member variables
Redist_options redistOptions; ///< Instantiate redistancing options.
grid_in_type &r_grid_in; ///< Define reference to input grid.
double h_max = get_biggest_spacing(g_temp); ///< Grid spacing in less resolved direction.
double h_min = get_smallest_spacing(g_temp); ///< Grid spacing in higher resolved direction.
/// Transform the half-bandwidth in no_of_grid_points into physical half-bandwidth kappa.
double kappa = ceil(redistOptions.width_NB_in_grid_points / 2.0) * h_max;
double kappa = ceil(redistOptions.width_NB_in_grid_points / 2.0) * get_biggest_spacing(g_temp);
/**@brief Artificial timestep for the redistancing iterations.
* @see get_time_step_CFL(g_temp_type &grid), get_time_step(), set_user_time_step()
*/
double time_step;
// Member functions
/** @brief Checks if narrow band thickness >= 4 grid points. Else, sets it to 4 grid points.
*
* @details Makes sure, that the narrow band within which the convergence criteria are checked during the
* redistancing, is thick enough.
*/
void assure_minimal_thickness_of_NB()
{
if (redistOptions.width_NB_in_grid_points < 8)
{
std::cout << "The narrow band thickness that you chose for the convergence check is very small. Consider "
"setting redist_options.width_NB_in_grid_points to at least 8" << std::endl;
} // check narrow band width of convergence check if defined SE_CLASS1
}
/** @brief Copies values from input grid to internal temporary grid and initializes ghost layer with minimum
* value of input grid.
*/
......@@ -270,44 +282,24 @@ private:
void init_temp_grid()
{
double min_value = get_min_val<Phi_0_in>(r_grid_in); // get minimum Phi_0 value on the input grid
init_grid_and_ghost<Phi_0_temp>(g_temp, min_value); // init. Phi_0_temp (incl. ghost) with min. Phi_0
init_grid_and_ghost<Phi_n_temp>(g_temp, min_value); // init. Phi_n_temp (incl. ghost) with min. Phi_0
init_grid_and_ghost<Phi_nplus1_temp>(g_temp, min_value); // init. Phi_nplus1_temp (incl. ghost) with min. Phi_0
copy_gridTogrid<Phi_0_in, Phi_0_temp>(r_grid_in, g_temp); // Copy Phi_0 from the input grid to Phi_0_temp
copy_gridTogrid<Phi_0_in, Phi_n_temp>(r_grid_in, g_temp); // Copy Phi_0 from the input grid to Phi_n_temp
}
/** @brief Checks if narrow band thickness >= 4 grid points. Else, sets it to 4 grid points.
*
* @details Makes sure, that the narrow band within which the convergence criteria are checked during the
* redistancing, is thick enough.
*/
void assure_minimal_thickness_of_NB()
{
if (redistOptions.width_NB_in_grid_points < 4)
{
redistOptions.width_NB_in_grid_points = 4;
} // overwrite kappa if set too small by user
}
/** @brief Run one timestep of re-distancing and compute Phi_n+1.
*
* @param phi_n Phi value on current node and current time.
* @param phi_n_magnOfGrad Gradient magnitude of current Phi from upwinding FD.
* @param dt Time step.
* @param sign_phi_n Sign of the current Phi, should be the smooth sign.
* @param sgn_phi_n Sign of the current Phi, should be the smooth sign.
*
* @return Phi_n+1 which is the Phi of the next time step on current node.
*
*/
double get_phi_nplus1(double phi_n, double phi_n_magnOfGrad, double dt, double sign_phi_n)
double get_phi_nplus1(double phi_n, double phi_n_magnOfGrad, double dt, double sgn_phi_n)
{
double step = dt * sign_phi_n * (1 - phi_n_magnOfGrad); // <- original Sussman
// if (step > 10)
// {
// std::cout << "phi_n_magnOfGrad = " << phi_n_magnOfGrad << ", step = " << step
// << ", skip to prevent exploding peaks." << std::endl;
// step = 0;
// }
return phi_n + step;
return phi_n + dt * sgn_phi_n * (1 - phi_n_magnOfGrad);
}
/** @brief Go one re-distancing time-step on the whole grid.
......@@ -316,13 +308,13 @@ private:
*/
void go_one_redistancing_step_whole_grid(g_temp_type &grid)
{
grid.template ghost_get<Phi_0_temp, Phi_nplus1_temp, Phi_grad_temp>();
grid.template ghost_get<Phi_n_temp, Phi_nplus1_temp, Phi_grad_temp>();
double spacing_x = grid.getSpacing()[0];
auto dom = grid.getDomainIterator();
while (dom.isNext())
{
auto key = dom.get();
const double phi_n = grid.template get<Phi_0_temp>(key);
const double phi_n = grid.template get<Phi_n_temp>(key);
const double phi_n_magnOfGrad = grid.template get<Phi_grad_temp>(key).norm();
double epsilon = phi_n_magnOfGrad * spacing_x;
grid.template get<Phi_nplus1_temp>(key) = get_phi_nplus1(phi_n, phi_n_magnOfGrad, time_step,
......@@ -337,8 +329,8 @@ private:
*/
void update_grid(g_temp_type &grid)
{
copy_gridTogrid<Phi_nplus1_temp, Phi_0_temp>(grid, grid); // Update Phi_0
get_upwind_gradient<Phi_0_temp, Phi_0_sign_temp, Phi_grad_temp>(grid, redistOptions.order_space_op, true);
copy_gridTogrid<Phi_nplus1_temp, Phi_n_temp>(grid, grid); // Update Phi_0
get_upwind_gradient<Phi_n_temp, Phi_0_sign_temp, Phi_grad_temp>(grid, redistOptions.order_space_op, true);
}
/** @brief Checks if a node lays within the narrow band around the interface.
......@@ -373,7 +365,7 @@ private:
total_points_in_nb += 1.0;
auto dphi_magn = grid.template get<Phi_grad_temp>(key).norm();
total_residual += abs(dphi_magn - 1);
total_change += abs(grid.template get<Phi_nplus1_temp>(key) - grid.template get<Phi_0_temp>(key));
total_change += abs(grid.template get<Phi_nplus1_temp>(key) - grid.template get<Phi_n_temp>(key));
}
++dom;
}
......@@ -451,49 +443,47 @@ private:
*/
void iterative_redistancing(g_temp_type &grid)
{
for (size_t i = 0; i <= redistOptions.max_iter; i++)
int i = 0;
while (i <= redistOptions.max_iter)
{
go_one_redistancing_step_whole_grid(grid);
if (i % redistOptions.interval_check_convergence == 0) // after some iterations check if steady state
// is reached in the narrow band
for (int j = 0; j < redistOptions.interval_check_convergence; j++)
{
if (redistOptions.print_current_iterChangeResidual)
{
print_out_iteration_change_residual(grid, i);
}
if (i >= redistOptions.min_iter)
go_one_redistancing_step_whole_grid(grid);
++i;
}
if (redistOptions.print_current_iterChangeResidual)
{
print_out_iteration_change_residual(grid, i);
}
if (i >= redistOptions.min_iter)
{
if (steady_state_NB(grid))
{
if (steady_state_NB(grid))
if (redistOptions.print_steadyState_iter)
{
if (redistOptions.print_steadyState_iter)
auto &v_cl = create_vcluster();
if (v_cl.rank() == 0)
{
auto &v_cl = create_vcluster();
if (v_cl.rank() == 0)
{
std::cout << "Steady state reached at iteration: " << i << std::endl;
std::cout << "Final_Iteration,Change,Residual" << std::endl;
}
print_out_iteration_change_residual(grid, i);
std::cout << "Steady state reached at iteration: " << i << std::endl;
std::cout << "Final_Iteration,Change,Residual" << std::endl;
}
update_grid(grid); // Update Phi
if (redistOptions.save_temp_grid)
{
g_temp.setPropNames({"Phi_0", "Phi_nplus1_temp", "Phi_grad_temp", "Phi_0_sign_temp"});
g_temp.save("g_temp_redistancing.hdf5"); // HDF5 file}
}
break;
print_out_iteration_change_residual(grid, i);
}
update_grid(grid); // Update Phi
if (redistOptions.save_temp_grid)
{
g_temp.setPropNames({"Phi_n", "Phi_nplus1_temp", "Phi_grad_temp", "Phi_0_sign_temp"});
g_temp.save("g_temp_redistancing.hdf5"); // HDF5 file}
}
break;
}
}
update_grid(grid);
}
// If save_temp_grid set true, save the temporary grid to an hdf5 file that can be reloaded onto a grid and
// reused
// If save_temp_grid set true, save the temporary grid as hdf5 that can be reloaded onto a grid
if (redistOptions.save_temp_grid)
{
g_temp.setPropNames({"Phi_0", "Phi_nplus1_temp", "Phi_grad_temp", "Phi_0_sign_temp"});
g_temp.setPropNames({"Phi_n", "Phi_nplus1_temp", "Phi_grad_temp", "Phi_0_sign_temp"});
g_temp.save("g_temp_redistancing.hdf5"); // HDF5 file}
}
}
......
......@@ -8,60 +8,81 @@
// Include redistancing files
#include "util/PathsAndFiles.hpp"
#include "level_set/redistancing_Sussman/RedistancingSussman.hpp"
#include "level_set/redistancing_Sussman/NarrowBand.hpp"
//#include "level_set/redistancing_Sussman/NarrowBand.hpp"
// Include header files for testing
#include "Draw/DrawSphere.hpp"
#include "l_norms/LNorms.hpp"
//#include "Draw/DrawSphere.hpp"
//#include "l_norms/LNorms.hpp"
#include "analytical_SDF/AnalyticalSDF.hpp"
BOOST_AUTO_TEST_SUITE(RedistancingSussmanTestSuite)
BOOST_AUTO_TEST_CASE(RedistancingSussman_unit_sphere)
{
const double EPSILON = std::numeric_limits<double>::epsilon();
const size_t grid_dim = 3;
// some indices
const size_t x = 0;
const size_t y = 1;
const size_t z = 2;
const size_t Phi_0_grid = 0;
const size_t SDF_sussman_grid = 1;
const size_t SDF_exact_grid = 2;
const size_t Error_grid = 3;
size_t N = 32;
const double dt = 0.000165334;
const size_t Phi_0_grid = 0;
size_t N = 4;
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});
Ghost<grid_dim, long int> ghost(0);
Box<grid_dim, double> box({0,0,0}, {1,1,1});
Ghost<grid_dim, long int> ghost(1);
typedef aggregate<double, double, double, double> props;
typedef grid_dist_id<grid_dim, double, props > grid_in_type;
grid_in_type g_dist(sz, box, ghost);
g_dist.setPropNames({"Phi_0", "SDF_sussman", "SDF_exact", "Relative error"});
std::cout << "---------------------------" << std::endl;
auto dom = g_dist.getDomainIterator();
std::cout << "demangle(typeid(decltype(dom)).name())" << demangle(typeid(decltype(dom)).name()) << std::endl;
auto key2 = dom.get();
std::cout << key2.to_string() << std::endl;
const double 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
while(dom.isNext())
{
std::cout << "start while loop" << std::endl;
auto key = dom.get();
std::cout << key.to_string() << std::endl;
if (key.getKeyRef().get(2) >= 34 || key.getKeyRef().get(2) <= 0)
{
std::cout << "ERROR at key" << key.getKeyRef().get(2) << std::endl;
abort();
}
g_dist.template get<Phi_0_grid>(key) = 1;
++dom;
}
#if 0
const double center[grid_dim] = {0.5*(box_lower + box_upper),
0.5*(box_lower + box_upper),
0.5*(box_lower + box_upper)};
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.min_iter = 1e4;
redist_options.max_iter = 1e4; // max. number of iterations you want to run the
// redistancing, even if steady state might not yet have been reached (default: 1e6)
redist_options.order_space_op = 5;
// define here which of the convergence criteria above should be used. If both are false, termination occurs
// when when iter > max_iter
redist_options.convTolChange.check = false;
redist_options.convTolResidual.check = false;
// 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
redist_options.convTolResidual.check = false; // (default: false)
redist_options.interval_check_convergence = 1; // 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.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)
redist_options.interval_check_convergence = 1000; // interval of #iterations at which
// convergence is checked.
redist_options.width_NB_in_grid_points = 8; // width of narrow band in number of grid points.
redist_options.print_current_iterChangeResidual = true; // if true, prints out every current iteration + corresponding change from the previous iteration + residual from SDF
redist_options.print_steadyState_iter = true; // if true, prints out the final iteration number when steady state was reached + final change + residual
redist_options.save_temp_grid = true;
RedistancingSussman<grid_in_type> redist_obj(g_dist, redist_options); // Instantiation of Sussman-redistancing class
redist_obj.set_user_time_step(dt);
......@@ -74,12 +95,11 @@ BOOST_AUTO_TEST_SUITE(RedistancingSussmanTestSuite)
// 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);
/////////////////////////////////////////////////////////////////////////////////////////////
// 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};
size_t bc[grid_dim] = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
typedef aggregate<double> props_nb;
typedef vector_dist<grid_dim, double, props_nb> vd_type;
Ghost<grid_dim, double> ghost_vd(0);
......@@ -93,20 +113,20 @@ BOOST_AUTO_TEST_SUITE(RedistancingSussmanTestSuite)
// Compute the L_2- and L_infinity-norm and save to file
L_norms lNorms_vd;
lNorms_vd = get_l_norms_vector<0>(vd_narrow_band);
std::cout << lNorms_vd.l2 << ", " << lNorms_vd.linf << std::endl;
std::cout << "L2:" << lNorms_vd.l2 << ", Linf:" << lNorms_vd.linf << std::endl;
// l-norms original implementation of 1st order upwinding:
// 32,0.033643,0.063065; dt = 0.000165334; 1e4 iterations
// 0.0336846, 0.0630651
// new impl
// N = 32
// dt = 0.000165334;
// 1e4 iterations
// order 1: 0.0336846, 0.0630651
// order 3: 0.02794, 0.0586704
// order 3: 0.02794, 0.0586704
// order 5: 0.0187199, 0.0367638
BOOST_CHECK(lNorms_vd.l2 < 0.03369 + EPSILON);
BOOST_CHECK(lNorms_vd.linf < 0.06307 + EPSILON);
// BOOST_CHECK(narrow_band_width > 0);
#endif
}
BOOST_AUTO_TEST_SUITE_END()
......
//
// Created by jstark on 12.07.21.
//
#define SE_CLASS1
#define BOOST_TEST_DYN_LINK
//#define BOOST_TEST_MAIN // in only one cpp file
#include <boost/test/unit_test.hpp>
// Include redistancing files
#include "util/PathsAndFiles.hpp"
#include "level_set/redistancing_Sussman/RedistancingSussman.hpp"
#include "level_set/redistancing_Sussman/NarrowBand.hpp"
// Include header files for testing
#include "Draw/DrawSphere.hpp"
#include "Draw/DrawDisk.hpp"
#include "l_norms/LNorms.hpp"
#include "analytical_SDF/AnalyticalSDF.hpp"
BOOST_AUTO_TEST_SUITE(RedistancingSussmanSinglePointConvergenceTestSuite)
#if 0
BOOST_AUTO_TEST_CASE(RedistancingSussmanSinglePoint_1D_test)
{
// 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
const size_t x = 0;
const size_t Phi_0_grid = 0;
const size_t SDF_sussman_grid = 1;
const size_t SDF_exact_grid = 2;
const size_t Error_grid = 3;
size_t N = 64;
double dt = 0.000503905;
double tf = 1e4 * dt;
dt *= 2.0;
// for (int i = 0; i < 3; i++)
int i = 0;
{
dt /= 2.0;
size_t max_iter = (size_t) std::round(tf / dt);
const size_t sz[grid_dim] = {N};
const double box_lower = -1.0;
const double box_upper = 1.0;
Box<grid_dim, double> 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;
grid_in_type g_dist(sz, box, ghost);
g_dist.setPropNames({"Phi_0", "SDF_sussman", "SDF_exact", "Relative error"});
// Now we initialize the grid with the indicator function and the analytical solution
auto dom = g_dist.getDomainIterator();
while(dom.isNext())
{
auto key = dom.get();
Point<grid_dim, double> 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;
}
g_dist.write("grid_1D_preRedistancing", FORMAT_BINARY); // Save the grid
Redist_options redist_options;
// redist_options.min_iter = max_iter;
// redist_options.max_iter = max_iter;
redist_options.order_space_op = 1;
redist_options.order_timestepper = 1;
// set both convergence criteria to false s.t. termination only when max_iterations reached
redist_options.convTolChange.check = true; // 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
redist_options.convTolChange.value = EPSILON;
redist_options.convTolResidual.check = false; // (default: false)
redist_options.interval_check_convergence = 100; // interval of #iterations at which
// convergence is checked (default: 100)
redist_options.width_NB_in_grid_points = 32; // 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)
redist_options.coord_checkpoint = 10;
RedistancingSussman<grid_in_type> redist_obj(g_dist,
redist_options); // Instantiation of Sussman-redistancing class
redist_obj.set_user_time_step(dt);
// std::cout << "CFL dt = " << redist_obj.get_time_step() << std::endl;
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>();
// 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_RKorder" + std::to_string(redist_options.order_timestepper) + "_i_" + std::to_string(i),
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);
// vd_type vd_narrow_band(0, box, bc, ghost_vd);
// 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 = 64;
// 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) + "_RKorder" + std::to_string(redist_options.order_timestepper),
// 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;
// lNorms_vd = get_l_norms_vector<Error_vd>(vd_narrow_band);