LNorms.hpp

Go to the documentation of this file.

//
// Created by jstark on 2020-05-17.
//
#ifndef ACCURACY_TESTS_LNORMS_HPP
#define ACCURACY_TESTS_LNORMS_HPP
 
#include <iostream>
#include <typeinfo>
#include <cmath>
#include <cstdio>
#include <stdlib.h>
 
// Include OpenFPM header files
#include "Vector/vector_dist.hpp"
#include "Grid/grid_dist_id.hpp"
#include "util/PathsAndFiles.hpp"
 
#include "level_set/redistancing_Sussman/HelpFunctions.hpp" // for printing to_string_with_precision
 
template <size_t PropNumeric, size_t PropAnalytic, size_t Error, typename gridtype>
void get_absolute_error(gridtype & grid)
{
    auto dom = grid.getDomainIterator();
    typedef typename std::remove_const_t<std::remove_reference_t<decltype(
    grid.template get<PropNumeric>(dom.get()))>> numeric_type;
    typedef typename std::remove_const_t<std::remove_reference_t<decltype(
    grid.template get<PropAnalytic>(dom.get()))>> analytic_type;
    
    if(!(std::is_same<numeric_type, analytic_type>::value))
    {
        std::cout << "Type of numerical and analytical solution must be the same! Aborting..." << std::endl;
        abort();
    }   while(dom.isNext())
    {
        auto key = dom.get();
        grid.template getProp<Error> (key) = abs(grid.template getProp<PropAnalytic> (key) - (grid.template getProp<PropNumeric> (key)));
        ++dom;
    }
}
template <size_t PropNumeric, size_t PropAnalytic, size_t Error, typename gridtype>
void get_relative_error(gridtype & grid)
{
    auto dom = grid.getDomainIterator();
    typedef typename std::remove_const_t<std::remove_reference_t<decltype(
    grid.template get<PropNumeric>(dom.get()))>> numeric_type;
    typedef typename std::remove_const_t<std::remove_reference_t<decltype(
    grid.template get<PropAnalytic>(dom.get()))>> analytic_type;
    
    if(!(std::is_same<numeric_type, analytic_type>::value))
    {
        std::cout << "Type of numerical and analytical solution must be the same! Aborting..." << std::endl;
        abort();
    }
    
    while(dom.isNext())
    {
        auto key = dom.get();
        grid.template getProp<Error> (key) =
                abs((grid.template getProp<PropNumeric> (key) - grid.template getProp<PropAnalytic> (key)) /
                    (grid.template getProp<PropAnalytic> (key) + std::numeric_limits<analytic_type>::epsilon()));
        ++dom;
    }
}
 
template <typename lnorm_type=double>
class LNorms
{
public:
    LNorms()
    {
        l2 = 0;
        linf = 0;
    }
    // Member variables
    lnorm_type l2; // L2 norm
    lnorm_type linf; // L_infinity norm
    
    // Member functions
    template <size_t Error, typename gridtype>
    void get_l_norms_grid(gridtype & grid)
    {
        auto dom = grid.getDomainIterator();
        typedef typename std::remove_const_t<std::remove_reference_t<decltype(
        grid.template get<Error>(dom.get()))>> error_type;
        
        error_type maxError = 0;
        error_type sumErrorSq = 0;
        
        while(dom.isNext())
        {
            auto key = dom.get();
            sumErrorSq += grid.template getProp<Error> (key) * grid.template getProp<Error> (key);
            if (grid.template getProp<Error> (key) > maxError) maxError = grid.template getProp<Error> (key); // update maxError
            ++dom;
        }
        auto &v_cl = create_vcluster();
        v_cl.sum(sumErrorSq);
        v_cl.max(maxError);
        v_cl.execute();
        l2 = (lnorm_type) sqrt( sumErrorSq / (error_type)grid.size());
        linf = (lnorm_type) maxError;
    }
    
    template <size_t Error, typename vectortype>
    void get_l_norms_vector(vectortype & vd)
    {
        auto dom = vd.getDomainIterator();
        typedef typename std::remove_const_t<std::remove_reference_t<decltype(
        vd.template getProp<Error>(dom.get()))>> error_type;
        
        error_type maxError = 0;
        error_type sumErrorSq = 0;
        int count_particles = 0;
        while(dom.isNext())
        {
            auto key = dom.get();
            sumErrorSq += vd.template getProp<Error> (key) * vd.template getProp<Error> (key);
            if (vd.template getProp<Error> (key) > maxError) maxError = vd.template getProp<Error> (key); // update maxError
            ++dom;
            ++count_particles;
        }
        auto &v_cl = create_vcluster();
        v_cl.sum(sumErrorSq);
        v_cl.sum(count_particles);
        v_cl.max(maxError);
        v_cl.execute();
        l2 = (lnorm_type) sqrt( sumErrorSq / (error_type)count_particles);
        linf = (lnorm_type) maxError;
    }
    
    template <size_t Error, typename grid_type>
    void get_l_norms(grid_type & grid)
    {
        auto dom = grid.getDomainIterator();
        typedef typename std::remove_const_t<std::remove_reference_t<decltype(
        grid.template getProp<Error>(dom.get()))>> error_type;
        
        error_type maxError = 0;
        error_type sumErrorSq = 0;
        int count_points = 0;
        while(dom.isNext())
        {
            auto key = dom.get();
            sumErrorSq += grid.template getProp<Error> (key) * grid.template getProp<Error> (key);
            if (grid.template getProp<Error> (key) > maxError) maxError = grid.template getProp<Error> (key); // update maxError
            ++dom;
            ++count_points;
        }
        auto &v_cl = create_vcluster();
        v_cl.sum(sumErrorSq);
        v_cl.sum(count_points);
        v_cl.max(maxError);
        v_cl.execute();
        l2 = (lnorm_type) sqrt( sumErrorSq / (error_type)count_points);
        linf = (lnorm_type) maxError;
    }
    
    
     template <typename T>
     void write_to_file(const T col_0,
                        const int precision,
                        const std::string & filename,
                        const std::string & path_output)
    {
        auto &v_cl = create_vcluster();
        if (v_cl.rank() == 0)
        {
            std::string path_output_lnorm = path_output + "/" + filename + ".csv";
            create_file_if_not_exist(path_output_lnorm);
            
            std::ofstream l_out;
            l_out.open(path_output_lnorm, std::ios_base::app); // append instead of overwrite
            l_out << std::to_string(col_0)
                    << ',' << to_string_with_precision(l2, precision)
                    << ',' << to_string_with_precision(linf, precision) << std::endl;
            l_out.close();
        }
    }
 
private:
 
 
};
 
 
 
 
 
 
 
 
 
#endif //ACCURACY_TESTS_LNORMS_HPP