diff --git a/src/DCPSE_op/DCPSE_Solver.hpp b/src/DCPSE_op/DCPSE_Solver.hpp
index 7df0f0810700bd8c94650f1375997eda6e25a290..64be04933d08198323908d32313f4b5366d2733b 100644
--- a/src/DCPSE_op/DCPSE_Solver.hpp
+++ b/src/DCPSE_op/DCPSE_Solver.hpp
@@ -16,6 +16,12 @@
template<unsigned int prp_id>
struct prop_id {};
+enum options_solver
+{
+ STANDARD,
+ LAGRANGE_MULTIPLIER
+};
+
//template<unsigned int prp_id> using prop_id = boost::mpl::int_<prp_id>;
template<typename Sys_eqs, typename particles_type>
@@ -58,11 +64,14 @@ class DCPSE_scheme: public MatrixAssembler
//! Total number of points
size_t tot;
+ //! solver options
+ options_solver opt;
+
/*! \brief Construct the gmap structure
*
*/
- void construct_pmap()
+ void construct_pmap(options_solver opt = options_solver::STANDARD)
{
Vcluster<> & v_cl = create_vcluster();
@@ -83,7 +92,11 @@ class DCPSE_scheme: public MatrixAssembler
v_cl.execute();
// resize b if needed
- b.resize(Sys_eqs::nvar * tot,Sys_eqs::nvar * sz);
+ if (opt == options_solver::STANDARD) {
+ b.resize(Sys_eqs::nvar * tot, Sys_eqs::nvar * sz);
+ } else{
+ b.resize(Sys_eqs::nvar * tot + 1, Sys_eqs::nvar * sz + 1);
+ }
// Calculate the starting point
@@ -259,7 +272,7 @@ public:
void solve(expr_type exp)
{
umfpack_solver<double> solver;
- auto x = solver.solve(getA(),getB());
+ auto x = solver.solve(getA(opt),getB(opt));
//petsc_solver<double> solver;
//solver.setSolver(KSPBCGS);
//solver.setMaxIter(1000);
@@ -272,9 +285,7 @@ public:
while (it.isNext())
{
auto p = it.get();
-
exp.value(p) = x(p.getKey());
-
++it;
}
}
@@ -289,12 +300,12 @@ public:
* \param b_g object grid that will store the solution
*
*/
- DCPSE_scheme(const typename Sys_eqs::stype r_cut, particles_type & part)
- :parts(part),p_map(0,part.getDecomposition().getDomain(),part.getDecomposition().periodicity(),Ghost<particles_type::dims,typename particles_type::stype>(r_cut)),row(0),row_b(0)
+ DCPSE_scheme(const typename Sys_eqs::stype r_cut, particles_type & part, options_solver opt = options_solver::STANDARD)
+ :parts(part),p_map(0,part.getDecomposition().getDomain(),part.getDecomposition().periodicity(),Ghost<particles_type::dims,typename particles_type::stype>(r_cut)),row(0),row_b(0),opt(opt)
{
p_map.resize(part.size_local());
- construct_pmap();
+ construct_pmap(opt);
}
@@ -357,15 +368,55 @@ public:
* \return the Sparse matrix produced
*
*/
- typename Sys_eqs::SparseMatrix_type & getA()
+ typename Sys_eqs::SparseMatrix_type & getA(options_solver opt = options_solver::STANDARD)
{
#ifdef SE_CLASS1
consistency();
#endif
- A.resize(tot*Sys_eqs::nvar,tot*Sys_eqs::nvar,
- p_map.size_local()*Sys_eqs::nvar,
- p_map.size_local()*Sys_eqs::nvar);
+ if (opt == options_solver::STANDARD) {
+ A.resize(tot * Sys_eqs::nvar, tot * Sys_eqs::nvar,
+ p_map.size_local() * Sys_eqs::nvar,
+ p_map.size_local() * Sys_eqs::nvar);
// std::cout << Eigen::MatrixXd(A) << std::endl;
+ } else
+ {
+ auto & v_cl = create_vcluster();
+ openfpm::vector<triplet> & trpl = A.getMatrixTriplets();
+
+ if (v_cl.rank() == v_cl.size() - 1)
+ {
+ A.resize(tot * Sys_eqs::nvar + 1, tot * Sys_eqs::nvar + 1,
+ p_map.size_local() * Sys_eqs::nvar + 1,
+ p_map.size_local() * Sys_eqs::nvar + 1);
+
+ for (int i = 0 ; i < tot * Sys_eqs::nvar ; i++)
+ {
+ triplet t1;
+ triplet t2;
+
+ t1.row() = tot * Sys_eqs::nvar;
+ t1.col() = i;
+ t1.value() = 1;
+
+ t2.row() = i;
+ t2.col() = tot * Sys_eqs::nvar;
+ t2.value() = 1;
+
+ trpl.add(t1);
+ trpl.add(t2);
+ }
+
+ row_b++;
+ row++;
+ }
+ else {
+ A.resize(tot * Sys_eqs::nvar + 1, tot * Sys_eqs::nvar + 1,
+ p_map.size_local() * Sys_eqs::nvar,
+ p_map.size_local() * Sys_eqs::nvar);
+ }
+
+
+ }
return A;
}
@@ -375,11 +426,20 @@ public:
* \return the vector produced
*
*/
- typename Sys_eqs::Vector_type & getB()
+ typename Sys_eqs::Vector_type & getB(options_solver opt = options_solver::STANDARD)
{
#ifdef SE_CLASS1
consistency();
#endif
+ if (opt == options_solver::LAGRANGE_MULTIPLIER)
+ {
+ auto & v_cl = create_vcluster();
+ if (v_cl.rank() == v_cl.size() - 1) {
+
+ b(tot * Sys_eqs::nvar) = 0;
+ }
+ }
+
return b;
}
diff --git a/src/DCPSE_op/DCPSE_copy/Dcpse.hpp b/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
index d24f97a21dcee46b9916b82963d0d1bfd430846a..74185aee1e19246926a888f5b944a617d6fc5bbb 100644
--- a/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
+++ b/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
@@ -86,7 +86,7 @@ public:
EMatrix<T, Eigen::Dynamic, 1> &a = localCoefficients[k];
Support<dim, T, part_type> support = localSupports[k];
auto eps = localEps[k];
-
+ int DrawKernelKounter=0;
size_t xpK = k;
Point<dim, T> xp = support.getReferencePoint();
for (auto &xqK : support.getKeys())
@@ -95,7 +95,9 @@ public:
Point<dim, T> normalizedArg = (xp - xq) / eps;
particles.template getProp<prp>(xqK) += computeKernel(normalizedArg, a);
+ DrawKernelKounter++;
}
+ std::cout<<"Number of Neighbours: "<<DrawKernelKounter<<std::endl;
}
template<unsigned int prp>
@@ -512,9 +514,6 @@ private:
res += coeff * mbValue * expFactor;
++counter;
}
-// if(fabs(res)<1e-8) {
-//res = 0;
-// }
return res;
}
diff --git a/src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp b/src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp
index 73309fe13de39478a6f1691b680a00400feffad3..f465a58bbece263d5646844ab7802179abcf5f89 100644
--- a/src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp
+++ b/src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp
@@ -93,7 +93,7 @@ size_t SupportBuilder<vector_type>::getNumElementsInSetOfCells(const std::set<gr
template<typename vector_type>
void SupportBuilder<vector_type>::enlargeSetOfCellsUntilSize(std::set<grid_key_dx<vector_type::dims>> &set, unsigned int requiredSize)
{
- while (getNumElementsInSetOfCells(set) < requiredSize)
+ while (getNumElementsInSetOfCells(set) < 5*requiredSize)
{
auto tmpSet = set;
for (const auto el : tmpSet)
diff --git a/src/DCPSE_op/DCPSE_op.hpp b/src/DCPSE_op/DCPSE_op.hpp
index 4c44bcd01513024158e37c1fee08a54459b3497e..f72b6dea9a31886d528ab19d2c543ccf676a7bc0 100644
--- a/src/DCPSE_op/DCPSE_op.hpp
+++ b/src/DCPSE_op/DCPSE_op.hpp
@@ -284,6 +284,69 @@ public:
};
+template <>
+class vector_dist_expression_op<void,void,VECT_COPY_N_TO_N>
+{
+ mutable int i;
+
+ //! expression 1
+ openfpm::vector<aggregate<int>> & l1;
+ openfpm::vector<aggregate<int>> & l2;
+
+public:
+
+
+ inline vector_dist_expression_op(openfpm::vector<aggregate<int>> & l1,openfpm::vector<aggregate<int>> & l2)
+ :l1(l1),l2(l2)
+ {}
+
+ template<typename pmap_type, typename unordered_map_type, typename coeff_type>
+ inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff) const
+ {
+ if (l1.template get<0>(i) != key.getKey())
+ {
+ std::cout << "ERROR" << std::endl;
+ }
+
+ cols[p_map. template getProp<0>(key)] += coeff;
+ std::cout << "L2: " << l2.template get<0>(i) << std::endl;
+ cols[p_map. template getProp<0>(l2.template get<0>(i))] -= coeff;
+
+ i++;
+ }
+};
+
+
+
+template <>
+class vector_dist_expression_op<void,void,VECT_COPY_1_TO_N>
+{
+ mutable int i = 0;
+
+ //! expression 1
+ openfpm::vector<aggregate<int>> & l1;
+ int l2_key;
+
+public:
+
+ inline vector_dist_expression_op(openfpm::vector<aggregate<int>> & l1, int l2_key)
+ :l1(l1),l2_key(l2_key)
+ {}
+
+ template<typename pmap_type, typename unordered_map_type, typename coeff_type>
+ inline void value_nz(pmap_type & p_map, const vect_dist_key_dx & key, unordered_map_type & cols, coeff_type & coeff) const
+ {
+ if (l1.template get<0>(i) != key.getKey())
+ {
+ std::cout << "ERROR" << std::endl;
+ }
+
+ cols[p_map. template getProp<0>(key)] += coeff;
+ cols[p_map. template getProp<0>(l2_key)] -= coeff;
+ i++;
+ }
+};
+
template <typename exp1,typename DCPSE_type>
@@ -747,6 +810,31 @@ public:
return vector_dist_expression_op<operand_type,dcpse_type,VECT_DCPSE_V_SUM>(arg,*(dcpse_type(*)[operand_type::vtype::dims])dcpse);
}
+
+
+ template<typename particles_type>
+ void checkMomenta(particles_type &particles)
+ {
+ Dcpse<particles_type::dims,particles_type> * dcpse_ptr = (Dcpse<particles_type::dims,particles_type> *)dcpse;
+
+ for (int i = 0 ; i < particles_type::dims ; i++)
+ {
+ dcpse_ptr[i].checkMomenta(particles);
+ }
+
+ }
+ template<unsigned int prp, typename particles_type>
+ void DrawKernel(particles_type &particles,int k)
+ {
+ Dcpse<particles_type::dims,particles_type> * dcpse_ptr = (Dcpse<particles_type::dims,particles_type> *)dcpse;
+
+ for (int i = 0 ; i < particles_type::dims ; i++)
+ {
+ dcpse_ptr[i].template DrawKernel<prp>(particles,k);
+ }
+
+ }
+
};
diff --git a/src/DCPSE_op/DCPSE_op_test.cpp b/src/DCPSE_op/DCPSE_op_test.cpp
index c84d6eaf29b4a925c9d69e84e00363540bbf8f65..4b7aa18f8cced56e533354e7a6dd07afeb1ce9c9 100644
--- a/src/DCPSE_op/DCPSE_op_test.cpp
+++ b/src/DCPSE_op/DCPSE_op_test.cpp
@@ -538,12 +538,13 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
BOOST_AUTO_TEST_CASE(dcpse_Lid_sf) {
- const size_t sz[2] = {81,81};
+ const size_t sz[2] = {8,8};
Box<2, double> box({0, 0}, {1,1});
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing = box.getHigh(0) / (sz[0] - 1);
Ghost<2, double> ghost(spacing * 3);
- double rCut = 2.0 * spacing;
+ double rCut = spacing*0.8 + spacing*1e-02;
+ std::cout<<spacing<<std::endl;
// sf W DW RHS Wnew V
vector_dist<2, double, aggregate<double,double,double,double,double,VectorS<2, double>>> Particles(0, box, bc, ghost);
@@ -567,93 +568,272 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
openfpm::vector<aggregate<int>> l_p;
openfpm::vector<aggregate<int>> r_p;
- auto Sf = getV<0>(Particles);
- auto W = getV<1>(Particles);
- auto dW = getV<2>(Particles);
- auto RHS = getV<3>(Particles);
- auto Wnew = getV<4>(Particles);
- auto V = getV<5>(Particles);
+ openfpm::vector<aggregate<int>> up_p1;
+ openfpm::vector<aggregate<int>> dw_p1;
+ openfpm::vector<aggregate<int>> l_p1;
+ openfpm::vector<aggregate<int>> r_p1;
+
// Here fill up the boxes for particle detection.
Box<2, double> up({box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0},
- {box.getHigh(0) - spacing / 2.0, box.getHigh(1) + spacing / 2.0});
+ {box.getHigh(0) - spacing / 2.0, box.getHigh(1) + spacing / 2.0});
- Box<2, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
- {box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0});
+ Box<2, double> up_d({box.getLow(0) + spacing / 2.0, box.getHigh(1) - 3 * spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0});
- Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
- {box.getLow(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
+ Box<2, double> down({box.getLow(0) + spacing / 2.0, box.getLow(1) - spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0});
+
+ Box<2, double> down_u({box.getLow(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getLow(1) + 3 * spacing / 2.0});
+
+ Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
+ {box.getLow(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
+
+ Box<2, double> left_r({box.getLow(0) + spacing / 2.0, box.getLow(1) + 3 * spacing / 2.0},
+ {box.getLow(0) + 3 * spacing / 2.0, box.getHigh(1) - 3 * spacing / 2.0});
+
+ Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
+
+ Box<2, double> right_l({box.getHigh(0) - 3 * spacing / 2.0, box.getLow(1) + 3 * spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getHigh(1) - 3 * spacing / 2.0});
+ Box<2, double> Bulk_box({box.getLow(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getHigh(1) -spacing / 2.0});
- Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
- {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
openfpm::vector<Box<2, double>> boxes;
boxes.add(up);
+ boxes.add(up_d);
boxes.add(down);
+ boxes.add(down_u);
boxes.add(left);
+ boxes.add(left_r);
boxes.add(right);
+ boxes.add(right_l);
+ boxes.add(Bulk_box);
VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
vtk_box.add(boxes);
- vtk_box.write("Re1000-1e-3-vtk_box.vtk");
+ vtk_box.write("box_sf.vtk");
+ // Particles.write_frame("Re1000-1e-3-Lid_sf",0);
auto it2 = Particles.getDomainIterator();
while (it2.isNext()) {
auto p = it2.get();
Point<2, double> xp = Particles.getPos(p);
+ Particles.getProp<0>(p) =0;//-M_PI*M_PI*sin(M_PI*xp[0])*sin(M_PI*xp[1]);
+ //Particles.getProp<1>(p) =0;//sin(M_PI*xp[0])*sin(M_PI*xp[1]);
+ Particles.getProp<2>(p) =0.0;
+ Particles.getProp<3>(p) =0.0;
+ Particles.getProp<4>(p) =0.0;
+ Particles.getProp<5>(p)[0] =0.0;
+ Particles.getProp<5>(p)[1] =0.0;
+
if (up.isInside(xp) == true) {
up_p.add();
up_p.last().get<0>() = p.getKey();
- Particles.getProp<1>(p) = -2.0/sz[0];
+ Particles. template getProp<1>(p) = 12.0/spacing;// -2.0*Particles.getProp<0>(p)/(spacing*spacing) - 12.0/spacing;
}
else if (down.isInside(xp) == true) {
dw_p.add();
dw_p.last().get<0>() = p.getKey();
- }else if (left.isInside(xp) == true) {
+ }
+ else if (left.isInside(xp) == true) {
l_p.add();
l_p.last().get<0>() = p.getKey();
} else if (right.isInside(xp) == true) {
r_p.add();
r_p.last().get<0>() = p.getKey();
}
- else {
+ else if (Bulk_box.isInside(xp) == true)
+ {
+ if (up_d.isInside(xp) == true) {
+ up_p1.add();
+ up_p1.last().get<0>() = p.getKey();
+ }
+ else if (down_u.isInside(xp) == true) {
+ dw_p1.add();
+ dw_p1.last().get<0>() = p.getKey();
+ }else if (left_r.isInside(xp) == true) {
+ l_p1.add();
+ l_p1.last().get<0>() = p.getKey();
+ } else if (right_l.isInside(xp) == true) {
+ r_p1.add();
+ r_p1.last().get<0>() = p.getKey();
+ }
bulk.add();
bulk.last().get<0>() = p.getKey();
}
++it2;
}
+/* //test for copy
+ for(int j=0;j<up_p1.size();j++)
+ {
+ auto p1=up_p1.get<0>(j);
+ Point<2, double> xp = Particles.getPos(p1);
+ Particles.getProp<3>(p1) = xp[0];
+ }
+ for(int j=0;j<l_p1.size();j++)
+ {
+ auto p1=l_p1.get<0>(j);
+ Point<2, double> xp = Particles.getPos(p1);
+ Particles.getProp<3>(p1) = xp[1];
+ }
+ for(int j=0;j<r_p1.size();j++)
+ {
+ auto p1=r_p1.get<0>(j);
+ Point<2, double> xp = Particles.getPos(p1);
+ Particles.getProp<3>(p1) = xp[1];
+ }
+ for(int j=0;j<dw_p1.size();j++)
+ {
+ auto p1=dw_p1.get<0>(j);
+ Point<2, double> xp = Particles.getPos(p1);
+ Particles.getProp<3>(p1)= xp[0];
+ }
+ //Copy Module
+ for(int j=0;j<up_p1.size();j++)
+ { auto p=up_p.get<0>(j);
+ auto p1=up_p1.get<0>(j);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ if(j==0)
+ { auto p=l_p.get<0>(l_p.size()-1);
+ auto p2=l_p.get<0>(l_p.size()-2);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2) = Particles.getProp<3>(p1);
+ }
+ if(j==up_p1.size()-1)
+ { auto p=r_p.get<0>(r_p.size()-1);
+ auto p2=r_p.get<0>(r_p.size()-2);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2) = Particles.getProp<3>(p1);
+ }
+ }
+
+ for(int j=0;j<l_p1.size();j++)
+ {
+ auto p=l_p.get<0>(j+2);
+ auto p1=l_p1.get<0>(j);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ }
+ for(int j=0;j<r_p1.size();j++)
+ {
+ auto p=r_p.get<0>(j+2);
+ auto p1=r_p1.get<0>(j);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ }
+ for(int j=0;j<dw_p1.size();j++)
+ { auto p=dw_p.get<0>(j);
+ auto p1=dw_p1.get<0>(j);
+ Particles.getProp<3>(p)= Particles.getProp<3>(p1);
+ if(j==0)
+ { auto p=l_p.get<0>(0);
+ auto p2=l_p.get<0>(1);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2) = Particles.getProp<3>(p1);
+ }
+ if(j==dw_p1.size()-1)
+ { auto p=r_p.get<0>(0);
+ auto p2=r_p.get<0>(1);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2) = Particles.getProp<3>(p1);
+ }
+ }*/
+
Particles.write_frame("Re1000-1e-3-Lid_sf",0);
+ auto Sf = getV<0>(Particles);
+ auto W = getV<1>(Particles);
+ auto dW = getV<2>(Particles);
+ auto RHS = getV<3>(Particles);
+ auto Wnew = getV<4>(Particles);
+ auto V = getV<5>(Particles);
+
+ for(int j=0;j<up_p.size();j++) {
+ auto p = up_p.get<0>(j);
+ Particles.getProp<1>(p) = -12.0/spacing;
+ }
+
+ Particles.write_frame("Re1000-1e-3-Lid_sf",1);
+ Derivative_x Dx(Particles, 2, rCut,1);
+ Derivative_y Dy(Particles, 2, rCut,1);
+ Gradient Grad(Particles, 2, rCut,1);
+ Laplacian Lap(Particles, 2, rCut, 1);
+ Curl2D Curl(Particles, 2, rCut, 1);
+ return;
- Derivative_x Dx(Particles, 1, rCut,3);
- Derivative_y Dy(Particles, 1, rCut,3);
- Gradient Grad(Particles, 1, rCut,3);
- Laplacian Lap(Particles, 1, rCut, 3);
- Curl2D Curl(Particles, 1, rCut, 3);
- DCPSE_scheme<equations,decltype(Particles)> Solver(2*rCut, Particles);
- auto Sf_Poisson = Lap(Sf);
+ /*DCPSE_scheme<equations,decltype(Particles)> Solver(2*rCut, Particles);
+ auto Sf_Poisson = -Lap(Sf);
Solver.impose(Sf_Poisson, bulk, prop_id<1>());
Solver.impose(Sf, up_p, 0);
Solver.impose(Sf, dw_p,0);
Solver.impose(Sf, l_p,0);
Solver.impose(Sf, r_p, 0);
Solver.solve(Sf);
- V=Curl(Sf);
-
- Particles.write("Initialization");
-
-
+ RHS=-Lap(Sf);*/
+ Particles.write_frame("Re1000-1e-3-Lid_sf",2);
+ return;
double dt=0.003;
double nu=0.01;
int n=150;
std::cout<<"Init Done"<<std::endl;
for(int i =1; i<=n ;i++)
- { dW=nu*Lap(W)+Dx(Sf)*Dy(W)-Dy(Sf)*Dx(W);
+ { dW=Dx(Sf)*Dy(W)-Dy(Sf)*Dx(W)+nu*Lap(W);
+ //Lap.DrawKernel<3>(Particles,837);
Wnew=W+dt*dW;
W=Wnew;
+ //Copy Module
+ for(int j=0;j<up_p1.size();j++)
+ { auto p=up_p.get<0>(j);
+ auto p1=up_p1.get<0>(j);
+ Particles.getProp<1>(p) = -2.0*Particles.getProp<0>(p1)/(spacing*spacing)-12.0/spacing;
+ if(j==0)
+ { auto p=l_p.get<0>(l_p.size()-1);
+ auto p2=l_p.get<0>(l_p.size()-2);
+ //Particles.getProp<1>(p) = 0;//-2.0*Particles.getProp<0>(p1)/(spacing*spacing);
+ Particles.getProp<1>(p2) = 0;//-2.0*Particles.getProp<0>(p1)/(spacing*spacing);
+ }
+ if(j==up_p1.size()-1)
+ { auto p=r_p.get<0>(r_p.size()-1);
+ auto p2=r_p.get<0>(r_p.size()-2);
+ //Particles.getProp<1>(p) = 0;//-2.0*Particles.getProp<0>(p1)/(spacing*spacing);
+ Particles.getProp<1>(p2) = 0;//-2.0*Particles.getProp<0>(p1)/(spacing*spacing);
+ }
+ }
+ for(int j=0;j<l_p1.size();j++)
+ {
+ auto p=l_p.get<0>(j+2);
+ auto p1=l_p1.get<0>(j);
+ Particles.getProp<1>(p) = -2.0*Particles.getProp<0>(p1)/(spacing*spacing);
+ }
+ for(int j=0;j<r_p1.size();j++)
+ {
+ auto p=r_p.get<0>(j+2);
+ auto p1=r_p1.get<0>(j);
+ Particles.getProp<1>(p) = -2.0*Particles.getProp<0>(p1)/(spacing*spacing);
+ }
+ for(int j=0;j<dw_p1.size();j++)
+ { auto p=dw_p.get<0>(j);
+ auto p1=dw_p1.get<0>(j);
+ Particles.getProp<1>(p) = -2.0*Particles.getProp<0>(p1)/(spacing*spacing);
+ if(j==0)
+ { auto p=l_p.get<0>(0);
+ auto p2=l_p.get<0>(1);
+ //Particles.getProp<1>(p) = 0;//-2.0*Particles.getProp<0>(p1)/(spacing*spacing);
+ Particles.getProp<1>(p2) = 0;//-2.0*Particles.getProp<0>(p1)/(spacing*spacing);
+ }
+ if(j==dw_p1.size()-1)
+ { auto p=r_p.get<0>(0);
+ auto p2=r_p.get<0>(1);
+ // Particles.getProp<1>(p) = 0;//-2.0*Particles.getProp<0>(p1)/(spacing*spacing);
+ Particles.getProp<1>(p2) = 0;//-2.0*Particles.getProp<0>(p1)/(spacing*spacing);
+ }
+ }
std::cout<<"W Done"<<std::endl;
+
DCPSE_scheme<equations,decltype(Particles)> Solver(2*rCut, Particles);
- auto Sf_Poisson = Lap(Sf);
+ auto Sf_Poisson = -Lap(Sf);
Solver.impose(Sf_Poisson, bulk, prop_id<1>());
Solver.impose(Sf, up_p, 0);
Solver.impose(Sf, dw_p,0);
@@ -662,18 +842,458 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Solver.solve(Sf);
V=Curl(Sf);
std::cout<<"Poisson Solved"<<std::endl;
- //if (i%10==0)
Particles.write_frame("Re1000-1e-3-Lid_sf",i);
+ //if (i%10==0)
std::cout<<i<<std::endl;
}
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ BOOST_AUTO_TEST_CASE(dcpse_poisson_Robin_anal) {
+// int rank;
+// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+ const size_t sz[2] = {81,81};
+ Box<2, double> box({0, 0}, {0.5, 0.5});
+ size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
+ double spacing = box.getHigh(0) / (sz[0] - 1);
+ Ghost<2, double> ghost(spacing * 3);
+ double rCut = 2.0 * spacing;
+ BOOST_TEST_MESSAGE("Init vector_dist...");
+
+ vector_dist<2, double, aggregate<double,double,double,double,double,double>> domain(0, box, bc, ghost);
+
+
+ //Init_DCPSE(domain)
+ BOOST_TEST_MESSAGE("Init domain...");
+
+ auto it = domain.getGridIterator(sz);
+ while (it.isNext()) {
+ domain.add();
+
+ auto key = it.get();
+ double x = key.get(0) * it.getSpacing(0);
+ domain.getLastPos()[0] = x;
+ double y = key.get(1) * it.getSpacing(1);
+ domain.getLastPos()[1] = y;
+
+ ++it;
+ }
+ BOOST_TEST_MESSAGE("Sync domain across processors...");
+
+ domain.map();
+ domain.ghost_get<0>();
+
+ Derivative_x Dx(domain, 2, rCut,2);
+ Derivative_y Dy(domain, 2, rCut,2);
+ //Gradient Grad(domain, 2, rCut);
+ Laplacian Lap(domain, 2, rCut, 2);
+ //Advection Adv(domain, 3, rCut, 3);
+ //Solver Sol_Lap(Lap),Sol_Dx(Dx);
+
+
+ openfpm::vector<aggregate<int>> bulk;
+ openfpm::vector<aggregate<int>> up_p;
+ openfpm::vector<aggregate<int>> dw_p;
+ openfpm::vector<aggregate<int>> l_p;
+ openfpm::vector<aggregate<int>> r_p;
+ openfpm::vector<aggregate<int>> ref_p;
+
+ auto v = getV<0>(domain);
+ auto RHS=getV<1>(domain);
+ auto sol = getV<2>(domain);
+ auto anasol = getV<3>(domain);
+ auto err = getV<4>(domain);
+ auto DCPSE_sol=getV<5>(domain);
+
+ // Here fill me
+
+ Box<2, double> up({box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
+
+ Box<2, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0});
+
+ Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
+ Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
+ openfpm::vector<Box<2, double>> boxes;
+ boxes.add(up);
+ boxes.add(down);
+ boxes.add(left);
+ boxes.add(right);
+
+ // Create a writer and write
+ VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
+ vtk_box.add(boxes);
+ vtk_box.write("vtk_box.vtk");
+
+
+ auto it2 = domain.getDomainIterator();
+
+ while (it2.isNext()) {
+ auto p = it2.get();
+ Point<2, double> xp = domain.getPos(p);
+ //domain.getProp<3>(p)=1+xp[0]*xp[0]+2*xp[1]*xp[1];
+ if (up.isInside(xp) == true) {
+ up_p.add();
+ up_p.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ } else if (down.isInside(xp) == true) {
+ dw_p.add();
+ dw_p.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+
+ } else if (left.isInside(xp) == true) {
+ l_p.add();
+ l_p.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+
+ } else if (right.isInside(xp) == true) {
+ r_p.add();
+ r_p.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+
+ } else {
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ }
+ ++it2;
+ }
+ DCPSE_scheme<equations,decltype(domain)> Solver(2 * rCut, domain);
+ auto Poisson = Lap(v);
+ auto D_x = Dx(v);
+ auto D_y = Dy(v);
+ Solver.impose(Poisson, bulk, prop_id<1>());
+ Solver.impose(D_y, up_p, 0);
+ Solver.impose(D_x, r_p, 0);
+ Solver.impose(v, dw_p, 0);
+ Solver.impose(v, l_p, 0);
+ Solver.solve(sol);
+ DCPSE_sol=Lap(sol);
+ double worst1 = 0.0;
+
+ v=abs(DCPSE_sol-RHS);
+
+ for(int j=0;j<bulk.size();j++)
+ { auto p=bulk.get<0>(j);
+ if (fabs(domain.getProp<3>(p) - domain.getProp<2>(p)) >= worst1) {
+ worst1 = fabs(domain.getProp<3>(p) - domain.getProp<2>(p));
+ }
+ domain.getProp<4>(p) = fabs(domain.getProp<3>(p) - domain.getProp<2>(p));
+
+ }
+ std::cout << "Maximum Analytic Error: " << worst1 << std::endl;
+
+ domain.write("Robin_anasol");
+ }
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ BOOST_AUTO_TEST_CASE(dcpse_poisson_Dirichlet_anal) {
+// int rank;
+// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+ const size_t sz[2] = {81,81};
+ Box<2, double> box({0, 0}, {1, 1});
+ size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
+ double spacing = box.getHigh(0) / (sz[0] - 1);
+ Ghost<2, double> ghost(spacing * 3);
+ double rCut = 2.0 * spacing;
+ BOOST_TEST_MESSAGE("Init vector_dist...");
+
+ vector_dist<2, double, aggregate<double,double,double,double,double,double>> domain(0, box, bc, ghost);
+
+
+ //Init_DCPSE(domain)
+ BOOST_TEST_MESSAGE("Init domain...");
+
+ auto it = domain.getGridIterator(sz);
+ while (it.isNext()) {
+ domain.add();
+
+ auto key = it.get();
+ double x = key.get(0) * it.getSpacing(0);
+ domain.getLastPos()[0] = x;
+ double y = key.get(1) * it.getSpacing(1);
+ domain.getLastPos()[1] = y;
+
+ ++it;
+ }
+ BOOST_TEST_MESSAGE("Sync domain across processors...");
+
+ domain.map();
+ domain.ghost_get<0>();
+
+ Derivative_x Dx(domain, 2, rCut,2);
+ Derivative_y Dy(domain, 2, rCut,2);
+ //Gradient Grad(domain, 2, rCut);
+ Laplacian Lap(domain, 2, rCut, 2);
+ //Advection Adv(domain, 3, rCut, 3);
+ //Solver Sol_Lap(Lap),Sol_Dx(Dx);
+
+ openfpm::vector<aggregate<int>> bulk;
+ openfpm::vector<aggregate<int>> up_p;
+ openfpm::vector<aggregate<int>> dw_p;
+ openfpm::vector<aggregate<int>> l_p;
+ openfpm::vector<aggregate<int>> r_p;
+ openfpm::vector<aggregate<int>> ref_p;
+
+ auto v = getV<0>(domain);
+ auto RHS=getV<1>(domain);
+ auto sol = getV<2>(domain);
+ auto anasol = getV<3>(domain);
+ auto err = getV<4>(domain);
+ auto DCPSE_sol=getV<5>(domain);
+
+ // Here fill me
+
+ Box<2, double> up({box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
+
+ Box<2, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0});
+
+ Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
+ Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
+ openfpm::vector<Box<2, double>> boxes;
+ boxes.add(up);
+ boxes.add(down);
+ boxes.add(left);
+ boxes.add(right);
+
+ // Create a writer and write
+ VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
+ vtk_box.add(boxes);
+ vtk_box.write("vtk_box.vtk");
+
+
+ auto it2 = domain.getDomainIterator();
+
+ while (it2.isNext()) {
+ auto p = it2.get();
+ Point<2, double> xp = domain.getPos(p);
+ //domain.getProp<3>(p)=1+xp[0]*xp[0]+2*xp[1]*xp[1];
+ if (up.isInside(xp) == true) {
+ up_p.add();
+ up_p.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ } else if (down.isInside(xp) == true) {
+ dw_p.add();
+ dw_p.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+
+ } else if (left.isInside(xp) == true) {
+ l_p.add();
+ l_p.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+
+ } else if (right.isInside(xp) == true) {
+ r_p.add();
+ r_p.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+
+ } else {
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ }
+ ++it2;
+ }
+ DCPSE_scheme<equations,decltype(domain)> Solver(2 * rCut, domain);
+ auto Poisson = Lap(v);
+ auto D_x = Dx(v);
+ auto D_y = Dy(v);
+ Solver.impose(Poisson, bulk, prop_id<1>());
+ Solver.impose(v, up_p, 0);
+ Solver.impose(v, r_p, 0);
+ Solver.impose(v, dw_p, 0);
+ Solver.impose(v, l_p, 0);
+ Solver.solve(sol);
+ DCPSE_sol=Lap(sol);
+ double worst1 = 0.0;
+
+ v=abs(DCPSE_sol-RHS);
+
+ for(int j=0;j<bulk.size();j++)
+ { auto p=bulk.get<0>(j);
+ if (fabs(domain.getProp<3>(p) - domain.getProp<2>(p)) >= worst1) {
+ worst1 = fabs(domain.getProp<3>(p) - domain.getProp<2>(p));
+ }
+ domain.getProp<4>(p) = fabs(domain.getProp<3>(p) - domain.getProp<2>(p));
+
+ }
+ std::cout << "Maximum Analytic Error: " << worst1 << std::endl;
+
+ domain.write("Dirichlet_anasol");
+ }
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
BOOST_AUTO_TEST_CASE(dcpse_poisson_Robin) {
+ //http://e6.ijs.si/medusa/wiki/index.php/Poisson%27s_equation#Full_Neumann_boundary_conditions
// int rank;
// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
- const size_t sz[2] = {100,100};
+ const size_t sz[2] = {81,81};
+ Box<2, double> box({0, 0}, {1.0, 1.0});
+ size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
+ double spacing = box.getHigh(0) / (sz[0] - 1);
+ Ghost<2, double> ghost(spacing * 3);
+ double rCut = 2.0 * spacing;
+ BOOST_TEST_MESSAGE("Init vector_dist...");
+
+ vector_dist<2, double, aggregate<double,double,double,double,double,VectorS<2, double>>> domain(0, box, bc, ghost);
+
+
+ //Init_DCPSE(domain)
+ BOOST_TEST_MESSAGE("Init domain...");
+
+ auto it = domain.getGridIterator(sz);
+ while (it.isNext()) {
+ domain.add();
+
+ auto key = it.get();
+ double x = key.get(0) * it.getSpacing(0);
+ domain.getLastPos()[0] = x;
+ double y = key.get(1) * it.getSpacing(1);
+ domain.getLastPos()[1] = y;
+
+ ++it;
+ }
+ BOOST_TEST_MESSAGE("Sync domain across processors...");
+
+ domain.map();
+ domain.ghost_get<0>();
+
+ //Derivative_x Dx(domain, 2, rCut,2);
+ Derivative_y Dy(domain, 2, rCut,2);
+ //Gradient Grad(domain, 2, rCut);
+ Laplacian Lap(domain, 3, rCut, 3);
+ //Advection Adv(domain, 3, rCut, 3);
+ //Solver Sol_Lap(Lap),Sol_Dx(Dx);
+ //DCPSE_scheme<equations,decltype(domain)> Solver(2 * rCut, domain);
+ DCPSE_scheme<equations,decltype(domain)> Solver(2 * rCut, domain);
+
+ openfpm::vector<aggregate<int>> bulk;
+ openfpm::vector<aggregate<int>> up_p;
+ openfpm::vector<aggregate<int>> dw_p;
+ openfpm::vector<aggregate<int>> l_p;
+ openfpm::vector<aggregate<int>> r_p;
+
+ auto v = getV<0>(domain);
+ //auto RHS=getV<1>(domain);
+ auto sol = getV<2>(domain);
+ auto anasol = getV<3>(domain);
+ auto err = getV<4>(domain);
+ auto u = getV<5>(domain);
+
+ // Here fill me
+
+ Box<2, double> up({box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
+
+ Box<2, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0});
+
+ Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
+ Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
+ openfpm::vector<Box<2, double>> boxes;
+ boxes.add(up);
+ boxes.add(down);
+ boxes.add(left);
+ boxes.add(right);
+
+ // Create a writer and write
+ VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
+ vtk_box.add(boxes);
+ vtk_box.write("vtk_box.vtk");
+
+
+ auto it2 = domain.getDomainIterator();
+
+ while (it2.isNext()) {
+ auto p = it2.get();
+ Point<2, double> xp = domain.getPos(p);
+ //domain.getProp<3>(p)=1+xp[0]*xp[0]+2*xp[1]*xp[1];
+ if (up.isInside(xp) == true) {
+ up_p.add();
+ up_p.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = sin(5.0*xp.get(0));
+ } else if (down.isInside(xp) == true) {
+ dw_p.add();
+ dw_p.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = sin(5.0*xp.get(0));
+ } else if (left.isInside(xp) == true) {
+ l_p.add();
+ l_p.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = sin(5.0*xp.get(0));
+ } else if (right.isInside(xp) == true) {
+ r_p.add();
+ r_p.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = sin(5.0*xp.get(0));
+ } else {
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = -10.0*exp(-((xp.get(0)-0.5)*(xp.get(0)-0.5)+(xp.get(1)-0.5)*(xp.get(1)-0.5))/0.02);
+ }
+
+ ++it2;
+ }
+
+
+ auto Poisson = Lap(v);
+ //auto D_x = Dx(v);
+ auto D_y = Dy(v);
+ Solver.impose(Poisson, bulk, prop_id<1>());
+ Solver.impose(D_y, up_p, prop_id<1>());
+ Solver.impose(-D_y, dw_p, prop_id<1>());
+ Solver.impose(v, l_p, 0);
+ Solver.impose(v, r_p, 0);
+ Solver.solve(sol);
+ anasol=Lap(sol);
+ double worst1 = 0.0;
+
+ for(int j=0;j<bulk.size();j++)
+ { auto p=bulk.get<0>(j);
+ if (fabs(domain.getProp<3>(p) - domain.getProp<1>(p)) >= worst1) {
+ worst1 = fabs(domain.getProp<3>(p) - domain.getProp<1>(p));
+ }
+ domain.getProp<4>(p) = fabs(domain.getProp<1>(p) - domain.getProp<3>(p));
+
+ }
+ std::cout << "Maximum Auto Error: " << worst1 << std::endl;
+
+ domain.write("Mixed");
+ }
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+
+ BOOST_AUTO_TEST_CASE(dcpse_poisson_Neumann) {
+ //https://fenicsproject.org/docs/dolfin/1.4.0/python/demo/documented/neumann-poisson/python/documentation.html
+// int rank;
+// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+ const size_t sz[2] = {81,81};
Box<2, double> box({0, 0}, {1.0, 1.0});
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing = box.getHigh(0) / (sz[0] - 1);
@@ -704,13 +1324,14 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
domain.map();
domain.ghost_get<0>();
- Derivative_x Dx(domain, 3, rCut,2);
- Derivative_y Dy(domain, 4, rCut,3);
+ Derivative_x Dx(domain, 2, rCut,2);
+ Derivative_y Dy(domain, 2, rCut,2);
//Gradient Grad(domain, 2, rCut);
Laplacian Lap(domain, 3, rCut, 2);
//Advection Adv(domain, 3, rCut, 3);
//Solver Sol_Lap(Lap),Sol_Dx(Dx);
- DCPSE_scheme<equations,decltype(domain)> Solver(2 * rCut, domain);
+ //DCPSE_scheme<equations,decltype(domain)> Solver(2 * rCut, domain);
+ DCPSE_scheme<equations,decltype(domain)> Solver(2 * rCut, domain,options_solver::LAGRANGE_MULTIPLIER);
openfpm::vector<aggregate<int>> bulk;
openfpm::vector<aggregate<int>> up_p;
@@ -775,7 +1396,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
} else {
bulk.add();
bulk.last().get<0>() = p.getKey();
- domain.getProp<1>(p) = 10*exp(-((xp.get(0)-0.5)*(xp.get(0)-0.5)+(xp.get(1)-0.5)*(xp.get(1)-0.5))/0.02);
+ domain.getProp<1>(p) = -10*exp(-((xp.get(0)-0.5)*(xp.get(0)-0.5)+(xp.get(1)-0.5)*(xp.get(1)-0.5))/0.02);
}
++it2;
@@ -788,10 +1409,8 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Solver.impose(Poisson, bulk, prop_id<1>());
Solver.impose(D_y, up_p, prop_id<1>());
Solver.impose(-D_y, dw_p, prop_id<1>());
-// Solver.impose(-D_x, l_p, prop_id<1>());
-// Solver.impose(D_x, r_p, prop_id<1>());
- Solver.impose(v, l_p, 0);
- Solver.impose(v, r_p, 0);
+ Solver.impose(-D_x, l_p, prop_id<1>());
+ Solver.impose(D_x, r_p, prop_id<1>());
Solver.solve(sol);
anasol=-Lap(sol);
double worst1 = 0.0;
@@ -804,9 +1423,9 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
domain.getProp<4>(p) = fabs(domain.getProp<1>(p) - domain.getProp<3>(p));
}
- std::cout << "Maximum Error: " << worst1 << std::endl;
+ std::cout << "Maximum Auto Error: " << worst1 << std::endl;
- domain.write("particles1");
+ domain.write("Neumann");
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -868,27 +1487,43 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Box<2, double> up({box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0},
{box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
+ Box<2, double> up_d({box.getLow(0) - spacing / 2.0, box.getHigh(1) - 8*spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - 6*spacing / 2.0});
+
Box<2, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
{box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0});
+ Box<2, double> down_u({box.getLow(0) - spacing / 2.0, box.getLow(1) + 3*spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getLow(1) + 4*spacing / 2.0});
+
Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
{box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+ Box<2, double> left_r({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
{box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+ Box<2, double> right_l({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
+
openfpm::vector<Box<2, double>> boxes;
boxes.add(up);
+ boxes.add(up_d);
boxes.add(down);
+ boxes.add(down_u);
boxes.add(left);
+ boxes.add(left_r);
boxes.add(right);
+ boxes.add(right_l);
// Create a writer and write
VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
vtk_box.add(boxes);
vtk_box.write("vtk_box.vtk");
-
auto it2 = domain.getDomainIterator();
while (it2.isNext()) {
@@ -916,6 +1551,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
bulk.last().get<0>() = p.getKey();
}
+
++it2;
}
@@ -953,9 +1589,6 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
std::cout << "Maximum Error: " << worst1 << std::endl;
-
-
-
domain.write("particles");
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -966,7 +1599,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
BOOST_AUTO_TEST_CASE(dcpse_op_vec) {
// int rank;
// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
- size_t edgeSemiSize = 40;
+ size_t edgeSemiSize = 80;
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};
diff --git a/src/DCPSE_op/DCPSE_op_test2.cpp b/src/DCPSE_op/DCPSE_op_test2.cpp
index 1a724b478d8600bb7fd52f6835364889d83e2bdc..b91366e34dedecdaa8e1f059383abaa6e5a1fb48 100644
--- a/src/DCPSE_op/DCPSE_op_test2.cpp
+++ b/src/DCPSE_op/DCPSE_op_test2.cpp
@@ -1,8 +1,8 @@
/*
* DCPSE_op_test.cpp
*
- * Created on: Jan 7, 2020
- * Author: Abhinav Singh, Pietro Incardona
+ * Created on: Feb 25, 2020
+ * Author: Abhinav Singh
*
*/
@@ -18,6 +18,8 @@
#include "Operators/Vector/vector_dist_operators.hpp"
#include "Vector/vector_dist_subset.hpp"
+//int vector_dist_expression_op<void,void,VECT_COPY_N_TO_N>::i = 0;
+//int vector_dist_expression_op<void,void,VECT_COPY_1_TO_N>::i = 0;
//! Specify the general characteristic of system to solve
struct equations2 {
@@ -48,16 +50,1803 @@ const bool equations2::boundary[] = {NON_PERIODIC, NON_PERIODIC};
//struct Debug;
BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
+ BOOST_AUTO_TEST_CASE(dcpse_Lid_mirror) {
+ const size_t sz[2] = {81,81};
+ Box<2, double> box_mirror({-0.3, -0.3}, {1.3,1.3});
+ Box<2, double> box_domain({0.0, 0.0}, {1.0,1.0});
+ size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
+ double spacing = box_domain.getHigh(0) / (sz[0] - 1);
+ Ghost<2, double> ghost(spacing * 3);
+ double rCut = 2.0 * spacing;
+ // P V Dv RHS Vtemp
+ vector_dist<2, double, aggregate<double,VectorS<2, double>,VectorS<2, double>,double,VectorS<2, double>,double>> Particles(0, box_mirror, bc, ghost);
+// vector_dist<2, double, aggregate<double,VectorS<2, double>,VectorS<2, double>,double,VectorS<2, double>,double>> Particles_V(0, box, bc, ghost);
+// vector_dist<2, double, aggregate<double,VectorS<2, double>,VectorS<2, double>,double,VectorS<2, double>,double>> Particles_P(0, box, bc, ghost);
+
+ auto it = Particles.getGridIterator(sz);
+ while (it.isNext()) {
+ Particles.add();
+ auto key = it.get();
+ double x = key.get(0) * spacing;
+ Particles.getLastPos()[0] = x;
+ double y = key.get(1) * spacing;
+ Particles.getLastPos()[1] = y;
+ ++it;
+ }
+ Particles.map();
+ Particles.ghost_get<0>();
+
+ // Here fill up the boxes for particle detection.
+
+ Box<2, double> up({box_domain.getLow(0) - spacing / 2.0, box_domain.getHigh(1) - spacing / 2.0},
+ {box_domain.getHigh(0) + spacing / 2.0, box_domain.getHigh(1) + spacing / 2.0});
+
+ Box<2, double> down({box_domain.getLow(0) - spacing / 2.0, box_domain.getLow(1) - spacing / 2.0},
+ {box_domain.getHigh(0) + spacing / 2.0, box_domain.getLow(1) + spacing / 2.0});
+
+ Box<2, double> left({box_domain.getLow(0) - spacing / 2.0, box_domain.getLow(1) + spacing / 2.0},
+ {box_domain.getLow(0) + spacing / 2.0, box_domain.getHigh(1) - spacing / 2.0});
+
+ Box<2, double> right({box_domain.getHigh(0) - spacing / 2.0, box_domain.getLow(1) + spacing / 2.0},
+ {box_domain.getHigh(0) + spacing / 2.0, box_domain.getHigh(1) - spacing / 2.0});
+
+ /*Box<2, double> up_m({box_domain.getLow(0) + 7*spacing / 2.0, box_domain.getHigh(1) - 5*spacing / 2.0},
+ {box_domain.getHigh(0) - 7* spacing / 2.0, box_domain.getHigh(1) + spacing / 2.0});
+
+ Box<2, double> down_m({box_domain.getLow(0) + 7*spacing / 2.0, box_domain.getLow(1) - spacing / 2.0},
+ {box_domain.getHigh(0) - 7* spacing / 2.0, box_domain.getLow(1) + 5*spacing / 2.0});
+
+ Box<2, double> left_m({box_domain.getLow(0) - spacing / 2.0, box_domain.getLow(1) + 7*spacing / 2.0},
+ {box_domain.getLow(0) + 5*spacing / 2.0, box.getHigh(1) - 7*spacing / 2.0});
+
+ Box<2, double> right_m({box_domain.getHigh(0) - 5*spacing / 2.0, box_domain.getLow(1) + 7*spacing / 2.0},
+ {box_domain.getHigh(0) + spacing / 2.0, box_domain.getHigh(1) - 7*spacing / 2.0});
+*/
+ Box<2, double> up_mI({box_domain.getLow(0) + spacing / 2.0, box_domain.getHigh(1) - 7*spacing / 2.0},
+ {box_domain.getHigh(0) - spacing / 2.0, box_domain.getHigh(1) - spacing / 2.0});
+
+ Box<2, double> down_mI({box_domain.getLow(0) + spacing / 2.0, box_domain.getLow(1) + spacing / 2.0},
+ {box_domain.getHigh(0) - spacing / 2.0, box_domain.getLow(1) + 7*spacing / 2.0});
+
+ Box<2, double> left_mI({box_domain.getLow(0) + spacing / 2.0, box_domain.getLow(1) + spacing / 2.0},
+ {box_domain.getLow(0) + 7*spacing / 2.0, box_domain.getHigh(1) - spacing / 2.0});
+
+ Box<2, double> right_mI({box_domain.getHigh(0) - 7*spacing / 2.0, box_domain.getLow(1) + spacing / 2.0},
+ {box_domain.getHigh(0) - spacing / 2.0, box_domain.getHigh(1) - spacing / 2.0});
+
+
+ Box<2, double> Bulk_box({box_domain.getLow(0) + spacing / 2.0, box_domain.getLow(1) + spacing / 2.0},
+ {box_domain.getHigh(0) - spacing / 2.0, box_domain.getHigh(1) -spacing / 2.0});
+ Box<2, double> bulk_box({box_domain.getLow(0) + 7*spacing / 2.0, box_domain.getLow(1) + 7*spacing / 2.0},
+ {box_domain.getHigh(0) - 7*spacing / 2.0, box_domain.getHigh(1) - 7*spacing / 2.0});
+
+ openfpm::vector<Box<2, double>> boxes;
+ boxes.add(up);
+ boxes.add(down);
+ boxes.add(left);
+ boxes.add(right);
+ /*boxes.add(up_m);
+ boxes.add(down_m);
+ boxes.add(left_m);
+ boxes.add(right_m);*/
+ boxes.add(up_mI);
+ boxes.add(down_mI);
+ boxes.add(left_mI);
+ boxes.add(right_mI);
+ boxes.add(Bulk_box) ;
+ boxes.add(bulk_box);
+ VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
+ vtk_box.add(boxes);
+ vtk_box.write("box.vtk");
+
+
+ openfpm::vector<aggregate<int>> Bulk;
+ openfpm::vector<aggregate<int>> bulk;
+ openfpm::vector<aggregate<int>> up_p;
+ openfpm::vector<aggregate<int>> dw_p;
+ openfpm::vector<aggregate<int>> l_p;
+ openfpm::vector<aggregate<int>> r_p;
+ openfpm::vector<aggregate<int>> up_pM;
+ openfpm::vector<aggregate<int>> dw_pM;
+ openfpm::vector<aggregate<int>> l_pM;
+ openfpm::vector<aggregate<int>> r_pM;
+ openfpm::vector<aggregate<int>> up_pMI;
+ openfpm::vector<aggregate<int>> dw_pMI;
+ openfpm::vector<aggregate<int>> l_pMI;
+ openfpm::vector<aggregate<int>> r_pMI;
+ openfpm::vector<aggregate<int>> ref_p;
+ openfpm::vector<aggregate<int>> V_p;
+ openfpm::vector<aggregate<int>> P_p;
+
+ openfpm::vector<aggregate<int>> corner_ul;
+ openfpm::vector<openfpm::vector<aggregate<int>>> corner_ref_ul;
+
+ openfpm::vector<aggregate<int>> corner_ur;
+ openfpm::vector<openfpm::vector<aggregate<int>>> corner_ref_ur;
+
+ openfpm::vector<aggregate<int>> corner_dr;
+ openfpm::vector<openfpm::vector<aggregate<int>>> corner_ref_dr;
+
+ openfpm::vector<aggregate<int>> corner_dl;
+ openfpm::vector<openfpm::vector<aggregate<int>>> corner_ref_dl;
+
+ bool bulk_ref = false;
+ bool Bulk_ref = false;
+
+ auto it2 = Particles.getDomainIterator();
+ while (it2.isNext()) {
+ auto p = it2.get();
+ Point<2, double> xp = Particles.getPos(p);
+ if (up.isInside(xp) == true){
+ up_p.add();
+ up_p.last().get<0>() = p.getKey();
+ Particles.getProp<1>(p)[0] = 1;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ else if (down.isInside(xp) == true) {
+ if (xp[0]==0 && xp[1]==0) {
+ ref_p.add();
+ ref_p.last().get<0>() = p.getKey();
+ }
+ else {
+ dw_p.add();
+ dw_p.last().get<0>() = p.getKey();
+ }
+ }
+ else if (left.isInside(xp) == true) {
+ l_p.add();
+ l_p.last().get<0>() = p.getKey();
+ }
+ else if (right.isInside(xp) == true) {
+ r_p.add();
+ r_p.last().get<0>() = p.getKey();
+ }
+/* else if (up_m.isInside(xp) == true){
+ up_pM.add();
+ up_pM.last().get<0>() = p.getKey();
+ }
+ else if (down_m.isInside(xp) == true) {
+ dw_pM.add();
+ dw_pM.last().get<0>() = p.getKey();
+ }
+ else if (left_m.isInside(xp) == true) {
+ l_pM.add();
+ l_pM.last().get<0>() = p.getKey();
+ }
+ else if (right_m.isInside(xp) == true) {
+ r_pM.add();
+ r_pM.last().get<0>() = p.getKey();
+ }*/
+
+
+ else if(Bulk_box.isInside(xp) == true) {
+ Bulk.add();
+ Bulk.last().get<0>() = p.getKey();
+ if (up_mI.isInside(xp) == true){
+ if (left_mI.isInside(xp) == true)
+ {
+ corner_ref_ul.add();
+
+ Particles.add();
+ Particles.getLastPos()[0] = Particles.getPos(p)[0];
+ Particles.getLastPos()[1] = 2 * box_domain.getHigh(1) - Particles.getPos(p)[1];
+
+ Particles.getLastProp<1>()[0] = 0.0;
+ Particles.getLastProp<1>()[1] = 0.0;
+
+ corner_ref_ul.last().add();
+ corner_ref_ul.last().last().get<0>() = Particles.size_local() - 1;
+
+ Particles.add();
+ Particles.getLastPos()[0] = box_domain.getLow(0) - Particles.getPos(p)[0];
+ Particles.getLastPos()[1] = Particles.getPos(p)[1];
+
+ Particles.getLastProp<1>()[0] = 0.0;
+ Particles.getLastProp<1>()[1] = 0.0;
+
+ corner_ref_ul.last().add();
+ corner_ref_ul.last().last().get<0>() = Particles.size_local() - 1;
+
+ corner_ul.add();
+ corner_ul.last().get<0>() = p.getKey();
+ }
+ else if (right_mI.isInside(xp) == true)
+ {
+ corner_ref_ur.add();
+
+ Particles.add();
+ Particles.getLastPos()[0] = Particles.getPos(p)[0];
+ Particles.getLastPos()[1] = 2 * box_domain.getHigh(1) - Particles.getPos(p)[1];
+
+ Particles.getLastProp<1>()[0] = 0.0;
+ Particles.getLastProp<1>()[1] = 0.0;
+
+ corner_ref_ur.last().add();
+ corner_ref_ur.last().last().get<0>() = Particles.size_local() - 1;
+
+ Particles.add();
+ Particles.getLastPos()[0] = 2 * box_domain.getHigh(0) - Particles.getPos(p)[0];
+ Particles.getLastPos()[1] = Particles.getPos(p)[1];
+
+ Particles.getLastProp<1>()[0] = 0.0;
+ Particles.getLastProp<1>()[1] = 0.0;
+
+ corner_ref_ur.last().add();
+ corner_ref_ur.last().last().get<0>() = Particles.size_local() - 1;
+
+ corner_ur.add();
+ corner_ur.last().get<0>() = p.getKey();
+ }
+ else
+ {
+ Particles.add();
+ Particles.getLastPos()[0] = Particles.getPos(p)[0];
+ Particles.getLastPos()[1] = 2 * box_domain.getHigh(1) - Particles.getPos(p)[1];
+
+ Particles.getLastProp<1>()[0] = 0.0;
+ Particles.getLastProp<1>()[1] = 0.0;
+
+ up_pM.add();
+ up_pM.last().get<0>() = Particles.size_local()-1;
+
+ up_pMI.add();
+ up_pMI.last().get<0>() = p.getKey();
+ }
+ }
+ if (down_mI.isInside(xp) == true) {
+ if (left_mI.isInside(xp) == true)
+ {
+ corner_ref_dl.add();
+
+ Particles.add();
+ Particles.getLastPos()[0] = Particles.getPos(p)[0];
+ Particles.getLastPos()[1] = box_domain.getLow(1) - Particles.getPos(p)[1];
+
+ Particles.getLastProp<1>()[0] = 0.0;
+ Particles.getLastProp<1>()[1] = 0.0;
+
+ corner_ref_dl.last().add();
+ corner_ref_dl.last().last().get<0>() = Particles.size_local() - 1;
+
+ Particles.add();
+ Particles.getLastPos()[0] = box_domain.getLow(0) - Particles.getPos(p)[0];
+ Particles.getLastPos()[1] = Particles.getPos(p)[1];
+
+ Particles.getLastProp<1>()[0] = 0.0;
+ Particles.getLastProp<1>()[1] = 0.0;
+
+ corner_ref_dl.last().add();
+ corner_ref_dl.last().last().get<0>() = Particles.size_local() - 1;
+
+ corner_dl.add();
+ corner_dl.last().get<0>() = p.getKey();
+ }
+ else if (right_mI.isInside(xp) == true)
+ {
+ corner_ref_dr.add();
+
+ Particles.add();
+ Particles.getLastPos()[0] = Particles.getPos(p)[0];
+ Particles.getLastPos()[1] = box_domain.getLow(1) - Particles.getPos(p)[1];
+
+ Particles.getLastProp<1>()[0] = 0.0;
+ Particles.getLastProp<1>()[1] = 0.0;
+
+ corner_ref_dr.last().add();
+ corner_ref_dr.last().last().get<0>() = Particles.size_local() - 1;
+
+ Particles.add();
+ Particles.getLastPos()[0] = 2 * box_domain.getHigh(0) - Particles.getPos(p)[0];
+ Particles.getLastPos()[1] = Particles.getPos(p)[1];
+
+ Particles.getLastProp<1>()[0] = 0.0;
+ Particles.getLastProp<1>()[1] = 0.0;
+
+ corner_ref_dr.last().add();
+ corner_ref_dr.last().last().get<0>() = Particles.size_local() - 1;
+
+ corner_dr.add();
+ corner_dr.last().get<0>() = p.getKey();
+ }
+ else {
+
+ Particles.add();
+ Particles.getLastPos()[0] = Particles.getPos(p)[0];
+ Particles.getLastPos()[1] = box_domain.getLow(1) - Particles.getPos(p)[1];
+
+ Particles.getLastProp<1>()[0] = 0.0;
+ Particles.getLastProp<1>()[1] = 0.0;
+
+ dw_pM.add();
+ dw_pM.last().get<0>() = Particles.size_local()-1;
+
+ dw_pMI.add();
+ dw_pMI.last().get<0>() = p.getKey();
+ }
+ }
+ if (left_mI.isInside(xp) == true) {
+ if (down_mI.isInside(xp) == true)
+ {
+/*
+ Particles.add();
+ Particles.getLastPos()[0] = box_domain.getLow(0) - Particles.getPos(p)[0];
+ Particles.getLastPos()[1] = Particles.getPos(p)[1];
+
+ Particles.getLastProp<1>()[0] = 0.0;
+ Particles.getLastProp<1>()[1] = 0.0;
+*/
+ }
+ else if (up_mI.isInside(xp) == true)
+ {
+
+ }
+ else {
+ Particles.add();
+ Particles.getLastPos()[0] = box_domain.getLow(0) - Particles.getPos(p)[0];
+ Particles.getLastPos()[1] = Particles.getPos(p)[1];
+
+ Particles.getLastProp<1>()[0] = 0.0;
+ Particles.getLastProp<1>()[1] = 0.0;
+
+ l_pM.add();
+ l_pM.last().get<0>() = Particles.size_local()-1;
+
+ l_pMI.add();
+ l_pMI.last().get<0>() = p.getKey();
+ }
+ }
+ if (right_mI.isInside(xp) == true) {
+ if (down_mI.isInside(xp) == true)
+ {
+
+ }
+ else if (up_mI.isInside(xp) == true)
+ {
+
+ }
+ else {
+
+ Particles.add();
+ Particles.getLastPos()[0] = 2 * box_domain.getHigh(0) - Particles.getPos(p)[0];
+ Particles.getLastPos()[1] = Particles.getPos(p)[1];
+
+ Particles.getLastProp<1>()[0] = 0.0;
+ Particles.getLastProp<1>()[1] = 0.0;
+
+ r_pM.add();
+ r_pM.last().get<0>() = Particles.size_local()-1;
+
+ r_pMI.add();
+ r_pMI.last().get<0>() = p.getKey();
+ }
+ }
+ if(bulk_box.isInside(xp)==true){
+/* if (bulk_ref == false)
+ {
+ bulk_ref = true;
+ ref_p.add();
+ ref_p.last().get<0>() = p.getKey();
+ }*/
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ }
+ }
+
+ ++it2;
+ }
+
+ Particles.write("Particles");
+ Derivative_x Dx(Particles, 2, rCut, 2);
+ Derivative_y Dy(Particles, 2, rCut, 2);
+
+ Gradient Grad(Particles, 2, rCut, 2);
+ Laplacian Lap(Particles, 2, rCut, 2);
+ Advection Adv(Particles, 2, rCut, 2);
+ Divergence Div(Particles, 2, rCut, 2);
+
+ auto P = getV<0>(Particles);
+ auto V = getV<1>(Particles);
+ auto dV = getV<2>(Particles);
+ auto RHS = getV<3>(Particles);
+ auto Vtemp = getV<4>(Particles);
+ auto divV = getV<5>(Particles);
+
+ double dt=0.005;
+ int n=3;
+ double nu=1e-2;
+
+ Particles.write_frame("Re1000-1e-4-Lid",0);
+ std::cout<<"Init Done"<<std::endl;
+
+ vector_dist_expression_op<void, void, VECT_COPY_N_TO_N> vcp_up(up_pM, up_pMI);
+ vector_dist_expression_op<void, void, VECT_COPY_N_TO_N> vcp_l(l_pM, l_pMI);
+ vector_dist_expression_op<void, void, VECT_COPY_N_TO_N> vcp_r(r_pM, r_pMI);
+ vector_dist_expression_op<void, void, VECT_COPY_N_TO_N> vcp_dw(dw_pM, dw_pMI);
+
+ for(int i =1; i<=n ;i++)
+ { dV=V+dt*(nu*Lap(V)-Adv(V,V));
+ std::cout<<"dV Done"<<std::endl;
+ RHS=1.0/dt*Div(dV);
+
+ std::cout<<"RHS Done"<<std::endl;
+ DCPSE_scheme<equations2,decltype(Particles)> Solver(2*rCut, Particles);
+ auto Pressure_Poisson = Lap(P);
+ auto D_x=Dx(P);
+ auto D_y=Dy(P);
+ Solver.impose(Pressure_Poisson, Bulk, prop_id<3>());
+ Solver.impose(D_y, up_p, prop_id<3>());
+ Solver.impose(-D_y, dw_p, prop_id<3>());
+ Solver.impose(-D_x, l_p, prop_id<3>());
+ Solver.impose(D_x, r_p, prop_id<3>());
+ Solver.impose(vcp_up, up_pM, 0);
+ Solver.impose(vcp_l, l_pM,0);
+ Solver.impose(vcp_r, r_pM, 0);
+ Solver.impose(vcp_dw, dw_pM, 0);
+
+ for (int i = 0 ; i < corner_ur.size() ; i++)
+ {
+ vector_dist_expression_op<void, void, VECT_COPY_1_TO_N> cor(corner_ref_ur.template get<0>(i),corner_ur.template get<0>(i));
+
+ Solver.impose(cor,corner_ref_ur.template get<0>(i), 0);
+ }
+
+ for (int i = 0 ; i < corner_ul.size() ; i++)
+ {
+ vector_dist_expression_op<void, void, VECT_COPY_1_TO_N> cor(corner_ref_ul.template get<0>(i),corner_ul.template get<0>(i));
+
+ Solver.impose(cor,corner_ref_ul.template get<0>(i), 0);
+ }
+
+ for (int i = 0 ; i < corner_dl.size() ; i++)
+ {
+ vector_dist_expression_op<void, void, VECT_COPY_1_TO_N> cor(corner_ref_dl.template get<0>(i),corner_dl.template get<0>(i));
+
+ Solver.impose(cor,corner_ref_dl.template get<0>(i), 0);
+ }
+
+ for (int i = 0 ; i < corner_dr.size() ; i++)
+ {
+ vector_dist_expression_op<void, void, VECT_COPY_1_TO_N> cor(corner_ref_dr.template get<0>(i),corner_dr.template get<0>(i));
+
+ Solver.impose(cor,corner_ref_dr.template get<0>(i), 0);
+ }
+ Solver.impose(P,ref_p,0);
+ Solver.solve(P);
+ std::cout<<"Poisson Solved"<<std::endl;
+ divV=P-Lap(P);
+ Particles.write_frame("Re1000-1e-4-Lid",i);
+ return;
+ Vtemp=dV-dt*Grad(P);
+ V=Vtemp;
+ for(int j=0;j<up_p.size();j++)
+ { auto p=up_p.get<0>(j);
+ Particles.getProp<1>(p)[0] = 1;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ for(int j=0;j<l_p.size();j++)
+ { auto p=l_p.get<0>(j);
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ for(int j=0;j<r_p.size();j++)
+ { auto p=r_p.get<0>(j);
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ for(int j=0;j<dw_p.size();j++)
+ { auto p=dw_p.get<0>(j);
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ Particles.getProp<1>(0)[0] = 0;
+ Particles.getProp<1>(0)[1] = 0;
+ std::cout<<"Velocity updated"<<std::endl;
+ Particles.write_frame("Re1000-1e-4-Lid",i);
+ std::cout<<i<<std::endl;
+ }
+
+ }
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ BOOST_AUTO_TEST_CASE(dcpse_Lid_mirr) {
+ const size_t sz[2] = {31,31};
+ Box<2, double> box({0, 0}, {1,1});
+ size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
+ double spacing = box.getHigh(0) / (sz[0] - 1);
+ Ghost<2, double> ghost(spacing * 3);
+ double rCut = 2.0 * spacing;
+ // P V Dv RHS Vtemp
+ vector_dist<2, double, aggregate<double,VectorS<2, double>,VectorS<2, double>,double,VectorS<2, double>,double>> Particles(0, box, bc, ghost);
+ vector_dist<2, double, aggregate<double,VectorS<2, double>,VectorS<2, double>,double,VectorS<2, double>,double>> Particles_V(0, box, bc, ghost);
+ vector_dist<2, double, aggregate<double,VectorS<2, double>,VectorS<2, double>,double,VectorS<2, double>,double>> Particles_P(0, box, bc, ghost);
+
+ auto it = Particles.getGridIterator(sz);
+ while (it.isNext()) {
+ Particles.add();
+ auto key = it.get();
+ double x = key.get(0) * it.getSpacing(0);
+ Particles.getLastPos()[0] = x;
+ double y = key.get(1) * it.getSpacing(1);
+ Particles.getLastPos()[1] = y;
+ ++it;
+ }
+ Particles.map();
+ Particles.ghost_get<0>();
+
+ openfpm::vector<aggregate<int>> Bulk;
+ openfpm::vector<aggregate<int>> bulk;
+ openfpm::vector<aggregate<int>> up_p;
+ openfpm::vector<aggregate<int>> dw_p;
+ openfpm::vector<aggregate<int>> l_p;
+ openfpm::vector<aggregate<int>> r_p;
+ openfpm::vector<aggregate<int>> up_pM;
+ openfpm::vector<aggregate<int>> dw_pM;
+ openfpm::vector<aggregate<int>> l_pM;
+ openfpm::vector<aggregate<int>> r_pM;
+ openfpm::vector<aggregate<int>> up_pMI;
+ openfpm::vector<aggregate<int>> dw_pMI;
+ openfpm::vector<aggregate<int>> l_pMI;
+ openfpm::vector<aggregate<int>> r_pMI;
+ openfpm::vector<aggregate<int>> ref_p;
+ openfpm::vector<aggregate<int>> V_p;
+ openfpm::vector<aggregate<int>> P_p;
+
+ // Here fill up the boxes for particle detection.
+
+ Box<2, double> up({box.getLow(0) + 5*spacing / 2.0, box.getHigh(1) - 7*spacing / 2.0},
+ {box.getHigh(0) - 5* spacing / 2.0, box.getHigh(1) - 5*spacing / 2.0});
+
+ Box<2, double> down({box.getLow(0) + 5*spacing / 2.0, box.getLow(1) + 5*spacing / 2.0},
+ {box.getHigh(0) - 5* spacing / 2.0, box.getLow(1) + 7*spacing / 2.0});
+
+ Box<2, double> left({box.getLow(0) + 5*spacing / 2.0, box.getLow(1) + 7*spacing / 2.0},
+ {box.getLow(0) + 7*spacing / 2.0, box.getHigh(1) - 7*spacing / 2.0});
+
+ Box<2, double> right({box.getHigh(0) - 7*spacing / 2.0, box.getLow(1) + 7*spacing / 2.0},
+ {box.getHigh(0) - 5*spacing / 2.0, box.getHigh(1) - 7*spacing / 2.0});
+
+ Box<2, double> up_m({box.getLow(0) + 7*spacing / 2.0, box.getHigh(1) - 5*spacing / 2.0},
+ {box.getHigh(0) - 7* spacing / 2.0, box.getHigh(1) + spacing / 2.0});
+
+ Box<2, double> down_m({box.getLow(0) + 7*spacing / 2.0, box.getLow(1) - spacing / 2.0},
+ {box.getHigh(0) - 7* spacing / 2.0, box.getLow(1) + 5*spacing / 2.0});
+
+ Box<2, double> left_m({box.getLow(0) - spacing / 2.0, box.getLow(1) + 7*spacing / 2.0},
+ {box.getLow(0) + 5*spacing / 2.0, box.getHigh(1) - 7*spacing / 2.0});
+
+ Box<2, double> right_m({box.getHigh(0) - 5*spacing / 2.0, box.getLow(1) + 7*spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - 7*spacing / 2.0});
+
+ Box<2, double> up_mI({box.getLow(0) + 7*spacing / 2.0, box.getHigh(1) - 13*spacing / 2.0},
+ {box.getHigh(0) - 7* spacing / 2.0, box.getHigh(1) - 7*spacing / 2.0});
+
+ Box<2, double> down_mI({box.getLow(0) + 7*spacing / 2.0, box.getLow(1) + 7*spacing / 2.0},
+ {box.getHigh(0) - 7* spacing / 2.0, box.getLow(1) + 13*spacing / 2.0});
+
+ Box<2, double> left_mI({box.getLow(0) + 7*spacing / 2.0, box.getLow(1) + 7*spacing / 2.0},
+ {box.getLow(0) + 13*spacing / 2.0, box.getHigh(1) - 7*spacing / 2.0});
+
+ Box<2, double> right_mI({box.getHigh(0) - 13*spacing / 2.0, box.getLow(1) + 7*spacing / 2.0},
+ {box.getHigh(0) - 7* spacing / 2.0, box.getHigh(1) - 7*spacing / 2.0});
+
+
+ Box<2, double> Bulk_box({box.getLow(0) + 7*spacing / 2.0, box.getLow(1) + 7*spacing / 2.0},
+ {box.getHigh(0) - 7*spacing / 2.0, box.getHigh(1) - 7*spacing / 2.0});
+ Box<2, double> bulk_box({box.getLow(0) + 13*spacing / 2.0, box.getLow(1) + 13*spacing / 2.0},
+ {box.getHigh(0) - 13*spacing / 2.0, box.getHigh(1) - 13*spacing / 2.0});
+
+ openfpm::vector<Box<2, double>> boxes;
+ boxes.add(up);
+ boxes.add(down);
+ boxes.add(left);
+ boxes.add(right);
+ boxes.add(up_m);
+ boxes.add(down_m);
+ boxes.add(left_m);
+ boxes.add(right_m);
+ boxes.add(up_mI);
+ boxes.add(down_mI);
+ boxes.add(left_mI);
+ boxes.add(right_mI);
+ boxes.add(Bulk_box) ;
+ boxes.add(bulk_box);
+ VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
+ vtk_box.add(boxes);
+ vtk_box.write("box.vtk");
+ auto it2 = Particles.getDomainIterator();
+ while (it2.isNext()) {
+ auto p = it2.get();
+ Point<2, double> xp = Particles.getPos(p);
+ if (up.isInside(xp) == true){
+ up_p.add();
+ up_p.last().get<0>() = p.getKey();
+ Particles.getProp<1>(p)[0] = 1;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ else if (down.isInside(xp) == true) {
+ dw_p.add();
+ dw_p.last().get<0>() = p.getKey();
+ }
+ else if (left.isInside(xp) == true) {
+ l_p.add();
+ l_p.last().get<0>() = p.getKey();
+ }
+ else if (right.isInside(xp) == true) {
+ r_p.add();
+ r_p.last().get<0>() = p.getKey();
+ }
+ else if (up_m.isInside(xp) == true){
+ up_pM.add();
+ up_pM.last().get<0>() = p.getKey();
+ }
+ else if (down_m.isInside(xp) == true) {
+ dw_pM.add();
+ dw_pM.last().get<0>() = p.getKey();
+ }
+ else if (left_m.isInside(xp) == true) {
+ l_pM.add();
+ l_pM.last().get<0>() = p.getKey();
+ }
+ else if (right_m.isInside(xp) == true) {
+ r_pM.add();
+ r_pM.last().get<0>() = p.getKey();
+ }
+ else if(Bulk_box.isInside(xp) == true){
+ Bulk.add();
+ Bulk.last().get<0>() = p.getKey();
+ if (up_mI.isInside(xp) == true){
+ up_pMI.add();
+ up_pMI.last().get<0>() = p.getKey();
+ }
+ if (down_mI.isInside(xp) == true) {
+ dw_pMI.add();
+ dw_pMI.last().get<0>() = p.getKey();
+ }
+ if (left_mI.isInside(xp) == true) {
+ l_pMI.add();
+ l_pMI.last().get<0>() = p.getKey();
+ }
+ if (right_m.isInside(xp) == true) {
+ r_pMI.add();
+ r_pMI.last().get<0>() = p.getKey();
+ }
+ if(bulk_box.isInside(xp)==true){
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ }
+ }
+ ++it2;
+ }
+ for (int i = 0 ; i < up_p.size() ; i++) {
+ Particles_V.add();
+ Particles_V.getLastPos()[0] = Particles.getPos(up_p.template get<0>(i))[0];
+ Particles_V.getLastPos()[1] = Particles.getPos(up_p.template get<0>(i))[1];
+ Particles_V.getLastProp<1>() = Particles.template getProp<1>(up_p.template get<0>(i));
+ V_p.add();
+ V_p.last().get<0>() = up_p.template get<0>(i);
+ }
+
+ for (int i = 0 ; i < dw_p.size() ; i++) {
+ Particles_V.add();
+ Particles_V.getLastPos()[0] = Particles.getPos(dw_p.template get<0>(i))[0];
+ Particles_V.getLastPos()[1] = Particles.getPos(dw_p.template get<0>(i))[1];
+ Particles_V.getLastProp<1>() = Particles.template getProp<1>(dw_p.template get<0>(i));
+ V_p.add();
+ V_p.last().get<0>() = dw_p.template get<0>(i);
+ }
+ for (int i = 0 ; i < l_p.size() ; i++) {
+ Particles_V.add();
+ Particles_V.getLastPos()[0] = Particles.getPos(l_p.template get<0>(i))[0];
+ Particles_V.getLastPos()[1] = Particles.getPos(l_p.template get<0>(i))[1];
+ Particles_V.getLastProp<1>() = Particles.template getProp<1>(l_p.template get<0>(i));
+ V_p.add();
+ V_p.last().get<0>() = l_p.template get<0>(i);
+ }
+ for (int i = 0 ; i < r_p.size() ; i++) {
+ Particles_V.add();
+ Particles_V.getLastPos()[0] = Particles.getPos(r_p.template get<0>(i))[0];
+ Particles_V.getLastPos()[1] = Particles.getPos(r_p.template get<0>(i))[1];
+ Particles_V.getLastProp<1>() = Particles.template getProp<1>(r_p.template get<0>(i));
+ V_p.add();
+ V_p.last().get<0>() = r_p.template get<0>(i);
+ }
+
+ for (int i = 0 ; i < Bulk.size() ; i++) {
+ Particles_V.add();
+ Particles_V.getLastPos()[0] = Particles.getPos(Bulk.template get<0>(i))[0];
+ Particles_V.getLastPos()[1] = Particles.getPos(Bulk.template get<0>(i))[1];
+ Particles_V.getLastProp<1>() = Particles.template getProp<1>(Bulk.template get<0>(i));
+ V_p.add();
+ V_p.last().get<0>() = Bulk.template get<0>(i);
+
+ Particles_P.add();
+ Particles_P.getLastPos()[0] = Particles.getPos(Bulk.template get<0>(i))[0];
+ Particles_P.getLastPos()[1] = Particles.getPos(Bulk.template get<0>(i))[1];
+ Particles_P.getLastProp<1>() = Particles.template getProp<1>(Bulk.template get<0>(i));
+ P_p.add();
+ P_p.last().get<0>() = Bulk.template get<0>(i);
+ }
+
+ for (int i = 0 ; i < up_pM.size() ; i++) {
+ Particles_P.add();
+ Particles_P.getLastPos()[0] = Particles.getPos(up_pM.template get<0>(i))[0];
+ Particles_P.getLastPos()[1] = Particles.getPos(up_pM.template get<0>(i))[1];
+ Particles_P.getLastProp<1>() = Particles.template getProp<1>(up_pM.template get<0>(i));
+ P_p.add();
+ P_p.last().get<0>() = Bulk.template get<0>(i);
+ }
+ for (int i = 0 ; i < dw_pM.size() ; i++) {
+ Particles_P.add();
+ Particles_P.getLastPos()[0] = Particles.getPos(dw_pM.template get<0>(i))[0];
+ Particles_P.getLastPos()[1] = Particles.getPos(dw_pM.template get<0>(i))[1];
+ Particles_P.getLastProp<1>() = Particles.template getProp<1>(dw_pM.template get<0>(i));
+ P_p.add();
+ P_p.last().get<0>() = Bulk.template get<0>(i);
+ }
+ for (int i = 0 ; i < l_pM.size() ; i++) {
+ Particles_P.add();
+ Particles_P.getLastPos()[0] = Particles.getPos(l_pM.template get<0>(i))[0];
+ Particles_P.getLastPos()[1] = Particles.getPos(l_pM.template get<0>(i))[1];
+ Particles_P.getLastProp<1>() = Particles.template getProp<1>(l_pM.template get<0>(i));
+ P_p.add();
+ P_p.last().get<0>() = Bulk.template get<0>(i);
+ }
+ for (int i = 0 ; i < r_pM.size() ; i++) {
+ Particles_P.add();
+ Particles_P.getLastPos()[0] = Particles.getPos(r_pM.template get<0>(i))[0];
+ Particles_P.getLastPos()[1] = Particles.getPos(r_pM.template get<0>(i))[1];
+ Particles_P.getLastProp<1>() = Particles.template getProp<1>(r_pM.template get<0>(i));
+ P_p.add();
+ P_p.last().get<0>() = Bulk.template get<0>(i);
+ }
+ Particles_V.map();
+ Particles_V.ghost_get<0>();
+ Particles_P.map();
+ Particles_P.ghost_get<0>();
+ Particles_V.write_frame("V",0);
+ Particles_P.write_frame("P",0);
+
+ Gradient Grad_P(Particles_P, 2, rCut, 2);
+ Laplacian Lap_P(Particles_P, 2, rCut, 2),Lap_V(Particles_V, 2, rCut, 2);
+ Advection Adv(Particles_V, 2, rCut, 2);
+ Divergence Div_V(Particles_V, 2, rCut, 2);
+
+ auto P = getV<0>(Particles_P);
+ auto V = getV<1>(Particles_V);
+ auto dV = getV<2>(Particles_V);
+ auto RHS = getV<3>(Particles_V);
+ auto Vtemp = getV<4>(Particles_V);
+ auto divV = getV<5>(Particles_V);
+
+ double dt=0.005;
+ int n=3;
+ double nu=1e-2;
+
+ vector_dist_expression_op<void, void, VECT_COPY_N_TO_N> vcp_up(up_pM, up_pMI);
+ vector_dist_expression_op<void, void, VECT_COPY_N_TO_N> vcp_l(l_pM, l_pMI);
+ vector_dist_expression_op<void, void, VECT_COPY_N_TO_N> vcp_r(r_pM, r_pMI);
+ vector_dist_expression_op<void, void, VECT_COPY_N_TO_N> vcp_dw(dw_pM, dw_pMI);
+
+ for (int i = 0 ; i < V_p.size() ; i++) {
+ Particles_V.getProp<1>(i) = Particles.template getProp<1>(V_p.template get<0>(i));
+ }
+ Particles_V.write_frame("V",1);
+ Particles.write_frame("Re1000-1e-4-Lid",0);
+ std::cout<<"Init Done"<<std::endl;
+ for(int i =1; i<=n ;i++)
+ { dV=V+dt*(nu*Lap_V(V)-Adv(V,V));
+ std::cout<<"dV Done"<<std::endl;
+ RHS=1.0/dt*Div_V(dV);
+ for (int i = 0 ; i < V_p.size() ; i++) {
+ Particles.getProp<3>(V_p.template get<0>(i)) = Particles_V.template getProp<3>(i);
+ }
+ for (int i = 0 ; i < P_p.size() ; i++) {
+ Particles_P.getProp<3>(i) = Particles.template getProp<3>(P_p.template get<0>(i));
+ }
+ std::cout<<"RHS Done"<<std::endl;
+ DCPSE_scheme<equations2,decltype(Particles_P)> Solver(2*rCut, Particles_P);
+ auto Pressure_Poisson = Lap_P(P);
+ Solver.impose(Pressure_Poisson, bulk, prop_id<3>());
+ Solver.impose(vcp_up, up_pM, 0);
+ Solver.impose(vcp_l, l_pM,0);
+ Solver.impose(vcp_r, r_pM, 0);
+ Solver.impose(vcp_dw, dw_pM, 0);
+ Solver.solve(P);
+ std::cout<<"Poisson Solved"<<std::endl;
+ for (int i = 0 ; i < P_p.size() ; i++) {
+ Particles.getProp<3>(i) = Particles_P.template getProp<3>(P_p.template get<0>(i));
+ }
+ Particles.write_frame("Re1000-1e-4-Lid",i);
+ return;
+/* k1=dt*(dV-Grad(P));
+ Vtemp=V+k1*0.5;
+ k2=dt*(nu*Lap(Vtemp)-Adv(Vtemp,Vtemp)-Grad(P));
+ Vtemp=V+k2*0.5;
+ k3=dt*(nu*Lap(Vtemp)-Adv(Vtemp,Vtemp)-Grad(P));
+ Vtemp=V+k3;
+ k4=dt*(nu*Lap(Vtemp)-Adv(Vtemp,Vtemp)-Grad(P));
+ Vtemp=V+0.16666666666*(k1+2*k2+2*k3+k4);*/
+ Vtemp=dV-dt*Grad_P(P);
+ V=Vtemp;
+ for(int j=0;j<up_p.size();j++)
+ { auto p=up_p.get<0>(j);
+ Particles.getProp<1>(p)[0] = 1;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ for(int j=0;j<l_p.size();j++)
+ { auto p=l_p.get<0>(j);
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ for(int j=0;j<r_p.size();j++)
+ { auto p=r_p.get<0>(j);
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ for(int j=0;j<dw_p.size();j++)
+ { auto p=dw_p.get<0>(j);
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ Particles.getProp<1>(0)[0] = 0;
+ Particles.getProp<1>(0)[1] = 0;
+ for (int i = 0 ; i < bulk.size() ; i++) {
+ Particles_P.getProp<1>(i) = Particles.template getProp<1>(bulk.template get<0>(i));
+ }
+// divV_bulk=Div_bulk(V_bulk);
+ for (int i = 0 ; i < bulk.size() ; i++)
+ {
+ Particles.template getProp<5>(bulk.template get<0>(i)) = Particles_P.getProp<0>(i);
+ }
+ //divV=Div(V);
+ std::cout<<"Velocity updated"<<std::endl;
+ Particles.write_frame("Re1000-1e-4-Lid",i);
+ std::cout<<i<<std::endl;
+ }
+ }
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+
+
+
+
+
+
+
+
+
+
+ BOOST_AUTO_TEST_CASE(dcpse_Lid2) {
+ const size_t sz[2] = {31, 31};
+ Box<2, double> box({0, 0}, {1, 1});
+ size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
+ double spacing = box.getHigh(0) / (sz[0] - 1);
+ Ghost<2, double> ghost(spacing * 3);
+ double rCut = 2.0 * spacing;
+ // P V Dv RHS Vtemp
+ vector_dist<2, double, aggregate<double, VectorS<2, double>, VectorS<2, double>, double, VectorS<2, double>, double>> Particles(
+ 0, box, bc, ghost);
+ auto it = Particles.getGridIterator(sz);
+ while (it.isNext()) {
+ Particles.add();
+ auto key = it.get();
+ double x = key.get(0) * it.getSpacing(0);
+ Particles.getLastPos()[0] = x;
+ double y = key.get(1) * it.getSpacing(1);
+ Particles.getLastPos()[1] = y;
+ ++it;
+ }
+
+ Particles.map();
+ Particles.ghost_get<0>();
+
+ openfpm::vector<aggregate<int>> bulk;
+ openfpm::vector<aggregate<int>> up_p;
+ openfpm::vector<aggregate<int>> dw_p;
+ openfpm::vector<aggregate<int>> l_p;
+ openfpm::vector<aggregate<int>> r_p;
+ openfpm::vector<aggregate<int>> up_p1;
+ openfpm::vector<aggregate<int>> dw_p1;
+ openfpm::vector<aggregate<int>> l_p1;
+ openfpm::vector<aggregate<int>> r_p1;
+ openfpm::vector<aggregate<int>> ref_p;
+ openfpm::vector<aggregate<int>> ref_p2;
+ openfpm::vector<aggregate<int>> FullBulk;
+ openfpm::vector<aggregate<int>> UL_p;
+ openfpm::vector<aggregate<int>> UL_p_ref;
+ openfpm::vector<aggregate<int>> UR_p;
+ openfpm::vector<aggregate<int>> UR_p_ref;
+ openfpm::vector<aggregate<int>> LL_p;
+ openfpm::vector<aggregate<int>> LL_p_ref;
+ openfpm::vector<aggregate<int>> LR_p;
+ openfpm::vector<aggregate<int>> LR_p_ref;
+
+ auto P = getV<0>(Particles);
+ auto V = getV<1>(Particles);
+ auto dV = getV<2>(Particles);
+ auto RHS = getV<3>(Particles);
+ auto Vtemp = getV<4>(Particles);
+ auto divV = getV<5>(Particles);
+
+
+ // Here fill up the boxes for particle detection.
+
+/*
+ Box<2, double> up({box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
+
+ Box<2, double> up_d({box.getLow(0) + spacing / 2.0, box.getHigh(1) - 3*spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
+ Box<2, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0});
+
+ Box<2, double> down_u({box.getLow(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getLow(1) + 3*spacing / 2.0});
+
+ Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
+ Box<2, double> left_r({box.getLow(0) + spacing / 2.0, box.getLow(1) + 3 *spacing / 2.0},
+ {box.getLow(0) + 3*spacing / 2.0, box.getHigh(1) - 3*spacing / 2.0});
+
+ Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
+ Box<2, double> right_l({box.getHigh(0) - 3*spacing / 2.0, box.getLow(1) + 3*spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getHigh(1) - 3*spacing / 2.0});
+*/
+ Box<2, double> up({box.getLow(0) + 3 * spacing / 2.0, box.getHigh(1) - spacing / 2.0},
+ {box.getHigh(0) - 3 * spacing / 2.0, box.getHigh(1) + spacing / 2.0});
+
+ Box<2, double> up_d({box.getLow(0) + 3 * spacing / 2.0, box.getHigh(1) - 3 * spacing / 2.0},
+ {box.getHigh(0) - 3 * spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
+ Box<2, double> down({box.getLow(0) + 3 * spacing / 2.0, box.getLow(1) - spacing / 2.0},
+ {box.getHigh(0) - 3 * spacing / 2.0, box.getLow(1) + spacing / 2.0});
+
+ Box<2, double> down_u({box.getLow(0) + 3 * spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getHigh(0) - 3 * spacing / 2.0, box.getLow(1) + 3 * spacing / 2.0});
+
+ Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + 3 * spacing / 2.0},
+ {box.getLow(0) + spacing / 2.0, box.getHigh(1) - 3 * spacing / 2.0});
+
+ Box<2, double> left_r({box.getLow(0) + spacing / 2.0, box.getLow(1) + 3 * spacing / 2.0},
+ {box.getLow(0) + 3 * spacing / 2.0, box.getHigh(1) - 3 * spacing / 2.0});
+
+ Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + 3 * spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - 3 * spacing / 2.0});
+
+ Box<2, double> right_l({box.getHigh(0) - 3 * spacing / 2.0, box.getLow(1) + 3 * spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getHigh(1) - 3 * spacing / 2.0});
+
+ Box<2, double> UL({box.getLow(0) - spacing / 2.0, box.getHigh(1) - 3 * spacing / 2.0},
+ {box.getLow(0) + 3 * spacing / 2.0, box.getHigh(1) + spacing / 2.0});
+
+ Box<2, double> UL_ref({box.getLow(0) + spacing / 2.0, box.getHigh(1) - 3 * spacing / 2.0},
+ {box.getLow(0) + 3 * spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
+
+ Box<2, double> UR({box.getHigh(0) - 3 * spacing / 2.0, box.getHigh(1) - 3 * spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
+
+ Box<2, double> UR_ref({box.getHigh(0) - 3 * spacing / 2.0, box.getHigh(1) - 3 * spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
+ Box<2, double> LL({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
+ {box.getLow(0) + 3 * spacing / 2.0, box.getLow(1) + 3 * spacing / 2.0});
+
+ Box<2, double> LL_ref({box.getLow(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getLow(0) + 3 * spacing / 2.0, box.getLow(1) + 3 * spacing / 2.0});
+
+ Box<2, double> LR({box.getHigh(0) - 3 * spacing / 2.0, box.getLow(1) - spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getLow(1) + 3 * spacing / 2.0});
+
+ Box<2, double> LR_ref({box.getHigh(0) - 3 * spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getLow(1) + 3 * spacing / 2.0});
+
+
+ openfpm::vector<Box<2, double>> boxes;
+ boxes.add(up);
+ boxes.add(up_d);
+ boxes.add(down);
+ boxes.add(down_u);
+ boxes.add(left);
+ boxes.add(left_r);
+ boxes.add(right);
+ boxes.add(right_l);
+ boxes.add(UL);
+ boxes.add(UL_ref);
+ boxes.add(UR_ref);
+ boxes.add(UR);
+ boxes.add(LL_ref);
+ boxes.add(LL);
+ boxes.add(LR_ref);
+ boxes.add(LR);
+ VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
+ vtk_box.add(boxes);
+ vtk_box.write("boxes.vtk");
+ auto it2 = Particles.getDomainIterator();
+
+ while (it2.isNext()) {
+ auto p = it2.get();
+ Point<2, double> xp = Particles.getPos(p);
+ if (xp[0] != 0 || xp[1] != 0) {
+ if (xp[0] != 0 + spacing || xp[1] != 0) {
+ FullBulk.add();
+ FullBulk.last().get<0>() = p.getKey();
+ }
+ }
+ if (up.isInside(xp) == true) {
+ up_p.add();
+ up_p.last().get<0>() = p.getKey();
+ Particles.getProp<1>(p)[0] = 1;
+ Particles.getProp<1>(p)[1] = 0;
+ } else if (down.isInside(xp) == true) {
+ dw_p.add();
+ dw_p.last().get<0>() = p.getKey();
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ } else if (left.isInside(xp) == true) {
+ l_p.add();
+ l_p.last().get<0>() = p.getKey();
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ } else if (right.isInside(xp) == true) {
+ r_p.add();
+ r_p.last().get<0>() = p.getKey();
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ } else if (up_d.isInside(xp) == true) {
+ up_p1.add();
+ up_p1.last().get<0>() = p.getKey();
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ } else if (down_u.isInside(xp) == true) {
+ dw_p1.add();
+ dw_p1.last().get<0>() = p.getKey();
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+
+ } else if (left_r.isInside(xp) == true) {
+ l_p1.add();
+ l_p1.last().get<0>() = p.getKey();
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ } else if (right_l.isInside(xp) == true) {
+ r_p1.add();
+ r_p1.last().get<0>() = p.getKey();
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ } else if (UL.isInside(xp) == true) {
+ if (UL_ref.isInside(xp) == true) {
+ UL_p_ref.add();
+ UL_p_ref.last().get<0>() = p.getKey();
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ } else {
+ UL_p.add();
+ UL_p.last().get<0>() = p.getKey();
+ }
+ } else if (UR.isInside(xp) == true) {
+ if (UR_ref.isInside(xp) == true) {
+ UR_p_ref.add();
+ UR_p_ref.last().get<0>() = p.getKey();
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ } else {
+ UR_p.add();
+ UR_p.last().get<0>() = p.getKey();
+ }
+ } else if (LR.isInside(xp) == true) {
+ if (LR_ref.isInside(xp) == true) {
+ LR_p_ref.add();
+ LR_p_ref.last().get<0>() = p.getKey();
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ } else {
+ LR_p.add();
+ LR_p.last().get<0>() = p.getKey();
+ }
+ } else if (LL.isInside(xp) == true) {
+ if (LL_ref.isInside(xp) == true) {
+ LL_p_ref.add();
+ LL_p_ref.last().get<0>() = p.getKey();
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ } else if (xp[0] == 0 && xp[1] == 0) {
+ ref_p.add();
+ ref_p.last().get<0>() = p.getKey();
+ } else {
+ LL_p.add();
+ LL_p.last().get<0>() = p.getKey();
+ }
+ } else {
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ Particles.getProp<0>(p) = 0.0;
+ }
+
+ ++it2;
+ }
+ Particles.getProp<1>(959)[1] = 0;
+ Particles.getProp<1>(931)[0] = 1;
+
+ /* for(int j=0;j<up_p1.size();j++)
+ {
+ auto p1=up_p1.get<0>(j);
+ Point<2, double> xp = Particles.getPos(p1);
+ Particles.getProp<3>(p1) = xp[0];
+ }
+ for(int j=0;j<l_p1.size();j++)
+ {
+ auto p1=l_p1.get<0>(j);
+ Point<2, double> xp = Particles.getPos(p1);
+ Particles.getProp<3>(p1) = xp[1];
+ }
+ for(int j=0;j<r_p1.size();j++)
+ {
+ auto p1=r_p1.get<0>(j);
+ Point<2, double> xp = Particles.getPos(p1);
+ Particles.getProp<3>(p1) = xp[1];
+ }
+ for(int j=0;j<dw_p1.size();j++)
+ {
+ auto p1=dw_p1.get<0>(j);
+ Point<2, double> xp = Particles.getPos(p1);
+ Particles.getProp<3>(p1)= xp[0];
+ }*/
+ Particles.write_frame("Re1000-1e-4-Lid", 0);
+/*
+ for(int j=0;j<up_p1.size();j++)
+ { if(j==0)
+ { auto p=up_p.get<0>(j);
+ auto p1=up_p1.get<0>(j);
+ auto p2=l_p.get<0>(l_p.size()-1);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2) = Particles.getProp<3>(p1);
+ }
+ auto p=up_p.get<0>(j+1);
+ auto p1=up_p1.get<0>(j);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ if(j==up_p1.size()-1)
+ { auto p=up_p.get<0>(j+2);
+ auto p1=up_p1.get<0>(j);
+ auto p2=r_p.get<0>(r_p.size()-1);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2) = Particles.getProp<3>(p1);
+ }
+ }
+ for(int j=0;j<l_p1.size();j++)
+ {
+ auto p=l_p.get<0>(j+1);
+ auto p1=l_p1.get<0>(j);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ }
+ for(int j=0;j<r_p1.size();j++)
+ {
+ auto p=r_p.get<0>(j+1);
+ auto p1=r_p1.get<0>(j);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ }
+ for(int j=0;j<dw_p1.size();j++)
+ { auto p=dw_p.get<0>(j);
+ auto p1=dw_p1.get<0>(j);
+ Particles.getProp<3>(p)= Particles.getProp<3>(p1);
+ if(j==0)
+ {
+ auto p = ref_p.get<0>(j);
+ auto p2=l_p.get<0>(0);
+ Particles.getProp<3>(p)= Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2)= Particles.getProp<3>(p1);
+ }
+ if(j==dw_p1.size()-1)
+ { auto p=dw_p.get<0>(j+1);
+ auto p2=r_p.get<0>(0);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2)= Particles.getProp<3>(p1);
+ }
+ }*/
+ Particles.write_frame("Re1000-1e-4-Lid", 1);
+ Derivative_x Dx(Particles, 2, rCut, 2);
+ Derivative_y Dy(Particles, 2, rCut, 2);
+ Gradient Grad(Particles, 2, rCut, 2);
+ Laplacian Lap(Particles, 2, rCut, 2);
+ Advection Adv(Particles, 2, rCut, 2);
+ Divergence Div(Particles, 2, rCut, 2);
+ double dt = 5e-4;
+ int n = 3;
+ double nu = 1e-2;
+/* divV=Div(V);
+ DCPSE_scheme<equations2,decltype(Particles)> Solver(2*rCut, Particles);
+ auto Pressure_Poisson = Lap(P);
+ auto D_y=Dy(P);
+ auto D_x=Dx(P);
+ Solver.impose(Pressure_Poisson,bulk,prop_id<5>());
+ Solver.impose(D_y, up_p,0);
+ Solver.impose(D_y, dw_p,0);
+ Solver.impose(D_x, l_p,0);
+ Solver.impose(D_x, r_p, 0);
+ Solver.impose(P, ref_p, 0);
+ Solver.solve(P);
+
+ std::cout<<"Poisson Solved"<<std::endl;
+ V=V-Grad(P);
+ Particles.write_frame("Re1000-1e-4-Lid",0);
+ std::cout<<"Init Done"<<std::endl;
+ return;*/
+ for (int i = 2; i <= n; i++) {
+ dV = dt * (nu * Lap(V) - Adv(V, V));
+ RHS = 1 / dt * Div(dV);
+ std::cout << "RHS Done" << std::endl;
+
+ //Neumann Copy
+/*
+ for(int j=0;j<up_p1.size();j++)
+ { if(j==0)
+ { auto p=up_p.get<0>(j);
+ auto p1=up_p1.get<0>(j);
+ auto p2=l_p.get<0>(l_p.size()-1);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2) = Particles.getProp<3>(p1);
+ }
+ auto p=up_p.get<0>(j+1);
+ auto p1=up_p1.get<0>(j);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ if(j==up_p1.size()-1)
+ { auto p=up_p.get<0>(j+2);
+ auto p1=up_p1.get<0>(j);
+ auto p2=r_p.get<0>(r_p.size()-1);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2) = Particles.getProp<3>(p1);
+ }
+ }
+ for(int j=0;j<l_p1.size();j++)
+ {
+ auto p=l_p.get<0>(j+1);
+ auto p1=l_p1.get<0>(j);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ }
+ for(int j=0;j<r_p1.size();j++)
+ {
+ auto p=r_p.get<0>(j+1);
+ auto p1=r_p1.get<0>(j);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ }
+ for(int j=0;j<dw_p1.size();j++)
+ { auto p=dw_p.get<0>(j);
+ auto p1=dw_p1.get<0>(j);
+ Particles.getProp<3>(p)= Particles.getProp<3>(p1);
+ if(j==0)
+ {
+ auto p = ref_p.get<0>(j);
+ auto p2=l_p.get<0>(0);
+ Particles.getProp<3>(p)= Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2)= Particles.getProp<3>(p1);
+ }
+ if(j==dw_p1.size()-1)
+ { auto p=dw_p.get<0>(j+1);
+ auto p2=r_p.get<0>(0);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2)= Particles.getProp<3>(p1);
+ }
+ }*/
+
+ vector_dist_expression_op<void, void, VECT_COPY_N_TO_N> vcp_up(up_p, up_p1);
+ vector_dist_expression_op<void, void, VECT_COPY_N_TO_N> vcp_l(l_p, l_p1);
+ vector_dist_expression_op<void, void, VECT_COPY_N_TO_N> vcp_r(r_p, r_p1);
+ vector_dist_expression_op<void, void, VECT_COPY_N_TO_N> vcp_dw(dw_p, dw_p1);
+
+ vector_dist_expression_op<void, void, VECT_COPY_1_TO_N> vcp_ul(UL_p, UL_p_ref.get<0>(0));
+ vector_dist_expression_op<void, void, VECT_COPY_1_TO_N> vcp_ur(UR_p, UR_p_ref.get<0>(0));
+ vector_dist_expression_op<void, void, VECT_COPY_1_TO_N> vcp_dl(LL_p, LL_p_ref.get<0>(0));
+ vector_dist_expression_op<void, void, VECT_COPY_1_TO_N> vcp_dr(LR_p, LR_p_ref.get<0>(0));
+
+ DCPSE_scheme<equations2, decltype(Particles)> Solver(2 * rCut, Particles);
+ auto Pressure_Poisson = Lap(P);
+ Solver.impose(Pressure_Poisson, bulk, prop_id<3>());
+ Solver.impose(vcp_up, up_p, 0);
+ Solver.impose(vcp_l, l_p, 0);
+ Solver.impose(vcp_r, r_p, 0);
+ Solver.impose(vcp_dw, dw_p, 0);
+ Solver.impose(vcp_ul, UL_p, 0);
+ Solver.impose(vcp_ur, UR_p, 0);
+ Solver.impose(vcp_dl, LL_p, 0);
+ Solver.impose(vcp_dr, LR_p, 0);
+/* Solver.impose(P, up_p, up_p1);
+ Solver.impose(P, dw_p,prop_id<3>());
+ Solver.impose(P, l_p,prop_id<3>());
+ Solver.impose(P, r_p, prop_id<3>());*/
+ Solver.impose(P, ref_p, 0);
+ Solver.solve(P);
+ std::cout << "Poisson Solved" << std::endl;
+ Vtemp = V + dV - dt * Grad(P);
+ V = Vtemp;
+ divV = Div(V);
+ for (int j = 0; j < up_p.size(); j++) {
+ auto p = up_p.get<0>(j);
+ Particles.getProp<1>(p)[0] = 1;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ for (int j = 0; j < l_p.size(); j++) {
+ auto p = l_p.get<0>(j);
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ for (int j = 0; j < r_p.size(); j++) {
+ auto p = r_p.get<0>(j);
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ for (int j = 0; j < dw_p.size(); j++) {
+ auto p = dw_p.get<0>(j);
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ Particles.getProp<1>(0)[0] = 0;
+ Particles.getProp<1>(0)[1] = 0;
+ Particles.getProp<1>(959)[1] = 0;
+ Particles.getProp<1>(931)[0] = 1;
+ std::cout << "Velocity updated" << std::endl;
+ Particles.write_frame("Re1000-1e-4-Lid", i);
+ std::cout << i << std::endl;
+ }
+ }
+/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+
+
+
+
+
+ BOOST_AUTO_TEST_CASE(dcpse_Lid_copy_RHS) {
+ const size_t sz[2] = {31,31};
+ Box<2, double> box({0, 0}, {1,1});
+ size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
+ double spacing = box.getHigh(0) / (sz[0] - 1);
+ Ghost<2, double> ghost(spacing * 3);
+ double rCut = 2.0 * spacing;
+ // P V Dv RHS Vtemp
+ vector_dist<2, double, aggregate<double,VectorS<2, double>,VectorS<2, double>,double,VectorS<2, double>,double>> Particles(0, box, bc, ghost);
+ auto it = Particles.getGridIterator(sz);
+ while (it.isNext()) {
+ Particles.add();
+ auto key = it.get();
+ double x = key.get(0) * it.getSpacing(0);
+ Particles.getLastPos()[0] = x;
+ double y = key.get(1) * it.getSpacing(1);
+ Particles.getLastPos()[1] = y;
+ ++it;
+ }
+
+ Particles.map();
+ Particles.ghost_get<0>();
+
+ openfpm::vector<aggregate<int>> bulk;
+ openfpm::vector<aggregate<int>> up_p;
+ openfpm::vector<aggregate<int>> dw_p;
+ openfpm::vector<aggregate<int>> l_p;
+ openfpm::vector<aggregate<int>> r_p;
+ openfpm::vector<aggregate<int>> up_p1;
+ openfpm::vector<aggregate<int>> dw_p1;
+ openfpm::vector<aggregate<int>> l_p1;
+ openfpm::vector<aggregate<int>> r_p1;
+ openfpm::vector<aggregate<int>> ref_p;
+ openfpm::vector<aggregate<int>> ref_p2;
+ openfpm::vector<aggregate<int>> FullBulk;
+ /*openfpm::vector<aggregate<int>> UL_p;
+ openfpm::vector<aggregate<int>> UL_p_ref;
+ openfpm::vector<aggregate<int>> UR_p;
+ openfpm::vector<aggregate<int>> UR_p_ref;
+ openfpm::vector<aggregate<int>> LL_p;
+ openfpm::vector<aggregate<int>> LL_p_ref;
+ openfpm::vector<aggregate<int>> LR_p;
+ openfpm::vector<aggregate<int>> LR_p_ref;*/
+
+ auto P = getV<0>(Particles);
+ auto V = getV<1>(Particles);
+ auto dV = getV<2>(Particles);
+ auto RHS = getV<3>(Particles);
+ auto Vtemp = getV<4>(Particles);
+ auto divV = getV<5>(Particles);
+
+
+ // Here fill up the boxes for particle detection.
+
+ Box<2, double> up({box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
+
+ Box<2, double> up_d({box.getLow(0) + spacing / 2.0, box.getHigh(1) - 3*spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
+ Box<2, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0});
+
+ Box<2, double> down_u({box.getLow(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getLow(1) + 3*spacing / 2.0});
+
+ Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
+ Box<2, double> left_r({box.getLow(0) + spacing / 2.0, box.getLow(1) + 3 *spacing / 2.0},
+ {box.getLow(0) + 3*spacing / 2.0, box.getHigh(1) - 3*spacing / 2.0});
+
+ Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
+ Box<2, double> right_l({box.getHigh(0) - 3*spacing / 2.0, box.getLow(1) + 3*spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getHigh(1) - 3*spacing / 2.0});
+
+
+ openfpm::vector<Box<2, double>> boxes;
+ boxes.add(up);
+ boxes.add(up_d);
+ boxes.add(down);
+ boxes.add(down_u);
+ boxes.add(left);
+ boxes.add(left_r);
+ boxes.add(right);
+ boxes.add(right_l);
+ VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
+ vtk_box.add(boxes);
+ vtk_box.write("boxes.vtk");
+ auto it2 = Particles.getDomainIterator();
+
+ while (it2.isNext()) {
+ auto p = it2.get();
+ Point<2, double> xp = Particles.getPos(p);
+ if (xp[0]!=0 || xp[1]!=0){
+ FullBulk.add();
+ FullBulk.last().get<0>() = p.getKey();
+ }
+ if (up.isInside(xp) == true) {
+ up_p.add();
+ up_p.last().get<0>() = p.getKey();
+ Particles.getProp<1>(p)[0] = 1;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ else if (down.isInside(xp) == true) {
+ if (xp[0]==0 && xp[1]==0) {
+ ref_p.add();
+ ref_p.last().get<0>() = p.getKey();
+ }
+ else {
+ dw_p.add();
+ dw_p.last().get<0>() = p.getKey();
+ }
+ }
+ else if (left.isInside(xp) == true) {
+ l_p.add();
+ l_p.last().get<0>() = p.getKey();
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ else if (right.isInside(xp) == true) {
+ r_p.add();
+ r_p.last().get<0>() = p.getKey();
+ }
+ else if (up_d.isInside(xp) == true) {
+ up_p1.add();
+ up_p1.last().get<0>() = p.getKey();
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ }
+ else if (down_u.isInside(xp) == true) {
+ dw_p1.add();
+ dw_p1.last().get<0>() = p.getKey();
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
- BOOST_AUTO_TEST_CASE(dcpse_Lid2) {
+ }
+ else if (left_r.isInside(xp) == true) {
+ l_p1.add();
+ l_p1.last().get<0>() = p.getKey();
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ }
+ else if (right_l.isInside(xp) == true) {
+ r_p1.add();
+ r_p1.last().get<0>() = p.getKey();
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ }
+ else {
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ Particles.getProp<0>(p) = 0.0;
+ }
+
+ ++it2;
+ }
+
+ for(int j=0;j<up_p1.size();j++)
+ {
+ auto p1=up_p1.get<0>(j);
+ Point<2, double> xp = Particles.getPos(p1);
+ Particles.getProp<3>(p1) = xp[0];
+ }
+ for(int j=0;j<l_p1.size();j++)
+ {
+ auto p1=l_p1.get<0>(j);
+ Point<2, double> xp = Particles.getPos(p1);
+ Particles.getProp<3>(p1) = xp[1];
+ }
+ for(int j=0;j<r_p1.size();j++)
+ {
+ auto p1=r_p1.get<0>(j);
+ Point<2, double> xp = Particles.getPos(p1);
+ Particles.getProp<3>(p1) = xp[1];
+ }
+ for(int j=0;j<dw_p1.size();j++)
+ {
+ auto p1=dw_p1.get<0>(j);
+ Point<2, double> xp = Particles.getPos(p1);
+ Particles.getProp<3>(p1)= xp[0];
+ }
+ Particles.write_frame("Re1000-1e-4-Lid",0);
+
+ for(int j=0;j<up_p1.size();j++)
+ { if(j==0)
+ { auto p=up_p.get<0>(j);
+ auto p1=up_p1.get<0>(j);
+ auto p2=l_p.get<0>(l_p.size()-1);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2) = Particles.getProp<3>(p1);
+ }
+ auto p=up_p.get<0>(j+1);
+ auto p1=up_p1.get<0>(j);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ if(j==up_p1.size()-1)
+ { auto p=up_p.get<0>(j+2);
+ auto p1=up_p1.get<0>(j);
+ auto p2=r_p.get<0>(r_p.size()-1);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2) = Particles.getProp<3>(p1);
+ }
+ }
+ for(int j=0;j<l_p1.size();j++)
+ {
+ auto p=l_p.get<0>(j+1);
+ auto p1=l_p1.get<0>(j);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ }
+ for(int j=0;j<r_p1.size();j++)
+ {
+ auto p=r_p.get<0>(j+1);
+ auto p1=r_p1.get<0>(j);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ }
+ for(int j=0;j<dw_p1.size();j++)
+ { auto p=dw_p.get<0>(j);
+ auto p1=dw_p1.get<0>(j);
+ Particles.getProp<3>(p)= Particles.getProp<3>(p1);
+ if(j==0)
+ {
+ auto p = ref_p.get<0>(j);
+ auto p2=l_p.get<0>(0);
+ Particles.getProp<3>(p)= Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2)= Particles.getProp<3>(p1);
+ }
+ if(j==dw_p1.size()-1)
+ { auto p=dw_p.get<0>(j+1);
+ auto p2=r_p.get<0>(0);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2)= Particles.getProp<3>(p1);
+ }
+ }
+ Particles.write_frame("Re1000-1e-4-Lid",1);
+ Derivative_x Dx(Particles, 2, rCut,2);
+ Derivative_y Dy(Particles, 2, rCut,2);
+ Gradient Grad(Particles, 2, rCut,2);
+ Laplacian Lap(Particles, 3, rCut, 3);
+ Advection Adv(Particles, 2, rCut, 2);
+ Divergence Div(Particles, 2, rCut, 2);
+ double dt=5e-4;
+ int n=3;
+ double nu=1e-2;
+/* divV=Div(V);
+ DCPSE_scheme<equations2,decltype(Particles)> Solver(2*rCut, Particles);
+ auto Pressure_Poisson = Lap(P);
+ auto D_y=Dy(P);
+ auto D_x=Dx(P);
+ Solver.impose(Pressure_Poisson,bulk,prop_id<5>());
+ Solver.impose(D_y, up_p,0);
+ Solver.impose(D_y, dw_p,0);
+ Solver.impose(D_x, l_p,0);
+ Solver.impose(D_x, r_p, 0);
+ Solver.impose(P, ref_p, 0);
+ Solver.solve(P);
+
+ std::cout<<"Poisson Solved"<<std::endl;
+ V=V-Grad(P);
+ Particles.write_frame("Re1000-1e-4-Lid",0);
+ std::cout<<"Init Done"<<std::endl;
+ return;*/
+ for(int i=2; i<=n ;i++)
+ { dV=V+dt*(nu*Lap(V)-Adv(V,V));
+ RHS=1/dt*Div(dV);
+ std::cout<<"RHS Done"<<std::endl;
+
+ /* //Neumann Copy
+
+ for(int j=0;j<up_p1.size();j++)
+ { if(j==0)
+ { auto p=up_p.get<0>(j);
+ auto p1=up_p1.get<0>(j);
+ auto p2=l_p.get<0>(l_p.size()-1);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2) = Particles.getProp<3>(p1);
+ }
+ auto p=up_p.get<0>(j+1);
+ auto p1=up_p1.get<0>(j);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ if(j==up_p1.size()-1)
+ { auto p=up_p.get<0>(j+2);
+ auto p1=up_p1.get<0>(j);
+ auto p2=r_p.get<0>(r_p.size()-1);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2) = Particles.getProp<3>(p1);
+ }
+ }
+ for(int j=0;j<l_p1.size();j++)
+ {
+ auto p=l_p.get<0>(j+1);
+ auto p1=l_p1.get<0>(j);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ }
+ for(int j=0;j<r_p1.size();j++)
+ {
+ auto p=r_p.get<0>(j+1);
+ auto p1=r_p1.get<0>(j);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ }
+ for(int j=0;j<dw_p1.size();j++)
+ { auto p=dw_p.get<0>(j);
+ auto p1=dw_p1.get<0>(j);
+ Particles.getProp<3>(p)= Particles.getProp<3>(p1);
+ if(j==0)
+ {
+ auto p = ref_p.get<0>(j);
+ auto p2=l_p.get<0>(0);
+ Particles.getProp<3>(p)= Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2)= Particles.getProp<3>(p1);
+ }
+ if(j==dw_p1.size()-1)
+ { auto p=dw_p.get<0>(j+1);
+ auto p2=r_p.get<0>(0);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2)= Particles.getProp<3>(p1);
+ }
+ } for(int j=0;j<up_p1.size();j++)
+ { if(j==0)
+ { auto p=up_p.get<0>(j);
+ auto p1=up_p1.get<0>(j);
+ auto p2=l_p.get<0>(l_p.size()-1);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2) = Particles.getProp<3>(p1);
+ }
+ auto p=up_p.get<0>(j+1);
+ auto p1=up_p1.get<0>(j);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ if(j==up_p1.size()-1)
+ { auto p=up_p.get<0>(j+2);
+ auto p1=up_p1.get<0>(j);
+ auto p2=r_p.get<0>(r_p.size()-1);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2) = Particles.getProp<3>(p1);
+ }
+ }
+ for(int j=0;j<l_p1.size();j++)
+ {
+ auto p=l_p.get<0>(j+1);
+ auto p1=l_p1.get<0>(j);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ }
+ for(int j=0;j<r_p1.size();j++)
+ {
+ auto p=r_p.get<0>(j+1);
+ auto p1=r_p1.get<0>(j);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ }
+ for(int j=0;j<dw_p1.size();j++)
+ { auto p=dw_p.get<0>(j);
+ auto p1=dw_p1.get<0>(j);
+ Particles.getProp<3>(p)= Particles.getProp<3>(p1);
+ if(j==0)
+ {
+ auto p = ref_p.get<0>(j);
+ auto p2=l_p.get<0>(0);
+ Particles.getProp<3>(p)= Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2)= Particles.getProp<3>(p1);
+ }
+ if(j==dw_p1.size()-1)
+ { auto p=dw_p.get<0>(j+1);
+ auto p2=r_p.get<0>(0);
+ Particles.getProp<3>(p) = Particles.getProp<3>(p1);
+ Particles.getProp<3>(p2)= Particles.getProp<3>(p1);
+ }
+ }
+*/
+ DCPSE_scheme<equations2,decltype(Particles)> Solver(2*rCut, Particles);
+ auto Pressure_Poisson = Lap(P);
+ auto D_y=Dy(P);
+ auto D_x=Dx(P);
+ Solver.impose(Pressure_Poisson,FullBulk,prop_id<3>());
+/* Solver.impose(P, up_p,0);
+ Solver.impose(P, dw_p,0);
+ Solver.impose(P, l_p,0);
+ Solver.impose(P, r_p, 0);*/
+ Solver.impose(P, ref_p, 0);
+ Solver.solve(P);
+ std::cout<<"Poisson Solved"<<std::endl;
+ Vtemp=dV-dt*Grad(P);
+ V=Vtemp;
+ for(int j=0;j<up_p.size();j++)
+ { auto p=up_p.get<0>(j);
+ Particles.getProp<1>(p)[0] = 1;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ for(int j=0;j<l_p.size();j++)
+ { auto p=l_p.get<0>(j);
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ for(int j=0;j<r_p.size();j++)
+ { auto p=r_p.get<0>(j);
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ for(int j=0;j<dw_p.size();j++)
+ { auto p=dw_p.get<0>(j);
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ Particles.getProp<1>(0)[0] = 0;
+ Particles.getProp<1>(0)[1] = 0;
+ divV=Div(V);
+ std::cout<<"Velocity updated"<<std::endl;
+ Particles.write_frame("Re1000-1e-4-Lid",i);
+ std::cout<<i<<std::endl;
+ }
+ }
+
+
+
+
+
+
+
+
+
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+
+
+
+
+
+ BOOST_AUTO_TEST_CASE(dcpse_Lid_normal) {
const size_t sz[2] = {31,31};
Box<2, double> box({0, 0}, {1,1});
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing = box.getHigh(0) / (sz[0] - 1);
Ghost<2, double> ghost(spacing * 3);
- double rCut = 2.0 * spacing;
- // P V Dv RHS Vtemp
- vector_dist<2, double, aggregate<double,VectorS<2, double>,VectorS<2, double>,double,VectorS<2, double>,double>> Particles(0, box, bc, ghost);
+ double rCut = 0.7 * spacing;
+ // P V Dv RHS Vtemp Proj_lap
+ vector_dist<2, double, aggregate<double,VectorS<2, double>,VectorS<2, double>,double,VectorS<2, double>,double,double>> Particles(0, box, bc, ghost);
auto it = Particles.getGridIterator(sz);
while (it.isNext()) {
Particles.add();
@@ -77,9 +1866,6 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
openfpm::vector<aggregate<int>> dw_p;
openfpm::vector<aggregate<int>> l_p;
openfpm::vector<aggregate<int>> r_p;
- openfpm::vector<aggregate<int>> ref_p;
- openfpm::vector<aggregate<int>> ref_p2;
- openfpm::vector<aggregate<int>> FullBulk;
auto P = getV<0>(Particles);
@@ -88,6 +1874,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
auto RHS = getV<3>(Particles);
auto Vtemp = getV<4>(Particles);
auto divV = getV<5>(Particles);
+ auto Proj_lap=getV<6>(Particles);
// Here fill up the boxes for particle detection.
@@ -98,10 +1885,10 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
Box<2, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
{box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0});
- Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
+ Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
{box.getLow(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
- Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
+ Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
{box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
openfpm::vector<Box<2, double>> boxes;
@@ -112,17 +1899,205 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
vtk_box.add(boxes);
vtk_box.write("boxes.vtk");
-
auto it2 = Particles.getDomainIterator();
while (it2.isNext()) {
auto p = it2.get();
Point<2, double> xp = Particles.getPos(p);
- if (xp[0]!=0 || xp[1]!=0) {
- if (xp[0]!=0+spacing || xp[1]!=0) {
- FullBulk.add();
- FullBulk.last().get<0>() = p.getKey();
+ Particles.getProp<3>(p) =2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ if (up.isInside(xp) == true) {
+ up_p.add();
+ up_p.last().get<0>() = p.getKey();
+ Particles.getProp<1>(p)[0] = 1;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ else if (down.isInside(xp) == true) {
+ dw_p.add();
+ dw_p.last().get<0>() = p.getKey();
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ else if (left.isInside(xp) == true) {
+ l_p.add();
+ l_p.last().get<0>() = p.getKey();
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ else if (right.isInside(xp) == true) {
+ r_p.add();
+ r_p.last().get<0>() = p.getKey();
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ else {
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ }
+ ++it2;
+ }
+
+ Derivative_x Dx(Particles, 2, rCut,1.9);
+ Derivative_y Dy(Particles, 2, rCut,1.9);
+ Gradient Grad(Particles, 2, rCut,1.9);
+ Laplacian Lap(Particles, 2, rCut, 1.9);
+ Advection Adv(Particles, 2, rCut, 1.9);
+ Divergence Div(Particles, 2, rCut, 1.9);
+ double dt=3e-3;
+ int n=50;
+ double nu=1e-2;
+ dV=dt*(nu*Lap(V)-Adv(V,V));
+ DCPSE_scheme<equations2,decltype(Particles)> Solver(2 * rCut, Particles,options_solver::LAGRANGE_MULTIPLIER);
+ auto Pressure_Poisson = Lap(P);
+ auto D_y=Dy(P);
+ auto D_x=Dx(P);
+ Solver.impose(Pressure_Poisson,bulk,prop_id<3>());
+ Solver.impose(D_y, up_p,0);
+ Solver.impose(D_x, r_p, 0);
+ Solver.impose(-D_y, dw_p,0);
+ Solver.impose(-D_x, l_p,0);
+ Solver.solve(P);
+ std::cout << "Poisson Solved" << std::endl;
+ Vtemp = V + (dV - dt*Grad(P));
+ V = Vtemp;
+ divV=Div(V);
+ Particles.write_frame("Re1000-1e-4-Lid",0);
+ for(int i=25; i<=n ;i++)
+ { dV=dt*(nu*Lap(V)-Adv(V,V));
+ //dV=Lap(V);
+ //dV=Adv(V,V);
+ // Vtemp=dV;
+ //dV=Vtemp - Adv(V,V);
+ //Vtemp=dV;
+ //dV=Vtemp;
+ //Lap.DrawKernel<5>(Particles,837);
+ //Lap.DrawKernel<5>(Particles,867);
+ if(i>0) {
+ RHS=1/dt*Div(dV);
+ std::cout<<"RHS Done"<<std::endl;
+ DCPSE_scheme<equations2,decltype(Particles)> Solver(2 * rCut, Particles,options_solver::LAGRANGE_MULTIPLIER);
+ auto Pressure_Poisson = Lap(P);
+ auto D_y=Dy(P);
+ auto D_x=Dx(P);
+ Solver.impose(Pressure_Poisson,bulk,prop_id<3>());
+ Solver.impose(D_y, up_p,0);
+ Solver.impose(D_x, r_p, 0);
+ Solver.impose(-D_y, dw_p,0);
+ Solver.impose(-D_x, l_p,0);
+ Solver.solve(P);
+ std::cout << "Poisson Solved" << std::endl;
+ Vtemp = V + (dV - dt*Grad(P));
+ V = Vtemp;
}
+ else {
+ Vtemp = V + (dV);
+ V = Vtemp;
+ }
+ divV = Div(V);
+ for(int j=0;j<up_p.size();j++)
+ { auto p=up_p.get<0>(j);
+ Particles.getProp<1>(p)[0] = 1;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ for(int j=0;j<l_p.size();j++)
+ { auto p=l_p.get<0>(j);
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ for(int j=0;j<r_p.size();j++)
+ { auto p=r_p.get<0>(j);
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ for(int j=0;j<dw_p.size();j++)
+ { auto p=dw_p.get<0>(j);
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ }
+ std::cout<<"Velocity updated"<<std::endl;
+ Particles.write_frame("Re1000-1e-4-Lid",i);
+ std::cout<<i<<std::endl;
+ }
+ }
+
+
+
+
+
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+
+
+
+
+/*
+
+ BOOST_AUTO_TEST_CASE(dcpse_Lid_RotProj) {
+ const size_t sz[2] = {31,31};
+ Box<2, double> box({0, 0}, {1,1});
+ size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
+ double spacing = box.getHigh(0) / (sz[0] - 1);
+ Ghost<2, double> ghost(spacing * 3);
+ double rCut = 2.0 * spacing;
+ // P V V1 V2 Dv RHS Vtemp Proj_lap
+ vector_dist<2, double, aggregate<double,VectorS<2, double>,VectorS<2, double>,VectorS<2, double>,VectorS<2, double>,double,VectorS<2, double>,double,double>> Particles(0, box, bc, ghost);
+ auto it = Particles.getGridIterator(sz);
+ while (it.isNext()) {
+ Particles.add();
+ auto key = it.get();
+ double x = key.get(0) * it.getSpacing(0);
+ Particles.getLastPos()[0] = x;
+ double y = key.get(1) * it.getSpacing(1);
+ Particles.getLastPos()[1] = y;
+ ++it;
+ }
+
+ Particles.map();
+ Particles.ghost_get<0>();
+
+ openfpm::vector<aggregate<int>> bulk;
+ openfpm::vector<aggregate<int>> up_p;
+ openfpm::vector<aggregate<int>> dw_p;
+ openfpm::vector<aggregate<int>> l_p;
+ openfpm::vector<aggregate<int>> r_p;
+ openfpm::vector<aggregate<int>> up_p1;
+ openfpm::vector<aggregate<int>> dw_p1;
+ openfpm::vector<aggregate<int>> l_p1;
+ openfpm::vector<aggregate<int>> r_p1;
+ openfpm::vector<aggregate<int>> ref_p;
+ openfpm::vector<aggregate<int>> ref_p2;
+ openfpm::vector<aggregate<int>> FullBulk;
+
+
+ // Here fill up the boxes for particle detection.
+
+ Box<2, double> up({box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
+
+ Box<2, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0});
+
+ Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
+ Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
+ openfpm::vector<Box<2, double>> boxes;
+ boxes.add(up);
+ boxes.add(down);
+ boxes.add(left);
+ boxes.add(right);
+ VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
+ vtk_box.add(boxes);
+ vtk_box.write("boxes.vtk");
+ auto it2 = Particles.getDomainIterator();
+
+ while (it2.isNext()) {
+ auto p = it2.get();
+ Point<2, double> xp = Particles.getPos(p);
+ if (xp[0]!=0 || xp[1]!=0){
+ FullBulk.add();
+ FullBulk.last().get<0>() = p.getKey();
}
if (up.isInside(xp) == true) {
up_p.add();
@@ -134,20 +2109,10 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
if (xp[0]==0 && xp[1]==0) {
ref_p.add();
ref_p.last().get<0>() = p.getKey();
- Particles.getProp<1>(p)[0] = 0;
- Particles.getProp<1>(p)[1] = 0;
- }
- else if (xp[0]==0+spacing && xp[1]==0) {
- ref_p.add();
- ref_p.last().get<0>() = p.getKey();
- Particles.getProp<1>(p)[0] = 0;
- Particles.getProp<1>(p)[1] = 0;
}
- else{
- dw_p.add();
- dw_p.last().get<0>() = p.getKey();
- Particles.getProp<1>(p)[0] = 0;
- Particles.getProp<1>(p)[1] = 0;
+ else {
+ dw_p.add();
+ dw_p.last().get<0>() = p.getKey();
}
}
else if (left.isInside(xp) == true) {
@@ -155,20 +2120,28 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
l_p.last().get<0>() = p.getKey();
Particles.getProp<1>(p)[0] = 0;
Particles.getProp<1>(p)[1] = 0;
- } else if (right.isInside(xp) == true) {
+ }
+ else if (right.isInside(xp) == true) {
r_p.add();
r_p.last().get<0>() = p.getKey();
- Particles.getProp<1>(p)[0] = 0;
- Particles.getProp<1>(p)[1] = 0;
}
else {
- bulk.add();
- bulk.last().get<0>() = p.getKey();
- Particles.getProp<0>(p) = 0.0;
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ Particles.getProp<0>(p) = 0.0;
}
++it2;
}
+ auto P = getV<0>(Particles);
+ auto V = getV<1>(Particles);
+ auto V1 = getV<2>(Particles);
+ auto V2 = getV<3>(Particles);
+ auto dV = getV<4>(Particles);
+ auto RHS = getV<5>(Particles);
+ auto Vtemp = getV<6>(Particles);
+ auto divV = getV<7>(Particles);
+ auto Proj_lap=getV<8>(Particles);
Derivative_x Dx(Particles, 2, rCut,2);
Derivative_y Dy(Particles, 2, rCut,2);
@@ -176,47 +2149,32 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
Laplacian Lap(Particles, 2, rCut, 2);
Advection Adv(Particles, 2, rCut, 2);
Divergence Div(Particles, 2, rCut, 2);
- double dt=5e-4;
- int n=30;
+ double dt=1e-3;
+ int n=3;
double nu=1e-2;
-/* divV=Div(V);
- DCPSE_scheme<equations2,decltype(Particles)> Solver(2*rCut, Particles);
- auto Pressure_Poisson = Lap(P);
- auto D_y=Dy(P);
- auto D_x=Dx(P);
- Solver.impose(Pressure_Poisson,bulk,prop_id<5>());
- Solver.impose(D_y, up_p,0);
- Solver.impose(D_y, dw_p,0);
- Solver.impose(D_x, l_p,0);
- Solver.impose(D_x, r_p, 0);
- Solver.impose(P, ref_p, 0);
- Solver.solve(P);
-
- std::cout<<"Poisson Solved"<<std::endl;
- V=V-Grad(P);
- Particles.write_frame("Re1000-1e-4-Lid",0);
- std::cout<<"Init Done"<<std::endl;
- return;*/
Particles.write_frame("Re1000-1e-4-Lid",0);
- for(int i =1; i<=n ;i++)
- { dV=dt*(nu*Lap(V)-Adv(V,V));
+ for(int i=1; i<=n ;i++)
+ { dV=2*dt*(nu*Lap(V)-Adv(V,V)-Grad(P));
+
RHS=1/dt*Div(dV);
std::cout<<"RHS Done"<<std::endl;
+
DCPSE_scheme<equations2,decltype(Particles)> Solver(2*rCut, Particles);
auto Pressure_Poisson = Lap(P);
- auto D_x = Dx(P);
- auto D_y = Dy(P);
- Solver.impose(Pressure_Poisson, bulk, prop_id<3>());
- Solver.impose(D_y, up_p, 0);
- Solver.impose(-D_y, dw_p,0);
- Solver.impose(-D_x, l_p,0);
- Solver.impose(D_x, r_p, 0);
+ auto D_y=Dy(P);
+ auto D_x=Dx(P);
+ Solver.impose(Pressure_Poisson,bulk,prop_id<3>());
+ Solver.impose(D_y, up_p,prop_id<3>());
+ Solver.impose(D_y, dw_p,prop_id<3>());
+ Solver.impose(D_x, l_p,prop_id<3>());
+ Solver.impose(D_x, r_p, prop_id<3>());
Solver.impose(P, ref_p, 0);
Solver.solve(P);
std::cout<<"Poisson Solved"<<std::endl;
- Vtemp=V+dV-dt*Grad(P);
+ Proj_lap=Dx(P)-RHS;
+ Vtemp=dV-dt*Grad(P);
V=Vtemp;
- divV=Div(V);
+ //divV=Div(V);
for(int j=0;j<up_p.size();j++)
{ auto p=up_p.get<0>(j);
Particles.getProp<1>(p)[0] = 1;
@@ -244,6 +2202,10 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
std::cout<<i<<std::endl;
}
}
+*/
+
+
+
BOOST_AUTO_TEST_SUITE_END()
diff --git a/src/Operators/Vector/vector_dist_operators.hpp b/src/Operators/Vector/vector_dist_operators.hpp
index 62ae5170fc2deb6950fb105d520d4081521da764..d1c33b0776094147028f1d39287a673ded80ce01 100644
--- a/src/Operators/Vector/vector_dist_operators.hpp
+++ b/src/Operators/Vector/vector_dist_operators.hpp
@@ -1,9 +1,9 @@
/*
- * vector_dist_operators.hpp
- *
- * Created on: Jun 11, 2016
- * Author: i-bird
- */
+* vector_dist_operators.hpp
+*
+* Created on: Jun 11, 2016
+* Author: i-bird
+*/
#ifndef OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_HPP_
#define OPENFPM_NUMERICS_SRC_OPERATORS_VECTOR_VECTOR_DIST_OPERATORS_HPP_
@@ -74,6 +74,11 @@
#define VECT_DCPSE_V_DOT 103
#define VECT_DCPSE_V_DIV 104
#define VECT_DCPSE_V_CURL2D 105
+#define VECT_COPY_1_TO_N 300
+#define VECT_COPY_N_TO_N 301
+#define VECT_COPY_N_TO_1 302
+#define VECT_PMUL 91
+#define VECT_SUB_UNI 92