From c351a4bb13d79c82b9676f92fbdfb6568bee0618 Mon Sep 17 00:00:00 2001
From: absingh <absingh@mpi-cbg.de>
Date: Mon, 11 May 2020 14:35:31 +0200
Subject: [PATCH] Active 2d and DCPSE update fix
---
src/DCPSE_op/DCPSEActiveGelTest.cpp | 613 +++++++++++++--------
src/DCPSE_op/DCPSE_copy/Dcpse.hpp | 16 +-
src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp | 2 +-
src/DCPSE_op/DCPSE_op.hpp | 16 +
src/DCPSE_op/DCPSE_op_test.cpp | 62 ++-
src/DCPSE_op/DCPSE_op_test3d.cpp | 294 ++++++++++
src/DCPSE_op/Fast_test.cpp | 44 +-
7 files changed, 790 insertions(+), 257 deletions(-)
diff --git a/src/DCPSE_op/DCPSEActiveGelTest.cpp b/src/DCPSE_op/DCPSEActiveGelTest.cpp
index f26dfc22..d6c58d7c 100644
--- a/src/DCPSE_op/DCPSEActiveGelTest.cpp
+++ b/src/DCPSE_op/DCPSEActiveGelTest.cpp
@@ -20,7 +20,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
BOOST_AUTO_TEST_CASE(Active2DPetsc) {
timer tt2;
tt2.start();
- const size_t sz[2] = {41, 41};
+ const size_t sz[2] = {17, 17};
Box<2, double> box({0, 0}, {10, 10});
double Lx = box.getHigh(0);
double Ly = box.getHigh(1);
@@ -213,9 +213,9 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
double V_err=1;
int n=0;
int ctr = 0;
- double dt = 3e-3;
+ double dt = 1e-3;
double tim = 0;
- double tf = 7.5;
+ double tf = 1e-1;
while (tim <= tf) {
tt.start();
Particles.ghost_get<Polarization>();
@@ -317,7 +317,6 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Solver.impose(V[y], l_p, 0, vy);
Solver.impose(V[x], r_p, 0, vx);
Solver.impose(V[y], r_p, 0, vy);
-
Solver.solve_with_solver(solverPetsc, V[x], V[y]);
//std::cout << "Stokes Solved in " << tt.getwct() << " seconds." << std::endl;
@@ -337,27 +336,32 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Particles.ghost_get<Velocity>();
V[x] = V[x] + Dx(H);
V[y] = V[y] + Dy(H);
+ P = P + (Dxx(H)+Dyy(H));
for (int j = 0; j < up_p.size(); j++) {
auto p = up_p.get<0>(j);
Particles.getProp<1>(p)[0] = 0;
Particles.getProp<1>(p)[1] = 0;
+ //Particles.getProp<4>(p) = 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<4>(p) = 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;
+ //Particles.getProp<4>(p) = 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;
+ // Particles.getProp<4>(p) = 0;
}
- P = P + 1.0*(Dxx(H)+Dyy(H));
Particles.ghost_get<Velocity>();
Particles.ghost_get<Pressure>();
sum = 0;
@@ -378,7 +382,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
n++;
}
tt.stop();
- std::cout << "Rel l2 cgs err in V = " << V_err << " and took " << tt.getwct() << " seconds."
+ std::cout << "Rel l2 cgs err in V = " << V_err << " and took " << tt.getwct() << " seconds with "<<n<<" iterations."
<< std::endl;
u[x][x] = Dx(V[x]);
u[x][y] = 0.5 * (Dx(V[y]) + Dy(V[x]));
@@ -494,24 +498,6 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
@@ -519,27 +505,24 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
timer tt2;
tt2.start();
const size_t sz[2] = {81, 81};
- Box<2, double> box({0, 0}, {10, 10});
+ Box<2, double> box({0, 0}, {10,10});
double Lx = box.getHigh(0);
double Ly = box.getHigh(1);
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing = box.getHigh(0) / (sz[0] - 1);
- /*double rCut = 2.8 * spacing;
- int ord = 2;
- double sampling_factor = 1.6;*/
- double alpha_V = 1;
- double alpha_H = 1;
double rCut = 3.1 * spacing;
- int ord = 2;
- double sampling_factor = 1.9;
double rCut2 = 3.1 * spacing;
+ int ord = 2;
int ord2 = 2;
+ double sampling_factor = 1.9;
double sampling_factor2 = 1.9;
- Ghost<2, double> ghost(rCut2);
-
-/* 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, double[2], double[2], double[2], double[2], double[2], double[2]>> Particles(
+ double alpha_V = 1.0;
+ double alpha_P = 1.0;
+ 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 H_t delmu */
+ 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, double[2], double[2], double[2], double[2], double[2], double[2], double, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>>> Particles(
0, box, bc, ghost);
+ vector_dist<2, double, aggregate<double,double,VectorS<2, double>>> Particles_subset(0, box, bc, ghost);
auto it = Particles.getGridIterator(sz);
while (it.isNext()) {
@@ -574,12 +557,14 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
constexpr int Strain_rate = 5;
constexpr int Stress = 6;
constexpr int MolField = 7;
+ auto Pos = getV<PROP_POS>(Particles);
auto Pol = getV<Polarization>(Particles);
auto V = getV<Velocity>(Particles);
auto W = getV<Vorticity>(Particles);
auto g = getV<ExtForce>(Particles);
auto P = getV<Pressure>(Particles);
+ auto P_bulk = getV<0>(Particles_subset); //Pressure only on inside
auto u = getV<Strain_rate>(Particles);
auto sigma = getV<Stress>(Particles);
auto h = getV<MolField>(Particles);
@@ -593,6 +578,9 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
auto f5 = getV<15>(Particles);
auto f6 = getV<16>(Particles);
auto H = getV<17>(Particles);
+ auto H_bulk = getV<1>(Particles_subset); //Pressure only on inside
+ auto Grad_bulk = getV<2>(Particles_subset);
+
auto V_t = getV<18>(Particles);
auto div = getV<19>(Particles);
auto H_t = getV<20>(Particles);
@@ -602,6 +590,11 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
auto Df4 = getV<24>(Particles);
auto Df5 = getV<25>(Particles);
auto Df6 = getV<26>(Particles);
+ auto delmu = getV<27>(Particles);
+ auto k1 = getV<28>(Particles);
+ auto k2 = getV<29>(Particles);
+ auto k3 = getV<30>(Particles);
+ auto k4 = getV<31>(Particles);
double eta = 1.0;
@@ -611,11 +604,12 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
double Ks = 1.0;
double Kb = 1.0;
double lambda = 0.1;
- double delmu = -1.0;
+ //double delmu = -1.0;
g = 0;
- V = 0;
+ delmu = -1.0;
P = 0;
-
+ P_bulk=0;
+ V = 0;
// Here fill up the boxes for particle boundary detection.
Box<2, double> up({box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0},
@@ -629,8 +623,12 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
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> mid({box.getHigh(0) / 2.0 - spacing, box.getHigh(1) / 2.0 - spacing / 2.0},
+ /*Box<2, double> mid({box.getHigh(0) / 2.0 - spacing, box.getHigh(1) / 2.0 - spacing / 2.0},
{box.getHigh(0) / 2.0, box.getHigh(1) / 2.0 + spacing / 2.0});
+*/
+ Box<2, double> mid({box.getLow(0) + 3.1 * spacing,box.getLow(1) + 3.1 * spacing},
+ {box.getHigh(0) - 3.1 * spacing, box.getHigh(1) - 3.1 * spacing});
+
openfpm::vector<Box<2, double>> boxes;
boxes.add(up);
@@ -666,15 +664,13 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
if (mid.isInside(xp) == true) {
ref_p.add();
ref_p.last().get<0>() = p.getKey();
- //Particles.getProp<4>(p) = 0;
+ Particles.getProp<4>(p) = 0;
} else {
bulkP.add();
bulkP.last().get<0>() = p.getKey();
}
bulk.add();
bulk.last().get<0>() = p.getKey();
- // Particles.getProp<0>(p)[x]= cos(2 * M_PI * (cos((2 * xp[x]- sz[x]) / sz[x]) - sin((2 * xp[y] - sz[y]) / sz[y])));
- // Particles.getProp<0>(p)[y] = sin(2 * M_PI * (cos((2 * xp[x]- sz[x]) / sz[x]) - sin((2 * xp[y] - sz[y]) / sz[y])));
}
@@ -682,212 +678,380 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
}
- Derivative_x Dx(Particles, ord, rCut, sampling_factor, support_options::RADIUS), Dx2(Particles, ord2, rCut2,
- sampling_factor2,
- support_options::RADIUS);
- Derivative_y Dy(Particles, ord, rCut, sampling_factor, support_options::RADIUS), Dy2(Particles, ord2, rCut2,
- sampling_factor2,
- support_options::RADIUS);
- Derivative_xy Dxy(Particles, ord2, rCut2, sampling_factor2, support_options::RADIUS);
+
+ for (int i = 0 ; i < bulk.size() ; i++) {
+ Particles_subset.add();
+ Particles_subset.getLastPos()[0] = Particles.getPos(bulk.template get<0>(i))[0];
+ Particles_subset.getLastPos()[1] = Particles.getPos(bulk.template get<0>(i))[1];
+ }
+
+
+ Particles_subset.map();
+ Particles_subset.ghost_get<0>();
+
+ Derivative_x Dx(Particles, ord, rCut, sampling_factor, support_options::RADIUS), Bulk_Dx(Particles_subset, ord2, rCut2, sampling_factor2, support_options::RADIUS);
+ Derivative_y Dy(Particles, ord, rCut, sampling_factor, support_options::RADIUS), Bulk_Dy(Particles_subset, ord2, rCut2, sampling_factor2, support_options::RADIUS);
+ Derivative_xy Dxy(Particles, ord, rCut, sampling_factor, support_options::RADIUS);
auto Dyx = Dxy;
- Derivative_xx Dxx(Particles, ord2, rCut2, sampling_factor2, support_options::RADIUS);
- Derivative_yy Dyy(Particles, ord2, rCut2, sampling_factor2, support_options::RADIUS);
- Gradient Grad(Particles, ord, rCut, sampling_factor, support_options::RADIUS);
- Laplacian Lap(Particles, ord2, rCut2, sampling_factor2, support_options::RADIUS);
- Advection Adv(Particles, ord, rCut, sampling_factor, support_options::RADIUS);
- Divergence Div(Particles, ord, rCut, sampling_factor, support_options::RADIUS);
+ Derivative_xx Dxx(Particles, ord, rCut, sampling_factor, support_options::RADIUS),Bulk_Dxx(Particles_subset, ord2, rCut2, sampling_factor2, support_options::RADIUS);
+ Derivative_yy Dyy(Particles, ord, rCut, sampling_factor, support_options::RADIUS),Bulk_Dyy(Particles_subset, ord2, rCut2, sampling_factor2, support_options::RADIUS);
- Particles.ghost_get<Polarization>();
- sigma[x][x] =
- -Ks * Dx(Pol[x]) * Dx(Pol[x]) - Kb * Dx(Pol[y]) * Dx(Pol[y]) + (Kb - Ks) * Dy(Pol[x]) * Dx(Pol[y]);
- sigma[x][y] =
- -Ks * Dy(Pol[y]) * Dx(Pol[y]) - Kb * Dy(Pol[x]) * Dx(Pol[x]) + (Kb - Ks) * Dx(Pol[y]) * Dx(Pol[x]);
- sigma[y][x] =
- -Ks * Dx(Pol[x]) * Dy(Pol[x]) - Kb * Dx(Pol[y]) * Dy(Pol[y]) + (Kb - Ks) * Dy(Pol[x]) * Dy(Pol[y]);
- sigma[y][y] =
- -Ks * Dy(Pol[y]) * Dy(Pol[y]) - Kb * Dy(Pol[x]) * Dy(Pol[x]) + (Kb - Ks) * Dx(Pol[y]) * Dy(Pol[x]);
- Particles.ghost_get<Stress>();
+ /* Derivative_x Dx(Particles, ord, rCut, sampling_factor), Bulk_Dx(Particles_subset, ord, rCut, sampling_factor2);
+ Derivative_y Dy(Particles, ord, rCut, sampling_factor), Bulk_Dy(Particles_subset, ord, rCut, sampling_factor2);
+ Derivative_xy Dxy(Particles, ord2, rCut2, sampling_factor2);
+ auto Dyx = Dxy;
+ Derivative_xx Dxx(Particles, ord2, rCut2, sampling_factor2),Bulk_Dxx(Particles_subset, ord2, rCut2, sampling_factor2);
+ Derivative_yy Dyy(Particles, ord2, rCut2, sampling_factor2),Bulk_Dyy(Particles_subset, ord2, rCut2, sampling_factor2);*/
+ petsc_solver<double> solverPetsc;
+ solverPetsc.setSolver(KSPGMRES);
+ //solverPetsc.setRestart(250);
+ //solverPetsc.setPreconditioner(PCJACOBI);
+ petsc_solver<double> solverPetsc2;
+ solverPetsc2.setSolver(KSPGMRES);
- h[y] = (Pol[x] * (Ks * Dyy(Pol[y]) + Kb * Dxx(Pol[y]) + (Ks - Kb) * Dxy(Pol[x])) -
- Pol[y] * (Ks * Dxx(Pol[x]) + Kb * Dyy(Pol[x]) + (Ks - Kb) * Dxy(Pol[y])));
- Particles.ghost_get<MolField>();
-/* for (int j = 0; j < bulk.size(); j++) {
- auto p = bulk.get<0>(j);
- double Tempo;
- float Temp;
- Tempo=Particles.getProp<7>(p)[y];
- Temp=Particles.getProp<7>(p)[y];
- std::cout << Particles.getProp<7>(p)[x] << std::endl;
- std::cout << Particles.getProp<7>(p)[y] << std::endl;
- }*/
+ eq_id vx, vy;
+ timer tt;
+ timer tt3;
+ vx.setId(0);
+ vy.setId(1);
+ double V_err_eps = 1e-5;
+ double V_err=1,V_err_old;
+ int n=0;
+ int nmax=75;
+ int ctr = 0,errctr;
+ double dt = 2e-7;
+ double tim = 0;
+ double tf = 2e-6;
+ div=0;
+ while (tim <= tf) {
+ tt.start();
+ Particles.ghost_get<Polarization>();
+ sigma[x][x] =
+ -Ks * Dx(Pol[x]) * Dx(Pol[x]) - Kb * Dx(Pol[y]) * Dx(Pol[y]) + (Kb - Ks) * Dy(Pol[x]) * Dx(Pol[y]);
+ sigma[x][y] =
+ -Ks * Dy(Pol[y]) * Dx(Pol[y]) - Kb * Dy(Pol[x]) * Dx(Pol[x]) + (Kb - Ks) * Dx(Pol[y]) * Dx(Pol[x]);
+ sigma[y][x] =
+ -Ks * Dx(Pol[x]) * Dy(Pol[x]) - Kb * Dx(Pol[y]) * Dy(Pol[y]) + (Kb - Ks) * Dy(Pol[x]) * Dy(Pol[y]);
+ sigma[y][y] =
+ -Ks * Dy(Pol[y]) * Dy(Pol[y]) - Kb * Dy(Pol[x]) * Dy(Pol[x]) + (Kb - Ks) * Dx(Pol[y]) * Dy(Pol[x]);
+ Particles.ghost_get<Stress>();
- f1 = gama * nu * Pol[x] * Pol[x] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) / (Pol[x] * Pol[x] + Pol[y] * Pol[y]);
- f2 = 2.0 * gama * nu * Pol[x] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]);
- f3 = gama * nu * Pol[y] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) / (Pol[x] * Pol[x] + Pol[y] * Pol[y]);
- f4 = 2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]);
- 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[y] * Pol[y] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]);
- Particles.ghost_get<11, 12, 13, 14, 15, 16>();
- Df1[x] = Dx(f1);
- Df2[x] = Dx(f2);
- Df3[x] = Dx(f3);
- Df4[x] = Dx(f4);
- Df5[x] = Dx(f5);
- Df6[x] = Dx(f6);
- Df1[y] = Dy(f1);
- Df2[y] = Dy(f2);
- Df3[y] = Dy(f3);
- Df4[y] = Dy(f4);
- Df5[y] = Dy(f5);
- Df6[y] = Dy(f6);
- Particles.ghost_get<21, 22, 23, 24, 25, 26>();
+ h[y] = (Pol[x] * (Ks * Dyy(Pol[y]) + Kb * Dxx(Pol[y]) + (Ks - Kb) * Dxy(Pol[x])) -
+ Pol[y] * (Ks * Dxx(Pol[x]) + Kb * Dyy(Pol[x]) + (Ks - Kb) * Dxy(Pol[y])));
+ Particles.ghost_get<MolField>();
+ f1 = gama * nu * Pol[x] * Pol[x] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
+ (Pol[x] * Pol[x] + Pol[y] * Pol[y]);
+ f2 = 2.0 * gama * nu * Pol[x] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
+ (Pol[x] * Pol[x] + Pol[y] * Pol[y]);
+ f3 = gama * nu * Pol[y] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
+ (Pol[x] * Pol[x] + Pol[y] * Pol[y]);
+ f4 = 2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]);
+ 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[y] * Pol[y] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]);
+ Particles.ghost_get<11, 12, 13, 14, 15, 16>();
+ Df1[x] = Dx(f1);
+ Df2[x] = Dx(f2);
+ Df3[x] = Dx(f3);
+ Df4[x] = Dx(f4);
+ Df5[x] = Dx(f5);
+ Df6[x] = Dx(f6);
- 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.0 * 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.0 * gama * lambda * delmu * (Pol[x] * Pol[y]));
+ Df1[y] = Dy(f1);
+ Df2[y] = Dy(f2);
+ Df3[y] = Dy(f3);
+ Df4[y] = Dy(f4);
+ Df5[y] = Dy(f5);
+ Df6[y] = Dy(f6);
+ Particles.ghost_get<21, 22, 23, 24, 25, 26>();
- 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.0 * 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.0 * gama * lambda * delmu * (Pol[x] * Pol[y]));
- Particles.ghost_get<9>();
+ 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.0 * 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.0 * gama * lambda * delmu * (Pol[x] * Pol[y]));
- Particles.write_frame("Polar10", 0);
- //Velocity Solution n iterations
- eq_id vx, vy;
- timer tt;
- vx.setId(0);
- vy.setId(1);
- double sum = 0, sum1 = 0;
- int n = 100;
-/* 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 Stokes1 = eta * Lap(V[x])
- + 0.5 * nu * (Df1[x] * Dx(V[x]) + f1 * Dxx(V[x]))
- + 0.5 * nu * (Df2[x] * 0.5 * (Dx(V[y]) + Dy(V[x])) + f2 * 0.5 * (Dxx(V[y]) + Dyx(V[x])))
- + 0.5 * nu * (Df3[x] * Dy(V[y]) + f3 * Dyx(V[y]))
- + 0.5 * nu * (Df4[y] * Dx(V[x]) + f4 * Dxy(V[x]))
- + 0.5 * nu * (Df5[y] * 0.5 * (Dx(V[y]) + Dy(V[x])) + f5 * 0.5 * (Dxy(V[y]) + Dyy(V[x])))
- + 0.5 * nu * (Df6[y] * Dy(V[y]) + f6 * Dyy(V[y]));
- auto Stokes2 = eta * Lap(V[y])
- - 0.5 * nu * (Df1[y] * Dx(V[x]) + f1 * Dxy(V[x]))
- - 0.5 * nu * (Df2[y] * 0.5 * (Dx(V[y]) + Dy(V[x])) + f2 * 0.5 * (Dxy(V[y]) + Dyy(V[x])))
- - 0.5 * nu * (Df3[y] * Dy(V[y]) + f3 * Dyy(V[y]))
- + 0.5 * nu * (Df4[x] * Dx(V[x]) + f4 * Dxx(V[x]))
- + 0.5 * nu * (Df5[x] * 0.5 * (Dx(V[y]) + Dy(V[x])) + f5 * 0.5 * (Dxx(V[y]) + Dyx(V[x])))
- + 0.5 * nu * (Df6[x] * Dy(V[y]) + f6 * Dyx(V[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.0 * 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.0 * gama * lambda * delmu * (Pol[x] * Pol[y]));
+ Particles.ghost_get<9>();
- auto Helmholtz = Lap(H);
- auto D_y = Dy(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<equations2d2E, decltype(Particles)> Solver(Particles);
- Solver.impose(Stokes1, bulk, RHS[0], vx);
- Solver.impose(Stokes2, bulk, RHS[1], vy);
- Solver.impose(V[x], up_p, 0, vx);
- Solver.impose(V[y], up_p, 0, vy);
- Solver.impose(V[x], dw_p, 0, vx);
- Solver.impose(V[y], dw_p, 0, vy);
- Solver.impose(V[x], l_p, 0, vx);
- Solver.impose(V[y], l_p, 0, vy);
- Solver.impose(V[x], r_p, 0, vx);
- Solver.impose(V[y], r_p, 0, vy);
+ Particles.write("PolarI");
+ //Velocity Solution n iterations
+ double sum, sum1;
+
+ auto Stokes1 = eta * (Dxx(V[x])+Dyy(V[x]))
+ + 0.5 * nu * (Df1[x] * Dx(V[x]) + f1 * Dxx(V[x]))
+ + 0.5 * nu * (Df2[x] * 0.5 * (Dx(V[y]) + Dy(V[x])) + f2 * 0.5 * (Dxx(V[y]) + Dyx(V[x])))
+ + 0.5 * nu * (Df3[x] * Dy(V[y]) + f3 * Dyx(V[y]))
+ + 0.5 * nu * (Df4[y] * Dx(V[x]) + f4 * Dxy(V[x]))
+ + 0.5 * nu * (Df5[y] * 0.5 * (Dx(V[y]) + Dy(V[x])) + f5 * 0.5 * (Dxy(V[y]) + Dyy(V[x])))
+ + 0.5 * nu * (Df6[y] * Dy(V[y]) + f6 * Dyy(V[y]));
+ auto Stokes2 = eta * (Dxx(V[y])+Dyy(V[y]))
+ - 0.5 * nu * (Df1[y] * Dx(V[x]) + f1 * Dxy(V[x]))
+ - 0.5 * nu * (Df2[y] * 0.5 * (Dx(V[y]) + Dy(V[x])) + f2 * 0.5 * (Dxy(V[y]) + Dyy(V[x])))
+ - 0.5 * nu * (Df3[y] * Dy(V[y]) + f3 * Dyy(V[y]))
+ + 0.5 * nu * (Df4[x] * Dx(V[x]) + f4 * Dxx(V[x]))
+ + 0.5 * nu * (Df5[x] * 0.5 * (Dx(V[y]) + Dy(V[x])) + f5 * 0.5 * (Dxx(V[y]) + Dyx(V[x])))
+ + 0.5 * nu * (Df6[x] * Dy(V[y]) + f6 * Dyx(V[y]));
+ tt.stop();
+ std::cout << "Init of Velocity took " << tt.getwct() << " seconds." << std::endl;
tt.start();
- Solver.solve(V[x], V[y]);
+ V_err=1;
+ n=0;
+ errctr=0;
+ while (V_err >= V_err_eps && n<=nmax) {
+ RHS[x] = dV[x];
+ RHS[y] = dV[y];
+ Grad_bulk[x]=Bulk_Dx(P_bulk);
+ Grad_bulk[y]=Bulk_Dy(P_bulk);
+ for (int i = 0 ; i < bulk.size() ; i++) {
+ Particles.template getProp<10>(bulk.template get<0>(i))[x] += Particles_subset.getProp<2>(i)[x];
+ Particles.template getProp<10>(bulk.template get<0>(i))[y] += Particles_subset.getProp<2>(i)[y];
+ }
+ Particles.ghost_get<10>();
+ DCPSE_scheme<equations2d2, decltype(Particles)> Solver(Particles);
+ Solver.impose(Stokes1, bulk, RHS[0], vx);
+ Solver.impose(Stokes2, bulk, RHS[1], vy);
+ Solver.impose(V[x], up_p, 0, vx);
+ Solver.impose(V[y], up_p, 0, vy);
+ Solver.impose(V[x], dw_p, 0, vx);
+ Solver.impose(V[y], dw_p, 0, vy);
+ Solver.impose(V[x], l_p, 0, vx);
+ Solver.impose(V[y], l_p, 0, vy);
+ Solver.impose(V[x], r_p, 0, vx);
+ Solver.impose(V[y], r_p, 0, vy);
+ //Solver.solve(V[x], V[y]);
+ Solver.solve_with_solver(solverPetsc, V[x], V[y]);
+
+ //std::cout << "Stokes Solved in " << tt.getwct() << " seconds." << std::endl;
+ Particles.ghost_get<Velocity>();
+ /*for (int i = 0 ; i < bulk.size() ; i++) {
+ Particles_subset.getProp<2>(i)[x] = Particles.template getProp<1>(bulk.template get<0>(i))[x];
+ Particles_subset.getProp<2>(i)[y] = Particles.template getProp<1>(bulk.template get<0>(i))[y];
+ }*/
+
+ //
+ //H_bulk = -(Bulk_Dx(Grad_bulk[x])+Bulk_Dy(Grad_bulk[y]));
+ /*for (int i = 0 ; i < bulk.size() ; i++) {
+ Particles.template getProp<19>(bulk.template get<0>(i)) = Particles_subset.getProp<1>(i);
+ }*/
+ div = -(Dx(V[x])+Dy(V[y]));
+
+ Particles.ghost_get<19>();
+ auto Helmholtz = Dxx(H)+Dyy(H);
+ DCPSE_scheme<equations2d1, decltype(Particles)> SolverH(Particles);
+ SolverH.impose(Helmholtz, bulk, prop_id<19>());
+ SolverH.impose(H, up_p, 0);
+ SolverH.impose(H, dw_p, 0);
+ SolverH.impose(H, l_p, 0);
+ SolverH.impose(H, r_p, 0);
+ //SolverH.solve(H);
+ SolverH.solve_with_solver(solverPetsc2, H);
+ //std::cout << "Helmholtz Solved in " << tt.getwct() << " seconds." << std::endl;
+ Particles.ghost_get<17>();
+ Particles.ghost_get<Velocity>();
+
+ for (int i = 0 ; i < bulk.size() ; i++) {
+ Particles_subset.getProp<1>(i) = Particles.template getProp<17>(bulk.template get<0>(i));
+ }
+ Grad_bulk[x]=alpha_V*Bulk_Dx(H_bulk);
+ Grad_bulk[y]=alpha_V*Bulk_Dy(H_bulk);
+ P=P+div;
+
+ for (int i = 0 ; i < bulk.size() ; i++) {
+ Particles.template getProp<1>(bulk.template get<0>(i))[x] += Particles_subset.getProp<2>(i)[x];
+ Particles.template getProp<1>(bulk.template get<0>(i))[y] += Particles_subset.getProp<2>(i)[y];
+ Particles_subset.getProp<0>(i) = Particles.template getProp<4>(bulk.template get<0>(i));
+ }
+ //V[x] = V[x] + Dx(H);
+ //V[y] = V[y] + Dy(H);
+ //P = P + 1.0*(Dxx(H)+Dyy(H));
+ //P = P -(Dx(V[x])+Dy(V[y]));
+ //V[x] = 2*V[x];
+ //V[y] = 2*V[y];
+ for (int j = 0; j < up_p.size(); j++) {
+ auto p = up_p.get<0>(j);
+ //Particles.getProp<1>(p)[0] = 0;
+ //Particles.getProp<1>(p)[1] = 0;
+ Particles.getProp<4>(p) = 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<4>(p) = 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;
+ Particles.getProp<4>(p) = 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;
+ Particles.getProp<4>(p) = 0;
+ }
+ Particles.ghost_get<Velocity>();
+ Particles.ghost_get<Pressure>();
+ sum = 0;
+ sum1 = 0;
+ 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];
+ }
+ sum = sqrt(sum);
+ sum1 = sqrt(sum1);
+ V_t = V;
+ V_err_old=V_err;
+ V_err=sum / sum1;
+ if (V_err>V_err_old)
+ {
+ errctr++;
+ //alpha_P -= 0.1;
+ }
+ else
+ {
+ errctr=0;
+ }
+ if (errctr>5)
+ break;
+ n++;
+ std::cout << "Rel l2 cgs err in V = " << V_err<< std::endl;
+ }
tt.stop();
- std::cout << "Stokes Solved in " << tt.getwct() << " seconds." << std::endl;
- Particles.ghost_get<Velocity>();
- div = -Div(V);
- Particles.ghost_get<19>();
- Particles.write_frame("PolarV10", 1);
- return;
- DCPSE_scheme<equations2d1E, decltype(Particles)> SolverH(
- Particles);//, options_solver::LAGRANGE_MULTIPLIER);
- SolverH.impose(Helmholtz, bulk, prop_id<19>());
- SolverH.impose(H, up_p, 0);
- SolverH.impose(H, dw_p, 0);
- SolverH.impose(H, l_p, 0);
- SolverH.impose(H, r_p, 0);
+ std::cout << "Rel l2 cgs err in V = " << V_err << " and took " << tt.getwct() << " seconds with "<<n<<" iterations."
+ << std::endl;
+
+ 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]);
+
+ W[x][x] = 0;
+ W[x][y] = 0.5 * (Dy(V[x]) - Dx(V[y]));
+ W[y][x] = 0.5 * (Dx(V[y]) - Dy(V[x]));
+ W[y][y] = 0;
+
+
+
+ Particles.write_frame("Polar_Eigen", ctr);
+ ctr++;
+ h[x] = -gama * lambda * delmu + u[x][x] * gama * nu * Pol[x] * Pol[x] / (Pol[x] * Pol[x] + Pol[y] * Pol[y])
+ + u[y][y] * gama * nu * Pol[y] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y])
+ + 2 * u[x][y] * gama * nu * Pol[x] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]);
+
+ Pos = Pos + dt * V;
+
+ Particles.map();
tt.start();
- SolverH.solve(H);
+ Dx.update(Particles);
+ Dy.update(Particles);
+ Dxy.update(Particles);
+ auto Dyx = Dxy;
+ Dxx.update(Particles);
+ Dyy.update(Particles);
tt.stop();
- std::cout << "Helmholtz Solved in " << tt.getwct() << " seconds." << std::endl;
- Particles.ghost_get<17>();
- Particles.ghost_get<Velocity>();
- V = V + alpha_V * Grad(H);
+ std::cout << "Updation of operators took " << tt.getwct() << " seconds." << std::endl;
+
+
+
+ k1[x] = ((h[x] * Pol[x] - h[y] * Pol[y]) / gama + lambda * delmu * Pol[x] -
+ nu * (u[x][x] * Pol[x] + u[x][y] * Pol[y]) + W[x][x] * Pol[x] + W[x][y] * Pol[y] );// -V[x]*Dx(Pol[x]) -V[y]*Dy(Pol[x]));
+ k1[y] = ((h[x] * Pol[y] - h[y] * Pol[x]) / gama + lambda * delmu * Pol[y] -
+ nu * (u[y][x] * Pol[x] + u[y][y] * Pol[y]) + W[y][x] * Pol[x] + W[y][y] * Pol[y] );// -V[x]*Dx(Pol[y]) -V[y]*Dy(Pol[y]));
+
+ k2[x] = ((h[x] * 0.5 * dt * k1[x] * Pol[x] - h[y] * 0.5 * dt * k1[y] * Pol[y]) / gama +
+ lambda * delmu * 0.5 * dt * k1[x] * Pol[x] -
+ nu * (u[x][x] * 0.5 * dt * k1[x] * Pol[x] + u[x][y] * 0.5 * dt * k1[y] * Pol[y]) +
+ W[x][x] * 0.5 * dt * k1[x] * Pol[x] + W[x][y] * 0.5 * dt * k1[y] * Pol[y] );// -V[x]*Dx(0.5 * dt * k1[x]*Pol[x]) -V[y]*Dy(0.5 * dt * k1[x]*Pol[x]));
+ k2[y] = ((h[x] * 0.5 * dt * k1[y] * Pol[y] - h[y] * 0.5 * dt * k1[x] * Pol[x]) / gama +
+ lambda * delmu * 0.5 * dt * k1[y] * Pol[y] -
+ nu * (u[y][x] * 0.5 * dt * k1[x] * Pol[x] + u[y][y] * 0.5 * dt * k1[y] * Pol[y]) +
+ W[y][x] * 0.5 * dt * k1[x] * Pol[x] + W[y][y] * 0.5 * dt * k1[y] * Pol[y] );// -V[x]*Dx(0.5 * dt * k1[y]*Pol[y]) -V[y]*Dy(0.5 * dt * k1[y]*Pol[y]));
+
+ k3[x] = ((h[x] * 0.5 * dt * k2[x] * Pol[x] - h[y] * 0.5 * dt * k2[y] * Pol[y]) / gama +
+ lambda * delmu * 0.5 * dt * k2[x] * Pol[x] -
+ nu * (u[x][x] * 0.5 * dt * k2[x] * Pol[x] + u[x][y] * 0.5 * dt * k2[y] * Pol[y]) +
+ W[x][x] * 0.5 * dt * k2[x] * Pol[x] + W[x][y] * 0.5 * dt * k2[y] * Pol[y] );// -V[x]*Dx(0.5 * dt * k2[x]*Pol[x]) -V[y]*Dy(0.5 * dt * k2[x]*Pol[x]));
+ k3[y] = ((h[x] * 0.5 * dt * k2[y] * Pol[y] - h[y] * 0.5 * dt * k2[x] * Pol[x]) / gama +
+ lambda * delmu * 0.5 * dt * k2[y] * Pol[y] -
+ nu * (u[y][x] * 0.5 * dt * k2[x] * Pol[x] + u[y][y] * 0.5 * dt * k2[y] * Pol[y]) +
+ W[y][x] * 0.5 * dt * k2[x] * Pol[x] + W[y][y] * 0.5 * dt * k2[y] * Pol[y] );// -V[x]*Dx(0.5 * dt * k2[y]*Pol[y]) -V[y]*Dy(0.5 * dt * k2[y]*Pol[y]));
+
+ k4[x] = ((h[x] * dt * k3[x] * Pol[x] - h[y] * dt * k3[y] * Pol[y]) / gama +
+ lambda * delmu * dt * k3[x] * Pol[x] -
+ nu * (u[x][x] * dt * k3[x] * Pol[x] + u[x][y] * dt * k3[y] * Pol[y]) +
+ W[x][x] * dt * k3[x] * Pol[x] + W[x][y] * dt * k3[y] * Pol[y] );// -V[x]*Dx( dt * k3[x]*Pol[x]) -V[y]*Dy( dt * k3[x]*Pol[x]));
+ k4[y] = ((h[x] * dt * k3[y] * Pol[y] - h[y] * dt * k3[x] * Pol[x]) / gama +
+ lambda * delmu * dt * k3[y] * Pol[y] -
+ nu * (u[y][x] * dt * k3[x] * Pol[x] + u[y][y] * dt * k3[y] * Pol[y]) +
+ W[y][x] * dt * k3[x] * Pol[x] + W[y][y] * dt * k3[y] * Pol[y] );// -V[x]*Dx( dt * k3[y]*Pol[y]) -V[y]*Dy( dt * k3[y]*Pol[y]));
+
+ Pol = Pol + (dt / 6.0) * (k1 + (2.0 * k2) + (2.0 * k3) + k4);
+
for (int j = 0; j < up_p.size(); j++) {
auto p = up_p.get<0>(j);
- Particles.getProp<1>(p)[0] = 0;
- Particles.getProp<1>(p)[1] = 0;
+ Particles.getProp<0>(p)[x] = sin(2 * M_PI * (cos((2 * Particles.getPos(p)[x] - Lx) / Lx) -
+ sin((2 * Particles.getPos(p)[y] - Ly) / Ly)));
+ Particles.getProp<0>(p)[y] = cos(2 * M_PI * (cos((2 * Particles.getPos(p)[x] - Lx) / Lx) -
+ sin((2 * Particles.getPos(p)[y] - Ly) / Ly)));
}
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<0>(p)[x] = sin(2 * M_PI * (cos((2 * Particles.getPos(p)[x] - Lx) / Lx) -
+ sin((2 * Particles.getPos(p)[y] - Ly) / Ly)));
+ Particles.getProp<0>(p)[y] = cos(2 * M_PI * (cos((2 * Particles.getPos(p)[x] - Lx) / Lx) -
+ sin((2 * Particles.getPos(p)[y] - Ly) / Ly)));
}
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;
+ Particles.getProp<0>(p)[x] = sin(2 * M_PI * (cos((2 * Particles.getPos(p)[x] - Lx) / Lx) -
+ sin((2 * Particles.getPos(p)[y] - Ly) / Ly)));
+ Particles.getProp<0>(p)[y] = cos(2 * M_PI * (cos((2 * Particles.getPos(p)[x] - Lx) / Lx) -
+ sin((2 * Particles.getPos(p)[y] - Ly) / Ly)));
}
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;
+ Particles.getProp<0>(p)[x] = sin(2 * M_PI * (cos((2 * Particles.getPos(p)[x] - Lx) / Lx) -
+ sin((2 * Particles.getPos(p)[y] - Ly) / Ly)));
+ Particles.getProp<0>(p)[y] = cos(2 * M_PI * (cos((2 * Particles.getPos(p)[x] - Lx) / Lx) -
+ sin((2 * Particles.getPos(p)[y] - Ly) / Ly)));
}
- P = P + alpha_H * Lap(H);
- Particles.ghost_get<Velocity>();
- Particles.ghost_get<Pressure>();
- sum = 0;
- sum1 = 0;
- for (int j = 0; j < bulk.size(); j++) {
+/* 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];
- }
- sum = sqrt(sum);
- sum1 = sqrt(sum1);
- V_t = V;
- //double alpha_V=1;
- //double alpha_H=1;
- std::cout << "Rel l2 cgs err in V at " << i << "= " << sum / sum1 << std::endl;
+ sum=sqrt(Particles.getProp<0>(p)[x]*Particles.getProp<0>(p)[x]+Particles.getProp<0>(p)[y]*Particles.getProp<0>(p)[y]);
+ Particles.getProp<0>(p)[x]/=sum;
+ Particles.getProp<0>(p)[y]/=sum;
+ }*/
+
+ std::cout<<"Time step " << ctr<<" : "<<tim << " over." << std::endl;
+ tim += dt;
std::cout << "----------------------------------------------------------" << std::endl;
- Particles.write_frame("pPolar", i);
-// return;
+
}
Particles.deleteGhost();
- Particles.write_frame("Polar", n + 1);
tt2.stop();
- std::cout << "The simulation took" << tt2.getwct() << "Seconds.";
-
+ std::cout << "The simulation took " << tt2.getwct() << "Seconds.";
}
@@ -1255,6 +1419,8 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
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, double[2], double[2], double[2], double[2], double[2], double[2], double, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>>> Particles(
0, box, bc, ghost);
+ vector_dist<2, double, aggregate<double,double,VectorS<2, double>>> Particles_subset(0, box, bc, ghost);
+
auto it = Particles.getGridIterator(sz);
while (it.isNext()) {
Particles.add();
@@ -1403,12 +1569,19 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
++it2;
}
- Derivative_x Dx(Particles, ord2, rCut2, sampling_factor2, support_options::RADIUS);
- Derivative_y Dy(Particles, ord2, rCut2, sampling_factor2, support_options::RADIUS);
+
+ for (int i = 0 ; i < bulk.size() ; i++) {
+ Particles_subset.add();
+ Particles_subset.getLastPos()[0] = Particles.getPos(bulk.template get<0>(i))[0];
+ Particles_subset.getLastPos()[1] = Particles.getPos(bulk.template get<0>(i))[1];
+ }
+
+ Derivative_x Dx(Particles, ord, rCut, sampling_factor, support_options::RADIUS), Bulk_Dx(Particles_subset, ord2, rCut2, sampling_factor2, support_options::RADIUS);
+ Derivative_y Dy(Particles, ord, rCut, sampling_factor, support_options::RADIUS), Bulk_Dy(Particles_subset, ord2, rCut2, sampling_factor2, support_options::RADIUS);
Derivative_xy Dxy(Particles, ord, rCut, sampling_factor, support_options::RADIUS);
auto Dyx = Dxy;
- Derivative_xx Dxx(Particles, ord, rCut, sampling_factor, support_options::RADIUS);
- Derivative_yy Dyy(Particles, ord, rCut, sampling_factor, support_options::RADIUS);
+ Derivative_xx Dxx(Particles, ord, rCut, sampling_factor, support_options::RADIUS),Bulk_Dxx(Particles_subset, ord2, rCut2, sampling_factor2, support_options::RADIUS);
+ Derivative_yy Dyy(Particles, ord, rCut, sampling_factor, support_options::RADIUS),Bulk_Dyy(Particles_subset, ord2, rCut2, sampling_factor2, support_options::RADIUS);
petsc_solver<double> solverPetsc;
solverPetsc.setSolver(KSPGMRES);
@@ -1516,6 +1689,9 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
auto D_y = Dy(P);
auto D_x = Dx(P);
div = -(Dx(dV[x])+Dy(dV[y]));
+ for (int i = 0 ; i < bulk.size() ; i++) {
+ Particles_subset.getProp<1>(i) = Particles.template getProp<19>(bulk.template get<0>(i));
+ }
Particles.ghost_get<19>();
DCPSE_scheme<equations2d1E, decltype(Particles)> SolverH(Particles);
SolverH.impose(Pressure_Poisson, bulkP, prop_id<19>());
@@ -1574,6 +1750,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
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, double[2], double[2], double[2], double[2], double[2], double[2], double, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>>> Particles(
0, box, bc, ghost);
+
auto it = Particles.getGridIterator(sz);
while (it.isNext()) {
Particles.add();
diff --git a/src/DCPSE_op/DCPSE_copy/Dcpse.hpp b/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
index eeb00923..eebbf2fb 100644
--- a/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
+++ b/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
@@ -58,9 +58,9 @@ private:
std::vector<T> localEps; // Each MPI rank has just access to the local ones
vector_type & particles;
- T rCut;
+ double rCut;
unsigned int convergenceOrder;
- T supportSizeFactor;
+ double supportSizeFactor;
support_options opt;
@@ -80,9 +80,6 @@ public:
monomialBasis(differentialSignature.asArray(), convergenceOrder),
opt(opt)
{
-
-
-
if (supportSizeFactor < 1) {
initializeAdaptive(particles, convergenceOrder, rCut);
} else {
@@ -405,10 +402,10 @@ public:
void initializeUpdate(vector_type &particles)
{
+ localSupports.clear();
localEps.clear();
localCoefficients.clear();
- SupportBuilder<vector_type>
- supportBuilder(particles, differentialSignature, rCut);
+ SupportBuilder<vector_type> supportBuilder(particles, differentialSignature, rCut);
unsigned int requiredSupportSize = monomialBasis.size() * supportSizeFactor;
auto it = particles.getDomainIterator();
@@ -426,7 +423,7 @@ public:
T eps = vandermonde.getEps();
- //localSupports.push_back(support);
+ localSupports.push_back(support);
localEps.push_back(eps);
// Compute the diagonal matrix E
DcpseDiagonalScalingMatrix<dim> diagonalScalingMatrix(monomialBasis);
@@ -516,7 +513,6 @@ private:
unsigned int convergenceOrder,
T rCut,
T supportSizeFactor) {
-
this->rCut=rCut;
this->supportSizeFactor=supportSizeFactor;
this->convergenceOrder=convergenceOrder;
@@ -561,6 +557,8 @@ private:
EMatrix<T, Eigen::Dynamic, 1> a(monomialBasis.size(), 1);
// ...solve the linear system...
a = A.colPivHouseholderQr().solve(b);
+ /*std::cout<<"Cond: "<<conditionNumber(V, 1e2)<<std::endl;
+ std::cout<<a.sum()<<std::endl;*/
// ...and store the solution for later reuse
localCoefficients.push_back(a);
//
diff --git a/src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp b/src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp
index ad34d8e0..daec4c58 100644
--- a/src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp
+++ b/src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp
@@ -180,7 +180,7 @@ std::vector<size_t> SupportBuilder<vector_type>::getPointsInSetOfCells(std::set<
{
for (int i = 0 ; i < rp.size() ; i++)
{
- if (rp.get(i).dist <= rCut)
+ if (rp.get(i).dist < rCut)
{
points.push_back(rp.get(i).offset);
}
diff --git a/src/DCPSE_op/DCPSE_op.hpp b/src/DCPSE_op/DCPSE_op.hpp
index b5c647dd..642b9e1e 100644
--- a/src/DCPSE_op/DCPSE_op.hpp
+++ b/src/DCPSE_op/DCPSE_op.hpp
@@ -1202,6 +1202,14 @@ public:
}
+ template<unsigned int prp, typename particles_type>
+ void DrawKernel(particles_type &particles,int k)
+ {
+ auto dcpse_temp = (Dcpse<particles_type::dims,particles_type> *)dcpse;
+ dcpse_temp->template DrawKernel<prp>(particles,k);
+
+ }
+
template<typename particles_type>
void update(particles_type &particles)
{
@@ -1247,6 +1255,14 @@ public:
}
+ template<unsigned int prp, typename particles_type>
+ void DrawKernel(particles_type &particles,int k)
+ {
+ auto dcpse_temp = (Dcpse<particles_type::dims,particles_type> *)dcpse;
+ dcpse_temp->template DrawKernel<prp>(particles,k);
+
+ }
+
template<typename particles_type>
void update(particles_type &particles)
{
diff --git a/src/DCPSE_op/DCPSE_op_test.cpp b/src/DCPSE_op/DCPSE_op_test.cpp
index 8503db7e..d622036b 100644
--- a/src/DCPSE_op/DCPSE_op_test.cpp
+++ b/src/DCPSE_op/DCPSE_op_test.cpp
@@ -1724,7 +1724,7 @@ 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] = {50,50};
+ 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);
@@ -1940,7 +1940,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
//solver.setRestart(500);
solver.setRestart(500);
- solver.log_monitor();
+ //solver.log_monitor();
// auto A = Solver.getA();
@@ -1985,7 +1985,11 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
domain.write("Robin_anasol");
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
+ //81 0.00131586
+ //161 0.000328664
+ //320 8.30297e-05
+ //520 3.12398e-05
+ //1024 8.08087e-06
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
BOOST_AUTO_TEST_CASE(dcpse_poisson_Dirichlet_anal) {
// int rank;
@@ -1994,8 +1998,8 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
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;
+ 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);
@@ -2021,14 +2025,15 @@ 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,support_options::RADIUS);
+ Derivative_y Dy(domain, 2, rCut,1.9,support_options::RADIUS);
//Gradient Grad(domain, 2, rCut);
- Laplacian Lap(domain, 2, rCut, 2);
+ 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>> bulkF;
openfpm::vector<aggregate<int>> up_p;
openfpm::vector<aggregate<int>> dw_p;
openfpm::vector<aggregate<int>> l_p;
@@ -2079,29 +2084,39 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
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));
+ bulkF.add();
+ bulkF.last().get<0>() = p.getKey();
} 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));
+ bulkF.add();
+ bulkF.last().get<0>() = p.getKey();
} 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));
+ bulkF.add();
+ bulkF.last().get<0>() = p.getKey();
} 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));
+ bulkF.add();
+ bulkF.last().get<0>() = p.getKey();
} 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));
+ bulkF.add();
+ bulkF.last().get<0>() = p.getKey();
}
++it2;
}
@@ -2110,18 +2125,37 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
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.impose(v, up_p, prop_id<1>());
+ Solver.impose(v, r_p, prop_id<1>());
+ Solver.impose(v, dw_p, prop_id<1>());
+ Solver.impose(v, l_p, prop_id<1>());
Solver.solve(sol);
DCPSE_sol=Lap(sol);
+
+
+ for (int j = 0; j < up_p.size(); j++) {
+ auto p = up_p.get<0>(j);
+ domain.getProp<5>(p) = 0;
+ }
+ for (int j = 0; j < dw_p.size(); j++) {
+ auto p = dw_p.get<0>(j);
+ domain.getProp<5>(p) = 0;
+ }
+ for (int j = 0; j < l_p.size(); j++) {
+ auto p = l_p.get<0>(j);
+ domain.getProp<5>(p) = 0;
+ }
+ for (int j = 0; j < r_p.size(); j++) {
+ auto p = r_p.get<0>(j);
+ domain.getProp<5>(p) = 0;
+ }
+
double worst1 = 0.0;
v=abs(DCPSE_sol-RHS);
- for(int j=0;j<bulk.size();j++)
- { auto p=bulk.get<0>(j);
+ for(int j=0;j<bulkF.size();j++)
+ { auto p=bulkF.get<0>(j);
if (fabs(domain.getProp<3>(p) - domain.getProp<2>(p)) >= worst1) {
worst1 = fabs(domain.getProp<3>(p) - domain.getProp<2>(p));
}
diff --git a/src/DCPSE_op/DCPSE_op_test3d.cpp b/src/DCPSE_op/DCPSE_op_test3d.cpp
index c22c2833..74c03486 100644
--- a/src/DCPSE_op/DCPSE_op_test3d.cpp
+++ b/src/DCPSE_op/DCPSE_op_test3d.cpp
@@ -24,6 +24,300 @@
BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests3)
+ BOOST_AUTO_TEST_CASE(stokes_3d_petscP) {
+ 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};
+ double spacing = box.getHigh(0) / (sz[0] - 1);
+ double rCut = 3.1 * spacing;
+ Ghost<3, double> ghost(rCut);
+ // P V v_star RHS V_t Helmholtz
+ vector_dist<3, double, aggregate<double,VectorS<3, double>,VectorS<3, double>,double,VectorS<3, double>,double, double,VectorS<3, double>>> Particles(0, box, bc, ghost);
+ vector_dist<3, double, aggregate<double,double,VectorS<3, double>>> Particles_subset(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;
+ double z = key.get(2) * it.getSpacing(1);
+ Particles.getLastPos()[2] = z;
+ ++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>> f_p;
+ openfpm::vector<aggregate<int>> b_p;
+
+
+ auto P = getV<0>(Particles);
+ auto V = getV<1>(Particles);
+ auto V_star = getV<2>(Particles);
+ V.setVarId(0);
+ auto RHS = getV<3>(Particles);
+ auto V_t = getV<4>(Particles);
+ auto H = getV<5>(Particles);
+ H.setVarId(0);
+ auto temp=getV<6>(Particles);
+ auto RHS2 = getV<7>(Particles);
+ auto P_bulk = getV<0>(Particles_subset);
+ auto H_bulk = getV<1>(Particles_subset); //Pressure only on inside
+ auto Grad_bulk = getV<2>(Particles_subset);
+
+ // Here fill up the boxes for particle detection.
+
+ Box<3, double> up({box.getLow(0) + spacing / 2.0, box.getHigh(1)- spacing / 2.0,box.getLow(2)+ spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getHigh(1) + spacing / 2.0,box.getHigh(2)- spacing / 2.0});
+
+ Box<3, double> down({box.getLow(0) + spacing / 2.0, box.getLow(1) - spacing / 2.0,box.getLow(2)+ spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0,box.getHigh(2)- spacing / 2.0});
+
+ Box<3, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0,box.getLow(2) - spacing / 2.0},
+ {box.getLow(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0,box.getHigh(2) + spacing / 2.0});
+
+ Box<3, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0,box.getLow(2)- spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0,box.getHigh(2)+ spacing / 2.0});
+
+ Box<3, double> front({box.getLow(0) + spacing / 2.0, box.getLow(1) - spacing / 2.0,box.getLow(2) - spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getHigh(1) + spacing / 2.0,box.getLow(2) + spacing / 2.0});
+
+ Box<3, double> back({box.getLow(0) + spacing / 2.0, box.getLow(1) - spacing / 2.0,box.getHigh(2) - spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getHigh(1) + spacing / 2.0,box.getHigh(2) + spacing / 2.0});
+
+ openfpm::vector<Box<3, double>> boxes;
+ boxes.add(up);
+ boxes.add(down);
+ boxes.add(left);
+ boxes.add(right);
+ boxes.add(front);
+ boxes.add(back);
+ VTKWriter<openfpm::vector<Box<3, double>>, VECTOR_BOX> vtk_box;
+ vtk_box.add(boxes);
+ vtk_box.write("boxes_3d.vtk");
+ auto it2 = Particles.getDomainIterator();
+ Particles.ghost_get<0,1,2,3,4,5>();
+ while (it2.isNext()) {
+ auto p = it2.get();
+ Point<3, double> xp = Particles.getPos(p);
+ Particles.getProp<0>(p) =0;
+ 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;
+ Particles.getProp<1>(p)[2] = 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;
+ Particles.getProp<1>(p)[2] = 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;
+ Particles.getProp<1>(p)[2] = 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;
+ Particles.getProp<1>(p)[2] = 0;
+ }
+ else if (front.isInside(xp) == true) {
+ f_p.add();
+ f_p.last().get<0>() = p.getKey();
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ Particles.getProp<1>(p)[2] = 0;
+ }
+ else if (back.isInside(xp) == true) {
+ b_p.add();
+ b_p.last().get<0>() = p.getKey();
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ Particles.getProp<1>(p)[2] = 0;
+ }
+ else {
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ Particles.getProp<1>(p)[0] = 0;
+ Particles.getProp<1>(p)[1] = 0;
+ Particles.getProp<1>(p)[2] = 0;
+ }
+ ++it2;
+ }
+
+ for (int i = 0 ; i < bulk.size() ; i++) {
+ Particles_subset.add();
+ Particles_subset.getLastPos()[0] = Particles.getPos(bulk.template get<0>(i))[0];
+ Particles_subset.getLastPos()[1] = Particles.getPos(bulk.template get<0>(i))[1];
+ Particles_subset.getLastPos()[2] = Particles.getPos(bulk.template get<0>(i))[2];
+ }
+ V_t=V;
+
+
+ eq_id vx,vy,vz;
+
+ vx.setId(0);
+ vy.setId(1);
+ vz.setId(2);
+
+ Derivative_x Dx(Particles, 2, rCut,1.9,support_options::RADIUS),B_Dx(Particles_subset, 2, rCut,1.9,support_options::RADIUS);
+ Derivative_y Dy(Particles, 2, rCut,1.9,support_options::RADIUS),B_Dy(Particles_subset, 2, rCut,1.9,support_options::RADIUS);
+ Derivative_z Dz(Particles, 2, rCut,1.9,support_options::RADIUS),B_Dz(Particles_subset, 2, rCut,1.9,support_options::RADIUS);
+ Gradient Grad(Particles, 2, rCut,1.9,support_options::RADIUS);
+ Laplacian Lap(Particles, 2, rCut,1.9,support_options::RADIUS),B_Lap(Particles_subset, 2, rCut,1.9,support_options::RADIUS);
+ Advection Adv(Particles, 2, rCut,1.9,support_options::RADIUS);
+ Divergence Div(Particles, 2, rCut,1.9,support_options::RADIUS);
+
+
+ petsc_solver<double> solverPetsc;
+ solverPetsc.setRestart(250);
+ petsc_solver<double> solverPetsc2;
+ solverPetsc2.setRestart(250);
+
+
+ double nu=1e-2;
+ Particles.ghost_get<0,1,2,3,4,5,6,7>();
+ double sum=0,sum2=0;
+ int n=50;
+ Particles.write_frame("Stokes3d",0);
+ V_t=V;
+ for(int i=1; i<=n ;i++)
+ {
+ Grad_bulk[0]=-B_Dx(P_bulk);
+ Grad_bulk[1]=-B_Dy(P_bulk);
+ for (int i = 0 ; i < bulk.size() ; i++) {
+ Particles.template getProp<7>(bulk.template get<0>(i))[0] = Particles_subset.getProp<2>(i)[0];
+ Particles.template getProp<7>(bulk.template get<0>(i))[1] = Particles_subset.getProp<2>(i)[1];
+ Particles.template getProp<7>(bulk.template get<0>(i))[2] = Particles_subset.getProp<2>(i)[2];
+ }
+ 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]);
+ Solver.impose(Stokes1,bulk,RHS2[0],vx);
+ Solver.impose(Stokes2,bulk,RHS2[1],vy);
+ Solver.impose(Stokes3,bulk,RHS2[2],vz);
+ Solver.impose(V_star[0], up_p,1.0,vx);
+ Solver.impose(V_star[1], up_p,0,vy);
+ Solver.impose(V_star[2], up_p,0,vz);
+ Solver.impose(V_star[0], r_p, 0,vx);
+ Solver.impose(V_star[1], r_p, 0,vy);
+ Solver.impose(V_star[2], r_p, 0,vz);
+ Solver.impose(V_star[0], dw_p,0,vx);
+ Solver.impose(V_star[1], dw_p,0,vy);
+ Solver.impose(V_star[2], dw_p,0,vz);
+ Solver.impose(V_star[0], l_p, 0,vx);
+ Solver.impose(V_star[1], l_p, 0,vy);
+ Solver.impose(V_star[2], l_p, 0,vz);
+ Solver.impose(V_star[0], f_p, 0,vx);
+ Solver.impose(V_star[1], f_p, 0,vy);
+ Solver.impose(V_star[2], f_p, 0,vz);
+ Solver.impose(V_star[0], b_p, 0,vx);
+ Solver.impose(V_star[1], b_p, 0,vy);
+ Solver.impose(V_star[2], b_p, 0,vz);
+
+ Solver.solve_with_solver(solverPetsc,V_star[0],V_star[1],V_star[2]);
+ std::cout << "Stokes Solved" << std::endl;
+ Particles.ghost_get<2>();
+ RHS=-Div(V_star);
+ DCPSE_scheme<equations3d1,decltype(Particles)> SolverH(Particles);
+ 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(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_with_solver(solverPetsc2,H);
+ Particles.ghost_get<5>();
+ for (int i = 0 ; i < bulk.size() ; i++) {
+ Particles_subset.getProp<1>(i) = Particles.template getProp<5>(bulk.template get<0>(i));
+ }
+ 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;
+ Particles.getProp<1>(p)[1] = 0;
+ Particles.getProp<1>(p)[2] = 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;
+ Particles.getProp<1>(p)[2] = 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;
+ Particles.getProp<1>(p)[2] = 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>(p)[2] = 0;
+ }
+ for(int j=0;j<f_p.size();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=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_bulk=P_bulk+B_Lap(H_bulk);//-0.5*B_Adv(V_t,H);
+ for (int i = 0 ; i < bulk.size() ; i++) {
+ //Particles.template getProp<1>(bulk.template get<0>(i))[x] += Particles_subset.getProp<2>(i)[x];
+ //Particles.template getProp<1>(bulk.template get<0>(i))[y] += Particles_subset.getProp<2>(i)[y];
+ Particles.template getProp<0>(bulk.template get<0>(i)) = Particles_subset.getProp<0>(i);
+ }
+ Particles.ghost_get<0,1>();
+ std::cout << "V,P Corrected" << std::endl;
+ sum=0;
+ for(int j=0;j<bulk.size();j++)
+ { auto p=bulk.get<0>(j);
+ sum+=(Particles.getProp<4>(p)[0]-Particles.getProp<1>(p)[0])*(Particles.getProp<4>(p)[0]- Particles.getProp<1>(p)[0])+(Particles.getProp<4>(p)[1]- Particles.getProp<1>(p)[1])*(Particles.getProp<4>(p)[1]- Particles.getProp<1>(p)[1])+(Particles.getProp<4>(p)[2]- Particles.getProp<1>(p)[2])*(Particles.getProp<4>(p)[2]- Particles.getProp<1>(p)[2]);
+ sum2+= Particles.getProp<1>(p)[0]*Particles.getProp<1>(p)[0]+Particles.getProp<1>(p)[1]*Particles.getProp<1>(p)[1]+Particles.getProp<1>(p)[2]*Particles.getProp<1>(p)[2];
+ }
+ sum=sqrt(sum);
+ sum2=sqrt(sum2);
+ V_t=V;
+ std::cout << "Relative l2 convergence error = " <<sum/sum2<< std::endl;
+ Particles.write_frame("Stokes3d",i);
+ }
+ }
+
BOOST_AUTO_TEST_CASE(stokes_3d_petsc) {
size_t grd_sz=31;
diff --git a/src/DCPSE_op/Fast_test.cpp b/src/DCPSE_op/Fast_test.cpp
index d8e0b1b7..03ac9db9 100644
--- a/src/DCPSE_op/Fast_test.cpp
+++ b/src/DCPSE_op/Fast_test.cpp
@@ -19,12 +19,12 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Box<2, double> box({0, 0}, {10,10});
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing = box.getHigh(0) / (sz[0] - 1);
- double rCut =spacing*(3.1);
+ double rCut = 3.1* spacing;
double ord = 2;
double sampling = 1.9;
double rCut2 = 3.1*spacing;
- double ord2 = 2;
- double sampling2 = 1.9;
+ double ord2 = 1;
+ double sampling2 = 1.2;
double sigma2 = spacing * spacing / (2 * 4);
std::mt19937 rng{7};
@@ -179,20 +179,30 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
- Laplacian Lap(Particles, ord2, rCut2,sampling,support_options::RADIUS);
- Derivative_xy Dxy(Particles,ord, rCut, sampling,support_options::RADIUS);
- Derivative_xy Dxy2(Particles,ord2, rCut2, sampling2,support_options::RADIUS);
+ //Laplacian Lap(Particles, ord2, rCut2,sampling,support_options::RADIUS);
+ Derivative_xx Dxx(Particles, ord2, rCut2,sampling,support_options::RADIUS);
+ Derivative_yy Dyy(Particles, ord2, rCut2,sampling,support_options::RADIUS);
+ Derivative_xy Dxy(Particles,ord2, rCut, sampling,support_options::RADIUS);
+ //Derivative_xy Dxy2(Particles,ord2, rCut2, sampling2,support_options::RADIUS);
Derivative_x Dx(Particles,ord, rCut, sampling,support_options::RADIUS);
Derivative_y Dy(Particles,ord, rCut, sampling,support_options::RADIUS);
+/* Laplacian Lap(Particles, ord2, rCut2,sampling);
+ Derivative_xy Dxy(Particles,ord, rCut, sampling);
+ Derivative_xy Dxy2(Particles,ord2, rCut2, sampling2);
+ Derivative_x Dx(Particles,ord, rCut, sampling);
+ Derivative_y Dy(Particles,ord, rCut, sampling);*/
+
std::cout<<"Dx"<<std::endl;
Dx.checkMomenta(Particles);
std::cout<<"Dy"<<std::endl;
Dy.checkMomenta(Particles);
std::cout<<"Dxy"<<std::endl;
Dxy.checkMomenta(Particles);
- std::cout<<"Lap"<<std::endl;
- Lap.checkMomenta(Particles);
+ std::cout<<"Dxx"<<std::endl;
+ Dxx.checkMomenta(Particles);
+ std::cout<<"Dyy"<<std::endl;
+ Dyy.checkMomenta(Particles);
auto its2 = Particles.getDomainIterator();
int ctr=0;
while (its2.isNext()) {
@@ -200,8 +210,8 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Dx.DrawKernel<0>(Particles, p.getKey());
Dy.DrawKernel<1>(Particles, p.getKey());
Dxy.DrawKernel<2>(Particles, p.getKey());
- Dxy2.DrawKernel<3>(Particles, p.getKey());
- Lap.DrawKernel<4>(Particles, p.getKey());
+ Dxx.DrawKernel<3>(Particles, p.getKey());
+ Dyy.DrawKernel<4>(Particles, p.getKey());
Particles.write_frame("Kernel_moved",ctr);
f1=0;
f2=0;
@@ -224,7 +234,8 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Dy.update(Particles);
Dxy.update(Particles);
auto Dyx = Dxy;
- Lap.update(Particles);
+ Dxx.update(Particles);
+ Dyy.update(Particles);
std::cout<<"Dx"<<std::endl;
Dx.checkMomenta(Particles);
@@ -232,8 +243,10 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Dy.checkMomenta(Particles);
std::cout<<"Dxy"<<std::endl;
Dxy.checkMomenta(Particles);
- std::cout<<"Lap"<<std::endl;
- Lap.checkMomenta(Particles);
+ std::cout<<"Dxx"<<std::endl;
+ Dxx.checkMomenta(Particles);
+ std::cout<<"Dyy"<<std::endl;
+ Dyy.checkMomenta(Particles);
auto its = Particles.getDomainIterator();
@@ -243,13 +256,14 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Dx.DrawKernel<0>(Particles, p.getKey());
Dy.DrawKernel<1>(Particles, p.getKey());
Dxy.DrawKernel<2>(Particles, p.getKey());
- Dxy2.DrawKernel<3>(Particles, p.getKey());
- Lap.DrawKernel<4>(Particles, p.getKey());
+ Dxx.DrawKernel<3>(Particles, p.getKey());
+ Dyy.DrawKernel<4>(Particles, p.getKey());
Particles.write_frame("Kernel_unmoved",ctr);
f1=0;
f2=0;
f3=0;
f4=0;
+ f5=0;
++its;
ctr++;
}
--
GitLab