From 47e20e280cb6340c409e154e97b28bc8fc2868bb Mon Sep 17 00:00:00 2001
From: absingh <absingh@mpi-cbg.de>
Date: Fri, 17 Jan 2020 14:55:17 +0100
Subject: [PATCH] Adding Advection, Laplacian and Vectorizing DCPSE Code
---
src/DCPSE_op/DCPSE_copy/Dcpse.hpp | 48 ++-
src/DCPSE_op/DCPSE_op.hpp | 334 +++++++++++++++++-
src/DCPSE_op/DCPSE_op_test.cpp | 252 ++++++++++++-
.../Vector/vector_dist_operators.hpp | 3 +
4 files changed, 610 insertions(+), 27 deletions(-)
diff --git a/src/DCPSE_op/DCPSE_copy/Dcpse.hpp b/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
index 62e74db5..24f07ae4 100644
--- a/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
+++ b/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
@@ -17,8 +17,11 @@
template<unsigned int dim, typename vector_type>
class Dcpse {
+public:
+
typedef typename vector_type::stype T;
typedef typename vector_type::value_type part_type;
+ typedef vector_type vtype;
// This works in this way:
// 1) User constructs this by giving a domain of points (where one of the properties is the value of our f),
@@ -35,6 +38,7 @@ private:
std::vector<T> localEps; // Each MPI rank has just access to the local ones
public:
+
// Here we require the first element of the aggregate to be:
// 1) the value of the function f on the point
Dcpse(vector_type &particles,
@@ -100,6 +104,26 @@ public:
}
+ template<bool cond>
+ struct is_scalar
+ {
+ template<typename op_type>
+ static auto analyze(const vect_dist_key_dx &key, op_type & o1) -> typename std::remove_reference<decltype(o1.value(key))>::type
+ {
+ return o1.value(key);
+ };
+ };
+
+ template<>
+ struct is_scalar<false>
+ {
+ template<typename op_type>
+ static auto analyze(const vect_dist_key_dx &key, op_type & o1) -> typename std::remove_reference<decltype(o1.value(key))>::type
+ {
+ return o1.value(key);
+ };
+ };
+
/**
* Computes the value of the differential operator on all the particles,
* using the f values stored at the fValuePos position in the aggregate
@@ -109,8 +133,14 @@ public:
* @param particles The set of particles to iterate over.
*/
template<typename op_type>
- T computeDifferentialOperator(const vect_dist_key_dx &key, op_type & o1) {
- char sign = 1;
+ auto computeDifferentialOperator(const vect_dist_key_dx &key, op_type & o1) -> decltype(is_scalar<std::is_fundamental<decltype(o1.value(key))>::value>::analyze(key,o1))
+ {
+
+ typedef decltype(is_scalar<std::is_fundamental<decltype(o1.value(key))>::value>::analyze(key,o1)) expr_type;
+
+ //typedef typename decltype(o1.value(key))::blabla blabla;
+
+ T sign = 1.0;
if (differentialOrder % 2 == 0) {
sign = -1;
}
@@ -119,19 +149,19 @@ public:
auto & particles = o1.getVector();
- T Dfxp = 0;
+ expr_type Dfxp = 0;
Support<dim, T, part_type> support = localSupports[key.getKey()];
size_t xpK = support.getReferencePointKey();
Point<dim, T> xp = support.getReferencePoint();
- T fxp = sign * o1.value(key);
+ expr_type fxp = sign * o1.value(key);
for (auto &xqK : support.getKeys()) {
Point<dim, T> xq = particles.getPos(xqK);
- T fxq = o1.value(vect_dist_key_dx(xqK));
+ expr_type fxq = o1.value(vect_dist_key_dx(xqK));
Point<dim, T> normalizedArg = (xp - xq) / eps;
EMatrix<T, Eigen::Dynamic, 1> &a = localCoefficients[key.getKey()];
- Dfxp += (fxq + fxp) * computeKernel(normalizedArg, a);
+ Dfxp = Dfxp + (fxq + fxp) * computeKernel(normalizedArg, a);
}
- Dfxp /= pow(eps, differentialOrder);
+ Dfxp = Dfxp / pow(eps, differentialOrder);
//
//T trueDfxp = particles.template getProp<2>(xpK);
// Store Dfxp in the right position
@@ -211,13 +241,13 @@ private:
// Get the points in the support of the DCPSE kernel and store the support for reuse
Support<dim, T, part_type> support = supportBuilder.getSupport(it, requiredSupportSize);
EMatrix<T, Eigen::Dynamic, Eigen::Dynamic> V(support.size(), monomialBasis.size());
-
+/* Some Debug code
if (it.get().getKey() == 5564)
{
int debug = 0;
debug++;
}
-
+*/
// Vandermonde matrix computation
Vandermonde<dim, T, EMatrix<T, Eigen::Dynamic, Eigen::Dynamic>>
vandermonde(support, monomialBasis);
diff --git a/src/DCPSE_op/DCPSE_op.hpp b/src/DCPSE_op/DCPSE_op.hpp
index 449a1ab7..85f2f621 100644
--- a/src/DCPSE_op/DCPSE_op.hpp
+++ b/src/DCPSE_op/DCPSE_op.hpp
@@ -10,6 +10,8 @@
#include "DCPSE_op/DCPSE_copy/Dcpse.hpp"
#include "Operators/Vector/vector_dist_operators.hpp"
+const double dcpse_scaling_factor = 2;
+
/*! \brief Subtraction operation
*
* \tparam exp1 expression1
@@ -31,6 +33,47 @@ public:
:o1(o1),dcp(dcp)
{}
+ /*! \brief This function must be called before value
+ *
+ * it initialize the expression if needed
+ *
+ */
+ inline void init() const
+ {
+ o1.init();
+ }
+
+ /*! \brief Evaluate the expression
+ *
+ * \param key where to evaluate the expression
+ *
+ * \return the result of the expression
+ *
+ */
+ template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx()))>::type > inline r_type value(const vect_dist_key_dx & key) const
+ {
+ return dcp.computeDifferentialOperator(key,o1);
+ }
+};
+
+template <typename exp1,typename DCPSE_type>
+class vector_dist_expression_op<exp1,DCPSE_type,VECT_DCPSE_V>
+{
+ //! expression 1
+ const exp1 o1;
+
+ DCPSE_type (& dcp)[DCPSE_type::vtype::dims];
+
+ static const int dims = DCPSE_type::vtype::dims;
+ typedef typename DCPSE_type::vtype::stype stype;
+
+public:
+
+ //! Costruct a subtraction expression out of two expressions
+ inline vector_dist_expression_op(const exp1 & o1, DCPSE_type (& dcp)[DCPSE_type::vtype::dims])
+ :o1(o1),dcp(dcp)
+ {}
+
/*! \brief This function must be called before value
*
* it initialize the expression if needed
@@ -48,9 +91,99 @@ public:
* \return the result of the expression
*
*/
- template<typename r_type=typename std::remove_reference<decltype(o1.value(vect_dist_key_dx()))>::type > inline r_type value(const vect_dist_key_dx & key) const
+ template<typename r_type=VectorS<dims,stype> > inline r_type value(const vect_dist_key_dx & key) const
{
- return dcp.computeDifferentialOperator(key,o1);
+ VectorS<dims,stype> v_grad;
+
+ for (int i = 0 ; i < dims ; i++)
+ {
+ v_grad.get(i) = dcp[i].computeDifferentialOperator(key,o1);
+ }
+
+ return v_grad;
+ }
+};
+
+template <typename exp1,typename DCPSE_type>
+class vector_dist_expression_op<exp1,DCPSE_type,VECT_DCPSE_V_SUM>
+{
+ //! expression 1
+ const exp1 o1;
+
+ DCPSE_type (& dcp)[DCPSE_type::vtype::dims];
+
+ static const int dims = DCPSE_type::vtype::dims;
+ typedef typename DCPSE_type::vtype::stype stype;
+
+public:
+
+ inline vector_dist_expression_op(const exp1 & o1, DCPSE_type (& dcp)[DCPSE_type::vtype::dims])
+ :o1(o1),dcp(dcp)
+ {}
+
+ inline void init() const
+ {
+ o1.init();
+ }
+
+ template<typename r_type=VectorS<dims,stype> > inline r_type value(const vect_dist_key_dx & key) const
+ {
+ //typedef typename std::remove_reference<decltype(o1.value(key))>::type::blabla blabla;
+
+ typename std::remove_reference<decltype(o1.value(key))>::type v_lap;
+ v_lap=0.0;
+
+
+ for (int i = 0 ; i < dims ; i++)
+ {
+ v_lap += dcp[i].computeDifferentialOperator(key,o1);
+ }
+
+ return v_lap;
+ }
+};
+
+
+
+template <typename exp1,typename exp2_pr>
+class vector_dist_expression_op<exp1,exp2_pr,VECT_DCPSE_V_DOT>
+{
+ typedef typename std::tuple_element<1,exp2_pr>::type DCPSE_type;
+ typedef typename std::tuple_element<0,exp2_pr>::type exp2;
+
+ //! expression 1
+ const exp1 o1;
+ const exp2 o2;
+
+ DCPSE_type (& dcp)[DCPSE_type::vtype::dims];
+
+ static const int dims = DCPSE_type::vtype::dims;
+ typedef typename DCPSE_type::vtype::stype stype;
+
+public:
+
+ inline vector_dist_expression_op(const exp1 & o1,const exp2 & o2, DCPSE_type (& dcp)[DCPSE_type::vtype::dims])
+ :o1(o1),o2(o2),dcp(dcp)
+ {}
+
+ inline void init() const
+ {
+ o1.init();
+ o2.init();
+ }
+
+ template<typename r_type=VectorS<dims,stype> > inline r_type value(const vect_dist_key_dx & key) const
+ {
+ //typedef typename std::remove_reference<decltype(o1.value(key))>::type::blabla blabla;
+ typename std::remove_reference<decltype(o1.value(key))>::type adv;
+ adv=0.0;
+ for (int i = 0 ; i < dims ; i++)
+ {
+ adv += o1.value(key)[i]*dcp[i].computeDifferentialOperator(key,o2);
+
+ }
+ return adv;
+
}
};
/*
@@ -75,13 +208,13 @@ class Derivative_x
public:
template<typename particles_type>
- Derivative_x(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut)
+ Derivative_x(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor)
{
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(0) = 1;
- dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,2.5);
+ dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor);
}
template<typename operand_type>
@@ -95,6 +228,199 @@ public:
};
+class Derivative_y
+{
+
+ void * dcpse;
+
+public:
+
+ template<typename particles_type>
+ Derivative_y(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor)
+ {
+ Point<particles_type::dims,unsigned int> p;
+ p.zero();
+ p.get(1) = 1;
+
+ dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor);
+ }
+
+ template<typename operand_type>
+
+ vector_dist_expression_op<operand_type,Dcpse<operand_type::vtype::dims,typename operand_type::vtype>,VECT_DCPSE> operator()(operand_type arg)
+ {
+ typedef Dcpse<operand_type::vtype::dims,typename operand_type::vtype> dcpse_type;
+
+ return vector_dist_expression_op<operand_type,dcpse_type,VECT_DCPSE>(arg,*(dcpse_type *)dcpse);
+ }
+};
+class Derivative_z
+{
+
+ void * dcpse;
+
+public:
+
+ template<typename particles_type>
+ Derivative_z(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor)
+ {
+ Point<particles_type::dims,unsigned int> p;
+ p.zero();
+ p.get(2) = 1;
+
+ dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor);
+ }
+
+ template<typename operand_type>
+
+ vector_dist_expression_op<operand_type,Dcpse<operand_type::vtype::dims,typename operand_type::vtype>,VECT_DCPSE> operator()(operand_type arg)
+ {
+ typedef Dcpse<operand_type::vtype::dims,typename operand_type::vtype> dcpse_type;
+
+ return vector_dist_expression_op<operand_type,dcpse_type,VECT_DCPSE>(arg,*(dcpse_type *)dcpse);
+ }
+};
+
+
+
+
+
+class Gradient
+{
+
+ void * dcpse;
+
+public:
+
+ template<typename particles_type>
+ Gradient(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor)
+ {
+ typedef Dcpse<particles_type::dims,particles_type> DCPSE_type;
+
+ dcpse = new unsigned char [particles_type::dims*sizeof(DCPSE_type)];
+
+ Dcpse<particles_type::dims,particles_type> * dcpse_ptr = (Dcpse<particles_type::dims,particles_type> *)dcpse;
+
+ for (int i = 0 ; i < particles_type::dims ; i++)
+ {
+ Point<particles_type::dims,unsigned int> p;
+ p.zero();
+ p.get(i) = 1;
+ new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor);
+ dcpse_ptr++;
+ }
+ }
+
+ template<typename operand_type>
+ vector_dist_expression_op<operand_type,Dcpse<operand_type::vtype::dims,typename operand_type::vtype>,VECT_DCPSE_V> operator()(operand_type arg)
+ {
+ typedef Dcpse<operand_type::vtype::dims,typename operand_type::vtype> dcpse_type;
+
+ return vector_dist_expression_op<operand_type,dcpse_type,VECT_DCPSE_V>(arg,*(dcpse_type(*)[operand_type::vtype::dims])dcpse);
+ }
+};
+
+class Laplacian
+{
+
+ void * dcpse;
+
+public:
+
+ template<typename particles_type>
+ Laplacian(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor)
+ {
+ typedef Dcpse<particles_type::dims,particles_type> DCPSE_type;
+
+ dcpse = new unsigned char [particles_type::dims*sizeof(DCPSE_type)];
+
+ Dcpse<particles_type::dims,particles_type> * dcpse_ptr = (Dcpse<particles_type::dims,particles_type> *)dcpse;
+
+ for (int i = 0 ; i < particles_type::dims ; i++)
+ {
+ Point<particles_type::dims,unsigned int> p;
+ p.zero();
+ p.get(i) = 2;
+ new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor);
+ dcpse_ptr++;
+ }
+ }
+
+ template<typename operand_type>
+ vector_dist_expression_op<operand_type,Dcpse<operand_type::vtype::dims,typename operand_type::vtype>,VECT_DCPSE_V_SUM> operator()(operand_type arg)
+ {
+ typedef Dcpse<operand_type::vtype::dims,typename operand_type::vtype> dcpse_type;
+
+ return vector_dist_expression_op<operand_type,dcpse_type,VECT_DCPSE_V_SUM>(arg,*(dcpse_type(*)[operand_type::vtype::dims])dcpse);
+ }
+};
+
+
+class Advection
+{
+
+ void * dcpse;
+
+public:
+
+ template<typename particles_type>
+ Advection(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor)
+ {
+ typedef Dcpse<particles_type::dims,particles_type> DCPSE_type;
+
+ dcpse = new unsigned char [particles_type::dims*sizeof(DCPSE_type)];
+
+ Dcpse<particles_type::dims,particles_type> * dcpse_ptr = (Dcpse<particles_type::dims,particles_type> *)dcpse;
+
+ for (int i = 0 ; i < particles_type::dims ; i++)
+ {
+ Point<particles_type::dims,unsigned int> p;
+ p.zero();
+ p.get(i) = 1;
+ new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor);
+ dcpse_ptr++;
+ }
+
+
+ }
+
+ template<typename operand_type1,typename operand_type2>
+ vector_dist_expression_op<operand_type1,std::pair<operand_type2,Dcpse<operand_type2::vtype::dims,typename operand_type2::vtype>>,VECT_DCPSE_V_DOT> operator()(operand_type1 arg, operand_type2 & arg2)
+ {
+ typedef Dcpse<operand_type2::vtype::dims,typename operand_type2::vtype> dcpse_type;
+
+ return vector_dist_expression_op<operand_type1,std::pair<operand_type2,dcpse_type>,VECT_DCPSE_V_DOT>(arg,arg2,*(dcpse_type(*)[operand_type2::vtype::dims])dcpse);
+ }
+};
+
+
+class Derivative_xy
+{
+
+ void * dcpse;
+
+public:
+
+ template<typename particles_type>
+ Derivative_xy(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor)
+ {
+ Point<particles_type::dims,unsigned int> p;
+ p.zero();
+ p.get(0) = 1;
+ p.get(1) = 1;
+
+ dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,dcpse_scaling_factor);
+ }
+
+ template<typename operand_type>
+
+ vector_dist_expression_op<operand_type,Dcpse<operand_type::vtype::dims,typename operand_type::vtype>,VECT_DCPSE> operator()(operand_type arg)
+ {
+ typedef Dcpse<operand_type::vtype::dims,typename operand_type::vtype> dcpse_type;
+
+ return vector_dist_expression_op<operand_type,dcpse_type,VECT_DCPSE>(arg,*(dcpse_type *)dcpse);
+ }
+};
template<typename operand_type1, typename operand_type2>
class plus
diff --git a/src/DCPSE_op/DCPSE_op_test.cpp b/src/DCPSE_op/DCPSE_op_test.cpp
index 67e73496..ce66e23c 100644
--- a/src/DCPSE_op/DCPSE_op_test.cpp
+++ b/src/DCPSE_op/DCPSE_op_test.cpp
@@ -3,7 +3,7 @@
*
* Created on: Jan 7, 2020
* Author: Abhinav Singh, Pietro Incardona
- * //adding labfolder hook
+ *
*/
#include "config.h"
@@ -19,11 +19,229 @@
BOOST_AUTO_TEST_SUITE( dcpse_op_suite_tests )
+
+BOOST_AUTO_TEST_CASE( dcpse_Navier)
+ {
+// int rank;
+// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+ size_t edgeSemiSize = 40;
+ const size_t sz[2] = {2 * edgeSemiSize, 2 * edgeSemiSize};
+ Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
+ size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
+ double spacing[2];
+ spacing[0] = 0.1 / (sz[0] - 1);
+ spacing[1] = 0.1 / (sz[1] - 1);
+ Ghost<2, double> ghost(spacing[0]*3);
+ double rCut = 2.0*spacing[0];
+ BOOST_TEST_MESSAGE("Init vector_dist...");
+ constexpr int Velocity = 0;
+ constexpr int Pressure = 1;
+ constexpr int bflag = 2;
+ double nu=1/3200.0;
+ // v P bflag
+ vector_dist<2, double, aggregate<VectorS<2,double>,double,int>> domain(0, box, bc, ghost);
+ //Init_DCPSE(domain)
+ BOOST_TEST_MESSAGE("Init domain...");
+ auto it = domain.getGridIterator(sz);
+ size_t pointId = 0;
+ size_t counter = 0;
+ while (it.isNext()) {
+ domain.add();
+ auto key = it.get();
+ mem_id k0 = key.get(0);
+ double x = k0 * spacing[0];
+ domain.getLastPos()[0] = x;
+ mem_id k1 = key.get(1);
+ double y = k1 * spacing[1];
+ domain.getLastPos()[1] = y;
+ // Initial Velocity field and Pressure
+ if (x == 0 || y==0 || x==0.1)
+ {
+ domain.template getLastProp<bflag>() = 1;
+ domain.template getLastProp<Velocity>()[0] = 0.0;
+ domain.template getLastProp<Velocity>()[1] = 0.0;
+ domain.template getLastProp<Pressure>() = 0.01;
+ if(x==0 && y==0)
+ {
+ domain.template getLastProp<bflag>() = 3;
+ domain.template getLastProp<Velocity>()[0] = 0.0;
+ domain.template getLastProp<Velocity>()[1] = 0.0;
+ domain.template getLastProp<Pressure>() = 0.0;
+ }
+ }
+ else{
+ domain.template getLastProp<bflag>() = 0;
+ domain.template getLastProp<Velocity>()[0] = 0.0;
+ domain.template getLastProp<Velocity>()[1] = 0.0;
+ domain.template getLastProp<Pressure>() = 0.01;
+ }
+ if(y==0.1)
+ {
+ domain.template getLastProp<bflag>() = 2;
+ domain.template getLastProp<Velocity>()[0] = 1.0;
+ domain.template getLastProp<Velocity>()[1] = 0.0;
+ domain.template getLastProp<Pressure>() = 0.01;
+ }
+ ++it;
+ }
+ BOOST_TEST_MESSAGE("Sync domain across processors...");
+
+ domain.map();
+ domain.ghost_get<0>();
+
+ Derivative_x Dx(domain,2,rCut);
+ Derivative_y Dy(domain,2,rCut);
+ Gradient Grad(domain,2,rCut);
+ Laplacian Lap(domain,2,rCut);
+ auto v = getV<0>(domain);
+ auto P = getV<1>(domain);
+ std::cout<<"RHS Initiated"<<std::endl;
+
+ v = -Grad(P);//-Grad(v)+Lap(v);
+
+ std::cout<<"RHS Computed"<<std::endl;
+ auto it2 = domain.getDomainIterator();
+ while (it2.isNext())
+ { auto p = it2.get();
+ if(domain.getProp<bflag>(p)==2)
+ {
+ domain.getProp<Velocity>(p)[0] = 1.0;
+ domain.getProp<Velocity>(p)[1] = 0.0;
+ }
+ ++it2;
+ }
+ domain.deleteGhost();
+ domain.write("Nav");
+
+ //std::cout << demangle(typeid(decltype(v)).name()) << "\n";
+
+ //Debug<decltype(expr)> a;
+
+ //typedef decltype(expr)::blabla blabla;
+
+ //auto err = Dx + Dx;
+ }
+
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+ BOOST_AUTO_TEST_CASE(dcpse_op_vec)
+ {
+// int rank;
+// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+ size_t edgeSemiSize = 40;
+ const size_t sz[2] = {2 * edgeSemiSize, 2 * edgeSemiSize};
+ Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
+ size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
+ double spacing[2];
+ spacing[0] = 2*M_PI / (sz[0] - 1);
+ spacing[1] = 2*M_PI / (sz[1] - 1);
+ Ghost<2, double> ghost(spacing[0]*3);
+ double rCut = 2.0*spacing[0];
+ BOOST_TEST_MESSAGE("Init vector_dist...");
+ double sigma2 = spacing[0] * spacing[1] / (2 * 4);
+
+ vector_dist<2, double, aggregate<double, VectorS<2,double>, VectorS<2,double>, VectorS<2,double>,VectorS<2,double>>> domain(0, box, bc, ghost);
+
+ //Init_DCPSE(domain)
+ BOOST_TEST_MESSAGE("Init domain...");
+// std::random_device rd{};
+// std::mt19937 rng{rd()};
+ std::mt19937 rng{6666666};
+
+ std::normal_distribution<> gaussian{0, sigma2};
+
+ auto it = domain.getGridIterator(sz);
+ size_t pointId = 0;
+ size_t counter = 0;
+ double minNormOne = 999;
+ while (it.isNext())
+ {
+ domain.add();
+ auto key = it.get();
+ mem_id k0 = key.get(0);
+ double x = k0 * spacing[0];
+ domain.getLastPos()[0] = x ;//+ gaussian(rng);
+ mem_id k1 = key.get(1);
+ double y = k1 * spacing[1];
+ domain.getLastPos()[1] = y ;//+gaussian(rng);
+ // Here fill the function value
+ domain.template getLastProp<1>()[0] = sin(domain.getLastPos()[0])+sin(domain.getLastPos()[1]);
+ domain.template getLastProp<1>()[1] = cos(domain.getLastPos()[0])+cos(domain.getLastPos()[1]);
+// domain.template getLastProp<0>() = x * x;
+// domain.template getLastProp<0>() = x;
+ // Here fill the validation value for Df/Dx
+ domain.template getLastProp<2>()[0] = 0;//cos(domain.getLastPos()[0]);//+cos(domain.getLastPos()[1]);
+ domain.template getLastProp<2>()[1] = 0;//-sin(domain.getLastPos()[0]);//+cos(domain.getLastPos()[1]);
+ domain.template getLastProp<4>()[0] = cos(domain.getLastPos()[0])*(sin(domain.getLastPos()[0])+sin(domain.getLastPos()[1]))+cos(domain.getLastPos()[1])*(cos(domain.getLastPos()[0])+cos(domain.getLastPos()[1]));
+ domain.template getLastProp<4>()[1] = -sin(domain.getLastPos()[0])*(sin(domain.getLastPos()[0])+sin(domain.getLastPos()[1]))-sin(domain.getLastPos()[1])*(cos(domain.getLastPos()[0])+cos(domain.getLastPos()[1]));
+
+
+// domain.template getLastProp<2>() = 2 * x;
+// domain.template getLastProp<2>() = 1;
+
+ ++counter;
+ ++it;
+ }
+ BOOST_TEST_MESSAGE("Sync domain across processors...");
+
+ domain.map();
+ domain.ghost_get<0>();
+
+ Derivative_x Dx(domain,2,rCut);
+ Derivative_y Dy(domain,2,rCut);
+ Gradient Grad(domain,2,rCut);
+ Laplacian Lap(domain,2,rCut,3);
+ Advection Adv(domain,2,rCut,3);
+ auto v = getV<1>(domain);
+ auto P = getV<0>(domain);
+ auto dv = getV<3>(domain);
+
+// typedef boost::mpl::int_<std::is_fundamental<point_expression_op<Point<2U, double>, point_expression<double>, Point<2U, double>, 3>>::value>::blabla blabla;
+
+// std::is_fundamental<decltype(o1.value(key))>
+
+ //vv=Lap(P);
+ //dv=Lap(v);//+Dy(P);
+ dv=Adv(v,v);//+Dy(P);
+ auto it2 = domain.getDomainIterator();
+
+ double worst1 = 0.0;
+
+ while (it2.isNext())
+ {
+ auto p = it2.get();
+
+ std::cout << "VALS: " << domain.getProp<3>(p)[0] << " " << domain.getProp<4>(p)[0] << std::endl;
+ std::cout << "VALS: " << domain.getProp<3>(p)[1] << " " << domain.getProp<4>(p)[1] << std::endl;
+
+ if (fabs(domain.getProp<3>(p)[0] - domain.getProp<4>(p)[0]) > worst1)
+ {
+ worst1 = fabs(domain.getProp<3>(p)[0] - domain.getProp<4>(p)[0]);
+ }
+
+ ++it2;
+ }
+
+ std::cout << "Maximum Error: " << worst1 << std::endl;
+
+ domain.deleteGhost();
+ domain.write("v");
+
+ //std::cout << demangle(typeid(decltype(v)).name()) << "\n";
+
+ //Debug<decltype(expr)> a;
+
+ //typedef decltype(expr)::blabla blabla;
+
+ //auto err = Dx + Dx;
+ }
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
BOOST_AUTO_TEST_CASE( dcpse_op_tests )
{
// int rank;
// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
- size_t edgeSemiSize = 160;
+ size_t edgeSemiSize = 40;
const size_t sz[2] = {2 * edgeSemiSize, 2 * edgeSemiSize};
Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
@@ -35,7 +253,7 @@ BOOST_AUTO_TEST_CASE( dcpse_op_tests )
BOOST_TEST_MESSAGE("Init vector_dist...");
double sigma2 = spacing[0] * spacing[1] / (2 * 4);
- vector_dist<2, double, aggregate<double, double, double>> domain(0, box, bc, ghost);
+ vector_dist<2, double, aggregate<double, double, double, VectorS<2,double>,VectorS<2,double>>> domain(0, box, bc, ghost);
//Init_DCPSE(domain)
BOOST_TEST_MESSAGE("Init domain...");
@@ -55,16 +273,20 @@ BOOST_AUTO_TEST_CASE( dcpse_op_tests )
auto key = it.get();
mem_id k0 = key.get(0);
double x = k0 * spacing[0];
- domain.getLastPos()[0] = x + gaussian(rng);
+ domain.getLastPos()[0] = x ;//+ gaussian(rng);
mem_id k1 = key.get(1);
double y = k1 * spacing[1];
- domain.getLastPos()[1] = y + gaussian(rng);
+ domain.getLastPos()[1] = y ;//+gaussian(rng);
// Here fill the function value
- domain.template getLastProp<0>() = sin(domain.getLastPos()[0]);
+ domain.template getLastProp<0>() = sin(domain.getLastPos()[0])+sin(domain.getLastPos()[1]);
// domain.template getLastProp<0>() = x * x;
// domain.template getLastProp<0>() = x;
// Here fill the validation value for Df/Dx
- domain.template getLastProp<2>() = cos(domain.getLastPos()[0]);
+ domain.template getLastProp<2>() = cos(domain.getLastPos()[0])+cos(domain.getLastPos()[1]);
+ domain.template getLastProp<4>()[0] = -sin(domain.getLastPos()[0]);
+ domain.template getLastProp<4>()[1] = -sin(domain.getLastPos()[1]);
+
+
// domain.template getLastProp<2>() = 2 * x;
// domain.template getLastProp<2>() = 1;
@@ -77,13 +299,15 @@ BOOST_AUTO_TEST_CASE( dcpse_op_tests )
domain.ghost_get<0>();
Derivative_x Dx(domain,2,rCut);
- //Dy Dy;
+ Derivative_y Dy(domain,2,rCut);
+ Gradient Grad(domain,2,rCut);
+ Laplacian Lap(domain,2,rCut);
auto v = getV<1>(domain);
auto P = getV<0>(domain);
+ auto vv = getV<3>(domain);
-
- v = Dx(P);//+ Dy(F);
-
+ vv=Lap(P);
+ v=Dx(P)+Dy(P);
auto it2 = domain.getDomainIterator();
double worst = 0.0;
@@ -100,12 +324,12 @@ BOOST_AUTO_TEST_CASE( dcpse_op_tests )
++it2;
}
- std::cout << "WORST " << worst << std::endl;
+ std::cout << "Maximum Error: " << worst << std::endl;
domain.deleteGhost();
- domain.write("Dx");
+ domain.write("v");
- std::cout << demangle(typeid(decltype(v)).name()) << "\n";
+ //std::cout << demangle(typeid(decltype(v)).name()) << "\n";
//Debug<decltype(expr)> a;
diff --git a/src/Operators/Vector/vector_dist_operators.hpp b/src/Operators/Vector/vector_dist_operators.hpp
index 0e49da99..6a06ea2a 100644
--- a/src/Operators/Vector/vector_dist_operators.hpp
+++ b/src/Operators/Vector/vector_dist_operators.hpp
@@ -63,6 +63,9 @@
#define VECT_NEARBYINT 89
#define VECT_RINT 90
#define VECT_DCPSE 100
+#define VECT_DCPSE_V 101
+#define VECT_DCPSE_V_SUM 102
+#define VECT_DCPSE_V_DOT 103
#define VECT_PMUL 91
#define VECT_SUB_UNI 92
--
GitLab