HelpFunctionsForGrid.hpp

Go to the documentation of this file.

//
// Created by jstark on 2020-05-12.
//
#ifndef REDISTANCING_SUSSMAN_HELPFUNCTIONSFORGRID_HPP
#define REDISTANCING_SUSSMAN_HELPFUNCTIONSFORGRID_HPP

#include <iostream>
#include <limits>

#include "HelpFunctions.hpp"
#include "VCluster/VCluster.hpp"
#include "Grid/grid_dist_id.hpp"

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

template <typename grid_type>
typename grid_type::stype get_time_step_CFL(grid_type & grid, typename grid_type::stype u [grid_type::dims], typename grid_type::stype C)
{
    typename grid_type::stype sum = 0;
    for (size_t d = 0; d < grid_type::dims; d++)
    {
        sum += u[d] / grid.spacing(d);
    }
    return C / sum;
}
template <typename grid_type>
typename grid_type::stype get_time_step_CFL(grid_type & grid, typename grid_type::stype u, typename grid_type::stype C)
{
    typename grid_type::stype sum = 0;
    for (size_t d = 0; d < grid_type::dims; d++)
    {
        sum += u / grid.spacing(d);
    }
    return C / sum;
}
#endif

template <size_t Prop, typename grid_type, typename T>
void init_grid_and_ghost(grid_type & grid, T init_value)
{
    auto dom_ghost = grid.getDomainGhostIterator();
    while(dom_ghost.isNext())
    {
        auto key = dom_ghost.get();
        grid.template get<Prop>(key) = init_value;
        ++dom_ghost;
    }
}

template <size_t attr_sc, size_t attr_ds, typename grid_source_type, typename grid_dest_type>
void copy_gridTogrid(const grid_source_type & grid_sc, grid_dest_type & grid_ds, bool include_ghost=false)
{
    assert(grid_source_type::dims == grid_dest_type::dims);
    assert(grid_sc.size() == grid_ds.size());
    if (include_ghost)
    {
        auto dom_sc = grid_sc.getDomainGhostIterator();
        auto dom_ds = grid_ds.getDomainGhostIterator();
        while (dom_sc.isNext())
        {
            auto key_sc = dom_sc.get();
            auto key_ds = dom_ds.get();
            grid_ds.template get<attr_ds>(key_ds) = grid_sc.template get<attr_sc>(key_sc);
            ++dom_sc;
            ++dom_ds;
        }
    }
    else
    {
        auto dom_sc = grid_sc.getDomainIterator();
        auto dom_ds = grid_ds.getDomainIterator();
        while (dom_sc.isNext())
        {
            auto key_sc = dom_sc.get();
            auto key_ds = dom_ds.get();
            grid_ds.template get<attr_ds>(key_ds) = grid_sc.template get<attr_sc>(key_sc);
            ++dom_sc;
            ++dom_ds;
        }
    }
}

template <typename grid_type>
typename grid_type::stype get_biggest_spacing(grid_type & grid)
{
    typename grid_type::stype h_max = 0;
    for (size_t d = 0; d < grid_type::dims; d++)
    {
        if (grid.spacing(d) > h_max) h_max = grid.spacing(d);
    }
    return h_max;
}

template <typename grid_type>
typename grid_type::stype get_smallest_spacing(grid_type & grid)
{
    typename grid_type::stype spacing [grid_type::dims];
    for (size_t d = 0; d < grid_type::dims; d++)
    {
        spacing[d] = grid.spacing(d);
    }
    std::sort(std::begin(spacing), std::end(spacing));
    return spacing[0];
}

template <size_t Prop1, size_t Prop2, typename grid_type>
auto average_difference(grid_type & grid)
{
    auto dom = grid.getDomainIterator();
    typedef typename std::decay_t<decltype(grid.template get<Prop1>(dom.get()))> prop_type;
    prop_type total_diff = 0;
    while (dom.isNext())
    {
        auto key = dom.get();
        total_diff += abs( grid.template get<Prop1>(key) - grid.template get<Prop2>(key) );
        ++dom;
    }
    return total_diff / grid.size();
}

template <size_t Prop, typename grid_type>
auto get_max_val(grid_type & grid)
{
    auto dom = grid.getDomainIterator();
    typedef typename std::decay_t<decltype(grid.template get<Prop>(dom.get()))> prop_type;
    prop_type max_value = std::numeric_limits<prop_type>::lowest();
    while(dom.isNext())
    {
        auto key = dom.get();
        if (grid.template get<Prop>(key) > max_value) max_value = grid.template get<Prop>(key);
        ++dom;
    }
    auto & v_cl = create_vcluster();
    v_cl.max(max_value);
    v_cl.execute();
    return max_value;
}

template <size_t Prop, typename grid_type>
auto get_min_val(grid_type & grid)
{
    auto dom = grid.getDomainIterator();
    typedef typename std::decay_t<decltype(grid.template get<Prop>(dom.get()))> prop_type;
    prop_type min_value = std::numeric_limits<prop_type>::max();
    while(dom.isNext())
    {
        auto key = dom.get();
        if (grid.template get<Prop>(key) < min_value) min_value = grid.template get<Prop>(key);
        ++dom;
    }
    auto & v_cl = create_vcluster();
    v_cl.min(min_value);
    v_cl.execute();
    return min_value;
}

template <size_t Prop_in, size_t Prop_out, typename grid_type>
void init_sign_prop(grid_type & grid)
{
    auto dom = grid.getDomainIterator();
    while(dom.isNext())
    {
        auto key = dom.get();
        grid.template get<Prop_out>(key) = sgn(grid.template get<Prop_in>(key));
        ++dom;
    }
}


template <size_t Vector_in, size_t Magnitude_out, typename magnitude_type, typename gridtype>
void get_vector_magnitude(gridtype & grid)
{
    grid.template ghost_get<Vector_in>();
    auto dom = grid.getDomainGhostIterator();
    while(dom.isNext())
    {
        magnitude_type sum = 0;
        auto key = dom.get();
        for(size_t d = 0; d < gridtype::dims; d++)
        {
            sum += grid.template get<Vector_in> (key)[d] * grid.template get<Vector_in> (key)[d];
        }
        grid.template get<Magnitude_out> (key) = sqrt(sum);
        ++dom;
    }
}

template <size_t Vector_in, typename key_type, typename gridtype>
auto get_vector_magnitude(gridtype & grid, key_type & key)
{
    typedef typename std::remove_const_t<std::remove_reference_t<decltype(grid.template get<Vector_in>(key)[0])>>
    prop_type;
    prop_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