RedistancingSussman.hpp

Go to the documentation of this file.

//
// Created by jstark on 2020-04-06.
//
#ifndef REDISTANCING_SUSSMAN_REDISTANCINGSUSSMAN_HPP
#define REDISTANCING_SUSSMAN_REDISTANCINGSUSSMAN_HPP

// Include standard library header files
#include <iostream>
#include <string>
#include <fstream>
// Include OpenFPM header files
#include "Vector/vector_dist.hpp"
#include "Grid/grid_dist_id.hpp"
#include "data_type/aggregate.hpp"
#include "Decomposition/CartDecomposition.hpp"

// Include other header files
#include "HelpFunctions.hpp"
#include "HelpFunctionsForGrid.hpp"
//#include "ComputeGradient.hpp"
#include "FiniteDifference/Upwind_gradient.hpp"

template <typename phi_type=double>
struct Conv_tol_change
{
    bool check = true;  
    phi_type value = 1e-15;     
};

template <typename phi_type=double>
struct Conv_tol_residual
{
    bool check = true; 
    phi_type value = 1e-3;      
};

template <typename phi_type=double>
struct Redist_options
{
    size_t min_iter = 1e5;
    size_t max_iter = 1e12;
    
    Conv_tol_change<phi_type> convTolChange;
    Conv_tol_residual<phi_type> convTolResidual;
    
    size_t interval_check_convergence = 100;
    size_t width_NB_in_grid_points = 8;
    bool print_current_iterChangeResidual = false;
    bool print_steadyState_iter = true;
    bool save_temp_grid = false;
};

template <typename phi_type=double>
struct DistFromSol
{
    phi_type change;    
    phi_type residual;  
    int count; 
};


template <typename grid_in_type, typename phi_type=double>
class RedistancingSussman
{
public:
    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(),
                                                                                       Ghost<grid_in_type::dims, long int>(3))
    {
        time_step = get_time_step_CFL(grid_in);
        order_upwind_gradient = 1;
#ifdef SE_CLASS1
        assure_minimal_thickness_of_NB();
#endif // SE_CLASS1
    }
    
    typedef aggregate<phi_type, phi_type[grid_in_type::dims], int>
            props_temp;
    typedef grid_dist_id<grid_in_type::dims, typename grid_in_type::stype, props_temp> g_temp_type;
    g_temp_type g_temp;
    
    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
        // 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
        copy_gridTogrid<Phi_n_temp, Phi_SDF_out>(g_temp, r_grid_in); // Copy resulting SDF function to input grid
    }
    template<typename T>
    void set_user_time_step(T dt)
    {
        time_step = dt;
    }
    auto get_time_step()
    {
        return time_step;
    }
    
    int get_finalIteration()
    {
        return final_iter;
    }
    
    auto get_finalChange()
    {
        return distFromSol.change;
    }
    
    auto get_finalResidual()
    {
        return distFromSol.residual;
    }
    
    int get_finalNumberNbPoints()
    {
        return distFromSol.count;
    }
    
private:
    //  Some indices for better readability
    static constexpr size_t Phi_n_temp          = 0; 
    static constexpr size_t Phi_grad_temp       = 1; 
    static constexpr size_t Phi_0_sign_temp     = 2; 
    
    //  Member variables
    Redist_options<phi_type> redistOptions; 
    grid_in_type &r_grid_in; 
    
    DistFromSol<phi_type> distFromSol; 
    // point in NB.
    int final_iter = 0; 
    
    phi_type kappa = ceil(redistOptions.width_NB_in_grid_points / 2.0) * get_biggest_spacing(g_temp);
    typename grid_in_type::stype time_step;
    int order_upwind_gradient;
    //  Member functions
#ifdef SE_CLASS1

    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
    }
