Commit 8980f2a5 by incardon

### Reverting to previous interpolation untill clarification

parent 0f84683a
 ... ... @@ -717,6 +717,23 @@ class petsc_solver KSPDestroy(&ksp); } /*! \brief Return the norm error of the solution * * \param x_ the solution * \param b_ the right-hand-side * * \return the solution error * */ static solError getSolNormError(const Vec & b_, const Vec & x_,KSP ksp) { Mat A_; Mat P_; KSPGetOperators(ksp,&A_,&P_); return getSolNormError(A_,b_,x_); } /*! \brief Return the norm error of the solution * * \param A_ the matrix that identity the linear system ... ... @@ -1259,6 +1276,20 @@ public: return getSolNormError(A.getMat(),b.getVec(),x.getVec()); } /*! \brief It return the resiual error * * \param A Sparse matrix * \param x solution * \param b right-hand-side * * \return the solution error norms * */ solError get_residual_error(const Vector & x, const Vector & b) { return getSolNormError(b.getVec(),x.getVec(),ksp); } /*! \brief Here we invert the matrix and solve the system * * \param A sparse matrix ... ...
 ... ... @@ -16,48 +16,18 @@ constexpr int inte_m2p = 0; constexpr int inte_p2m = 1; template struct agg_arr { T ele[n_ele]; }; /*! \brief multiply the src by coeff for several types T * * \tparam type T * */ template 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 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)[N1], const T & coeff, const T (& src)[N1]) { for (size_t i = 0 ; i < N1 ; i++) ... ... @@ -65,21 +35,9 @@ struct mul_inte } }; /*! \brief multiply the src by coeff for several types T * * \tparam type T * */ template 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)[N1][N2], const T & coeff, const T (& src)[N1][N2]) { for (size_t i = 0 ; i < N1 ; i++) ... ... @@ -88,158 +46,34 @@ struct mul_inte } }; /*! \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 template 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 inline static void value(grid & gd, template inline static void value(grid & gd, vector & vd, const grid_dist_lin_dx & k_dist, const grid_dist_key_dx & k_dist, iterator & key_p, typename vector::stype (& a_int)[openfpm::math::pow(vector::dims,np)], const size_t & key) const grid_cpu> & a_int, const key_type & key) { mul_inte(k_dist))>::type>::value(gd.template get(k_dist),a_int[key],vd.template getProp(key_p)); mul_inte(k_dist))>::type>::value(gd.template get(k_dist),a_int.template get<0>(key),vd.template getProp(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 struct inte_template template struct inte_template { /*! \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 inline static void value(grid & gd, template inline static void value(grid & gd, vector & vd, const grid_dist_lin_dx & k_dist, const grid_dist_key_dx & k_dist, iterator & key_p, typename vector::stype (& a_int)[openfpm::math::pow(vector::dims,np)], const size_t & key) { mul_inte(k_dist))>::type>::value(vd.template getProp(key_p),a_int[key],gd.template get(k_dist)); } }; /*! \brief Calculate aint * * */ template 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 gs(sz); grid_key_dx_iterator 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 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]) const grid_cpu> & a_int, const key_type & key) { 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++; } } mul_inte(k_dist))>::type>::value(vd.template getProp(key_p),a_int.template get<0>(key),gd.template get(k_dist)); } }; /*! \brief Calculate aint * * */ template 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 * ... ... @@ -271,133 +105,113 @@ inline size_t getSub(Point & p, return (size_t)-1; } template struct inte_calc_impl template inline void inte_calc(iterator & it, vector & vd, const Box & domain, int (& ip)[vector::dims][kernel::np], grid & gd, const typename vector::stype (& dx)[vector::dims], typename vector::stype (& xp)[vector::dims], const grid_inte & a_int, typename vector::stype (& a)[vector::dims][kernel::np], typename vector::stype (& x)[vector::dims][kernel::np], size_t (& sz)[vector::dims], const CellList,shift> & geo_cell) { //! Type of the calculations typedef typename vector::stype arr_type; template static inline void inte_calc(iterator & it, vector & vd, const Box & 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,shift> & geo_cell, openfpm::vector> & offsets) { auto key_p = it.get(); auto key_p = it.get(); Point p = vd.getPos(key_p); Point p = vd.getPos(key_p); // On which sub-domain we interpolate the particle size_t sub = getSub(p,geo_cell,gd); // On which sub-domain we interpolate the particle size_t sub = getSub(p,geo_cell,gd); typename vector::stype x0[vector::dims]; 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]; // 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 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 base; // convert the global grid position into local grid position grid_key_dx 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++) base.set_d(i,ip[i][0] - gd.getLocalGridsInfo().get(sub).origin.get(i) - (long int)kernel::np/2 + 1); ip[i][j+1] = (int)ip[i][j]+1; } // convenient grid of points for (size_t i = 0 ; i < vector::dims ; i++) xp[i] = x0[i] - ip[i][0]; 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 (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++) xp[i] = x0[i] - ip[i][0]; a[i][j] = kernel::value(x[i][j],j); } 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); } grid_sm gs(sz); grid_key_dx_iterator kit(gs); 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); } while (kit.isNext()) { auto key = kit.get(); calculate_aint::value(sz,a_int,a); a_int.template get<0>(key) = 1; grid_dist_lin_dx k_dist_lin; k_dist_lin.setSub(sub); for (size_t i = 0 ; i < vector::dims ; i++) a_int.template get<0>(key) *= a[i][key.get(i)]; size_t k = 0; auto & gs_info = gd.get_loc_grid(sub).getGrid(); ++kit; } size_t lin_base = gs_info.LinId(base); grid_key_dx_iterator kit2(gs); grid_dist_key_dx k_dist; k_dist.setSub(sub); #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; while (kit2.isNext()) { auto key = kit2.get(); inte_template::value(gd,vd,k_dist_lin,key_p,a_int,k); for (size_t i = 0 ; i < vector::dims ; i++) k_dist.getKeyRef().set_d(i,key.get(i) + base.get(i)); k++; } inte_template::value(gd,vd,k_dist,key_p,a_int,key); ++kit2; } }; } /*! \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 class interpolate { //! Cell list used to convert particles position to sub-domain CellList,shift> geo_cell; /*! Structure to order boxes by volume * * * */ struct Box_vol { //! Box Box 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) void operator<(const Box_vol & bv) { return vol < bv.vol; } ... ... @@ -409,53 +223,9 @@ class interpolate //! 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> & 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 & gs = loc_grid.getGrid(); grid_sm gs2(sz); grid_key_dx_iterator 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) { ... ... @@ -496,16 +266,6 @@ public: } }; /*! \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 void p2m(vector & vd, grid & gd) { #ifdef SE_CLASS1 ... ... @@ -529,16 +289,6 @@ public: 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> offsets; calculate_the_offsets(offsets,sz); // point position typename vector::stype xp[vector::dims]; ... ... @@ -546,30 +296,24 @@ public: typename vector::stype x[vector::dims][kernel::np]; typename vector::stype a[vector::dims][kernel::np]; // grid_cpu> a_int(sz); // a_int.setMemory(); typename vector::stype a_int[openfpm::math::pow(vector::dims,kernel::np)]; size_t sz[vector::dims]; for (size_t i = 0 ; i < vector::dims ; i++) sz[i] = kernel::np; grid_cpu> a_int(sz); a_int.setMemory();