Skip to content
Snippets Groups Projects
Commit 173b92ab authored by Pietro Incardona's avatar Pietro Incardona
Browse files

Test working

parent 6a195bfe
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 761 additions and 87079 deletions
src/FiniteDifference/sum.hpp
+ 62
− 37
View file @ 173b92ab
......@@ -13,6 +13,11 @@
#include <unordered_map>
#include "util/for_each_ref.hpp"
/*! \brief sum functor value
*
* \param v_expr vector expression
*
*/
template<typename v_expr>
struct sum_functor_value
{
......@@ -28,10 +33,10 @@
//! Grid info
const grid_sm<last::dims,void> & gs;
// grid mapping
//! grid mapping
const typename stub_or_real<last,last::dims,typename last::stype,typename last::b_grid::decomposition::extended_type>::type & g_map;
// grid position
//! grid position
grid_dist_key_dx<last::dims> & kmap;
//! coefficent
......@@ -42,13 +47,29 @@
/*! \brief constructor
*
* \param g_map Grid mapping, it convert grid position into vector/Matrix row
* \param kmap grid position
* \param gs grid information
* \param spacing grid spacing
* \param cols unordered map contain the map colum -> value
* \param coeff it contain an additional actual coefficients in front of the values
*
*/
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)
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),coeff(coeff),spacing(spacing)
{};
//! It call this function for every expression in the sum
/*! \brief It call this function for every expression operand in the sum
*
* \param t expression operand id
*
*/
template<typename T>
void operator()(T& t) const
{
......@@ -70,27 +91,35 @@
template<typename ... expr>
struct sum
{
// Transform from variadic template to boost mpl vector
//! Transform from variadic template to boost mpl vector
typedef boost::mpl::vector<expr... > v_expr;
// size of v_expr
//! size of v_expr
typedef typename boost::mpl::size<v_expr>::type v_sz;
//! struct that specify the syste, of equations
typedef typename boost::mpl::at<v_expr, boost::mpl::int_<v_sz::type::value - 1> >::type Sys_eqs;
/*! \brief Calculate which colums of the Matrix are non zero
*
* \param pos position where the derivative is calculated
* \param gs Grid info
* \param cols non-zero colums calculated by the function
* \param coeff coefficent (constant in front of the derivative)
* \param g_map Grid mapping, it convert grid position into vector/Matrix row
* \param kmap grid position
* \param gs grid information
* \param spacing grid spacing
* \param cols unordered map contain the map colum -> value
* \param coeff it contain an additional actual coefficients in front of the values
*
* ### Example
*
* \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::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)
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);
......@@ -99,23 +128,15 @@ struct sum
boost::mpl::for_each_ref< boost::mpl::range_c<int,0,v_sz::type::value - 1> >(sm);
}
/*! \brief Calculate the position in the cell where the sum operator is performed
*
* it is done for the first element, the others are supposed to be in the same position
* it just return the position of the staggered property in the first expression
*
* \param position where we are calculating the derivative
* \param gs Grid info
* \param s_pos staggered position of the properties
*
*/
inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
{
return boost::mpl::at<v_expr, boost::mpl::int_<0> >::type::position(pos,gs,s_pos);
}
};
};
/*! \brief It ancapsulate the minus operation
*
* \tparam arg
* \tparam Sys_eqs system of equation
*
*/
template<typename arg, typename Sys_eqs>
struct minus
{
......@@ -123,23 +144,27 @@ struct minus
*
* \tparam ord
*
* \snipper FDScheme_unit_tests.hpp minus example
* \snippet FDScheme_unit_tests.hpp minus example
*
* \param g_map Grid mapping, it convert grid position into vector/Matrix row
* \param kmap grid position
* \param gs grid information
* \param spacing grid spacing
* \param cols unordered map contain the map colum -> value
* \param coeff it contain an additional actual coefficients in front of the values
*
*/
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)
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);
}
/*! \brief Calculate the position where the minus is calculated
*
* it just return the position of the staggered property at first expression
*
*/
inline static grid_key_dx<Sys_eqs::dims> position(grid_key_dx<Sys_eqs::dims> & pos, const grid_sm<Sys_eqs::dims,void> & gs, const comb<Sys_eqs::dims> (& s_pos)[Sys_eqs::nvar])
{
std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " only CENTRAL, FORWARD, BACKWARD derivative are defined";
}
};
#endif /* OPENFPM_NUMERICS_SRC_FINITEDIFFERENCE_SUM_HPP_ */
src/Solvers/petsc_solver.hpp
+ 37
− 5
View file @ 173b92ab
......@@ -875,7 +875,15 @@ public:
*/
void setRelTol(PetscReal rtol_)
{
PetscOptionsSetValue("-ksp_rtol",std::to_string(rtol_).c_str());
PetscReal rtol;
PetscReal abstol;
PetscReal dtol;
PetscInt maxits;
PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol_,abstol,dtol,maxits));
// PetscOptionsSetValue("-ksp_rtol",std::to_string(rtol_).c_str());
}
/*! \brief Set the absolute tolerance as stop criteria
......@@ -887,7 +895,15 @@ public:
*/
void setAbsTol(PetscReal abstol_)
{
PetscOptionsSetValue("-ksp_atol",std::to_string(abstol_).c_str());
PetscReal rtol;
PetscReal abstol;
PetscReal dtol;
PetscInt maxits;
PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol_,dtol,maxits));
// PetscOptionsSetValue("-ksp_atol",std::to_string(abstol_).c_str());
}
/*! \brief Set the divergence tolerance
......@@ -899,7 +915,15 @@ public:
*/
void setDivTol(PetscReal dtol_)
{
PetscOptionsSetValue("-ksp_dtol",std::to_string(dtol_).c_str());
PetscReal rtol;
PetscReal abstol;
PetscReal dtol;
PetscInt maxits;
PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol,dtol_,maxits));
// PetscOptionsSetValue("-ksp_dtol",std::to_string(dtol_).c_str());
}
/*! \brief Set the maximum number of iteration for Krylov solvers
......@@ -909,8 +933,16 @@ public:
*/
void setMaxIter(PetscInt n)
{
PetscOptionsSetValue("-ksp_max_it",std::to_string(n).c_str());
maxits = n;
PetscReal rtol;
PetscReal abstol;
PetscReal dtol;
PetscInt maxits;
PETSC_SAFE_CALL(KSPGetTolerances(ksp,&rtol,&abstol,&dtol,&maxits));
PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol,dtol,n));
// PetscOptionsSetValue("-ksp_max_it",std::to_string(n).c_str());
this->maxits = n;
}
/*! For the BiCGStab(L) it define the number of search directions
......
src/interpolation/interpolation.hpp
+ 480
− 101
View file @ 173b92ab
This diff is collapsed.
src/interpolation/interpolation_broken.hpp deleted 100644 → 0
+ 0
− 629
View file @ 6a195bfe
/*
* mp4_interpolation.hpp
*
* Created on: May 4, 2017
* Author: i-bird
*/
#ifndef OPENFPM_NUMERICS_SRC_INTERPOLATION_INTERPOLATION_HPP_
#define OPENFPM_NUMERICS_SRC_INTERPOLATION_INTERPOLATION_HPP_
#include "NN/CellList/MemFast.hpp"
#include "NN/CellList/CellList.hpp"
#define INTERPOLATION_ERROR_OBJECT std::runtime_error("Runtime interpolation error");
constexpr int inte_m2p = 0;
constexpr int inte_p2m = 1;
template<unsigned int n_ele, typename T>
struct agg_arr
{
T ele[n_ele];
};
/*! \brief multiply the src by coeff for several types T
*
* \tparam type T
*
*/
template<typename T>
struct mul_inte
{
/*! \brief multiply the src by coeff for several types T
*
* \param result the result of the multiplication
* \param coeff coefficent to use for of the multiplication
* \param src source
*
*/
inline static void value(T & result, const T & coeff, const T & src)
{
result += coeff * src;
}
};
/*! \brief multiply the src by coeff for several types T
*
* \tparam type T
*
*/
template<typename T, unsigned int N1>
struct mul_inte<T[N1]>
{
/*! \brief multiply the src by coeff for several types T
*
* \param result the result of the multiplication
* \param coeff coefficent to use for of the multiplication
* \param src source
*
*/
inline static void value(T (& result)[N1], const T & coeff, const T (& src)[N1])
{
for (size_t i = 0 ; i < N1 ; i++)
result[i] += coeff * src[i];
}
};
/*! \brief multiply the src by coeff for several types T
*
* \tparam type T
*
*/
template<typename T, unsigned int N1, unsigned int N2>
struct mul_inte<T[N1][N2]>
{
/*! \brief multiply the src by coeff for several types T
*
* \param result the result of the multiplication
* \param coeff coefficent to use for of the multiplication
* \param src source
*
*/
inline static void value(T (& result)[N1][N2], const T & coeff, const T (& src)[N1][N2])
{
for (size_t i = 0 ; i < N1 ; i++)
for (size_t j = 0 ; j < N2 ; j++)
result[i][j] += coeff * src[i][j];
}
};
/*! \brief Class that select the operation to do differently if we are doing Mesh to particle (m2p) or particle to mesh (p2m)
*
* \tparam prp_g property of the grid to interpolate
* \tparam prp_v property of the vector to interpolate
* \param M2P or P2M
*
*/
template<unsigned int np, unsigned int prp_g, unsigned int prp_v,unsigned int m2p_or_p2m>
struct inte_template
{
/*! \brief Evaluate the interpolation
*
* \tparam np number of kernel points in one direction
* \tparam grid type of grid
* \tparam vector type of vector
* \tparam iterator type of the iterator
*
* \param gd grid for interpolation
* \param vd vector for interpolation
* \param k_dist grid key grid point for interpolation
* \param key_p particle for interpolation
* \param a_int interpolation coefficients pre-calculated
* \param key indicate which pre-calculated coefficient we have to use
*
*/
template<typename grid, typename vector, typename iterator> inline static void value(grid & gd,
vector & vd,
const grid_dist_lin_dx & k_dist,
iterator & key_p,
typename vector::stype (& a_int)[openfpm::math::pow(vector::dims,np)],
const size_t & key)
{
mul_inte<typename std::remove_reference<decltype(gd.template get<prp_g>(k_dist))>::type>::value(gd.template get<prp_g>(k_dist),a_int[key],vd.template getProp<prp_v>(key_p));
}
};
/*! \brief Class that select the operation to do differently if we are doing Mesh to particle (m2p) or particle to mesh (p2m)
*
* \tparam prp_g property of the grid to interpolate
* \tparam prp_v property of the vector to interpolate
* \param M2P or P2M
*
*/
template<unsigned int np, unsigned int prp_g, unsigned int prp_v>
struct inte_template<np,prp_g,prp_v,inte_m2p>
{
/*! \brief Evaluate the interpolation
*
* \tparam grid type of grid
* \tparam vector type of vector
* \tparam iterator type of the iterator
* \tparam key_type type of the key
*
* \param gd grid for interpolation
* \param vd vector for interpolation
* \param k_dist grid key grid point for interpolation
* \param key_p particle for interpolation
* \param a_int interpolation coefficients pre-calculated
* \param key indicate which pre-calculated coefficient we have to use
*
*/
template<typename grid, typename vector, typename iterator> inline static void value(grid & gd,
vector & vd,
const grid_dist_lin_dx & k_dist,
iterator & key_p,
typename vector::stype (& a_int)[openfpm::math::pow(vector::dims,np)],
const size_t & key)
{
mul_inte<typename std::remove_reference<decltype(gd.template get<prp_g>(k_dist))>::type>::value(vd.template getProp<prp_v>(key_p),a_int[key],gd.template get<prp_g>(k_dist));
}
};
/*! \brief Calculate aint
*
*
*/
template<unsigned int dims, typename vector, unsigned int np>
struct calculate_aint
{
static inline void value(size_t (& sz)[vector::dims],
typename vector::stype a_int[openfpm::math::pow(vector::dims,np)],
typename vector::stype (& a)[vector::dims][np])
{
grid_sm<vector::dims,void> gs(sz);
grid_key_dx_iterator<vector::dims> kit(gs);
size_t s = 0;
while (kit.isNext())
{
auto key = kit.get();
a_int[s] = 1;
for (size_t i = 0 ; i < vector::dims ; i++)
a_int[s] *= a[i][key.get(i)];
s++;
++kit;
}
}
};
/*! \brief Calculate aint 2D
*
*
*/
template<typename vector, unsigned int np>
struct calculate_aint<2,vector,np>
{
static inline void value(size_t (& sz)[vector::dims],
typename vector::stype a_int[openfpm::math::pow(vector::dims,np)],
typename vector::stype (& a)[vector::dims][np])
{
size_t s = 0;
for (size_t i = 0 ; i < np ; i++)
{
for (size_t j = 0 ; j < np ; j++)
{
a_int[s] = a[0][j]*a[1][i];
s++;
}
}
}
};
/*! \brief Calculate aint
*
*
*/
template<typename vector, unsigned int np>
struct calculate_aint<3,vector,np>
{
static inline void value(size_t (& sz)[vector::dims],
typename vector::stype a_int[openfpm::math::pow(vector::dims,np)],
typename vector::stype (& a)[vector::dims][np])
{
size_t s = 0;
for (size_t i = 0 ; i < np ; i++)
{
for (size_t j = 0 ; j < np ; j++)
{
for (size_t k = 0 ; k < np ; k++)
{
a_int[s] = a[0][k]*a[1][j]*a[2][i];
s++;
}
}
}
}
};
/*! \brief return the sub-domain where this particle must be interpolated
*
* \param p particle position
*
* \return the sub-domain id
*
*/
template<typename vector, typename grid>
inline size_t getSub(Point<vector::dims,typename vector::stype> & p,
const CellList<vector::dims,typename vector::stype,Mem_fast<vector::dims,typename vector::stype>,shift<vector::dims,typename vector::stype>> & geo_cell,
grid & gd)
{
size_t cell = geo_cell.getCell(p);
for (size_t i = 0 ; i < geo_cell.getNelements(cell) ; i++)
{
size_t ns = geo_cell.get(cell,i);
if (gd.getDecomposition().getSubDomain(ns).isInside(p))
return ns;
}
#ifdef SE_CLASS1
std::cerr << __FILE__ << ":" << __LINE__ << " Error: " << std::endl;
ACTION_ON_ERROR(INTERPOLATION_ERROR_OBJECT)
#endif
return (size_t)-1;
}
template<typename vector,typename kernel>
struct inte_calc_impl
{
//! Type of the calculations
typedef typename vector::stype arr_type;
template<unsigned int prp_g, unsigned int prp_v, unsigned int m2p_or_p2m, typename iterator, typename grid>
static inline void inte_calc(iterator & it,
vector & vd,
const Box<vector::dims,typename vector::stype> & domain,
int (& ip)[vector::dims][kernel::np],
grid & gd,
const typename vector::stype (& dx)[vector::dims],
typename vector::stype (& xp)[vector::dims],
typename vector::stype (& a_int)[openfpm::math::pow(vector::dims,kernel::np)],
typename vector::stype (& a)[vector::dims][kernel::np],
typename vector::stype (& x)[vector::dims][kernel::np],
size_t (& sz)[vector::dims],
const CellList<vector::dims,typename vector::stype,Mem_fast<vector::dims,typename vector::stype>,shift<vector::dims,typename vector::stype>> & geo_cell,
openfpm::vector<agg_arr<openfpm::math::pow(vector::dims,kernel::np),typename vector::stype>> & offsets)
{
auto key_p = it.get();
Point<vector::dims,typename vector::stype> p = vd.getPos(key_p);
// On which sub-domain we interpolate the particle
size_t sub = getSub<vector>(p,geo_cell,gd);
typename vector::stype x0[vector::dims];
// calculate the position of the particle in the grid
// coordinated
for (size_t i = 0 ; i < vector::dims ; i++)
x0[i] = (p.get(i)-domain.getLow(i))*dx[i];
// convert into integer
for (size_t i = 0 ; i < vector::dims ; i++)
ip[i][0] = (int)x0[i];
// convert the global grid position into local grid position
grid_key_dx<vector::dims> base;
for (size_t i = 0 ; i < vector::dims ; i++)
base.set_d(i,ip[i][0] - gd.getLocalGridsInfo().get(sub).origin.get(i) - (long int)kernel::np/2 + 1);
// convenient grid of points
for (size_t j = 0 ; j < kernel::np-1 ; j++)
{
for (size_t i = 0 ; i < vector::dims ; i++)
ip[i][j+1] = (int)ip[i][j]+1;
}
for (size_t i = 0 ; i < vector::dims ; i++)
xp[i] = x0[i] - ip[i][0];
for (long int j = 0 ; j < kernel::np ; j++)
{
for (size_t i = 0 ; i < vector::dims ; i++)
x[i][j] = - xp[i] + typename vector::stype((long int)j - (long int)kernel::np/2 + 1);
}
for (size_t j = 0 ; j < kernel::np ; j++)
{
for (size_t i = 0 ; i < vector::dims ; i++)
a[i][j] = kernel::value(x[i][j],j);
}
calculate_aint<vector::dims,vector,kernel::np>::value(sz,a_int,a);
grid_dist_lin_dx k_dist_lin;
k_dist_lin.setSub(sub);
size_t k = 0;
auto & gs_info = gd.get_loc_grid(sub).getGrid();
size_t lin_base = gs_info.LinId(base);
#pragma omp simd
for (size_t i = 0 ; i < openfpm::math::pow(vector::dims,kernel::np) ; i++)
{
size_t lin = offsets.get(sub).ele[k] + lin_base;
k_dist_lin.getKeyRef() = lin;
inte_template<kernel::np,prp_g,prp_v,m2p_or_p2m>::value(gd,vd,k_dist_lin,key_p,a_int,k);
k++;
}
}
};
/*! \brief Main class for interpolation Particle to mest p2m and Mesh to particle m2p
*
* This function is the main class to interpolate from particle to mesh and mesh to particle
*
* \tparam vector type of vector for interpolation
* \tparam grid type of grid for interpolation
* \tparam interpolation kernel
*
*/
template<typename vector,typename grid, typename kernel>
class interpolate
{
//! Cell list used to convert particles position to sub-domain
CellList<vector::dims,typename vector::stype,Mem_fast<vector::dims,typename vector::stype>,shift<vector::dims,typename vector::stype>> geo_cell;
/*! Structure to order boxes by volume
*
*
*
*/
struct Box_vol
{
//! Box
Box<vector::dims,size_t> bv;
//! Calculated volume
size_t vol;
/*! \brief operator to reorder
*
* \param bv box to compare with
*
* \return true if bv has volume bigger volume
*
*/
bool operator<(const Box_vol & bv)
{
return vol < bv.vol;
}
};
//! particles
vector & vd;
//! grid
grid & gd;
//! Type of the calculations
typedef typename vector::stype arr_type;
/*! \brief It calculate the interpolation stencil offsets
*
* \param array where to store offsets
* \param sz kernel stencil points in each direction
*
*/
void calculate_the_offsets(openfpm::vector<agg_arr<openfpm::math::pow(vector::dims,kernel::np),typename vector::stype>> & offsets, size_t (& sz)[vector::dims])
{
offsets.resize(gd.getN_loc_grid());
for (size_t i = 0 ; i < offsets.size() ; i++)
{
auto & loc_grid = gd.get_loc_grid(i);
const grid_sm<vector::dims,void> & gs = loc_grid.getGrid();
grid_sm<vector::dims,void> gs2(sz);
grid_key_dx_iterator<vector::dims> kit2(gs2);
size_t k = 0;
while (kit2.isNext())
{
auto key = kit2.get();
long int lin = gs.LinId(key);
offsets.get(i).ele[k] = lin;
++k;
++kit2;
}
}
}
public:
/*! \brief construct an interpolation object between a grid and a vector
*
* When possible and easy to do we suggest to retain the object
*
* \param vd interpolation vector
* \param gd interpolation grid
*
*/
interpolate(vector & vd, grid & gd)
:vd(vd),gd(gd)
{
// get the processor bounding box in grid units
Box<vector::dims,typename vector::stype> bb = gd.getDecomposition().getProcessorBounds();
Box<vector::dims,typename vector::stype> bunit = gd.getDecomposition().getCellDecomposer().getCellBox();
size_t div[vector::dims];
for (size_t i = 0 ; i < vector::dims ; i++)
div[i] = (bb.getHigh(i) - bb.getLow(i)) / bunit.getHigh(i);
geo_cell.Initialize(bb,div);
// Now draw the domain into the cell list
auto & dec = gd.getDecomposition();
for (size_t i = 0 ; i < dec.getNSubDomain() ; i++)
{
const Box<vector::dims,typename vector::stype> & bx = dec.getSubDomain(i);
// get the cells this box span
const grid_key_dx<vector::dims> p1 = geo_cell.getCellGrid(bx.getP1());
const grid_key_dx<vector::dims> p2 = geo_cell.getCellGrid(bx.getP2());
// Get the grid and the sub-iterator
auto & gi = geo_cell.getGrid();
grid_key_dx_iterator_sub<vector::dims> g_sub(gi,p1,p2);
// add the box-id to the cell list
while (g_sub.isNext())
{
auto key = g_sub.get();
geo_cell.addCell(gi.LinId(key),i);
++g_sub;
}
}
};
/*! \brief Interpolate particles to mesh
*
* Most of the time the particle set and the mesh are the same
* as the one used in the constructor. They also can be different
* as soon as they have the same decomposition
*
* \param gd grid or mesh
* \param vd particle set
*
*/
template<unsigned int prp_v, unsigned int prp_g> void p2m(vector & vd, grid & gd)
{
#ifdef SE_CLASS1
if (!vd.getDecomposition().is_equal_ng(gd.getDecomposition()) )
{
std::cerr << __FILE__ << ":" << __LINE__ << " Error: the distribution of the vector of particles" <<
" and the grid is different. In order to interpolate the two data structure must have the" <<
" same decomposition" << std::endl;
ACTION_ON_ERROR(INTERPOLATION_ERROR_OBJECT)
}
#endif
Box<vector::dims,typename vector::stype> domain = vd.getDecomposition().getDomain();
// grid spacing
typename vector::stype dx[vector::dims];
for (size_t i = 0 ; i < vector::dims ; i++)
dx[i] = 1.0/gd.spacing(i);
size_t sz[vector::dims];
for (size_t i = 0 ; i < vector::dims ; i++)
sz[i] = kernel::np;
// Precalculate the offset for each sub-sub-domain
openfpm::vector<agg_arr<openfpm::math::pow(vector::dims,kernel::np),typename vector::stype>> offsets;
calculate_the_offsets(offsets,sz);
// point position
typename vector::stype xp[vector::dims];
int ip[vector::dims][kernel::np];
typename vector::stype x[vector::dims][kernel::np];
typename vector::stype a[vector::dims][kernel::np];
// grid_cpu<vector::dims,aggregate<typename vector::stype>> a_int(sz);
// a_int.setMemory();
typename vector::stype a_int[openfpm::math::pow(vector::dims,kernel::np)];
auto it = vd.getDomainIterator();
while (it.isNext() == true)
{
inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_p2m>(it,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
++it;
}
}
/*! \brief Interpolate mesh to particle
*
* Most of the time the particle set and the mesh are the same
* as the one used in the constructor. They also can be different
* as soon as they have the same decomposition
*
* \param gd grid or mesh
* \param vd particle set
*
*/
template<unsigned int prp_g, unsigned int prp_v> void m2p(grid & gd, vector & vd)
{
#ifdef SE_CLASS1
if (!vd.getDecomposition().is_equal_ng(gd.getDecomposition()) )
{
std::cerr << __FILE__ << ":" << __LINE__ << " Error: the distribution of the vector of particles" <<
" and the grid is different. In order to interpolate the two data structure must have the" <<
" same decomposition" << std::endl;
ACTION_ON_ERROR(INTERPOLATION_ERROR_OBJECT)
}
#endif
Box<vector::dims,typename vector::stype> domain = vd.getDecomposition().getDomain();
// grid spacing
typename vector::stype dx[vector::dims];
for (size_t i = 0 ; i < vector::dims ; i++)
dx[i] = 1.0/gd.spacing(i);
// point position
typename vector::stype xp[vector::dims];
int ip[vector::dims][kernel::np];
typename vector::stype x[vector::dims][kernel::np];
typename vector::stype a[vector::dims][kernel::np];
size_t sz[vector::dims];
for (size_t i = 0 ; i < vector::dims ; i++)
sz[i] = kernel::np;
// Precalculate the offset for each sub-sub-domain
openfpm::vector<agg_arr<openfpm::math::pow(vector::dims,kernel::np),typename vector::stype>> offsets;
calculate_the_offsets(offsets,sz);
// grid_cpu<vector::dims,aggregate<typename vector::stype>> a_int(sz);
// a_int.setMemory();
typename vector::stype a_int[openfpm::math::pow(vector::dims,kernel::np)];
auto it = vd.getDomainIterator();
while (it.isNext() == true)
{
inte_calc_impl<vector,kernel>::template inte_calc<prp_g,prp_v,inte_m2p>(it,vd,domain,ip,gd,dx,xp,a_int,a,x,sz,geo_cell,offsets);
++it;
}
}
};
#endif /* OPENFPM_NUMERICS_SRC_INTERPOLATION_INTERPOLATION_HPP_ */
src/interpolation/interpolation_unit_tests.hpp
+ 182
− 3
View file @ 173b92ab
......@@ -83,7 +83,7 @@ template<typename vector, unsigned int mom_p> void momenta_vector(vector & vd,ty
}
BOOST_AUTO_TEST_CASE( interpolation_full_test )
BOOST_AUTO_TEST_CASE( interpolation_full_test_2D )
{
Box<2,float> domain({0.0,0.0},{1.0,1.0});
size_t sz[2] = {64,64};
......@@ -107,9 +107,8 @@ BOOST_AUTO_TEST_CASE( interpolation_full_test )
vd.getPos(p)[0] = (double)rand()/RAND_MAX;
vd.getPos(p)[1] = (double)rand()/RAND_MAX;
vd.getPos(p)[2] = (double)rand()/RAND_MAX;
vd.getProp<0>(p) = 5.0/*(double)rand()/RAND_MAX*/;
vd.getProp<0>(p) = 5.0;
++it;
}
......@@ -249,6 +248,186 @@ BOOST_AUTO_TEST_CASE( interpolation_full_test )
}
}
BOOST_AUTO_TEST_CASE( interpolation_full_test_3D )
{
Box<3,double> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
size_t sz[3] = {64,64,64};
Ghost<3,long int> gg(2);
Ghost<3,double> gv(0.01);
size_t bc_v[3] = {PERIODIC,PERIODIC,PERIODIC};
{
vector_dist<3,double,aggregate<double>> vd(65536,domain,bc_v,gv);
grid_dist_id<3,double,aggregate<double>> gd(vd.getDecomposition(),sz,gg);
// set one particle on vd
auto it = vd.getDomainIterator();
while (it.isNext())
{
auto p = it.get();
vd.getPos(p)[0] = (double)rand()/RAND_MAX;
vd.getPos(p)[1] = (double)rand()/RAND_MAX;
vd.getPos(p)[2] = (double)rand()/RAND_MAX;
vd.getProp<0>(p) = 5.0;
++it;
}
vd.map();
// Reset the grid
auto it2 = gd.getDomainGhostIterator();
while (it2.isNext())
{
auto key = it2.get();
gd.template get<0>(key) = 0.0;
++it2;
}
interpolate<decltype(vd),decltype(gd),mp4_kernel<double>> inte(vd,gd);
inte.p2m<0,0>(vd,gd);
double mg[3];
double mv[3];
momenta_grid<decltype(gd),0>(gd,mg);
momenta_vector<decltype(vd),0>(vd,mv);
BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
BOOST_REQUIRE_CLOSE(mg[2],mv[2],0.001);
momenta_grid<decltype(gd),1>(gd,mg);
momenta_vector<decltype(vd),1>(vd,mv);
BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
BOOST_REQUIRE_CLOSE(mg[2],mv[2],0.001);
momenta_grid<decltype(gd),2>(gd,mg);
momenta_vector<decltype(vd),2>(vd,mv);
BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
BOOST_REQUIRE_CLOSE(mg[2],mv[2],0.001);
auto & v_cl = create_vcluster();
// We have to do a ghost get before interpolating m2p
// Before doing mesh to particle particle must be arranged
// into a grid like
vd.clear();
auto it4 = vd.getGridIterator(sz);
while(it4.isNext())
{
auto key = it4.get();
vd.add();
vd.getLastPos()[0] = key.get(0) * it4.getSpacing(0) + domain.getLow(0) + 0.1*it4.getSpacing(0);
vd.getLastPos()[1] = key.get(1) * it4.getSpacing(1) + domain.getLow(1) + 0.1*it4.getSpacing(1);
vd.getLastPos()[2] = key.get(2) * it4.getSpacing(2) + domain.getLow(2) + 0.1*it4.getSpacing(2);
vd.getLastProp<0>() = 0.0;
++it4;
}
// Reset also the grid
auto it5 = gd.getDomainGhostIterator();
while(it5.isNext())
{
auto key = it5.get();
gd.get<0>(key) = 0.0;
++it5;
}
gd.ghost_get<0>();
grid_key_dx<3> start({3,3,3});
grid_key_dx<3> stop({(long int)gd.size(0) - 4,(long int)gd.size(1) - 4,(long int)gd.size(2) - 4});
auto it6 = gd.getSubDomainIterator(start,stop);
while(it6.isNext())
{
auto key = it6.get();
gd.get<0>(key) = 5.0;
++it6;
}
gd.ghost_get<0>();
vd.map();
gd.ghost_get<0>();
inte.m2p<0,0>(gd,vd);
momenta_grid_domain<decltype(gd),0>(gd,mg);
momenta_vector<decltype(vd),0>(vd,mv);
v_cl.sum(mg[0]);
v_cl.sum(mg[1]);
v_cl.sum(mg[2]);
v_cl.sum(mv[0]);
v_cl.sum(mv[1]);
v_cl.sum(mv[2]);
v_cl.execute();
BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
BOOST_REQUIRE_CLOSE(mg[2],mv[2],0.001);
momenta_grid_domain<decltype(gd),1>(gd,mg);
momenta_vector<decltype(vd),1>(vd,mv);
v_cl.sum(mg[0]);
v_cl.sum(mg[1]);
v_cl.sum(mg[2]);
v_cl.sum(mv[0]);
v_cl.sum(mv[1]);
v_cl.sum(mv[2]);
v_cl.execute();
BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
BOOST_REQUIRE_CLOSE(mg[2],mv[2],0.001);
momenta_grid_domain<decltype(gd),2>(gd,mg);
momenta_vector<decltype(vd),2>(vd,mv);
v_cl.sum(mg[0]);
v_cl.sum(mg[1]);
v_cl.sum(mg[2]);
v_cl.sum(mv[0]);
v_cl.sum(mv[1]);
v_cl.sum(mv[2]);
v_cl.execute();
BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
BOOST_REQUIRE_CLOSE(mg[2],mv[2],0.001);
}
}
BOOST_AUTO_TEST_CASE( int_kernel_test )
{
mp4_kernel<float> mp4;
......
test/umfpack_solver<double>_lid_driven_cavity_p2_grid_0_test.vtk deleted 100644 → 0
+ 0
− 43712
View file @ 6a195bfe
This diff is collapsed.
test/umfpack_solver<double>_lid_driven_cavity_p2_grid_1_test.vtk deleted 100644 → 0
+ 0
− 42592
View file @ 6a195bfe
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment