Skip to content

GitLab

O
openfpm_numerics
Commit 64d29477 authored by Pietro Incardona's avatar Pietro Incardona
Browse files
Options

Staggered to normal for grid

parent dfb8b8f7
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 586 additions and 385 deletions
src/FiniteDifference/Average.hpp
View file @ 64d29477
......@@ -65,24 +65,24 @@ class Avg<d,arg,Sys_eqs,CENTRAL>
*
*
*/
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>::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())
{
Avg<d,arg,Sys_eqs,BACKWARD>::value(g_map,kmap,gs,cols,coeff);
Avg<d,arg,Sys_eqs,BACKWARD>::value(g_map,kmap,gs,spacing,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);
arg::value(g_map,kmap,gs,spacing,cols,coeff/2);
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);
arg::value(g_map,kmap,gs,spacing,cols,coeff/2);
kmap.getKeyRef().set_d(d,old_val);
}
......@@ -135,16 +135,16 @@ class Avg<d,arg,Sys_eqs,FORWARD>
*
*
*/
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>::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);
arg::value(g_map,kmap,gs,cols,coeff);
arg::value(g_map,kmap,gs,spacing,cols,coeff/2);
kmap.getKeyRef().set_d(d,old_val);
// backward
arg::value(g_map,kmap,gs,cols,coeff);
arg::value(g_map,kmap,gs,spacing,cols,coeff/2);
}
......@@ -177,15 +177,15 @@ class Avg<d,arg,Sys_eqs,BACKWARD>
*
*
*/
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>::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);
arg::value(g_map,kmap,gs,cols,coeff);
arg::value(g_map,kmap,gs,spacing,cols,coeff/2);
kmap.getKeyRef().set_d(d,old_val);
// forward
arg::value(g_map,kmap,gs,cols,coeff);
arg::value(g_map,kmap,gs,spacing,cols,coeff/2);
}
......
src/FiniteDifference/Derivative.hpp
View file @ 64d29477
......@@ -64,24 +64,24 @@ class D<d,arg,Sys_eqs,CENTRAL>
*
*
*/
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>::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())
{
D<d,arg,Sys_eqs,BACKWARD>::value(g_map,kmap,gs,cols,coeff);
D<d,arg,Sys_eqs,BACKWARD>::value(g_map,kmap,gs,spacing,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);
arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[d]/2.0 );
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);
arg::value(g_map,kmap,gs,spacing,cols,-coeff/spacing[d]/2.0 );
kmap.getKeyRef().set_d(d,old_val);
}
......@@ -135,7 +135,7 @@ public:
*
*
*/
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)
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)
{
#ifdef SE_CLASS1
if (Sys_eqs::boundary[d] == PERIODIC)
......@@ -146,42 +146,42 @@ public:
if (pos.get(d) == (long int)gs.size(d)-1 )
{
arg::value(g_map,kmap,gs,cols,1.5*coeff);
arg::value(g_map,kmap,gs,spacing,cols,1.5*coeff/spacing[d]);
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,-2.0*coeff);
arg::value(g_map,kmap,gs,spacing,cols,-2.0*coeff/spacing[d]);
kmap.getKeyRef().set_d(d,old_val);
old_val = kmap.getKeyRef().get(d);
kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) - 2);
arg::value(g_map,kmap,gs,cols,0.5*coeff);
arg::value(g_map,kmap,gs,spacing,cols,0.5*coeff/spacing[d]);
kmap.getKeyRef().set_d(d,old_val);
}
else if (pos.get(d) == 0)
{
arg::value(g_map,kmap,gs,cols,-1.5*coeff);
arg::value(g_map,kmap,gs,spacing,cols,-1.5*coeff/spacing[d]);
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,2.0*coeff);
arg::value(g_map,kmap,gs,spacing,cols,2.0*coeff/spacing[d]);
kmap.getKeyRef().set_d(d,old_val);
old_val = kmap.getKeyRef().get(d);
kmap.getKeyRef().set_d(d, kmap.getKeyRef().get(d) + 2);
arg::value(g_map,kmap,gs,cols,-0.5*coeff);
arg::value(g_map,kmap,gs,spacing,cols,-0.5*coeff/spacing[d]);
kmap.getKeyRef().set_d(d,old_val);
}
else
{
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);
arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[d]);
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);
arg::value(g_map,kmap,gs,spacing,cols,-coeff/spacing[d]);
kmap.getKeyRef().set_d(d,old_val);
}
}
......@@ -231,16 +231,16 @@ class D<d,arg,Sys_eqs,FORWARD>
*
*
*/
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>::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);
arg::value(g_map,kmap,gs,cols,coeff);
arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[d]);
kmap.getKeyRef().set_d(d,old_val);
// backward
arg::value(g_map,kmap,gs,cols,-coeff);
arg::value(g_map,kmap,gs,spacing,cols,-coeff/spacing[d]);
}
......@@ -273,16 +273,16 @@ class D<d,arg,Sys_eqs,BACKWARD>
*
*
*/
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>::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);
arg::value(g_map,kmap,gs,cols,-coeff);
arg::value(g_map,kmap,gs,spacing,cols,-coeff/spacing[d]);
kmap.getKeyRef().set_d(d,old_val);
// forward
arg::value(g_map,kmap,gs,cols,coeff);
arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[d]);
}
......
src/FiniteDifference/FDScheme.hpp
View file @ 64d29477
......@@ -13,6 +13,7 @@
#include "Grid/grid_dist_id_iterator_sub.hpp"
#include "eq.hpp"
#include "data_type/scalar.hpp"
#include "NN/CellList/CellDecomposer.hpp"
/*! \brief Finite Differences
*
......@@ -73,6 +74,15 @@
template<typename Sys_eqs>
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 Sys_eqs Sys_eqs_typ;
private:
// Padding
Padding<Sys_eqs::dims> pd;
......@@ -86,7 +96,8 @@ 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;
// Get the grid spacing
typename Sys_eqs::stype spacing[Sys_eqs::dims];
// mapping grid
g_map_type g_map;
......@@ -154,14 +165,14 @@ class FDScheme
for (size_t i = 0 ; i < nz_rows.size() ; i++)
{
if (nz_rows.get(i) == false)
std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix not all the rows are filled\n";
std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix row " << i << " is not filled\n";
}
// all the colums must have a non zero element
for (size_t i = 0 ; i < nz_cols.size() ; i++)
{
if (nz_cols.get(i) == false)
std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix not all the colums are filled\n";
std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix colum " << i << " is not filled\n";
}
}
......@@ -177,6 +188,26 @@ class FDScheme
public:
/*! \brief Get the grid padding
*
* \return the grid padding
*
*/
const Padding<Sys_eqs::dims> & getPadding()
{
return pd;
}
/*! \brief Return the map between the grid index position and the position in the distributed vector
*
* \return the map
*
*/
const g_map_type & getMap()
{
return g_map;
}
/*! \brief Constructor
*
* \param pd Padding
......@@ -219,27 +250,69 @@ public:
// sync the ghost
g_map.template ghost_get<0>();
// Create a CellDecomposer and calculate the spacing
size_t sz_g[Sys_eqs::dims];
for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
sz_g[i] = gs.getSize()[i] - 1;
CellDecomposer_sm<Sys_eqs::dims,typename Sys_eqs::stype> cd(domain,sz_g,0);
for (size_t i = 0 ; i < Sys_eqs::dims ; i++)
spacing[i] = cd.getCellBox().getHigh(i);
}
/*! \brief Impose an operator
*
* This function impose an operator on a box region to produce the system
*
* Ax = b
*
* ## Stokes equation, lid driven cavity with one splipping wall
*
* \param op Operator to impose (A term)
* \param num right hand side of the term (b term)
* \param id Equation id in the system that we are imposing
* \param start starting point of the box
* \param stop stop point of the box
*
*/
template<typename T> void imposeA(const T & op , grid_dist_iterator_sub<Sys_eqs::dims,typename Sys_eqs::b_grid::d_grid> it_d, bool skip_first = false)
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)
{
// 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();
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();
auto it = g_map.getSubDomainIterator(start,stop);
auto it = g_map.getSubDomainIterator(start_k,stop_k);
impose(op,num,id,it,skip_first);
}
/*! \brief Impose an operator
*
* This function impose an operator on a particular grid region to produce the system
*
* Ax = b
*
* ## Stokes equation, lid driven cavity with one splipping wall
*
* \param op Operator to impose (A term)
* \param num right hand side of the term (b term)
* \param id Equation id in the system that we are imposing
* \param it_d iterator that define where you want to impose
*
*/
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)
{
auto it = it_d;
grid_sm<Sys_eqs::dims,void> gs = g_map.getGridInfoVoid();
std::unordered_map<long int,float> cols;
// resize b if needed
b.resize(Sys_eqs::nvar * gs.size());
bool is_first = skip_first;
// iterate all the grid points
......@@ -254,68 +327,29 @@ public:
// get the position
auto key = it.get();
grid_key_dx<2> gkey = g_map.getGKey(key);
// Calculate the non-zero colums
T::value(g_map,key,gs,cols,1.0);
T::value(g_map,key,gs,spacing,cols,1.0);
// create the triplet
for ( auto it = cols.begin(); it != cols.end(); ++it )
{
trpl.add();
trpl.last().row() = row;
trpl.last().row() = Sys_eqs::nvar * gs.LinId(gkey) + id;
trpl.last().col() = it->first;
trpl.last().value() = it->second;
std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
}
b.get(Sys_eqs::nvar * gs.LinId(gkey) + id) = num;
cols.clear();
std::cout << "\n";
++row;
++it;
}
}
/*! \brief Impose an operator
*
*
*
*/
void imposeB(typename Sys_eqs::stype num , grid_dist_iterator_sub<Sys_eqs::dims,typename Sys_eqs::b_grid::d_grid> it_d, bool skip_first = false)
{
// 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;
bool is_first = skip_first;
// iterate all the grid points
while (it.isNext())
{
if (is_first == true)
{
++it;
is_first = false;
continue;
}
// get the position
auto key = it.get();
b.add(num);
cols.clear();
++row_b;
++it;
}
......
src/FiniteDifference/Laplacian.hpp
View file @ 64d29477
......@@ -53,23 +53,23 @@ public:
*
*
*/
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>::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++)
{
long int old_val = kmap.getKeyRef().get(i);
kmap.getKeyRef().set_d(i, kmap.getKeyRef().get(i) + 1);
arg::value(g_map,kmap,gs,cols,coeff);
arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[i]/spacing[i]);
kmap.getKeyRef().set_d(i,old_val);
old_val = kmap.getKeyRef().get(i);
kmap.getKeyRef().set_d(i, kmap.getKeyRef().get(i) - 1);
arg::value(g_map,kmap,gs,cols,coeff);
arg::value(g_map,kmap,gs,spacing,cols,coeff/spacing[i]/spacing[i]);
kmap.getKeyRef().set_d(i,old_val);
}
arg::value(g_map,kmap,gs,cols, - 2.0 * Sys_eqs::dims * coeff);
arg::value(g_map,kmap,gs,spacing,cols, - 2.0 * coeff/spacing[i]/spacing[i]);
}
}
......
src/FiniteDifference/eq.hpp
View file @ 64d29477
......@@ -92,7 +92,7 @@ public:
*
*
*/
static void value(const map_grid & 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)
static void value(const map_grid & 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 (Sys_eqs::ord == EQS_FIELD)
cols[g_map.template get<0>(kmap)*Sys_eqs::nvar + f] += coeff;
......
src/FiniteDifference/eq_unit_test.hpp
View file @ 64d29477
......@@ -75,18 +75,20 @@ typedef Field<P,lid_nn> Prs;
typedef mul<eta,Lap<v_x,lid_nn>,lid_nn> eta_lap_vx;
typedef D<x,Prs,lid_nn> p_x;
typedef sum<eta_lap_vx,p_x,lid_nn> vx_eq;
typedef minus<p_x,lid_nn> m_p_x;
typedef sum<eta_lap_vx,m_p_x,lid_nn> vx_eq;
// Eq2 V_y
typedef mul<eta,Lap<v_y,lid_nn>,lid_nn> eta_lap_vy;
typedef D<x,Prs,lid_nn> p_x;
typedef sum<eta_lap_vy,p_x,lid_nn> vy_eq;
typedef D<y,Prs,lid_nn> p_y;
typedef minus<p_y,lid_nn> m_p_y;
typedef sum<eta_lap_vy,m_p_y,lid_nn> vy_eq;
// Eq3 Incompressibility
typedef D<x,v_x,lid_nn> dx_vx;
typedef D<y,v_y,lid_nn> dy_vy;
typedef D<x,v_x,lid_nn,FORWARD> dx_vx;
typedef D<y,v_y,lid_nn,FORWARD> dy_vy;
typedef sum<dx_vx,dy_vy,lid_nn> ic_eq;
......@@ -94,6 +96,13 @@ typedef sum<dx_vx,dy_vy,lid_nn> ic_eq;
typedef Avg<x,v_y,lid_nn> avg_vy;
typedef Avg<y,v_x,lid_nn> avg_vx;
typedef Avg<x,v_y,lid_nn,FORWARD> avg_vy_f;
typedef Avg<y,v_x,lid_nn,FORWARD> avg_vx_f;
#define EQ_1 0
#define EQ_2 1
#define EQ_3 2
BOOST_AUTO_TEST_SUITE( eq_test_suite )
// Lid driven cavity, uncompressible fluid
......@@ -106,7 +115,10 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
// Ghost
Ghost<2,float> g(0.01);
size_t sz[] = {8,8};
long int sz[] = {8,8};
size_t szu[2];
szu[0] = (size_t)sz[0];
szu[1] = (size_t)sz[1];
Padding<2> pd({1,1},{0,0});
......@@ -114,70 +126,52 @@ BOOST_AUTO_TEST_CASE(lid_driven_cavity)
init_global_v_cluster(&boost::unit_test::framework::master_test_suite().argc,&boost::unit_test::framework::master_test_suite().argv);
// Initialize openfpm
grid_dist_id<2,float,aggregate<float[2],float>,CartDecomposition<2,float>> g_dist(sz,domain,g);
grid_dist_id<2,float,aggregate<float[2],float>,CartDecomposition<2,float>> g_dist(szu,domain,g);
// Distributed grid
FDScheme<lid_nn> fd(pd,domain,g_dist.getGridInfo(),g_dist.getDecomposition(),g_dist.getVC());
// start and end of the bulk
fd.imposeA(ic_eq(), g_dist.getSubDomainIterator(grid_key_dx<2>(0,0),grid_key_dx<2>(sz[1]-2,sz[0]-2)),true);
fd.imposeB(0.0,g_dist.getSubDomainIterator(grid_key_dx<2>(0,0),grid_key_dx<2>(sz[1]-2,sz[0]-2)),true);
fd.imposeA(Prs(), g_dist.getSubDomainIterator(grid_key_dx<2>(0,0),grid_key_dx<2>(0,0)));
fd.imposeB(0.0, g_dist.getSubDomainIterator(grid_key_dx<2>(0,0),grid_key_dx<2>(0,0)));
fd.imposeA(vx_eq(), g_dist.getSubDomainIterator(grid_key_dx<2>(0,1),grid_key_dx<2>(sz[1]-2,sz[0]-2)));
fd.imposeB(0.0, g_dist.getSubDomainIterator(grid_key_dx<2>(0,1),grid_key_dx<2>(sz[1]-2,sz[0]-2)));
fd.imposeA(vy_eq(), g_dist.getSubDomainIterator(grid_key_dx<2>(1,0),grid_key_dx<2>(sz[1]-2,sz[0]-2)));
fd.imposeB(0.0, g_dist.getSubDomainIterator(grid_key_dx<2>(1,0),grid_key_dx<2>(sz[1]-2,sz[0]-2)));
// v_x R L T B
fd.impose(ic_eq(),0.0, EQ_3, {0,0},{sz[0]-2,sz[1]-2},true);
fd.impose(Prs(), 0.0, EQ_3, {0,0},{0,0});
fd.impose(vx_eq(),0.0, EQ_1, {1,0},{sz[0]-2,sz[1]-2});
fd.impose(vy_eq(),0.0, EQ_2, {0,1},{sz[0]-2,sz[1]-2});
// v_x
// R L
fd.imposeA(v_x(), g_dist.getSubDomainIterator(grid_key_dx<2>(0,0),grid_key_dx<2>(sz[1]-2,0)));
fd.imposeB(0.0, g_dist.getSubDomainIterator(grid_key_dx<2>(0,0),grid_key_dx<2>(sz[1]-2,0)));
fd.imposeA(v_x(), g_dist.getSubDomainIterator(grid_key_dx<2>(0,sz[0]-1),grid_key_dx<2>(sz[1]-2,sz[0]-1)));
fd.imposeB(1.0, g_dist.getSubDomainIterator(grid_key_dx<2>(0,sz[0]-1),grid_key_dx<2>(sz[1]-2,sz[0]-1)));
fd.impose(v_x(),0.0, EQ_1, {0,0},{0,sz[1]-2});
fd.impose(v_x(),0.0, EQ_1, {sz[0]-1,0},{sz[0]-1,sz[1]-2});
// T B
fd.imposeA(avg_vx(), g_dist.getSubDomainIterator(grid_key_dx<2>(0,0),grid_key_dx<2>(0,sz[0]-1)));
fd.imposeB(0.0, g_dist.getSubDomainIterator(grid_key_dx<2>(0,0),grid_key_dx<2>(0,sz[0]-1)));
fd.imposeA(avg_vx(), g_dist.getSubDomainIterator(grid_key_dx<2>(sz[1]-1,0),grid_key_dx<2>(sz[1]-1,sz[0]-1)));
fd.imposeB(0.0, g_dist.getSubDomainIterator(grid_key_dx<2>(sz[1]-1,0),grid_key_dx<2>(sz[1]-1,sz[0]-1)));
// v_y R L T B
fd.impose(avg_vx_f(),0.0, EQ_1, {0,-1},{sz[0]-1,-1});
fd.impose(avg_vx(),0.0, EQ_1, {0,sz[1]-1},{sz[0]-1,sz[1]-1});
// v_y
// R L
fd.imposeA(avg_vy(), g_dist.getSubDomainIterator(grid_key_dx<2>(0,0),grid_key_dx<2>(sz[1]-1,0)));
fd.imposeB(0.0, g_dist.getSubDomainIterator(grid_key_dx<2>(0,0),grid_key_dx<2>(sz[1]-1,0)));
fd.imposeA(avg_vy(), g_dist.getSubDomainIterator(grid_key_dx<2>(0,sz[0]-1),grid_key_dx<2>(sz[1]-1,sz[0]-1)));
fd.imposeB(0.0, g_dist.getSubDomainIterator(grid_key_dx<2>(0,sz[0]-1),grid_key_dx<2>(sz[1]-1,sz[0]-1)));
fd.impose(avg_vy_f(),0.0, EQ_2 , {-1,0},{-1,sz[1]-1});
fd.impose(avg_vy(),1.0, EQ_2, {sz[0]-1,0},{sz[0]-1,sz[1]-1});
// T B
fd.imposeA(v_y(), g_dist.getSubDomainIterator(grid_key_dx<2>(0,0),grid_key_dx<2>(0,sz[0]-2)));
fd.imposeB(0.0, g_dist.getSubDomainIterator(grid_key_dx<2>(0,0),grid_key_dx<2>(0,sz[0]-2)));
fd.imposeA(v_y(), g_dist.getSubDomainIterator(