#endif // SE_CLASS1

    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
        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
    }

    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);
    }
    
    void go_one_redistancing_step_whole_grid(g_temp_type &grid)
    {
        get_upwind_gradient<Phi_n_temp, Phi_0_sign_temp, Phi_grad_temp>(grid, order_upwind_gradient, true);
        grid.template ghost_get<Phi_n_temp, Phi_grad_temp>();
        auto dom = grid.getDomainIterator();
        while (dom.isNext())
        {
            auto key = dom.get();
            const phi_type phi_n = grid.template get<Phi_n_temp>(key);
            const phi_type phi_n_magnOfGrad = get_vector_magnitude<Phi_grad_temp>(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));
            ++dom;
        }
    }

    bool lays_inside_NB(phi_type Phi)
    {
        return (abs(Phi) <= kappa);
    }
    
    void update_distFromSol(g_temp_type &grid)
    {
        phi_type max_residual = 0;
        phi_type max_change = 0;
        int count = 0;
        auto dom = grid.getDomainIterator();
        while (dom.isNext())
        {
            auto key = dom.get();
            if (lays_inside_NB(grid.template get<Phi_n_temp>(key)))
            {
                count++;
                phi_type phi_n_magnOfGrad = get_vector_magnitude<Phi_grad_temp>(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));
                
                if (abs(phi_nplus1 - grid.template get<Phi_n_temp>(key)) > max_change)
                {
                    max_change = abs(phi_nplus1 - grid.template get<Phi_n_temp>(key));
                }
                
                if (abs(phi_n_magnOfGrad - 1) > max_residual) { max_residual = abs(phi_n_magnOfGrad - 1); }
            }
            ++dom;
        }
        auto &v_cl = create_vcluster();
        v_cl.max(max_change);
        v_cl.max(max_residual);
        v_cl.sum(count);
        v_cl.execute();
        
        // Update member variable distFromSol
        distFromSol.change   = max_change;
        distFromSol.residual = max_residual;
        distFromSol.count    = count;
    }
    
    void print_out_iteration_change_residual(g_temp_type &grid, size_t iter)
    {
        update_distFromSol(grid);
        auto &v_cl = create_vcluster();
        if (v_cl.rank() == 0)
        {
            if (iter == 0)
            {
                std::cout << "Iteration,MaxChange,MaxResidual,NumberOfNarrowBandPoints" << std::endl;
            }
            std::cout << iter
            << "," << to_string_with_precision(distFromSol.change, 15)
            << "," << to_string_with_precision(distFromSol.residual, 15)
            << "," << distFromSol.count << std::endl;
        }
    }
    
    bool steady_state_NB(g_temp_type &grid)
    {
        bool steady_state = false;
        update_distFromSol(grid);
        if (redistOptions.convTolChange.check && redistOptions.convTolResidual.check)
        {
            steady_state = (
                    distFromSol.change <= redistOptions.convTolChange.value &&
                    distFromSol.residual <= redistOptions.convTolResidual.value &&
                    distFromSol.count > 0
                    );
        }
        else
        {
            if (redistOptions.convTolChange.check)
            {
                steady_state = (distFromSol.change <= redistOptions.convTolChange.value && distFromSol.count > 0);
            }       // Use the normalized total change between two iterations in the narrow bands steady-state criterion
            if (redistOptions.convTolResidual.check)
            {
                steady_state = (distFromSol.residual <= redistOptions.convTolResidual.value && distFromSol.count > 0);
            } // Use the normalized total residual of phi compared to SDF in the narrow bands steady-state criterion
        }
        return steady_state;
    }
    
    void iterative_redistancing(g_temp_type &grid)
    {
        int i = 0;
        while (i < redistOptions.max_iter)
        {
            for (int j = 0; j < redistOptions.interval_check_convergence; j++)
            {
                go_one_redistancing_step_whole_grid(grid);
                ++i;
            }
            if (redistOptions.print_current_iterChangeResidual)
            {
                print_out_iteration_change_residual(grid, i);
            }
            if (redistOptions.save_temp_grid)
            {
                get_upwind_gradient<Phi_n_temp, Phi_0_sign_temp, Phi_grad_temp>(g_temp, order_upwind_gradient, true);
                g_temp.setPropNames({"Phi_Sussman_Out", "Phi_upwind_gradient", "Phi_0_sign_temp"});
                g_temp.save("g_temp_redistancing_iteration_" + std::to_string(i) + ".hdf5"); // HDF5 file
                // g_temp.write_frame("g_temp_redistancing_iteration", i, FORMAT_BINARY); // VTK file
            }
            if (i >= redistOptions.min_iter)
            {
                if (steady_state_NB(grid))
                {
                    if (redistOptions.print_steadyState_iter)
                    {
                        auto &v_cl = create_vcluster();
                        if (v_cl.rank() == 0)
                        {
                            std::cout << "Steady state criterion reached at iteration: " << i << std::endl;
                            std::cout << "FinalIteration,MaxChange,MaxResidual" << std::endl;
                        }
                        print_out_iteration_change_residual(grid, i);
                    }
                    break;
                }
            }
        }
        update_distFromSol(grid);
        final_iter = i;
        // If save_temp_grid set true, save the temporary grid as hdf5 that can be reloaded onto a grid
        if (redistOptions.save_temp_grid)
        {
            get_upwind_gradient<Phi_n_temp, Phi_0_sign_temp, Phi_grad_temp>(g_temp, order_upwind_gradient, true);
            g_temp.setPropNames({"Phi_Sussman_Out", "Phi_upwind_gradient", "Phi_0_sign_temp"});
            g_temp.save("g_temp_redistancing_final.hdf5"); // HDF5 file
            // g_temp.write("g_temp_redistancing_final", FORMAT_BINARY); // VTK file
        }
    }
};

#endif //REDISTANCING_SUSSMAN_REDISTANCINGSUSSMAN_HPP