From 3c4e69a5642a36e60538a996d40a015c00a0a999 Mon Sep 17 00:00:00 2001
From: absingh <absingh@mpi-cbg.de>
Date: Tue, 7 Apr 2020 14:09:25 +0200
Subject: [PATCH] Slice Solver test
---
mpi_include | 2 +-
mpi_libs | 2 +-
src/DCPSE_op/DCPSE_op.hpp | 42 +++---
src/DCPSE_op/DCPSE_op_test.cpp | 227 +++++++++++++++++++++++++++++++-
src/DCPSE_op/DCPSE_op_test2.cpp | 36 +++--
5 files changed, 260 insertions(+), 49 deletions(-)
diff --git a/mpi_include b/mpi_include
index 4439b1ca..16981421 100644
--- a/mpi_include
+++ b/mpi_include
@@ -1 +1 @@
--I/home/i-bird/MPI/include
\ No newline at end of file
+-I
\ No newline at end of file
diff --git a/mpi_libs b/mpi_libs
index 3a37a5cd..0519ecba 100644
--- a/mpi_libs
+++ b/mpi_libs
@@ -1 +1 @@
--Wl,-rpath -Wl,/home/i-bird/MPI/lib -Wl,--enable-new-dtags -pthread /home/i-bird/MPI/lib/libmpi.so
\ No newline at end of file
+
\ No newline at end of file
diff --git a/src/DCPSE_op/DCPSE_op.hpp b/src/DCPSE_op/DCPSE_op.hpp
index e864655e..bdc10901 100644
--- a/src/DCPSE_op/DCPSE_op.hpp
+++ b/src/DCPSE_op/DCPSE_op.hpp
@@ -662,13 +662,13 @@ class Derivative_x
public:
template<typename particles_type>
- Derivative_x(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor)
+ Derivative_x(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor, double rcut_verlet = -1.0)
{
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(0) = 1;
- dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor);
+ dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor, rcut_verlet);
}
template<typename operand_type>
@@ -690,13 +690,13 @@ class Derivative_y
public:
template<typename particles_type>
- Derivative_y(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor)
+ Derivative_y(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor, double rcut_verlet = -1.0)
{
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(1) = 1;
- dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor);
+ dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor, rcut_verlet);
}
template<typename operand_type>
@@ -716,13 +716,13 @@ class Derivative_z
public:
template<typename particles_type>
- Derivative_z(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor)
+ Derivative_z(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor, double rcut_verlet = -1.0)
{
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(2) = 1;
- dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor);
+ dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor, rcut_verlet);
}
template<typename operand_type>
@@ -747,7 +747,7 @@ class Gradient
public:
template<typename particles_type>
- Gradient(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor)
+ Gradient(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor, double rcut_verlet = -1.0)
{
typedef Dcpse<particles_type::dims,particles_type> DCPSE_type;
@@ -760,7 +760,7 @@ public:
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(i) = 1;
- new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor);
+ new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor, rcut_verlet);
dcpse_ptr++;
}
}
@@ -781,7 +781,7 @@ class Curl2D
public:
template<typename particles_type>
- Curl2D(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor)
+ Curl2D(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor, double rcut_verlet = -1.0)
{
typedef Dcpse<particles_type::dims,particles_type> DCPSE_type;
@@ -791,11 +791,11 @@ public:
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(1) = 1;
- new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor);
+ new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor, rcut_verlet);
dcpse_ptr++;
p.zero();
p.get(0) = 1;
- new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor);
+ new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor, rcut_verlet);
dcpse_ptr++;
}
@@ -877,7 +877,7 @@ class Divergence
public:
template<typename particles_type>
- Divergence(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor)
+ Divergence(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor, double rcut_verlet = -1.0)
{
typedef Dcpse<particles_type::dims,particles_type> DCPSE_type;
@@ -890,7 +890,7 @@ public:
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(i) = 1;
- new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor);
+ new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor, rcut_verlet);
dcpse_ptr++;
}
}
@@ -913,7 +913,7 @@ class Advection
public:
template<typename particles_type>
- Advection(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor)
+ Advection(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor, double rcut_verlet = -1.0)
{
typedef Dcpse<particles_type::dims,particles_type> DCPSE_type;
@@ -926,7 +926,7 @@ public:
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(i) = 1;
- new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor);
+ new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,scaling_factor, rcut_verlet);
dcpse_ptr++;
}
@@ -976,14 +976,14 @@ class Derivative_xy
public:
template<typename particles_type>
- Derivative_xy(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor)
+ Derivative_xy(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor, double rcut_verlet = -1.0)
{
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(0) = 1;
p.get(1) = 1;
- dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,dcpse_scaling_factor);
+ dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,dcpse_scaling_factor, rcut_verlet);
}
template<typename operand_type>
@@ -1005,14 +1005,14 @@ class Derivative_xx
public:
template<typename particles_type>
- Derivative_xx(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor)
+ Derivative_xx(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor, double rcut_verlet = -1.0)
{
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(0) = 2;
p.get(1) = 0;
- dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,dcpse_scaling_factor);
+ dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,dcpse_scaling_factor, rcut_verlet);
}
template<typename operand_type>
@@ -1034,14 +1034,14 @@ class Derivative_yy
public:
template<typename particles_type>
- Derivative_yy(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor)
+ Derivative_yy(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double scaling_factor=dcpse_scaling_factor, double rcut_verlet = -1.0)
{
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(0) = 0;
p.get(1) = 2;
- dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,dcpse_scaling_factor);
+ dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,dcpse_scaling_factor, rcut_verlet);
}
template<typename operand_type>
diff --git a/src/DCPSE_op/DCPSE_op_test.cpp b/src/DCPSE_op/DCPSE_op_test.cpp
index 380c80d3..66ea2e27 100644
--- a/src/DCPSE_op/DCPSE_op_test.cpp
+++ b/src/DCPSE_op/DCPSE_op_test.cpp
@@ -42,6 +42,29 @@ struct equations {
typedef Vector<double> Vector_type;
};
+struct equations2d1 {
+ //! dimensionaly of the equation ( 3D problem ...)
+ static const unsigned int dims = 2;
+ //! number of fields in the system
+ static const unsigned int nvar = 2;
+
+ //! boundary at X and Y
+ static const bool boundary[];
+
+ //! type of space float, double, ...
+ typedef double stype;
+
+ //! type of base particles
+ typedef vector_dist<dims, double, aggregate<double>> b_part;
+
+ //! type of SparseMatrix for the linear solver
+ typedef SparseMatrix<double, int, EIGEN_BASE> SparseMatrix_type;
+
+ //! type of Vector for the linear solver
+ typedef Vector<double> Vector_type;
+};
+const bool equations2d1::boundary[] = {NON_PERIODIC, NON_PERIODIC};
+
const bool equations::boundary[] = {NON_PERIODIC, NON_PERIODIC};
//template<typename T>
@@ -552,12 +575,12 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
BOOST_AUTO_TEST_CASE(dcpse_Lid_sf) {
- const size_t sz[2] = {16,16};
+ 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 = spacing + spacing*1e-02;
+ double rCut = 2.5*spacing + 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);
@@ -780,7 +803,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
// Derivative_x Dx(Particles, 2, rCut,1);
// Derivative_y Dy(Particles, 2, rCut,1);
// Gradient Grad(Particles, 2, rCut,1);
- Laplacian Lap(Particles, 1, rCut, 3.9,3.0*spacing + 1e-6);
+ Laplacian Lap(Particles, 2, rCut, 1.9,3.1*spacing);
// Curl2D Curl(Particles, 2, rCut, 1);
auto its = Particles.getDomainIterator();
@@ -879,7 +902,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing = box.getHigh(0) / (sz[0] - 1);
Ghost<2, double> ghost(spacing * 3);
- double rCut = 0.5 * spacing;
+ double rCut = 2 * spacing;
BOOST_TEST_MESSAGE("Init vector_dist...");
vector_dist<2, double, aggregate<double,double,double,double,double,double>> domain(0, box, bc, ghost);
@@ -905,10 +928,10 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
domain.map();
domain.ghost_get<0>();
- Derivative_x Dx(domain, 2, rCut,2);
- Derivative_y Dy(domain, 2, rCut,2);
+ Derivative_x Dx(domain, 2, rCut,1.9,3.1*spacing);
+ Derivative_y Dy(domain, 2, rCut,1.9,3.1*spacing);
//Gradient Grad(domain, 2, rCut);
- Laplacian Lap(domain, 2, rCut, 2);
+ Laplacian Lap(domain, 2, rCut, 1.9,3.1*spacing);
//Advection Adv(domain, 3, rCut, 3);
//Solver Sol_Lap(Lap),Sol_Dx(Dx);
@@ -2263,6 +2286,196 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Particles.write("test_out");
}
+
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ BOOST_AUTO_TEST_CASE(dcpse_slice_solver) {
+// int rank;
+// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+ const size_t sz[2] = {31,31};
+ 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 * spacing;
+ BOOST_TEST_MESSAGE("Init vector_dist...");
+
+ vector_dist<2, double, aggregate<VectorS<2, double>,VectorS<2, double>,VectorS<2, double>,VectorS<2, double>,VectorS<2, double>,VectorS<2, double>>> domain(0, box, bc, ghost);
+
+
+ //Init_DCPSE(domain)
+ BOOST_TEST_MESSAGE("Init domain...");
+
+ 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,1.9,3.1*spacing);
+ Derivative_y Dy(domain, 2, rCut,1.9,3.1*spacing);
+ //Gradient Grad(domain, 2, rCut);
+ Laplacian Lap(domain, 2, rCut, 1.9,3.1*spacing);
+ //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);
+ if (up.isInside(xp) == true) {
+ up_p.add();
+ up_p.last().get<0>() = p.getKey();
+ domain.getProp<1>(p)[0] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p)[0] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+
+ domain.getProp<1>(p)[1] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p)[1] = 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)[0] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p)[0] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+
+ domain.getProp<1>(p)[1] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p)[1] = 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)[0] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p)[0] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+
+ domain.getProp<1>(p)[1] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p)[1] = 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)[0] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p)[0] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+
+ domain.getProp<1>(p)[1] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p)[1] = 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)[0] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p)[0] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+
+ domain.getProp<1>(p)[1] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ domain.getProp<3>(p)[1] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
+ }
+ ++it2;
+ }
+
+ eq_id vx,vy;
+
+ vx.setId(0);
+ vy.setId(1);
+
+ DCPSE_scheme<equations2d1,decltype(domain)> Solver(2 * rCut, domain);
+ auto Poisson0 = Lap(v[0]);
+ auto Poisson1 = Lap(v[1]);
+ //auto D_x = Dx(v[1]);
+ //auto D_y = Dy(v[1]);
+ Solver.impose(Poisson0, bulk, RHS[0],vx);
+ Solver.impose(Poisson1, bulk, RHS[1],vy);
+ Solver.impose(v[0], up_p, RHS[0],vx);
+ Solver.impose(v[1], up_p, RHS[1],vy);
+ Solver.impose(v[0], r_p, RHS[0],vx);
+ Solver.impose(v[1], r_p, RHS[1],vy);
+ Solver.impose(v[0], dw_p, RHS[0],vx);
+ Solver.impose(v[1], dw_p, RHS[1],vy);
+ Solver.impose(v[0], l_p, RHS[0],vx);
+ Solver.impose(v[1], l_p, RHS[1],vy);
+ Solver.solve(sol);
+ DCPSE_sol=Lap(sol);
+ double worst1 = 0.0;
+ double worst2 = 0.0;
+
+
+ v=sol-RHS;
+
+ for(int j=0;j<bulk.size();j++)
+ { auto p=bulk.get<0>(j);
+ if (fabs(domain.getProp<3>(p)[0] - domain.getProp<2>(p)[0]) >= worst1) {
+ worst1 = fabs(domain.getProp<3>(p)[0] - domain.getProp<2>(p)[0]);
+ }
+ domain.getProp<4>(p)[0] = fabs(domain.getProp<3>(p)[0] - domain.getProp<2>(p)[0]);
+
+ }
+ for(int j=0;j<bulk.size();j++)
+ { auto p=bulk.get<0>(j);
+ if (fabs(domain.getProp<3>(p)[1] - domain.getProp<2>(p)[1]) >= worst1) {
+ worst1 = fabs(domain.getProp<3>(p)[1] - domain.getProp<2>(p)[1]);
+ }
+ domain.getProp<4>(p)[1] = fabs(domain.getProp<3>(p)[1] - domain.getProp<2>(p)[1]);
+
+ }
+ std::cout << "Maximum Analytic Error in slice x: " << worst1 << std::endl;
+ std::cout << "Maximum Analytic Error in slice y: " << worst2 << std::endl;
+
+ domain.write("Slice_anasol");
+ }
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
BOOST_AUTO_TEST_SUITE_END()
diff --git a/src/DCPSE_op/DCPSE_op_test2.cpp b/src/DCPSE_op/DCPSE_op_test2.cpp
index 9ff87917..e92c2c1e 100644
--- a/src/DCPSE_op/DCPSE_op_test2.cpp
+++ b/src/DCPSE_op/DCPSE_op_test2.cpp
@@ -76,12 +76,12 @@ const bool equations1d::boundary[] = {NON_PERIODIC, NON_PERIODIC};
BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
BOOST_AUTO_TEST_CASE(dcpse_Lid_Stokes) {
- const size_t sz[2] = {81,81};
+ 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 = 0.5 * spacing;
+ double rCut = 2.5 * spacing;
// P V v_star RHS V_BC Helmholtz
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);
@@ -104,7 +104,6 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
openfpm::vector<aggregate<int>> l_p;
openfpm::vector<aggregate<int>> r_p;
-
auto P = getV<0>(Particles);
auto V = getV<1>(Particles);
auto V_star = getV<2>(Particles);
@@ -180,12 +179,12 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
vx.setId(0);
vy.setId(1);
- 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);
+ Derivative_x Dx(Particles, 2, rCut,1.9,3.1*spacing );
+ Derivative_y Dy(Particles, 2, rCut,1.9,3.1*spacing);
+ Gradient Grad(Particles, 2, rCut,1.9,3.1*spacing );
+ Laplacian Lap(Particles, 1, rCut, 4.1,3.1*spacing);
+ Advection Adv(Particles, 2, rCut, 1.9,3.1*spacing);
+ Divergence Div(Particles, 2, rCut, 1.9,3.1*spacing);
//starting the simulation at a nice *continuous* place
@@ -205,7 +204,6 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
V_star = V + (V_t - 1e-3*Grad(P));
V = V_star;
-
int n=10;
double nu=1e-2;
Particles.write_frame("Stokes",0);
@@ -259,12 +257,12 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
BOOST_AUTO_TEST_CASE(dcpse_Lid_normal) {
- const size_t sz[2] = {81,81};
+ 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 = 0.5 * spacing;
+ double rCut = 2.5 * 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);
@@ -356,13 +354,13 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
++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=5e-4;
+ Derivative_x Dx(Particles, 2, rCut,1.9,3.1*spacing);
+ Derivative_y Dy(Particles, 2, rCut,1.9,3.1*spacing);
+ Gradient Grad(Particles, 2, rCut,1.9,3.1*spacing);
+ Laplacian Lap(Particles, 2, rCut, 1.9,3.1*spacing);
+ Advection Adv(Particles, 2, rCut, 1.9,3.1*spacing);
+ Divergence Div(Particles, 2, rCut, 1.9,3.1*spacing);
+ double dt=1e-3;
int n=50;
double nu=1e-2;
dV=dt*(nu*Lap(V)-Adv(V,V));
--
GitLab