Commit 6c6b3e43 authored by Pietro Incardona's avatar Pietro Incardona

Finite differences it scale

parent 834b07c4
/*
* eq_unit_tests.hpp
*
* Created on: Sep 18, 2015
* Author: i-bird
*/
#ifndef OPENFPM_NUMERICS_SRC_EQUATIONS_EQ_UNIT_TESTS_HPP_
#define OPENFPM_NUMERICS_SRC_EQUATIONS_EQ_UNIT_TESTS_HPP_
#include "FiniteDifference/FDScheme.hpp"
BOOST_AUTO_TEST_SUITE( eq_test )
struct Sys_p
{
static const unsigned int dims = 3;
typedef float stype;
};
// define coordinates id
constexpr unsigned int x = 0;
constexpr unsigned int y = 1;
// define constants
constexpr unsigned int eta = 0;
constexpr unsigned int zeta = 1;
constexpr unsigned int gamma = 2;
constexpr unsigned int nu = 3;
constexpr unsigned int lambda = 4;
constexpr unsigned int dmu = 5;
constexpr unsigned int sigma[2][2] = {{6,7},{8,9}};
constexpr unsigned int p[] = {10,11};
constexpr unsigned int del_nu = 12;
constexpr unsigned int g_ext = 13;
// define variables fields
// velocity
constexpr unsigned int V[2] = {0,1};
// pressure
constexpr unsigned int P = 2;
// model Eq1 active Gel
//// First level expression
/*
*
* p_x^2
*
*/
//typedef mul<p[x],p[x],Sys_p> px2;
/*
*
* p_y^2
*
*/
//typedef mul<p[y],p[y],Sys_p> py2;
/*
* p_x * p_y
*
*/
//typedef mul<p[x],p[y],Sys_p> px_py;
BOOST_AUTO_TEST_CASE( eq_test_use)
{
// create an equation object
}
BOOST_AUTO_TEST_SUITE_END()
#endif /* OPENFPM_NUMERICS_SRC_EQUATIONS_EQ_UNIT_TESTS_HPP_ */
......@@ -87,7 +87,7 @@ class Avg<d,arg,Sys_eqs,CENTRAL>
* \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
*
*/
inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
// if the system is staggered the CENTRAL derivative is equivalent to a forward derivative
if (is_grid_staggered<Sys_eqs>::value())
......@@ -170,7 +170,7 @@ class Avg<d,arg,Sys_eqs,FORWARD>
* \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
*
*/
inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
long int old_val = kmap.getKeyRef().get(d);
......@@ -230,7 +230,7 @@ class Avg<d,arg,Sys_eqs,BACKWARD>
* \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
*
*/
inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims], std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims], std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
long int old_val = kmap.getKeyRef().get(d);
kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 1);
......
......@@ -90,7 +90,7 @@ class D<d,arg,Sys_eqs,CENTRAL>
* \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
*
*/
inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
// if the system is staggered the CENTRAL derivative is equivalent to a forward derivative
if (is_grid_staggered<Sys_eqs>::value())
......@@ -205,7 +205,7 @@ public:
* \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
*
*/
static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
#ifdef SE_CLASS1
if (Sys_eqs::boundary[d] == PERIODIC)
......@@ -318,7 +318,7 @@ class D<d,arg,Sys_eqs,FORWARD>
* \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
*
*/
inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] ,std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] ,std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
long int old_val = kmap.getKeyRef().get(d);
......@@ -374,7 +374,7 @@ class D<d,arg,Sys_eqs,BACKWARD>
* \snippet FDScheme_unit_tests.hpp Usage of stencil derivative
*
*/
inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims], std::unordered_map<long int,typename Sys_eqs::stype > & cols , typename Sys_eqs::stype coeff)
inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims], std::unordered_map<long int,typename Sys_eqs::stype > & cols , typename Sys_eqs::stype coeff)
{
long int old_val = kmap.getKeyRef().get(d);
......
......@@ -124,7 +124,7 @@ class FDScheme
public:
// Distributed grid map
typedef grid_dist_id<Sys_eqs::dims,typename Sys_eqs::stype,scalar<size_t>,typename Sys_eqs::b_grid::decomposition> g_map_type;
typedef grid_dist_id<Sys_eqs::dims,typename Sys_eqs::stype,scalar<size_t>,typename Sys_eqs::b_grid::decomposition::extended_type> g_map_type;
typedef Sys_eqs Sys_eqs_typ;
......@@ -135,7 +135,8 @@ private:
typedef typename Sys_eqs::SparseMatrix_type::triplet_type triplet;
openfpm::vector<typename Sys_eqs::stype> b;
// Vector b
typename Sys_eqs::Vector_type b;
// Domain Grid informations
const grid_sm<Sys_eqs::dims,void> & gs;
......@@ -325,13 +326,14 @@ public:
* \param pd Padding, how many points out of boundary are present
* \param domain extension of the domain
* \param gs grid infos where Finite differences work
* \param stencil maximum extension of the stencil on each directions
* \param dec Decomposition of the domain
*
*/
FDScheme(Padding<Sys_eqs::dims> & pd, const Box<Sys_eqs::dims,typename Sys_eqs::stype> & domain, const grid_sm<Sys_eqs::dims,void> & gs, const typename Sys_eqs::b_grid::decomposition & dec)
:pd(pd),gs(gs),g_map(dec,padding(gs.getSize(),pd),domain,Ghost<Sys_eqs::dims,typename Sys_eqs::stype>(1)),row(0),row_b(0)
FDScheme(Padding<Sys_eqs::dims> & pd, const Ghost<Sys_eqs::dims,long int> & stencil, const Box<Sys_eqs::dims,typename Sys_eqs::stype> & domain, const grid_sm<Sys_eqs::dims,void> & gs, const typename Sys_eqs::b_grid & b_g)
:pd(pd),gs(gs),g_map(b_g,stencil,pd),row(0),row_b(0)
{
Vcluster & v_cl = dec.getVC();
Vcluster & v_cl = b_g.getDecomposition().getVC();
// Calculate the size of the local domain
size_t sz = g_map.getLocalDomainSize();
......@@ -364,10 +366,7 @@ public:
}
// sync the ghost
g_map.write("g_map_before_ghost.vtk");
g_map.template ghost_get<0>();
g_map.write("g_map_after_ghost.vtk");
// Create a CellDecomposer and calculate the spacing
......@@ -399,8 +398,8 @@ public:
template<typename T> void impose(const T & op , typename Sys_eqs::stype num ,long int id ,const long int (& start)[Sys_eqs::dims], const long int (& stop)[Sys_eqs::dims], bool skip_first = false)
{
// add padding to start and stop
grid_key_dx<Sys_eqs::dims> start_k = grid_key_dx<Sys_eqs::dims>(start) + pd.getKP1();
grid_key_dx<Sys_eqs::dims> stop_k = grid_key_dx<Sys_eqs::dims>(stop) + pd.getKP1();
grid_key_dx<Sys_eqs::dims> start_k = grid_key_dx<Sys_eqs::dims>(start);
grid_key_dx<Sys_eqs::dims> stop_k = grid_key_dx<Sys_eqs::dims>(stop);
auto it = g_map.getSubDomainIterator(start_k,stop_k);
......@@ -423,6 +422,8 @@ public:
*/
template<typename T> void impose(const T & op , typename Sys_eqs::stype num ,long int id ,grid_dist_iterator_sub<Sys_eqs::dims,typename g_map_type::d_grid> it_d, bool skip_first = false)
{
Vcluster & v_cl = *global_v_cluster;
openfpm::vector<triplet> & trpl = A.getMatrixTriplets();
auto it = it_d;
......@@ -431,19 +432,21 @@ public:
std::unordered_map<long int,float> cols;
// resize b if needed
b.resize(Sys_eqs::nvar * gs.size());
b.resize(Sys_eqs::nvar * g_map.size());
bool is_first = skip_first;
// iterate all the grid points
while (it.isNext())
{
if (is_first == true)
if (is_first == true && v_cl.getProcessUnitID() == 0)
{
++it;
is_first = false;
continue;
}
else
is_first = false;
// get the position
auto key = it.get();
......@@ -463,7 +466,7 @@ public:
// std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
}
b.get(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num;
b(g_map.template get<0>(key)*Sys_eqs::nvar + id) = num;
cols.clear();
// std::cout << "\n";
......@@ -491,14 +494,12 @@ public:
#ifdef SE_CLASS1
consistency();
#endif
A.resize(row,row);
A.resize(g_map.size()*Sys_eqs::nvar,g_map.size()*Sys_eqs::nvar);
return A;
}
typename Sys_eqs::Vector_type B;
/*! \brief produce the B vector
*
* \return the vector produced
......@@ -510,13 +511,7 @@ public:
consistency();
#endif
B.resize(row_b);
// copy the vector
for (size_t i = 0; i < row_b; i++)
B.insert(i,b.get(i));
return B;
return b;
}
};
......
......@@ -14,6 +14,7 @@
#include "util/grid_dist_testing.hpp"
#include "FiniteDifference/Average.hpp"
#include "FiniteDifference/sum.hpp"
#include "eq.hpp"
constexpr unsigned int x = 0;
constexpr unsigned int y = 1;
......
......@@ -33,7 +33,7 @@ class Lap
* \snippet FDScheme_unit_tests.hpp Laplacian usage
*
*/
inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
}
......@@ -88,7 +88,7 @@ public:
* \snippet FDScheme_unit_tests.hpp Laplacian usage
*
*/
inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap , const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
// for each dimension
for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
......
......@@ -81,7 +81,7 @@ struct pos_val
template<unsigned int f, typename Sys_eqs>
class Field
{
typedef typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type map_grid;
typedef typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type::extended_type>::type map_grid;
public:
......
......@@ -18,6 +18,7 @@
#include "Vector/Vector.hpp"
#include "Solvers/umfpack_solver.hpp"
#include "data_type/aggregate.hpp"
#include "FiniteDifference/FDScheme.hpp"
BOOST_AUTO_TEST_SUITE( eq_test_suite )
......@@ -139,16 +140,18 @@ typedef Avg<y,v_x,lid_nn,FORWARD> avg_vx_f;
BOOST_AUTO_TEST_CASE(lid_driven_cavity)
{
Vcluster & v_cl = *global_v_cluster;
//! [lid-driven cavity 2D]
// Domain, a rectangle
Box<2,float> domain({0.0,0.0},{1.0,1.0});
Box<2,float> domain({0.0,0.0},{3.0,1.0});
// Ghost (Not important in this case but required)
Ghost<2,float> g(0.01);
// Grid points on x=256 and y=64
long int sz[] = {8,8};
long int sz[] = {256,64};
size_t szu[2];
szu[0] = (size_t)sz[0];
szu[1] = (size_t)sz[1];
......@@ -165,8 +168,11 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
// Distributed grid that store the solution
grid_dist_id<2,float,aggregate<float[2],float>,CartDecomposition<2,float>> g_dist(szu,domain,g);
// Ghost stencil
Ghost<2,long int> stencil_max(1);
// Finite difference scheme
FDScheme<lid_nn> fd(pd,domain,g_dist.getGridInfo(),g_dist.getDecomposition());
FDScheme<lid_nn> fd(pd, stencil_max, domain, g_dist.getGridInfo(), g_dist);
// Here we impose the equation, we start from the incompressibility Eq imposed in the bulk with the
// exception of the first point {0,0} and than we set P = 0 in {0,0}, why we are doing this is again
......@@ -218,12 +224,35 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
//! [lid-driven cavity 2D]
g_dist.write("lid_driven_cavity");
// Check that match
bool test = compare("lid_driven_cavity_grid_0_test.vtk","lid_driven_cavity_grid_0.vtk");
BOOST_REQUIRE_EQUAL(test,true);
g_dist.write("lid_driven_cavity_p" + std::to_string(v_cl.getProcessingUnits()));
if (v_cl.getProcessUnitID() == 0)
{
if (v_cl.getProcessingUnits() == 1)
{
// Check that match
bool test = compare("lid_driven_cavity_p1_grid_0_test.vtk","lid_driven_cavity_grid_0.vtk");
BOOST_REQUIRE_EQUAL(test,true);
}
else if (v_cl.getProcessingUnits() == 2)
{
// Check that match
bool test = compare("lid_driven_cavity_p2_grid_0_test.vtk","lid_driven_cavity_p2_grid_0.vtk");
BOOST_REQUIRE_EQUAL(test,true);
test = compare("lid_driven_cavity_p2_grid_1_test.vtk","lid_driven_cavity_p2_grid_1.vtk");
BOOST_REQUIRE_EQUAL(test,true);
}
else if (v_cl.getProcessingUnits() == 3)
{
// Check that match
bool test = compare("lid_driven_cavity_p3_grid_0_test.vtk","lid_driven_cavity_p3_grid_0.vtk");
BOOST_REQUIRE_EQUAL(test,true);
test = compare("lid_driven_cavity_p3_grid_1_test.vtk","lid_driven_cavity_p3_grid_1.vtk");
BOOST_REQUIRE_EQUAL(test,true);
test = compare("lid_driven_cavity_p3_grid_2_test.vtk","lid_driven_cavity_p3_grid_2.vtk");
BOOST_REQUIRE_EQUAL(test,true);
}
}
}
BOOST_AUTO_TEST_SUITE_END()
......
......@@ -130,6 +130,8 @@ typedef Avg<x,v_z,lid_nn_3d,FORWARD> avg_x_vz_f;
BOOST_AUTO_TEST_CASE(lid_driven_cavity)
{
Vcluster & v_cl = *global_v_cluster;
// Domain
Box<3,float> domain({0.0,0.0,0.0},{3.0,1.0,1.0});
......@@ -154,8 +156,11 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
// Initialize openfpm
grid_dist_id<3,float,aggregate<float[3],float>,CartDecomposition<3,float>> g_dist(szu,domain,g);
// Ghost stencil
Ghost<3,long int> stencil_max(1);
// Distributed grid
FDScheme<lid_nn_3d> fd(pd,domain,g_dist.getGridInfo(),g_dist.getDecomposition());
FDScheme<lid_nn_3d> fd(pd,stencil_max,domain,g_dist.getGridInfo(),g_dist);
// start and end of the bulk
......@@ -228,11 +233,35 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
// Bring the solution to grid
x.copy<FDScheme<lid_nn_3d>,decltype(g_dist),velocity,pressure>(fd,{0,0},{sz[0]-1,sz[1]-1,sz[2]-1},g_dist);
g_dist.write("lid_driven_cavity_3d");
// Check that match
bool test = compare("lid_driven_cavity_3d_grid_0.vtk","lid_driven_cavity_3d_grid_0_test.vtk");
BOOST_REQUIRE_EQUAL(test,true);
g_dist.write("lid_driven_cavity_3d_p" + std::to_string(v_cl.getProcessingUnits()));
if (v_cl.getProcessUnitID() == 0)
{
if (v_cl.getProcessingUnits() == 1)
{
// Check that match
bool test = compare("lid_driven_cavity_3d_p1_grid_0_test.vtk","lid_driven_cavity_3d_p1_grid_0.vtk");
BOOST_REQUIRE_EQUAL(test,true);
}
else if (v_cl.getProcessingUnits() == 2)
{
// Check that match
bool test = compare("lid_driven_cavity_p2_grid_0_test.vtk","lid_driven_cavity_p2_grid_0.vtk");
BOOST_REQUIRE_EQUAL(test,true);
test = compare("lid_driven_cavity_p2_grid_1_test.vtk","lid_driven_cavity_p2_grid_1.vtk");
BOOST_REQUIRE_EQUAL(test,true);
}
else if (v_cl.getProcessingUnits() == 3)
{
// Check that match
bool test = compare("lid_driven_cavity_p3_grid_0_test.vtk","lid_driven_cavity_p3_grid_0.vtk");
BOOST_REQUIRE_EQUAL(test,true);
test = compare("lid_driven_cavity_p3_grid_1_test.vtk","lid_driven_cavity_p3_grid_1.vtk");
BOOST_REQUIRE_EQUAL(test,true);
test = compare("lid_driven_cavity_p3_grid_2_test.vtk","lid_driven_cavity_p3_grid_2.vtk");
BOOST_REQUIRE_EQUAL(test,true);
}
}
}
BOOST_AUTO_TEST_SUITE_END()
......
......@@ -51,7 +51,7 @@ struct const_mul_functor_value
const grid_sm<last::dims,void> & gs;
// grid mapping
const grid_dist_id<last::dims,typename last::stype,scalar<size_t>,typename last::b_grid::decomposition> & g_map;
const grid_dist_id<last::dims,typename last::stype,scalar<size_t>,typename last::b_grid::decomposition::extended_type> & g_map;
// grid position
grid_dist_key_dx<last::dims> & kmap;
......@@ -65,7 +65,7 @@ struct const_mul_functor_value
/*! \brief constructor
*
*/
const_mul_functor_value(const grid_dist_id<last::dims,typename last::stype,scalar<size_t>,typename last::b_grid::decomposition> & g_map, grid_dist_key_dx<last::dims> & kmap, const grid_sm<last::dims,void> & gs, typename last::stype (& spacing)[last::dims], std::unordered_map<long int,typename last::stype> & cols, typename last::stype coeff)
const_mul_functor_value(const grid_dist_id<last::dims,typename last::stype,scalar<size_t>,typename last::b_grid::decomposition::extended_type> & g_map, grid_dist_key_dx<last::dims> & kmap, const grid_sm<last::dims,void> & gs, typename last::stype (& spacing)[last::dims], std::unordered_map<long int,typename last::stype> & cols, typename last::stype coeff)
:cols(cols),gs(gs),g_map(g_map),kmap(kmap),coeff(coeff),spacing(spacing)
{};
......@@ -113,7 +113,7 @@ struct mul
* \param coeff coefficent (constant in front of the derivative)
*
*/
inline static void value(const grid_dist_id<Sys_eqs::dims,typename Sys_eqs::stype,scalar<size_t>,typename Sys_eqs::b_grid::decomposition> & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap, const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
inline static void value(const grid_dist_id<Sys_eqs::dims,typename Sys_eqs::stype,scalar<size_t>,typename Sys_eqs::b_grid::decomposition::extended_type> & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap, const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
const_mul_functor_value<v_expr> mfv(g_map,kmap,gs,spacing,cols,coeff);
......
......@@ -29,7 +29,7 @@ struct sum_functor_value
const grid_sm<last::dims,void> & gs;
// grid mapping
const typename stub_or_real<last,last::dims,typename last::stype,typename last::b_grid::decomposition>::type & g_map;
const typename stub_or_real<last,last::dims,typename last::stype,typename last::b_grid::decomposition::extended_type>::type & g_map;
// grid position
grid_dist_key_dx<last::dims> & kmap;
......@@ -47,7 +47,7 @@ struct sum_functor_value
/*! \brief constructor
*
*/
sum_functor_value(const typename stub_or_real<last,last::dims,typename last::stype,typename last::b_grid::decomposition>::type & g_map, grid_dist_key_dx<last::dims> & kmap, const grid_sm<last::dims,void> & gs, typename last::stype (& spacing)[last::dims], std::unordered_map<long int,typename last::stype> & cols, typename last::stype coeff)
sum_functor_value(const typename stub_or_real<last,last::dims,typename last::stype,typename last::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<last::dims> & kmap, const grid_sm<last::dims,void> & gs, typename last::stype (& spacing)[last::dims], std::unordered_map<long int,typename last::stype> & cols, typename last::stype coeff)
:cols(cols),gs(gs),g_map(g_map),kmap(kmap),key(key),coeff(coeff),spacing(spacing)
{};
......@@ -93,7 +93,7 @@ struct sum
* \snippet FDScheme_unit_tests.hpp sum example
*
*/
inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap, const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap, const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
// Sum functor
sum_functor_value<v_expr> sm(g_map,kmap,gs,spacing,cols,coeff);
......@@ -129,7 +129,7 @@ struct minus
* \snipper FDScheme_unit_tests.hpp minus example
*
*/
inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap, const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
inline static void value(const typename stub_or_real<Sys_eqs,Sys_eqs::dims,typename Sys_eqs::stype,typename Sys_eqs::b_grid::decomposition::extended_type>::type & g_map, grid_dist_key_dx<Sys_eqs::dims> & kmap, const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::stype (& spacing )[Sys_eqs::dims] , std::unordered_map<long int,typename Sys_eqs::stype > & cols, typename Sys_eqs::stype coeff)
{
arg::value(g_map,kmap,gs,spacing,cols,-coeff);
}
......
......@@ -214,6 +214,103 @@ public:
{
return mat.coeffRef(i,j);
}
/*! \brief Save this object into file
*
* \param file filename
*
* \return true if succed
*
*/
bool save(const std::string & file) const
{
std::vector<size_t> pap_prp;
Packer<decltype(trpl),HeapMemory>::packRequest(trpl,pap_prp);
// Calculate how much preallocated memory we need to pack all the objects
size_t req = ExtPreAlloc<HeapMemory>::calculateMem(pap_prp);
// allocate the memory
HeapMemory pmem;
pmem.allocate(req);
ExtPreAlloc<HeapMemory> mem(pap_prp,pmem);
//Packing
Pack_stat sts;
Packer<decltype(trpl),HeapMemory>::pack(mem,trpl,sts);
// Save into a binary file
std::ofstream dump (file, std::ios::out | std::ios::binary);
if (dump.is_open() == false)
return false;
dump.write ((const char *)pmem.getPointer(), pmem.size());
return true;
}
/*! \brief Load this object from file
*
* \param file filename
*
* \return true if succed
*
*/
bool load(const std::string & file)
{
std::ifstream fs (file, std::ios::in | std::ios::binary | std::ios::ate );
if (fs.is_open() == false)
return false;
// take the size of the file
size_t sz = fs.tellg();
fs.close();
// reopen the file without ios::ate to read
std::ifstream input (file, std::ios::in | std::ios::binary );
if (input.is_open() == false)
return false;
// Create the HeapMemory and the ExtPreAlloc memory
std::vector<size_t> pap_prp;
pap_prp.push_back(sz);
HeapMemory pmem;
ExtPreAlloc<HeapMemory> mem(pap_prp,pmem);
// read
input.read((char *)pmem.getPointer(), sz);
//close the file
input.close();
//Unpacking
Unpack_stat ps;
Unpacker<decltype(trpl),HeapMemory>::unpack(mem,trpl,ps);
return true;
}
/*! \brief Get the value from triplet
*
* \warning It is extremly slow because it do a full search across the triplets elements
*
* \param r row
* \param c colum
*
*/
T getValue(size_t r, size_t c)
{
for (size_t i = 0 ; i < trpl.size() ; i++)
{
if (r == (size_t)trpl.get(i).row() && c == (size_t)trpl.get(i).col())
return trpl.get(i).value();
}
return 0;
}
};
......
......@@ -12,7 +12,7 @@
#include "util/mul_array_extents.hpp"
#include <fstream>
#include "Vector_eigen_util.hpp"
#include "Grid/staggered_dist_grid_util.hpp"
#include "Grid/staggered_dist_grid.hpp"
#include "Space/Ghost.hpp"
#include "FiniteDifference/util/common.hpp"
#include <boost/mpl/vector_c.hpp>
......@@ -88,146 +88,45 @@ class Vector<T,Eigen::Matrix<T, Eigen::Dynamic, 1>>
//size of each chunk
mutable openfpm::vector<size_t> sz;
/*! \brief Check that the size of the iterators match
*
* It check the the boxes that the sub iterator defines has same dimensions, for example
* if the first sub-iterator, iterate from (1,1) to (5,3) and the second from (2,2) to (6,4)
* they match (2,2) to (4,6) they do not match
*
* \tparam Grid_map type of the map grid
* \tparam Grid_dst type of the destination grid
*
* \param it1 Iterator1
* \param it2 Iterator2
*
* \return true if they match
*
*/
template<typename Eqs_sys, typename it1_type, typename it2_type> bool checkIterator(const it1_type & it1, const it2_type & it2)
{
#ifdef SE_CLASS1
grid_key_dx<Eqs_sys::dims> it1_k = it1.getStop() - it1.getStart();
grid_key_dx<Eqs_sys::dims> it2_k = it2.getStop() - it2.getStart();
for (size_t i = 0 ; i < Eqs_sys::dims ; i++)
{
if (it1_k.get(i) != it2_k.get(i))
{
std::cerr << __FILE__ << ":" << __LINE__ << " error src iterator and destination iterator does not match in size\n";
return false;
}
}
return true;
#else
return true;
#endif
}