Commit f2dc4cea authored by Pietro Incardona's avatar Pietro Incardona

Fixing numerics for Finite Differences

parent 3385f02c
/*
* eq_1_ag_unit_test.hpp
*
* This file model the Eq1 in 2D of Active Gel
*
* R. Ramaswamy, G. Bourantas, F. Jülicher, and I. F. Sbalzarini.
* A hybrid particle-mesh method for incompressible active polar
* viscous gels. J. Comput. Phys., 291:334–361, 2015.
*
* Created on: Sep 18, 2015
* Author: Pietro Incardona
*/
#ifndef OPENFPM_NUMERICS_SRC_EQUATIONS_EQ_1_AG_UNIT_TEST_HPP_
#define OPENFPM_NUMERICS_SRC_EQUATIONS_EQ_1_AG_UNIT_TEST_HPP_
////////////////////////////////////// Left part eq 1 //////////////////////////////
// Creating square parentesis term left part equation 1
// Derivation and multiplication of t1 ... t6
namespace eq1
{
typedef mul<etad2,D<x,t1>> t1_l;
typedef mul<eta,D<x,t2>> t2_l;
typedef mul<etad2,D<x,t3>> t3_l;
typedef mul<eta,D<y,t4>> t4_l;
typedef mul<eta2,D<y,t5>> t5_l;
typedef mul<eta,D<y,t6>> t6_l;
typedef mul<eta,Lap<v>> meLv;
typedef D<x,P> DxP;
// Assemble the left part
typedef add<meLV,minus<DxP>,t1_l,t2_l,t3_l,t4_l,t5_l,t6_l> left;
////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////// Right part eq 1 /////////////////////////////
typedef div<2,D<x,hp>> t1_r;
typedef mul<zeta,D<x,mul<del_nu,px2>>> t2_r;
typedef mul<zeta,D<y,mul<del_nu,px_py>>> t3_r;
typedef mul<zeta,D<x,mul<del_nu,dev<2,px2_p_py2>>>> t4_r;
typedef mul<nud2,D<x,mul<-2,hp,px_py>> t5_r;
typedef mul<etad2,D<y,mul<hp,px2_m_py2>>> t6_r;
typedef D<x,sigma[x][x]> t7_r;
typedef D<y,sigma[x][y]> t8_r;
typedef add<t1_r,t2_r,t3_r,t4_r,t5_r,t6_r,t7_r,t8_r,g_ext> eq1r;
typedef Eq<eq1l,eq1r> equation;
}
/////////////////////////////////////////////////////////////////////////////////////
#endif /* OPENFPM_NUMERICS_SRC_EQUATIONS_EQ_1_AG_UNIT_TEST_HPP_ */
/*
* eq_2_ag_unit_test.hpp
*
* Created on: Sep 18, 2015
* Author: i-bird
*/
#ifndef OPENFPM_NUMERICS_SRC_EQUATIONS_EQ_2_AG_UNIT_TEST_HPP_
#define OPENFPM_NUMERICS_SRC_EQUATIONS_EQ_2_AG_UNIT_TEST_HPP_
// Left part equation 2
namespace eq2
{
typedef mul<etad2,D<y,t1>> t1_l;
typedef mul<eta,D<y,t2>> t2_l;
typedef mul<etad2,D<y,t3>> t3_l;
typedef mul<eta,D<x,t4>> t4_l;
typedef mul<eta2,D<x,t5>> t5_l;
typedef mul<eta,D<x,t6>> t6_l;
typedef mul<eta,Lap<v>> meLv;
typedef D<y,P> DyP;
// left part equation 2
typedef add<meLv,DyP,Dt1,Dt2,Dt3,Dt4,Dt6> left;
typedef div<-2,D<x,hp>> t1_r;
typedef mul<zeta,D<y,mul<del_nu,py2>>> t2_r;
typedef mul<zeta,D<x,mul<del_nu,px_py>>> t3_r;
typedef mul<zeta,D<y,mul<del_nu,dev<2,px2_p_py2>>>> t4_r;
typedef mul<zeta,D<y,mul<del_nu,dev<2,px2_p_py2>>>> t5_r;
typedef mul<nud2,D<x,mul<hp,px2_m_py2>>> t6_r;
typedef D<x,sigma[x][y]> t7_r;
typedef D<y,sigma[y][x]> t8_r;
typedef mul<nud2,D<y,mul<gm,lambda,del_nu,px2_m_py2>> t9_r
typedef add<t1_r,t2_r,t3_r,t4_r,t5_r,t6_r,t7_r,t8_r,t9_r> right;
}
#endif /* OPENFPM_NUMERICS_SRC_EQUATIONS_EQ_2_AG_UNIT_TEST_HPP_ */
/*
* eq_base_unit_test.hpp
*
* Created on: Sep 18, 2015
* Author: i-bird
*/
#ifndef OPENFPM_NUMERICS_SRC_EQUATIONS_EQ_BASE_UNIT_TEST_HPP_
#define OPENFPM_NUMERICS_SRC_EQUATIONS_EQ_BASE_UNIT_TEST_HPP_
// quadratic and cubic term in p
typedef mul<p[x],p[x]> px2;
typedef mul<p[y],p[y]> py2;
typedef mul<p[x],p[y]> px_py;
typedef add<px2,py2> px2_p_py2;
typedef add<px2,py2> px2_s_py2;
typedef mul<px2,px> px3;
typedef mul<py2,py> py3;
typedef mul<gm,nu> gm_nu;
// quadratic term in p
typedef mul<px3,py> px3_m_py;
typedef mul<px2,py2> px2_m_py2;
typedef mul<px,py3> px_m_py3;
typedef mul<px2,px2_s_py2> px2_r;
typedef mul<px_py,px2_s_py2> px_py_r;
typedef mul<py2,px2_s_py2> py2_r;
// Division of the quadratic term by px2_s_py2
typedef div<px3_m_py,px2_p_py2> p_c1;
typedef div<px2_m_py2,px2_p_py2> p_c2;
typedef div<px_m_py3,px2_p_py2> p_c3;
typedef div<px2_r,px2_p_py2> p_c4;
typedef div<px_py_r,px2_p_py2> p_c5;
typedef div<py2_r,px2_p_py2> p_c6;
// Create u_alpha_beta equation
template <unsigned int i, unsigned int j> using u = div<2,add<D<i,v[j]>,D<j,v[i]>>>;
// terms
typedef mul<gm,p_c1,u[x][x]> t1;
typedef mul<gm,p_c2,u[x][y]> t2;
typedef mul<gm,p_c3,u[y][y]> t3;
typedef mul<gm,p_c4,u[x][x]> t4;
typedef mul<gm,p_c5,u[x][y]> t5;
typedef mul<gm,p_c6,u[y][y]> t6;
#endif /* OPENFPM_NUMERICS_SRC_EQUATIONS_EQ_BASE_UNIT_TEST_HPP_ */
/*
* 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_ */
/*
* Average.hpp
*
* Created on: Nov 18, 2015
* Author: Pietro Incardona
*/
#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_AVERAGE_HPP_
#define OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_AVERAGE_HPP_
#define CENTRAL 0
#define CENTRAL_B_ONE_SIDE 1
#define FORWARD 2
#include "util/mathutil.hpp"
#include "Vector/map_vector.hpp"
#include "Grid/comb.hpp"
#include "FiniteDifference/util/common.hpp"
#include "util/util_num.hpp"
/*! \brief Average
*
* \tparam d on which dimension average
* \tparam Field which field average
* \tparam impl which implementation
*
*/
template<unsigned int d, typename Field, typename Sys_eqs, unsigned int impl=CENTRAL>
class Avg
{
/*! \brief Create the row of the Matrix
*
* \tparam ord
*
*/
inline static void value(const grid_key_dx<Sys_eqs::dims> & pos, 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";
}
/*! \brief Calculate the position where the derivative is calculated
*
* In case on non staggered case this function just return pos, in case of staggered,
* it calculate where the operator is calculated on a staggered grid
*
*/
inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs)
{
std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
}
};
/*! \brief Average on direction i
*
*
*/
template<unsigned int d, typename arg, typename Sys_eqs>
class Avg<d,arg,Sys_eqs,CENTRAL>
{
public:
/*! \brief fill the row
*
*
*/
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)
{
// if the system is staggered the CENTRAL derivative is equivalent to a forward derivative
if (is_grid_staggered<Sys_eqs>::value() == true)
{
Avg<d,arg,Sys_eqs,FORWARD>::value(g_map,kmap,gs,cols,coeff);
return;
}
long int old_val = kmap.getKeyRef().get(d);
kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 1);
arg::value(g_map,kmap,gs,cols,coeff);
kmap.getKeyRef().set_d(d,old_val);
old_val = kmap.getKeyRef().get(d);
kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 1);
arg::value(g_map,kmap,gs,cols,coeff);
kmap.getKeyRef().set_d(d,old_val);
}
/*! \brief Calculate the position where the derivative is calculated
*
* In case on non staggered case this function just return pos, in case of staggered,
* it calculate where the operator is calculated on a staggered grid
*
* \param pos from the position
* \param fld Field we are deriving, if not provided the function just return pos
* \param s_pos position of the properties in the staggered grid
*
*/
inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, long int fld = -1, const openfpm::vector<comb<Sys_eqs::dims>> & s_pos = openfpm::vector<comb<Sys_eqs::dims>>())
{
if (is_grid_staggered<Sys_eqs>::value())
{
if (fld == -1)
return pos;
if (s_pos[fld][d] == 1)
{
grid_key_dx<Sys_eqs::dims> ret = pos;
ret.set_d(d,0);
return pos;
}
else
{
grid_key_dx<Sys_eqs::dims> ret = pos;
ret.set_d(d,1);
return pos;
}
}
return pos;
}
};
/*! \brief Average FORWARD on direction i
*
*g
*/
template<unsigned int d, typename arg, typename Sys_eqs>
class Avg<d,arg,Sys_eqs,FORWARD>
{
public:
/*! \brief fill the row
*
*
*/
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)
{
long int old_val = kmap.getKeyRef().get(d);
kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 1);
arg::value(g_map,kmap,gs,cols,coeff);
kmap.getKeyRef().set_d(d,old_val);
// backward
arg::value(g_map,kmap,gs,cols,coeff);
}
/*! \brief Calculate the position where the derivative is calculated
*
* In case on non staggered case this function just return pos, in case of staggered,
* it calculate where the operator is calculated on a staggered grid
*
* \param pos from the position
* \param fld Field we are deriving, if not provided the function just return pos
* \param s_pos position of the properties in the staggered grid
*
*/
inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, long int fld = -1, const openfpm::vector<comb<Sys_eqs::dims>> & s_pos = openfpm::vector<comb<Sys_eqs::dims>>())
{
return pos;
}
};
#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_AVERAGE_HPP_ */
......@@ -2,7 +2,7 @@
* Derivative.hpp
*
* Created on: Oct 5, 2015
* Author: i-bird
* Author: Pietro Incardona
*/
#ifndef OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_DERIVATIVE_HPP_
......
......@@ -84,8 +84,10 @@ class FDScheme
// Domain Grid informations
const grid_sm<Sys_eqs::dims,void> & gs;
typedef grid_dist_id<Sys_eqs::dims,typename Sys_eqs::stype,scalar<size_t>,typename Sys_eqs::b_grid::decomposition> g_map_type;
// mapping grid
grid_dist_id<Sys_eqs::dims,typename Sys_eqs::stype,scalar<size_t>,typename Sys_eqs::b_grid::decomposition> g_map;
g_map_type g_map;
// row of the matrix
size_t row;
......@@ -102,6 +104,23 @@ class FDScheme
// processor
size_t s_pnt;
/* \brief calculate the mapping grid size with padding
*
* \param gs original grid size
*
* \return padded grid size
*
*/
inline const std::vector<size_t> padding( const size_t (& sz)[Sys_eqs::dims], Padding<Sys_eqs::dims> & pd)
{
std::vector<size_t> g_sz_pad(Sys_eqs::dims);
for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
g_sz_pad[i] = sz[i] + pd.getLow(i) + pd.getHigh(i);
return g_sz_pad;
}
/*! \brief Check if the Matrix is consistent
*
*/
......@@ -142,14 +161,14 @@ class FDScheme
}
}
/*! \brief Convert from integer ghost to continuos
/* \brief Convert discrete ghost into continous ghost
*
* \return the continuos version of the ghost
* \param Ghost
*
*/
Ghost<dim,Sys_eqs::stype> convert_into_cg()
inline Ghost<Sys_eqs::dims,typename Sys_eqs::stype> convert_dg_cg(const Ghost<Sys_eqs::dims,typename Sys_eqs::stype> & g)
{
return g_map_type::convert_ghost(g,g_map.getCellDecomposer());
}
public:
......@@ -160,18 +179,20 @@ public:
* \param gs grid infos where Finite differences work
*
*/
FDScheme(Padding<Sys_eqs::dims> & pd, const Box<Sys_eqs::dims,typename Sys_eqs::stype> & domain, const grid_sm<Sys_eqs::dims,void> & gs, typename Sys_eqs::b_grid::decomposition & dec, Vcluster & v_cl)
:pd(pd),gs(gs),g_map(dec.duplicate(Ghost<Sys_eqs::dims,long int>(1)),gs.getSize(),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, Vcluster & v_cl)
: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)
{
// Calculate the size of the local domain
size_t sz = g_map.getLocalDomainSize();
// Get the total size of the local grids on each processors
v_cl.allGather(sz,pnt);
v_cl.execute();
s_pnt = 0;
// calculate the starting point for this processor
for (size_t i = 0 ; i < v_cl.getProcessUnitID() ; i++)
s_pnt += pnt.get(0);
s_pnt += pnt.get(i);
// Calculate the starting point
......@@ -204,8 +225,18 @@ public:
*
*
*/
template<typename T> void imposeA(const T & op , grid_dist_iterator_sub<Sys_eqs::dims,typename Sys_eqs::b_grid::d_grid> it)
template<typename T> void imposeA(const T & op , grid_dist_iterator_sub<Sys_eqs::dims,typename Sys_eqs::b_grid::d_grid> it_d)
{
// key one
grid_key_dx<Sys_eqs::dims> gk_one;
gk_one.one();
// add padding to start and stop
grid_key_dx<Sys_eqs::dims> start = it_d.getStart() + pd.getKP1();
grid_key_dx<Sys_eqs::dims> stop = it_d.getStop() + pd.getKP1();
auto it = g_map.getSubDomainIterator(start,stop);
std::unordered_map<long int,float> cols;
// iterate all the grid points
......
......@@ -12,6 +12,7 @@
#include "FiniteDifference/Laplacian.hpp"
#include "Decomposition/CartDecomposition.hpp"
#include "util/grid_dist_testing.hpp"
#include "FiniteDifference/Average.hpp"
constexpr unsigned int x = 0;
constexpr unsigned int y = 1;
......@@ -113,7 +114,7 @@ const bool syss_pp::boundary[] = {PERIODIC,PERIODIC};
BOOST_AUTO_TEST_SUITE( fd_test )
BOOST_AUTO_TEST_CASE( fd_test_use_non_periodic)
BOOST_AUTO_TEST_CASE( fd_test_use_der_non_periodic)
{
// grid size
size_t sz[2]={16,16};
......@@ -239,8 +240,78 @@ BOOST_AUTO_TEST_CASE( fd_test_use_non_periodic)
BOOST_REQUIRE_EQUAL(cols_y[13*16+15],0.5);
}
BOOST_AUTO_TEST_CASE( fd_test_use_avg_non_periodic)
{
// grid size
size_t sz[2]={16,16};
// grid_dist_testing
grid_dist_testing<2> g_map(sz);
// grid_sm
grid_sm<2,void> ginfo(sz);
// Create several keys
grid_dist_key_dx<2> key11(0,grid_key_dx<2>(1,1));
grid_dist_key_dx<2> key00(0,grid_key_dx<2>(0,0));
grid_dist_key_dx<2> key22(0,grid_key_dx<2>(2,2));
grid_dist_key_dx<2> key1515(0,grid_key_dx<2>(15,15));
// filled colums
std::unordered_map<long int,float> cols_x;
std::unordered_map<long int,float> cols_y;
Avg<x,Field<V,sys_nn>,sys_nn>::value(g_map,key11,ginfo,cols_x,1);
Avg<y,Field<V,sys_nn>,sys_nn>::value(g_map,key11,ginfo,cols_y,1);
BOOST_REQUIRE_EQUAL(cols_x.size(),2);
BOOST_REQUIRE_EQUAL(cols_y.size(),2);
BOOST_REQUIRE_EQUAL(cols_x[17+1],1);
BOOST_REQUIRE_EQUAL(cols_x[17-1],1);
BOOST_REQUIRE_EQUAL(cols_y[17+16],1);
BOOST_REQUIRE_EQUAL(cols_y[17-16],1);
// filled colums
std::unordered_map<long int,float> cols_xx;
std::unordered_map<long int,float> cols_xy;
std::unordered_map<long int,float> cols_yx;
std::unordered_map<long int,float> cols_yy;
// Composed derivative
Avg<x,Avg<x,Field<V,sys_nn>,sys_nn>,sys_nn>::value(g_map,key22,ginfo,cols_xx,1);
Avg<x,Avg<y,Field<V,sys_nn>,sys_nn>,sys_nn>::value(g_map,key22,ginfo,cols_xy,1);
Avg<y,Avg<x,Field<V,sys_nn>,sys_nn>,sys_nn>::value(g_map,key22,ginfo,cols_yx,1);
Avg<y,Avg<y,Field<V,sys_nn>,sys_nn>,sys_nn>::value(g_map,key22,ginfo,cols_yy,1);
BOOST_REQUIRE_EQUAL(cols_xx.size(),3);
BOOST_REQUIRE_EQUAL(cols_xy.size(),4);
BOOST_REQUIRE_EQUAL(cols_yx.size(),4);
BOOST_REQUIRE_EQUAL(cols_yy.size(),3);
BOOST_REQUIRE_EQUAL(cols_xx[32],1);
BOOST_REQUIRE_EQUAL(cols_xx[34],2);
BOOST_REQUIRE_EQUAL(cols_xx[36],1);
BOOST_REQUIRE_EQUAL(cols_xy[17],1);
BOOST_REQUIRE_EQUAL(cols_xy[19],1);
BOOST_REQUIRE_EQUAL(cols_xy[49],1);
BOOST_REQUIRE_EQUAL(cols_xy[51],1);
BOOST_REQUIRE_EQUAL(cols_yx[17],1);
BOOST_REQUIRE_EQUAL(cols_yx[19],1);
BOOST_REQUIRE_EQUAL(cols_yx[49],1);
BOOST_REQUIRE_EQUAL(cols_xy[51],1);
BOOST_REQUIRE_EQUAL(cols_yy[2],1);
BOOST_REQUIRE_EQUAL(cols_yy[34],2);
BOOST_REQUIRE_EQUAL(cols_yy[66],1);
}
BOOST_AUTO_TEST_CASE( fd_test_use_staggered_non_periodic)
BOOST_AUTO_TEST_CASE( fd_test_use_staggered_der_non_periodic)
{