diff --git a/mpi_include b/mpi_include
index 4439b1ca3ff6f7e0872a08610a3d8f99ee1db3a9..16981421116cb2ce7f9dcf3dba32df467827856c 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 3a37a5cd2100d56757bcde4ef6ae508b3f7c0c68..0519ecba6ea913e21689ec692e81e9e4973fbf73 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_test.cpp b/src/DCPSE_op/DCPSE_op_test.cpp
index 1c71fc527a35bdbbe73f7b48e7f721c1dcdb23ec..ea07731b133d219e86c02bf49de278609b08eab3 100644
--- a/src/DCPSE_op/DCPSE_op_test.cpp
+++ b/src/DCPSE_op/DCPSE_op_test.cpp
@@ -121,7 +121,7 @@ const bool equationsp::boundary[] = {PERIODIC, NON_PERIODIC};
BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
BOOST_AUTO_TEST_CASE(dcpse_Active2D) {
- const size_t sz[2] = {23, 23};
+ const size_t sz[2] = {31, 31};
Box<2, double> box({0, 0}, {10,10});
double Lx=box.getHigh(0);
double Ly=box.getHigh(1);
@@ -130,8 +130,8 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
double rCut = 3.1 * spacing;
Ghost<2, double> ghost(rCut);
-/* pol V vort Ext Press strain stress Mfield, dPol dV RHS f1 f2 f3 f4 f5 f6 H V_t div */
- vector_dist<2, double, aggregate<VectorS<2, double>,VectorS<2, double>, double[2][2], VectorS<2, double>, double, double[2][2], double[2][2], VectorS<2, double>,VectorS<2, double>, VectorS<2, double> , VectorS<2, double>, double,double,double,double,double,double,double,VectorS<2, double>,double>> Particles(0, box, bc, ghost);
+/* pol V vort Ext Press strain stress Mfield, dPol dV RHS f1 f2 f3 f4 f5 f6 H V_t div H_t */
+ vector_dist<2, double, aggregate<VectorS<2, double>,VectorS<2, double>, double[2][2], VectorS<2, double>, double, double[2][2], double[2][2], VectorS<2, double>,VectorS<2, double>, VectorS<2, double> , VectorS<2, double>, double,double,double,double,double,double,double,VectorS<2, double>,double,double>> Particles(0, box, bc, ghost);
auto it = Particles.getGridIterator(sz);
while (it.isNext()) {
@@ -187,7 +187,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
auto H = getV<17>(Particles);
auto V_t = getV<18>(Particles);
auto div = getV<19>(Particles);
-
+ auto H_t = getV<20>(Particles);
double eta = 1.0;
double nu = -0.5;
@@ -197,6 +197,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
double Kb = 1;
double lambda = 0.1;
double delmu = -1;
+ g = 0;
// Here fill up the boxes for particle boundary detection.
@@ -289,6 +290,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Derivative_x Dx(Particles, 2, rCut,1.9,support_options::RADIUS);
Derivative_y Dy(Particles, 2, rCut,1.9,support_options::RADIUS);
Derivative_xy Dxy(Particles, 2, rCut,1.9,support_options::RADIUS);
+ auto Dyx=Dxy;
Derivative_xx Dxx(Particles, 2, rCut,1.9,support_options::RADIUS);
Derivative_yy Dyy(Particles, 2, rCut,1.9,support_options::RADIUS);
Gradient Grad(Particles, 2, rCut,1.9,support_options::RADIUS);
@@ -311,12 +313,12 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
f5 = 4.0*gama*nu*Pol[x]*Pol[x]*Pol[y]*Pol[y] / (Pol[x]*Pol[x] + Pol[y]*Pol[y]);
f6 = 2.0*gama*nu*Pol[x]*Pol[x]*Pol[x]*Pol[y] / (Pol[x]*Pol[x] + Pol[y]*Pol[y]);
- dV[x] = 0.5*Dy(h[y]) + zeta*Dx(delmu*Pol[x]*Pol[x]) + zeta*Dy(delmu*Pol[x]*Pol[y]) - zeta*Dx(0.5*delmu*(Pol[x]*Pol[x]+Pol[y]*Pol[y])) - 0.5*nu*Dx(-2*h[y]*Pol[x]*Pol[y])
+ dV[x] = -0.5*Dy(h[y]) + zeta*Dx(delmu*Pol[x]*Pol[x]) + zeta*Dy(delmu*Pol[x]*Pol[y]) - zeta*Dx(0.5*delmu*(Pol[x]*Pol[x]+Pol[y]*Pol[y])) - 0.5*nu*Dx(-2*h[y]*Pol[x]*Pol[y])
- 0.5*nu*Dy(h[y]*(Pol[x]*Pol[x] - Pol[y]*Pol[y])) - Dx(sigma[x][x]) - Dy(sigma[x][y]) - g[x]
- 0.5*nu*Dx(-gama*lambda*delmu*(Pol[x]*Pol[x] - Pol[y]*Pol[y])) - 0.5*Dy(-2*gama*lambda*delmu*(Pol[x]*Pol[y]));
- dV[y] = 0.5*Dx(-h[y]) + zeta*Dy(delmu*Pol[y]*Pol[y]) + zeta*Dx(delmu*Pol[x]*Pol[y]) - zeta*Dy(0.5*delmu*(Pol[x]*Pol[x]+Pol[y]*Pol[y])) - 0.5*nu*Dy(-2*h[y]*Pol[x]*Pol[y])
+ dV[y] = -0.5*Dx(-h[y]) + zeta*Dy(delmu*Pol[y]*Pol[y]) + zeta*Dx(delmu*Pol[x]*Pol[y]) - zeta*Dy(0.5*delmu*(Pol[x]*Pol[x]+Pol[y]*Pol[y])) - 0.5*nu*Dy(-2*h[y]*Pol[x]*Pol[y])
- 0.5*nu*Dx(h[y]*(Pol[x]*Pol[x] - Pol[y]*Pol[y])) - Dx(sigma[y][x]) - Dy(sigma[y][y]) - g[y]
- 0.5*nu*Dy(gama*lambda*delmu*(Pol[x]*Pol[x] - Pol[y]*Pol[y])) - 0.5*Dx(-2*gama*lambda*delmu*(Pol[x]*Pol[y]));
@@ -327,17 +329,19 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
vy.setId(1);
double sum=0,sum1=0;
int n=30;
- auto Stokes1 = nu*Lap(V[x]) + 0.5*nu*(Dx(f1*Dx(V[x]))+Dx(f2*0.5*(Dx(V[y])+Dy(V[x])))+Dx(f3*Dy(V[y]))+Dy(f4*Dx(V[x]))+Dy(f5*0.5*(Dx(V[y])+Dy(V[x])))+Dy(f6*Dy(V[y])));
- auto Stokes2 = nu*Lap(V[y]) + 0.5*nu*(Dy(f1*Dx(V[x]))+Dy(f2*0.5*(Dx(V[y])+Dy(V[x])))+Dy(f3*Dy(V[y]))+Dx(f4*Dx(V[x]))+Dx(f5*0.5*(Dx(V[y])+Dy(V[x])))+Dx(f6*Dy(V[y])));
+/* auto Stokes1 = nu*Lap(V[x]) + 0.5*nu*(f1*Dxx(V[x])+Dx(f1)*Dx(V[x])) + (Dx(f2)*0.5*(Dx(V[y])+Dy(V[x])) + f2*0.5*(Dxx(V[y])+Dyx(V[x]))) + (Dx(f3)*Dy(V[y])+ f3*Dyx(V[y])) + (Dy(f4)*Dx(V[x])+f4*Dxy(V[x])) + (Dy(f5)*0.5*(Dx(V[y])+Dy(V[x]))+f5*0.5*(Dxy(V[y])+Dyy(V[x]))) + ( Dy(f6)*Dy(V[y])+f6*Dyy(V[y]) );
+ auto Stokes2 = nu*Lap(V[y]) + 0.5*nu*(f1*Dxy(V[x])+Dy(f1)*Dx(V[x])) + (Dy(f2)*0.5*(Dx(V[y])+Dy(V[x])) + f2*0.5*(Dxy(V[y])+Dyy(V[x]))) + (Dy(f3)*Dy(V[y])+ f3*Dyy(V[y])) + (Dx(f4)*Dx(V[x])+f4*Dxx(V[x])) + (Dx(f5)*0.5*(Dx(V[y])+Dy(V[x]))+f5*0.5*(Dxx(V[y])+Dyx(V[x]))) + ( Dx(f6)*Dy(V[y])+f6*Dyx(V[y]) );
auto Helmholtz = Lap(H);
auto D_y=Dy(H);
- auto D_x=Dx(H);
+ auto D_x=Dx(H);*/
for(int i=1; i<=n ;i++)
{ RHS[x]=Dx(P)+dV[x];
RHS[y]=Dy(P)+dV[y];
Particles.ghost_get<10>();
DCPSE_scheme<equations2dp,decltype(Particles)> Solver( Particles);
+ auto Stokes1 = nu*Lap(V[x]) + 0.5*nu*(f1*Dxx(V[x])+Dx(f1)*Dx(V[x])) + (Dx(f2)*0.5*(Dx(V[y])+Dy(V[x])) + f2*0.5*(Dxx(V[y])+Dyx(V[x]))) + (Dx(f3)*Dy(V[y])+ f3*Dyx(V[y])) + (Dy(f4)*Dx(V[x])+f4*Dxy(V[x])) + (Dy(f5)*0.5*(Dx(V[y])+Dy(V[x]))+f5*0.5*(Dxy(V[y])+Dyy(V[x]))) + ( Dy(f6)*Dy(V[y])+f6*Dyy(V[y]) );
+ auto Stokes2 = nu*Lap(V[y]) + 0.5*nu*(f1*Dxy(V[x])+Dy(f1)*Dx(V[x])) + (Dy(f2)*0.5*(Dx(V[y])+Dy(V[x])) + f2*0.5*(Dxy(V[y])+Dyy(V[x]))) + (Dy(f3)*Dy(V[y])+ f3*Dyy(V[y])) + (Dx(f4)*Dx(V[x])+f4*Dxx(V[x])) + (Dx(f5)*0.5*(Dx(V[y])+Dy(V[x]))+f5*0.5*(Dxx(V[y])+Dyx(V[x]))) + ( Dx(f6)*Dy(V[y])+f6*Dyx(V[y]) );
Solver.impose(Stokes1,bulk,RHS[0],vx);
Solver.impose(Stokes2,bulk,RHS[1],vy);
Solver.impose(V[x], up_p,0,vx);
@@ -349,12 +353,15 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Particles.ghost_get<Velocity>();
div=Div(V);
Particles.ghost_get<19>();
- DCPSE_scheme<equationsp,decltype(Particles)> SolverH(Particles,options_solver::LAGRANGE_MULTIPLIER);
+ DCPSE_scheme<equationsp,decltype(Particles)> SolverH(Particles,options_solver::LAGRANGE_MULTIPLIER);//,
+ auto Helmholtz = Lap(H);
+ auto D_y=Dy(H);
SolverH.impose(Helmholtz,bulk,prop_id<19>());
SolverH.impose(D_y, up_p,0);
SolverH.impose(-D_y, dw_p,0);
SolverH.solve(H);
Particles.ghost_get<17>();
+ Particles.ghost_get<Velocity>();
//std::cout << "Helmholtz Solved" << std::endl;
V=V-Grad(H);
for(int j=0;j<up_p.size();j++)
@@ -376,12 +383,12 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
for(int j=0;j<bulk.size();j++)
{ auto p=bulk.get<0>(j);
sum+=(Particles.getProp<18>(p)[0]-Particles.getProp<1>(p)[0])*(Particles.getProp<18>(p)[0]- Particles.getProp<1>(p)[0])+(Particles.getProp<18>(p)[1]- Particles.getProp<1>(p)[1])*(Particles.getProp<18>(p)[1]- Particles.getProp<1>(p)[1]);
- sum1+= Particles.getProp<1>(p)[0]*Particles.getProp<1>(p)[0]+Particles.getProp<1>(p)[1]*Particles.getProp<1>(p)[1];
+ sum1+=Particles.getProp<1>(p)[0]*Particles.getProp<1>(p)[0]+Particles.getProp<1>(p)[1]*Particles.getProp<1>(p)[1];
}
sum=sqrt(sum);
sum1=sqrt(sum1);
V_t=V;
- std::cout << "Relative l2 convergence error in velocity = " <<sum/sum1<< std::endl;
+ std::cout << "Rel l2 cgs err in V at "<<i<<"= " <<sum/sum1<< std::endl;
Particles.write_frame("Polar",i);
}
@@ -389,10 +396,10 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Particles.write_frame("Polar",n+1);
- u[x][x] = Dx(V[x]);
- u[x][y] = 0.5*(Dx(V[y])+Dy(V[x]));
- u[y][x] = 0.5*(Dy(V[x])+Dx(V[y]));
- u[y][y] = Dy(V[y]);
+ //u[x][x] = Dx(V[x]);
+ //u[x][y] = 0.5*(Dx(V[y])+Dy(V[x]));
+ //u[y][x] = 0.5*(Dy(V[x])+Dx(V[y]));
+ //u[y][y] = Dy(V[y]);
/*
@@ -1055,12 +1062,12 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
BOOST_AUTO_TEST_CASE(dcpse_poisson_Robin_anal) {
// int rank;
// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
- const size_t sz[2] = {81,81};
+ const size_t sz[2] = {161,161};
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;
+ Ghost<2, double> ghost(spacing * 3.1);
+ double rCut = 3.1 * spacing;
BOOST_TEST_MESSAGE("Init vector_dist...");
vector_dist<2, double, aggregate<double,double,double,double,double,double>> domain(0, box, bc, ghost);
@@ -1928,7 +1935,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 = 80;
+ size_t edgeSemiSize = 40;
const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
Box<2, double> box({0, 0}, {1,1});
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 4b185ee486206a45f9ae2c7b0f9de2dd846c60dc..244136baa334e7d3396808415ddf3b91ec0133d5 100644
--- a/src/DCPSE_op/DCPSE_op_test2.cpp
+++ b/src/DCPSE_op/DCPSE_op_test2.cpp
@@ -280,7 +280,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
BOOST_AUTO_TEST_CASE(dcpse_Lid_Stokes) {
- const size_t sz[2] = {31,31};
+ const size_t sz[2] = {161,161};
Box<2, double> box({0, 0}, {1,1});
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing = box.getHigh(0) / (sz[0] - 1);
@@ -393,7 +393,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
double nu=1e-2;
double sum=0,sum2=0;
- int n=15;
+ int n=100;
Particles.write_frame("Stokes",0);
V_t=V;
for(int i=1; i<=n ;i++)
@@ -413,19 +413,19 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
Solver.impose(V_star[1], l_p, 0,vy);
Solver.solve(V_star[0],V_star[1]);
//std::cout << "Stokes Solved" << std::endl;
- RHS=Div(V_star);
+ RHS=-Div(V_star);
DCPSE_scheme<equations1d,decltype(Particles)> SolverH( Particles,options_solver::LAGRANGE_MULTIPLIER);
auto Helmholtz = Lap(H);
auto D_y=Dy(H);
auto D_x=Dx(H);
SolverH.impose(Helmholtz,bulk,prop_id<3>());
- SolverH.impose(D_y, up_p,0);
- SolverH.impose(D_x, r_p, 0);
- SolverH.impose(-D_y, dw_p,0);
- SolverH.impose(-D_x, l_p,0);
+ SolverH.impose(Dy(H), up_p,0);
+ SolverH.impose(Dx(H), r_p, 0);
+ SolverH.impose(Dy(H), dw_p,0);
+ SolverH.impose(Dx(H), l_p,0);
SolverH.solve(H);
//std::cout << "Helmholtz Solved" << std::endl;
- V=V_star-Grad(H);
+ V=V_star+Grad(H);
for(int j=0;j<up_p.size();j++)
{ auto p=up_p.get<0>(j);
Particles.getProp<1>(p)[0] = 1;
@@ -446,7 +446,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
Particles.getProp<1>(p)[0] = 0;
Particles.getProp<1>(p)[1] = 0;
}
- P=P-Lap(H)+0.5*Adv(V_t,H);
+ P=P+Lap(H)-0.5*Adv(V_t,H);
//std::cout << "V,P Corrected" << std::endl;
sum=0;
for(int j=0;j<bulk.size();j++)
@@ -457,12 +457,12 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
sum=sqrt(sum);
sum2=sqrt(sum2);
V_t=V;
- std::cout << "Relative l2 convergence error = " <<sum/sum2<< std::endl;
+ std::cout << "Rel l2 cgs err at "<<i<<"= " <<sum/sum2<< std::endl;
+ if(i%5==0)
Particles.write_frame("Stokes",i);
-
}
- }
+ }
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
diff --git a/src/DCPSE_op/DCPSE_op_test3d.cpp b/src/DCPSE_op/DCPSE_op_test3d.cpp
index d62d07b136c3fcfa5e5c76490fd7d8e97db4bc7d..a90cb53bfbafbe416a3d201a550fab25df802d0f 100644
--- a/src/DCPSE_op/DCPSE_op_test3d.cpp
+++ b/src/DCPSE_op/DCPSE_op_test3d.cpp
@@ -76,7 +76,7 @@ const bool equations3d3::boundary[] = {NON_PERIODIC, NON_PERIODIC};
BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests3)
BOOST_AUTO_TEST_CASE(stokes_3d) {
- size_t grd_sz=8;
+ size_t grd_sz=21;
const size_t sz[3] = {grd_sz,grd_sz,grd_sz};
Box<3, double> box({0, 0,0}, {1,1,1});
size_t bc[3] = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
@@ -154,7 +154,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests3)
while (it2.isNext()) {
auto p = it2.get();
Point<3, double> xp = Particles.getPos(p);
- Particles.getProp<3>(p) =0;
+ Particles.getProp<0>(p) =0;
if (up.isInside(xp) == true) {
up_p.add();
up_p.last().get<0>() = p.getKey();
@@ -214,7 +214,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests3)
vx.setId(0);
vy.setId(1);
- vy.setId(2);
+ vz.setId(2);
Derivative_x Dx(Particles, 2, rCut,1.9,support_options::RADIUS );
Derivative_y Dy(Particles, 2, rCut,1.9,support_options::RADIUS);
@@ -227,12 +227,12 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests3)
double nu=1e-2;
double sum=0,sum2=0;
- int n=15;
+ int n=50;
Particles.write_frame("Stokes3d",0);
V_t=V;
for(int i=1; i<=n ;i++)
{ RHS2=-Grad(P);
- DCPSE_scheme<equations3d3,decltype(Particles)> Solver( Particles);
+ DCPSE_scheme<equations3d3,decltype(Particles)> Solver(Particles);
auto Stokes1 = Adv(V[0],V_star[0])-nu*Lap(V_star[0]);
auto Stokes2 = Adv(V[1],V_star[1])-nu*Lap(V_star[1]);
auto Stokes3 = Adv(V[2],V_star[2])-nu*Lap(V_star[2]);
@@ -258,23 +258,23 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests3)
Solver.impose(V_star[1], b_p, 0,vy);
Solver.impose(V_star[2], b_p, 0,vz);
Solver.solve(V_star[0],V_star[1],V_star[2]);
- //std::cout << "Stokes Solved" << std::endl;
- RHS=Div(V_star);
+ std::cout << "Stokes Solved" << std::endl;
+ RHS=-Div(V_star);
DCPSE_scheme<equations3d,decltype(Particles)> SolverH(Particles,options_solver::LAGRANGE_MULTIPLIER);
auto Helmholtz = Lap(H);
auto D_x=Dx(H);
auto D_y=Dy(H);
auto D_z=Dz(H);
SolverH.impose(Helmholtz,bulk,prop_id<3>());
- SolverH.impose(D_y, up_p,0);
- SolverH.impose(D_x, r_p, 0);
- SolverH.impose(-D_y, dw_p,0);
- SolverH.impose(-D_x, l_p,0);
- SolverH.impose(D_z, f_p,0);
- SolverH.impose(-D_z, b_p,0);
+ SolverH.impose(H, up_p,0);
+ SolverH.impose(H, r_p, 0);
+ SolverH.impose(H, dw_p,0);
+ SolverH.impose(H, l_p,0);
+ SolverH.impose(H, f_p,0);
+ SolverH.impose(H, b_p,0);
SolverH.solve(H);
- //std::cout << "Helmholtz Solved" << std::endl;
- V=V_star-Grad(H);
+ std::cout << "Helmholtz Solved" << std::endl;
+ V=V_star+Grad(H);
for(int j=0;j<up_p.size();j++)
{ auto p=up_p.get<0>(j);
Particles.getProp<1>(p)[0] = 1;
@@ -300,18 +300,18 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests3)
Particles.getProp<1>(p)[2] = 0;
}
for(int j=0;j<f_p.size();j++)
- { auto p=dw_p.get<0>(j);
+ { auto p=f_p.get<0>(j);
Particles.getProp<1>(p)[0] = 0;
Particles.getProp<1>(p)[1] = 0;
Particles.getProp<1>(p)[2] = 0;
}
for(int j=0;j<b_p.size();j++)
- { auto p=dw_p.get<0>(j);
+ { auto p=b_p.get<0>(j);
Particles.getProp<1>(p)[0] = 0;
Particles.getProp<1>(p)[1] = 0;
Particles.getProp<1>(p)[2] = 0;
}
- P=P-Lap(H)+0.5*Adv(V_t,H);
+ P=P+Lap(H)-0.5*Adv(V_t,H);
std::cout << "V,P Corrected" << std::endl;
sum=0;
for(int j=0;j<bulk.size();j++)
@@ -328,6 +328,325 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests3)
}
}
+
+ BOOST_AUTO_TEST_CASE(dcpse_op_vec3d) {
+// int rank;
+// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+ size_t edgeSemiSize = 20;
+ const size_t sz[3] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1,2 * edgeSemiSize+1};
+ Box<3, double> box({0, 0,0}, {1,1,1});
+ size_t bc[3] = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
+ double spacing = box.getHigh(0) / (sz[0] - 1);
+ double rCut = 3.1 * spacing;
+ Ghost<3, double> ghost(rCut);
+ BOOST_TEST_MESSAGE("Init vector_dist...");
+ double sigma2 = spacing * spacing/ (2 * 4);
+
+ vector_dist<3, double, aggregate<double, VectorS<3, double>, VectorS<3, double>, VectorS<3, double>, VectorS<3, double>,double,double>> domain(
+ 0, box, bc, ghost);
+
+ //Init_DCPSE(domain)
+ BOOST_TEST_MESSAGE("Init domain...");
+// std::random_device rd{};
+// std::mt19937 rng{rd()};
+// std::mt19937 rng{6666666};
+
+// std::normal_distribution<> gaussian{0, sigma2};
+
+ auto it = domain.getGridIterator(sz);
+ size_t pointId = 0;
+ size_t counter = 0;
+ double minNormOne = 999;
+ while (it.isNext()) {
+ domain.add();
+ auto key = it.get();
+ mem_id k0 = key.get(0);
+ double x = k0 * spacing;
+ domain.getLastPos()[0] = x;//+ gaussian(rng);
+ mem_id k1 = key.get(1);
+ double y = k1 * spacing;
+ domain.getLastPos()[1] = y;//+gaussian(rng);
+ mem_id k2 = key.get(2);
+ double z = k2 * spacing;
+ domain.getLastPos()[1] = z;//+gaussian(rng);
+ // Here fill the function value
+ domain.template getLastProp<0>() = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]) + sin(domain.getLastPos()[2]) ;
+ domain.template getLastProp<1>()[0] = cos(domain.getLastPos()[0]);
+ domain.template getLastProp<1>()[1] = cos(domain.getLastPos()[1]) ;
+ domain.template getLastProp<1>()[2] = cos(domain.getLastPos()[2]);
+// domain.template getLastProp<0>() = x * x;
+// domain.template getLastProp<0>() = x;
+ // Here fill the validation value for Df/Dx
+ domain.template getLastProp<2>()[0] = 0;//cos(domain.getLastPos()[0]);//+cos(domain.getLastPos()[1]);
+ domain.template getLastProp<2>()[1] = 0;//-sin(domain.getLastPos()[0]);//+cos(domain.getLastPos()[1]);
+ domain.template getLastProp<3>()[0] = 0;//cos(domain.getLastPos()[0]);//+cos(domain.getLastPos()[1]);
+ domain.template getLastProp<3>()[1] = 0;//-sin(domain.getLastPos()[0]);//+cos(domain.getLastPos()[1]);
+ domain.template getLastProp<3>()[2] = 0;
+
+ domain.template getLastProp<4>()[0] = cos(domain.getLastPos()[0]) * (sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1])) +
+ cos(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));
+ domain.template getLastProp<4>()[1] = -sin(domain.getLastPos()[0]) * (sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1])) -
+ sin(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));
+ domain.template getLastProp<4>()[2] = -sin(domain.getLastPos()[0]) * (sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1])) -
+ sin(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));
+
+ domain.template getLastProp<5>() = -(sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]) + sin(domain.getLastPos()[2])) ;
+
+
+
+
+// domain.template getLastProp<2>() = 2 * x;
+// domain.template getLastProp<2>() = 1;
+
+ ++counter;
+ ++it;
+ }
+ BOOST_TEST_MESSAGE("Sync domain across processors...");
+
+ domain.map();
+ domain.ghost_get<0>();
+
+ //Derivative_x Dx(domain, 2, rCut);
+ //Derivative_y Dy(domain, 2, rCut);
+ //Gradient Grad(domain, 2, rCut);
+ //Laplacian Lap(domain, 2, rCut, 3);
+ Advection Adv(domain, 2, rCut, 1.9,support_options::RADIUS);
+ auto v = getV<1>(domain);
+ auto P = getV<0>(domain);
+ auto dv = getV<3>(domain);
+ auto dP = getV<6>(domain);
+
+
+// typedef boost::mpl::int_<std::is_fundamental<point_expression_op<Point<2U, double>, point_expression<double>, Point<2U, double>, 3>>::value>::blabla blabla;
+
+// std::is_fundamental<decltype(o1.value(key))>
+
+ //vv=Lap(P);
+ //dv=Lap(v);//+Dy(P);
+ dv = Adv(v, v);//+Dy(P);
+ auto it2 = domain.getDomainIterator();
+
+ double worst1 = 0.0;
+
+ while (it2.isNext()) {
+ auto p = it2.get();
+
+ //std::cout << "VALS: " << domain.getProp<3>(p)[0] << " " << domain.getProp<4>(p)[0] << std::endl;
+ //std::cout << "VALS: " << domain.getProp<3>(p)[1] << " " << domain.getProp<4>(p)[1] << std::endl;
+
+ //domain.getProp<0>(p)=std::sqrt((domain.getProp<3>(p)[0] - domain.getProp<4>(p)[0])*(domain.getProp<3>(p)[0] - domain.getProp<4>(p)[0])+(domain.getProp<3>(p)[1] - domain.getProp<4>(p)[1])*(domain.getProp<3>(p)[1] - domain.getProp<4>(p)[1]));
+
+ if (fabs(domain.getProp<3>(p)[1] - domain.getProp<4>(p)[1]) > worst1) {
+ worst1 = fabs(domain.getProp<3>(p)[1] - domain.getProp<4>(p)[1]);
+
+ }
+
+ ++it2;
+ }
+
+ std::cout << "Maximum Error in component 2: " << worst1 << std::endl;
+
+ //Adv.checkMomenta(domain);
+ //Adv.DrawKernel<2>(domain,0);
+
+ //domain.deleteGhost();
+ domain.write("v1");
+
+ dP = Adv(v, P);//+Dy(P);
+ auto it3 = domain.getDomainIterator();
+
+ double worst2 = 0.0;
+
+ while (it3.isNext()) {
+ auto p = it3.get();
+
+ //std::cout << "VALS: " << domain.getProp<3>(p)[0] << " " << domain.getProp<4>(p)[0] << std::endl;
+ //std::cout << "VALS: " << domain.getProp<3>(p)[1] << " " << domain.getProp<4>(p)[1] << std::endl;
+
+ //domain.getProp<0>(p)=std::sqrt((domain.getProp<3>(p)[0] - domain.getProp<4>(p)[0])*(domain.getProp<3>(p)[0] - domain.getProp<4>(p)[0])+(domain.getProp<3>(p)[1] - domain.getProp<4>(p)[1])*(domain.getProp<3>(p)[1] - domain.getProp<4>(p)[1]));
+
+ if (fabs(domain.getProp<6>(p) - domain.getProp<5>(p)) > worst2) {
+ worst2 = fabs(domain.getProp<6>(p) - domain.getProp<5>(p));
+
+ }
+
+ ++it3;
+ }
+
+ std::cout << "Maximum Error: " << worst2 << std::endl;
+
+ //Adv.checkMomenta(domain);
+ //Adv.DrawKernel<2>(domain,0);
+
+ domain.deleteGhost();
+ domain.write("v2");
+
+
+
+ //std::cout << demangle(typeid(decltype(v)).name()) << "\n";
+
+ //Debug<decltype(expr)> a;
+
+ //typedef decltype(expr)::blabla blabla;
+
+ //auto err = Dx + Dx;
+ }
+
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ BOOST_AUTO_TEST_CASE(dcpse_poisson_Robin_anal) {
+// int rank;
+// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+ const size_t sz[2] = {161,161};
+ 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.1);
+ double rCut = 3.1 * 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,1.9,support_options::RADIUS);
+ Derivative_y Dy(domain, 2, rCut,1.9,support_options::RADIUS);
+ //Gradient Grad(domain, 2, rCut);
+ Laplacian Lap(domain, 2, rCut,1.9,support_options::RADIUS);
+ //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<equations3d,decltype(domain)> Solver( 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_sphere) {
const size_t sz[3] = {31,31,31};
Box<3, double> box({0, 0}, {1,1});
@@ -365,8 +684,8 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests3)
Particles.ghost_get<0>();
Particles.write_frame("Sphere",i);\
- }
-*/
+ }*/
+
BOOST_AUTO_TEST_SUITE_END()