Skip to content
Snippets Groups Projects
Select Git revision
  • 9d4f38b4c5de279cfd14a349dab13e323de4afc4
  • master default
  • cgnsUnstructured
  • partitioning
  • poppler
  • HighOrderBLCurving
  • gmsh_3_0_4
  • gmsh_3_0_3
  • gmsh_3_0_2
  • gmsh_3_0_1
  • gmsh_3_0_0
  • gmsh_2_16_0
  • gmsh_2_15_0
  • gmsh_2_14_1
  • gmsh_2_14_0
  • gmsh_2_13_2
  • gmsh_2_13_1
  • gmsh_2_12_0
  • gmsh_2_11_0
  • gmsh_2_10_1
  • gmsh_2_10_0
  • gmsh_2_9_3
  • gmsh_2_9_2
  • gmsh_2_9_1
  • gmsh_2_9_0
  • gmsh_2_8_6
26 results

Callbacks.cpp

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    • Christophe Geuzaine's avatar
      9d4f38b4
      · 9d4f38b4
      Christophe Geuzaine authored
      - removed the sort view list stuff: it's a bad idea since we re-sort
      the list in many places when we search/replace through it
      
      - add a simple "set visibility->all on/off" switch
      9d4f38b4
      History
      Christophe Geuzaine authored
      - removed the sort view list stuff: it's a bad idea since we re-sort
      the list in many places when we search/replace through it
      
      - add a simple "set visibility->all on/off" switch
    finite_difference_baseline.cpp 5.86 KiB
    #include <iostream>
    #include <fstream>
    #include <sstream>
    #include <sys/time.h>
    #include <vector>
    #include <cmath>
    #include <chrono>
    
    #include <silo.h>
    
    template<typename T>
    class fd_grid
    {
        std::vector<T>  m_data;
        size_t m_rows, m_cols;
    
    public:
        fd_grid()
        {}
    
        fd_grid(size_t rows, size_t cols) :
            m_data(rows*cols), m_rows(rows), m_cols(cols)
        {}
    
        fd_grid(fd_grid&& other)
        {
            std::swap(m_data, other.m_data);
        }
    
        fd_grid& operator=(fd_grid&& other)
        {
            std::swap(m_data, other.m_data);
            return *this;
        }
    
        T operator()(size_t i, size_t j) const
        {
            return m_data[i*m_cols + j];
        }
    
        T& operator()(size_t i, size_t j)
        {
            return m_data[i*m_cols + j];
        }
    
        size_t rows() const { return m_rows; }
        size_t cols() const { return m_cols; }
    
        T* data() { return m_data.data(); }
        const T* data() const { return m_data.data(); }
    };
    
    int
    visit_dump(const fd_grid<double>& g, const std::string& fn)
    {
        DBfile *db = nullptr;
        db = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL, "Kokkos test", DB_HDF5);
        if (!db)
        {
            std::cout << "Cannot write simulation output" << std::endl;
            return -1;
        }
    
        auto size_x = g.cols();
        auto size_y = g.rows();
        
        std::vector<double> x(size_x);
        std::vector<double> y(size_y);
        
        for (size_t i = 0; i < size_x; i++)
            x.at(i) = double(i)/(size_x-1);
            
        for (size_t i = 0; i < size_y; i++)
            y.at(i) = double(i)/(size_y-1);
            
        int dims[] = { int(size_y), int(size_y) };
        int ndims = 2;
        double *coords[] = {x.data(), y.data()};
        
        DBPutQuadmesh(db, "mesh", NULL, coords, dims, ndims, DB_DOUBLE, DB_COLLINEAR, NULL);
            
        DBPutQuadvar1(db, "solution", "mesh", g.data(), dims, ndims, NULL, 0, DB_DOUBLE, DB_NODECENT, NULL);
        
        DBClose(db);
        return 0;
    }
    
    #define HIGH_ORDER
    
    #ifdef HIGH_ORDER
        #define HALO_SIZE 4
    #else
        #define HALO_SIZE 1
    #endif
    
    int main( int argc, char* argv[] )
    {
        using namespace std::chrono;
    
        int     Nx          = 500;
        int     Ny          = 500;
        int     halo        = HALO_SIZE;
    
        int     rows        = Ny + 2*halo;
        int     cols        = Nx + 2*halo;
        int     timesteps   = 5000;
        double  dx          = 1./(Nx-1);
        double  dy          = 1./(Ny-1);
        double  dt          = 0.001;
    
        double  velocity    = 1;
        double  damping     = 0.1;
    
        auto c2 = velocity*velocity;
        auto dx2 = dx*dx;
        auto dy2 = dy*dy;
        auto dt2 = dt*dt;
    
        auto D1 = 1./(1. + 0.5*damping*dt);
        auto D2 = 0.5*damping*dt - 1;
    
        typedef fd_grid<double>  matrix_view_type;
    
        matrix_view_type u_prev(rows, cols);
        matrix_view_type u_curr(rows, cols);
        matrix_view_type u_next(rows, cols);
    
        for (int i = 0; i < rows; i++)
        {
            for (int j = 0; j < cols; j++)
            {
                if ( i < halo or i > rows-halo or j < halo or j > cols-halo)
                {
                    u_prev(i,j) = 0.0;
                    u_curr(i,j) = 0.0;
                }
                else
                {
                    double x = dx*i - 0.5;
                    double y = dy*j - 0.1;
                    u_prev(i,j) = -std::exp(-2400*(x*x + y*y));
                    u_curr(i,j) = 2*dt*u_prev(i,j);
                }
            }
        }
    
        double iter_time = 0.0;
    
        double w0 = -205.0/72.0;
        double w1 = 8.0/5.0;
        double w2 = -1.0/5.0;
        double w3 = 8.0/315.0;
        double w4 = -1.0/560.0;
        double w[9] = { w4, w3, w2, w1, w0, w1, w2, w3, w4 };
    
        for (int ts = 0; ts < timesteps; ts++)
        {
            auto t_begin = high_resolution_clock::now();
    
            for (size_t i = halo; i < rows - halo; i++)
            {
                for (size_t j = halo; j < cols - halo; j++)
                {
                    /* 2nd order in space */
    #ifndef HIGH_ORDER
                    u_next(i,j) = D1*(2.0*u_curr(i,j) + D2*u_prev(i,j) +
                            + (c2*dt2/dy2)*(u_curr(i+1,j) - 2.0*u_curr(i,j) + u_curr(i-1,j))
                            + (c2*dt2/dx2)*(u_curr(i,j+1) - 2.0*u_curr(i,j) + u_curr(i,j-1)));
    
                    if (j == 1)
                        u_next(i,j-1) = 0;//u_next(i,j);
    
                    if (j == cols-2)
                        u_next(i,j+1) = 0;//u_next(i,j);
    
                    if (i == 1)
                        u_next(i-1,j) = 0;//u_next(i,j);
    
                    if (i == rows-2)
                        u_next(i+1,j) = 0;//u_next(i,j);
    #else
    
                    /* 8th order in space */
                    double deriv_x = 0.0;
                    double deriv_y = 0.0;
                    for (int k = -4; k <= 4; k++)
                    {
                        deriv_x += w[k+4]*u_curr(i+k,j);
                        deriv_y += w[k+4]*u_curr(i,j+k);
                    }
                    u_next(i,j) = D1*(2.0*u_curr(i,j) + D2*u_prev(i,j) +
                            + (c2*dt2/dy2)*(deriv_y)
                            + (c2*dt2/dx2)*(deriv_x));
    
                    if (j == halo)
                        for (size_t k = 0; k < halo; k++)
                            u_next(i,k) = 0;
    
                    if (j == cols-(halo+1))
                        for (size_t k = 0; k < halo; k++)
                            u_next(i,j+k+1) = 0;
    
                    if (i == 1)
                        for (size_t k = 0; k < halo; k++)
                            u_next(k,j) = 0;
    
                    if (i == rows-(halo+1))
                        for (size_t k = 0; k < halo; k++)
                            u_next(i+k+1,j) = 0;
    #endif
                }
            }
    
            std::swap(u_prev, u_curr);
            std::swap(u_curr, u_next);
    
            auto t_end = high_resolution_clock::now();
    
            duration<double> time_span = duration_cast<duration<double>>(t_end - t_begin);
            iter_time += time_span.count();
    
            /*
            if ( (ts % 100) == 0 )
            {
                std::stringstream ss;
                ss << "solution_" << ts << ".silo";
                visit_dump(u_curr, ss.str());
            }
            */
        }
    
        std::cout << "Average iteration time: " << iter_time/timesteps << std::endl;
    
        return 0;
    }