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/CMakeLists.txt b/src/CMakeLists.txt
index 0dbcbcb851d0ae51e6c46bc6e4f9e512592bd931..30840bd283c7afd162bb23519babbca83501c02c 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.8 FATAL_ERROR)
########################### Executables
-add_executable(numerics DCPSE_op/DCPSE_op_test3d.cpp DCPSE_op/DCPSE_op_test2.cpp DCPSE_op/DCPSE_op_test.cpp main.cpp Matrix/SparseMatrix_unit_tests.cpp interpolation/interpolation_unit_tests.cpp Vector/Vector_unit_tests.cpp Solvers/petsc_solver_unit_tests.cpp FiniteDifference/FDScheme_unit_tests.cpp FiniteDifference/eq_unit_test_3d.cpp FiniteDifference/eq_unit_test.cpp Operators/Vector/vector_dist_operators_unit_tests.cpp ../../openfpm_vcluster/src/VCluster/VCluster.cpp ../../openfpm_devices/src/memory/CudaMemory.cu ../../openfpm_devices/src/memory/HeapMemory.cpp ../../openfpm_devices/src/memory/PtrMemory.cpp ../../openfpm_devices/src/Memleak_check.cpp DCPSE_op/DCPSEActiveGelTest.cpp)
+add_executable(numerics DCPSE_op/DCPSE_op_test3d.cpp DCPSE_op/DCPSE_op_test2.cpp DCPSE_op/DCPSE_op_test.cpp main.cpp Matrix/SparseMatrix_unit_tests.cpp interpolation/interpolation_unit_tests.cpp Vector/Vector_unit_tests.cpp Solvers/petsc_solver_unit_tests.cpp FiniteDifference/FDScheme_unit_tests.cpp FiniteDifference/eq_unit_test_3d.cpp FiniteDifference/eq_unit_test.cpp Operators/Vector/vector_dist_operators_unit_tests.cpp ../../openfpm_vcluster/src/VCluster/VCluster.cpp ../../openfpm_devices/src/memory/CudaMemory.cu ../../openfpm_devices/src/memory/HeapMemory.cpp ../../openfpm_devices/src/memory/PtrMemory.cpp ../../openfpm_devices/src/Memleak_check.cpp DCPSE_op/DCPSEActiveGelTest.cpp DCPSE_op/Fast_test.cpp)
###########################
diff --git a/src/DCPSE_op/DCPSEActiveGelTest.cpp b/src/DCPSE_op/DCPSEActiveGelTest.cpp
index c54b2b06afe8a0fe5c8befa670b1ddf9fb47e7a2..130cd58cdad1de50ed73fcb6f37ee6bd2f9d0273 100644
--- a/src/DCPSE_op/DCPSEActiveGelTest.cpp
+++ b/src/DCPSE_op/DCPSEActiveGelTest.cpp
@@ -165,16 +165,19 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
int ord = 2;
double sampling_factor = 1.9;
- Derivative_x Dx(Particles, ord, rCut, sampling_factor, support_options::RADIUS);
- Derivative_y Dy(Particles, ord, rCut, sampling_factor, support_options::RADIUS);
+ int ord2 = 2;
+ double sampling_factor2 = 1.6;
+ double rCut2 = 2.8;
+ Derivative_x Dx(Particles, ord2, rCut2, sampling_factor2, support_options::RADIUS);
+ Derivative_y Dy(Particles, 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);
- Gradient Grad(Particles, ord, rCut, sampling_factor, support_options::RADIUS);
+ Gradient Grad(Particles, ord2, rCut2, sampling_factor2, support_options::RADIUS);
Laplacian Lap(Particles, ord, rCut, sampling_factor, support_options::RADIUS);
- Advection Adv(Particles, ord, rCut, sampling_factor, support_options::RADIUS);
- Divergence Div(Particles, ord, rCut, sampling_factor, support_options::RADIUS);
+ Advection Adv(Particles, ord2, rCut2, sampling_factor2, support_options::RADIUS);
+ Divergence Div(Particles, ord2, rCut2, sampling_factor2, support_options::RADIUS);
Particles.ghost_get<Polarization>();
sigma[x][x] =
@@ -364,21 +367,24 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
BOOST_AUTO_TEST_CASE(Active2DEigen) {
- const size_t sz[2] = {21, 21};
- Box<2, double> box({0, 0}, {10, 10});
+ const size_t sz[2] = {17, 17};
+ Box<2, double> box({0, 0}, {1,1});
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 = 4.1 * spacing;
- int ord = 3;
- double sampling_factor = 1.2;
-
+ /*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 ord2 = 2;
double sampling_factor2 = 1.9;
-
- Ghost<2, double> ghost(rCut);
+ 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(
@@ -491,8 +497,8 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
while (it2.isNext()) {
auto p = it2.get();
Point<2, double> xp = Particles.getPos(p);
- Particles.getProp<0>(p)[x] = cos(2 * M_PI * (cos((2 * xp[x] - Lx) / Lx) - sin((2 * xp[y] - Ly) / Ly)));
- Particles.getProp<0>(p)[y] = sin(2 * M_PI * (cos((2 * xp[x] - Lx) / Lx) - sin((2 * xp[y] - Ly) / Ly)));
+ Particles.getProp<0>(p)[x] = sin(2 * M_PI * (cos((2 * xp[x] - Lx) / Lx) - sin((2 * xp[y] - Ly) / Ly)));
+ Particles.getProp<0>(p)[y] = cos(2 * M_PI * (cos((2 * xp[x] - Lx) / Lx) - sin((2 * xp[y] - Ly) / Ly)));
if (up.isInside(xp) == true) {
up_p.add();
up_p.last().get<0>() = p.getKey();
@@ -544,7 +550,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
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[y]) * Dx(Pol[x]) + (Kb - Ks) * Dx(Pol[x]) * Dx(Pol[y]);
+ -Ks * Dy(Pol[y]) * Dx(Pol[y]) - Kb * Dy(Pol[x]) * Dx(Pol[x]) + (Kb - Ks) * Dx(Pol[x]) * Dx(Pol[y]);
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] =
@@ -564,24 +570,18 @@ 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]);
Particles.ghost_get<11, 12, 13, 14, 15, 16>();
- Df1[x] = Dx(gama * nu * Pol[x] * Pol[x] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df2[x] = Dx(2.0 * gama * nu * Pol[x] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df3[x] = Dx(gama * nu * Pol[y] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df4[x] = Dx(2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df5[x] = Dx(4.0 * gama * nu * Pol[x] * Pol[x] * Pol[y] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df6[x] = Dx(2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df1[y] = Dy(gama * nu * Pol[x] * Pol[x] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df2[y] = Dy(2.0 * gama * nu * Pol[x] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df3[y] = Dy(gama * nu * Pol[y] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df4[y] = Dy(2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df5[y] = Dy(4.0 * gama * nu * Pol[x] * Pol[x] * Pol[y] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df6[y] = Dy(2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
+ 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>();
@@ -633,8 +633,8 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
(Df6[x] * Dy(V[y]) + f6 * Dyx(V[y]));
auto Helmholtz = Lap(H);
- auto D_y = Dy2(H);
- auto D_x = Dx2(H);
+ auto D_y = Dy(H);
+ auto D_x = Dx(H);
for (int i = 1; i <= n; i++) {
RHS[x] = Dx(P) + dV[x];
@@ -658,19 +658,19 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Particles.ghost_get<Velocity>();
div = -Div(V);
Particles.ghost_get<19>();
- DCPSE_scheme<equations2d1E, decltype(Particles)> SolverH(Particles, options_solver::LAGRANGE_MULTIPLIER);
+ DCPSE_scheme<equations2d1E, decltype(Particles)> SolverH(Particles);//, options_solver::LAGRANGE_MULTIPLIER);
SolverH.impose(Helmholtz, bulk, prop_id<19>());
- SolverH.impose(D_y, up_p, 0);
- SolverH.impose(D_y, dw_p, 0);
- SolverH.impose(D_x, l_p, 0);
- SolverH.impose(D_x, r_p, 0);
+ SolverH.impose(H, up_p, 0);
+ SolverH.impose(H, dw_p, 0);
+ SolverH.impose(H, l_p, 0);
+ SolverH.impose(H, r_p, 0);
tt.start();
SolverH.solve(H);
tt.stop();
std::cout << "Helmholtz Solved in " << tt.getwct() << " seconds." << std::endl;
Particles.ghost_get<17>();
Particles.ghost_get<Velocity>();
- V = V + Grad(H);
+ V = V + alpha_V*Grad(H);
for (int j = 0; j < up_p.size(); j++) {
auto p = up_p.get<0>(j);
Particles.getProp<1>(p)[0] = 0;
@@ -691,7 +691,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Particles.getProp<1>(p)[0] = 0;
Particles.getProp<1>(p)[1] = 0;
}
- P = P + Lap(H);
+ P = P + alpha_H*Lap(H);
Particles.ghost_get<Velocity>();
Particles.ghost_get<Pressure>();
sum = 0;
@@ -708,6 +708,8 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
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;
std::cout << "----------------------------------------------------------" << std::endl;
Particles.write_frame("Polar", i);
@@ -722,7 +724,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
BOOST_AUTO_TEST_CASE(Active2DEigenP) {
- const size_t sz[2] = {31, 31};
+ const size_t sz[2] = {21,21};
Box<2, double> box({0, 0}, {10, 10});
double Lx = box.getHigh(0);
double Ly = box.getHigh(1);
@@ -779,9 +781,11 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
auto Pol = getV<Polarization>(Particles);
auto V = getV<Velocity>(Particles);
+ V.setVarId(0);
auto W = getV<Vorticity>(Particles);
auto g = getV<ExtForce>(Particles);
auto P = getV<Pressure>(Particles);
+ P.setVarId(0);
auto u = getV<Strain_rate>(Particles);
auto sigma = getV<Stress>(Particles);
auto h = getV<MolField>(Particles);
@@ -850,8 +854,8 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
while (it2.isNext()) {
auto p = it2.get();
Point<2, double> xp = Particles.getPos(p);
- Particles.getProp<0>(p)[x] = cos(2 * M_PI * (cos((2 * xp[x] - Lx) / Lx) - sin((2 * xp[y] - Ly) / Ly)));
- Particles.getProp<0>(p)[y] = sin(2 * M_PI * (cos((2 * xp[x] - Lx) / Lx) - sin((2 * xp[y] - Ly) / Ly)));
+ Particles.getProp<0>(p)[x] = sin(2 * M_PI * (cos((2 * xp[x] - Lx) / Lx) - sin((2 * xp[y] - Ly) / Ly)));
+ Particles.getProp<0>(p)[y] = cos(2 * M_PI * (cos((2 * xp[x] - Lx) / Lx) - sin((2 * xp[y] - Ly) / Ly)));
if (up.isInside(xp) == true) {
up_p.add();
up_p.last().get<0>() = p.getKey();
@@ -905,7 +909,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
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[y]) * Dx(Pol[x]) + (Kb - Ks) * Dx(Pol[x]) * Dx(Pol[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] =
@@ -926,41 +930,35 @@ 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]);
Particles.ghost_get<11, 12, 13, 14, 15, 16>();
- Df1[x] = Dx(gama * nu * Pol[x] * Pol[x] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df2[x] = Dx(2.0 * gama * nu * Pol[x] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df3[x] = Dx(gama * nu * Pol[y] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df4[x] = Dx(2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df5[x] = Dx(4.0 * gama * nu * Pol[x] * Pol[x] * Pol[y] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df6[x] = Dx(2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df1[y] = Dy(gama * nu * Pol[x] * Pol[x] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df2[y] = Dy(2.0 * gama * nu * Pol[x] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df3[y] = Dy(gama * nu * Pol[y] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df4[y] = Dy(2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df5[y] = Dy(4.0 * gama * nu * Pol[x] * Pol[x] * Pol[y] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df6[y] = Dy(2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
+ 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>();
- 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]));
+ - 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])
- 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]));
+ - 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]));
Particles.ghost_get<9>();
@@ -1018,8 +1016,8 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Particles.ghost_get<19>();
DCPSE_scheme<equations2d1pE, decltype(Particles)> SolverH(Particles, options_solver::LAGRANGE_MULTIPLIER);
SolverH.impose(Helmholtz, bulk, prop_id<19>());
- SolverH.impose(D_y, up_p, 0);
- SolverH.impose(D_y, dw_p, 0);
+ SolverH.impose(H, up_p, 0);
+ SolverH.impose(H, dw_p, 0);
tt.start();
SolverH.solve(H);
tt.stop();
@@ -1102,6 +1100,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
openfpm::vector<aggregate<int>> bulk;
openfpm::vector<aggregate<int>> bulkP;
+ openfpm::vector<aggregate<int>> bulkF;
openfpm::vector<aggregate<int>> up_p;
openfpm::vector<aggregate<int>> dw_p;
openfpm::vector<aggregate<int>> l_p;
@@ -1206,19 +1205,27 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
l_p.last().get<0>() = p.getKey();
bulk.add();
bulk.last().get<0>() = p.getKey();
+ bulkP.add();
+ bulkP.last().get<0>() = p.getKey();
} else if (right.isInside(xp) == true) {
r_p.add();
r_p.last().get<0>() = p.getKey();
bulk.add();
bulk.last().get<0>() = p.getKey();
+ bulkP.add();
+ bulkP.last().get<0>() = p.getKey();
} else {
if (mid.isInside(xp) == true) {
ref_p.add();
ref_p.last().get<0>() = p.getKey();
+ bulkP.add();
+ bulkP.last().get<0>() = p.getKey();
Particles.getProp<4>(p) = 0;
} else {
bulkP.add();
bulkP.last().get<0>() = p.getKey();
+ bulkF.add();
+ bulkF.last().get<0>() = p.getKey();
}
bulk.add();
bulk.last().get<0>() = p.getKey();
@@ -1248,7 +1255,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
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[y]) * Dx(Pol[x]) + (Kb - Ks) * Dx(Pol[x]) * Dx(Pol[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] =
@@ -1256,6 +1263,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Particles.ghost_get<Stress>();
+
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>();
@@ -1409,9 +1417,9 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
double Ly = box.getHigh(1);
size_t bc[2] = {PERIODIC, NON_PERIODIC};
double spacing = box.getHigh(0) / (sz[0] - 1);
- int ord = 2;
- double rCut = 3.1 * spacing;
- double sampling_factor = 1.9;
+ int ord = 3;
+ double rCut = 4.1 * spacing;
+ double sampling_factor = 1.3;
int ord2 = 2;
@@ -1582,18 +1590,20 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Advection Adv(Particles, ord, rCut, sampling_factor, support_options::RADIUS);
Divergence Div(Particles, ord, rCut, sampling_factor, 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[y]) * Dx(Pol[x]) + (Kb - Ks) * Dx(Pol[x]) * Dx(Pol[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) * Dy(Pol[x]) * Dx(Pol[y]);
+ H_t=-Ks * Dy(Pol[y]) * Dy(Pol[y]) - Kb * Dy(Pol[x]) * Dy(Pol[x]) + (Kb - Ks) * Dy(Pol[x]) * Dx(Pol[y]);
Particles.ghost_get<Stress>();
+
+
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>();
@@ -1607,24 +1617,18 @@ 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]);
Particles.ghost_get<11, 12, 13, 14, 15, 16>();
- Df1[x] = Dx(gama * nu * Pol[x] * Pol[x] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df2[x] = Dx(2.0 * gama * nu * Pol[x] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df3[x] = Dx(gama * nu * Pol[y] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df4[x] = Dx(2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df5[x] = Dx(4.0 * gama * nu * Pol[x] * Pol[x] * Pol[y] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df6[x] = Dx(2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df1[y] = Dy(gama * nu * Pol[x] * Pol[x] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df2[y] = Dy(2.0 * gama * nu * Pol[x] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df3[y] = Dy(gama * nu * Pol[y] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df4[y] = Dy(2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df5[y] = Dy(4.0 * gama * nu * Pol[x] * Pol[x] * Pol[y] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df6[y] = Dy(2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
+ 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>();
@@ -1732,28 +1736,44 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
}
- BOOST_AUTO_TEST_CASE(Active2DEigenP_saddle) {
- const size_t sz[2] = {21, 21};
+
+
+
+
+ //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+
+ BOOST_AUTO_TEST_CASE(Active2DEigen_saddle)
+ {
+ 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);
- size_t bc[2] = {PERIODIC, NON_PERIODIC};
+ size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing = box.getHigh(0) / (sz[0] - 1);
int ord = 2;
- double rCut = 3.1 * spacing;
- double sampling_factor = 1;
+ double rCut = 2.8* spacing;
+ double sampling_factor = 1.6;
int ord2 = 2;
double rCut2 = 3.1 * spacing;
- double sampling_factor2 = 1;
+ double sampling_factor2 = 1.9;
- Ghost<2, double> ghost(rCut);
+ double eta = 1.0;
+ double nu = -0.5;
+ double gama = 0.1;
+ double zeta = 0.07;
+ double Ks = 1.0;
+ double Kb = 1.0;
+ double lambda = 0.1;
+ double delmu = -1.0;
+
+ 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<2I , 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(
+ 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(
0, box, bc, ghost);
-
auto it = Particles.getGridIterator(sz);
while (it.isNext()) {
Particles.add();
@@ -1761,7 +1781,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
double x = key.get(0) * it.getSpacing(0);
Particles.getLastPos()[0] = x;
double y = key.get(1) * it.getSpacing(1);
- Particles.getLastPos()[1] = y * 0.99999;
+ Particles.getLastPos()[1] = y ;
++it;
}
@@ -1821,14 +1841,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
auto Df6 = getV<26>(Particles);
- double eta = 1.0;
- double nu = -0.5;
- double gama = 0.1;
- double zeta = 0.07;
- double Ks = 1.0;
- double Kb = 1.0;
- double lambda = 0.1;
- double delmu = -1.0;
+
g = 0;
V = 0;
P = 0;
@@ -1865,52 +1878,52 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
while (it2.isNext()) {
auto p = it2.get();
Point<2, double> xp = Particles.getPos(p);
- Particles.getProp<0>(p)[x] = cos(2 * M_PI * (cos((2 * xp[x] - Lx) / Lx) - sin((2 * xp[y] - Ly) / Ly)));
- Particles.getProp<0>(p)[y] = sin(2 * M_PI * (cos((2 * xp[x] - Lx) / Lx) - sin((2 * xp[y] - Ly) / Ly)));
+ Particles.getProp<0>(p)[x] = sin(2 * M_PI * (cos((2 * xp[x] - Lx) / Lx) - sin((2 * xp[y] - Ly) / Ly)));
+ Particles.getProp<0>(p)[y] = cos(2 * M_PI * (cos((2 * xp[x] - Lx) / Lx) - sin((2 * xp[y] - Ly) / Ly)));
if (up.isInside(xp) == true) {
up_p.add();
up_p.last().get<0>() = p.getKey();
bulkF.add();
bulkF.last().get<0>() = p.getKey();
+
} else if (down.isInside(xp) == true) {
dw_p.add();
dw_p.last().get<0>() = p.getKey();
bulkF.add();
bulkF.last().get<0>() = p.getKey();
+
} else if (left.isInside(xp) == true) {
l_p.add();
l_p.last().get<0>() = p.getKey();
- bulk.add();
- bulk.last().get<0>() = p.getKey();
bulkF.add();
bulkF.last().get<0>() = p.getKey();
+
} else if (right.isInside(xp) == true) {
r_p.add();
r_p.last().get<0>() = p.getKey();
- bulk.add();
- bulk.last().get<0>() = p.getKey();
bulkF.add();
bulkF.last().get<0>() = p.getKey();
+
} else {
if (mid.isInside(xp) == true) {
ref_p.add();
ref_p.last().get<0>() = p.getKey();
Particles.getProp<4>(p) = 0;
} else {
- bulkP.add();
- bulkP.last().get<0>() = p.getKey();
+
bulkF.add();
bulkF.last().get<0>() = p.getKey();
}
bulk.add();
bulk.last().get<0>() = p.getKey();
+ bulkP.add();
+ bulkP.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])));
}
++it2;
}
-
Derivative_x Dx(Particles, ord, rCut, sampling_factor, support_options::RADIUS), Dx2(Particles, ord2, rCut2,
sampling_factor2,
support_options::RADIUS);
@@ -1930,14 +1943,17 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
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[y]) * Dx(Pol[x]) + (Kb - Ks) * Dx(Pol[x]) * Dx(Pol[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) * Dy(Pol[x]) * Dx(Pol[y]);
+
+ H_t=-Ks * Dy(Pol[y]) * Dy(Pol[y]) - Kb * Dy(Pol[x]) * Dy(Pol[x]) + (Kb - Ks) * Dy(Pol[x]) * Dx(Pol[y]);
Particles.ghost_get<Stress>();
+
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>();
@@ -1951,24 +1967,18 @@ 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]);
Particles.ghost_get<11, 12, 13, 14, 15, 16>();
- Df1[x] = Dx(gama * nu * Pol[x] * Pol[x] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df2[x] = Dx(2.0 * gama * nu * Pol[x] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df3[x] = Dx(gama * nu * Pol[y] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df4[x] = Dx(2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df5[x] = Dx(4.0 * gama * nu * Pol[x] * Pol[x] * Pol[y] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df6[x] = Dx(2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df1[y] = Dy(gama * nu * Pol[x] * Pol[x] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df2[y] = Dy(2.0 * gama * nu * Pol[x] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df3[y] = Dy(gama * nu * Pol[y] * Pol[y] * (Pol[x] * Pol[x] - Pol[y] * Pol[y]) /
- (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df4[y] = Dy(2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df5[y] = Dy(4.0 * gama * nu * Pol[x] * Pol[x] * Pol[y] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
- Df6[y] = Dy(2.0 * gama * nu * Pol[x] * Pol[x] * Pol[x] * Pol[y] / (Pol[x] * Pol[x] + Pol[y] * Pol[y]));
+ 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>();
@@ -2009,24 +2019,23 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
(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 * (f1 * Dxx(V[x]) + Df1[x] * Dx(V[x])) +
+ auto Stokes1 = - Dx(P) + eta * Lap(V[x]) + 0.5 * nu * (f1 * Dxx(V[x]) + Df1[x] * Dx(V[x])) +
(Df2[x] * 0.5 * (Dx(V[y]) + Dy(V[x])) + f2 * 0.5 * (Dxx(V[y]) + Dyx(V[x]))) +
(Df3[x] * Dy(V[y]) + f3 * Dyx(V[y])) + (Df4[y] * Dx(V[x]) + f4 * Dxy(V[x])) +
(Df5[y] * 0.5 * (Dx(V[y]) + Dy(V[x])) + f5 * 0.5 * (Dxy(V[y]) + Dyy(V[x]))) +
- (Df6[y] * Dy(V[y]) + f6 * Dyy(V[y])) - Dx(P);
- auto Stokes2 = eta * Lap(V[y]) + 0.5 * nu * (f1 * Dxy(V[x]) + Df1[y] * Dx(V[x])) +
+ (Df6[y] * Dy(V[y]) + f6 * Dyy(V[y])) ;
+ auto Stokes2 = - Dy(P) + eta * Lap(V[y]) + 0.5 * nu * (f1 * Dxy(V[x]) + Df1[y] * Dx(V[x])) +
(Df2[y] * 0.5 * (Dx(V[y]) + Dy(V[x])) + f2 * 0.5 * (Dxy(V[y]) + Dyy(V[x]))) +
(Df3[y] * Dy(V[y]) + f3 * Dyy(V[y])) + (Df4[x] * Dx(V[x]) + f4 * Dxx(V[x])) +
(Df5[x] * 0.5 * (Dx(V[y]) + Dy(V[x])) + f5 * 0.5 * (Dxx(V[y]) + Dyx(V[x]))) +
- (Df6[x] * Dy(V[y]) + f6 * Dyx(V[y])) - Dy(P);
+ (Df6[x] * Dy(V[y]) + f6 * Dyx(V[y])) ;
auto continuity = Dx(V[x]) + Dy(V[y]);
Particles.ghost_get<Pressure>();
RHS[x] = dV[x];
RHS[y] = dV[y];
Particles.ghost_get<10>();
- DCPSE_scheme<equations2d3pE, decltype(Particles)> Solver(Particles);
- //Solver.setEquationStructure({eq_struct::VECTOR,eq_struct::SCALAR});
+ DCPSE_scheme<equations2d3E, decltype(Particles)> Solver(Particles);
Solver.impose(Stokes1, bulk, RHS[0], vx);
Solver.impose(Stokes2, bulk, RHS[1], vy);
Solver.impose(continuity, bulkF, 0, ic);
@@ -2034,7 +2043,12 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
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.impose(P, ref_p, 0, ic);
+
tt.start();
BOOST_TEST_MESSAGE("Solver Imposed...");
/*//auto A = Solver.getA();
diff --git a/src/DCPSE_op/DCPSE_Solver.hpp b/src/DCPSE_op/DCPSE_Solver.hpp
index c234d42a226f837fbe18470aa1c797b2e5c06854..6f16dc78790205733ca929b2be6ecf98bb408ab5 100644
--- a/src/DCPSE_op/DCPSE_Solver.hpp
+++ b/src/DCPSE_op/DCPSE_Solver.hpp
@@ -18,23 +18,23 @@ struct prop_id {};
class eq_id
{
- int id;
+ int id;
public:
- eq_id()
- :id(0)
- {}
+ eq_id()
+ :id(0)
+ {}
- int getId()
- {
- return id;
- }
+ int getId()
+ {
+ return id;
+ }
- void setId(int id)
- {
- this->id = id;
- }
+ void setId(int id)
+ {
+ this->id = id;
+ }
};
enum options_solver
@@ -128,14 +128,14 @@ class DCPSE_scheme: public MatrixAssembler
}
else
{
- if (v_cl.rank() == v_cl.size()-1)
- {
- b.resize(Sys_eqs::nvar * tot + 1, Sys_eqs::nvar * sz + 1);
- }
- else
- {
- b.resize(Sys_eqs::nvar * tot + 1, Sys_eqs::nvar * sz);
- }
+ if (v_cl.rank() == v_cl.size()-1)
+ {
+ b.resize(Sys_eqs::nvar * tot + 1, Sys_eqs::nvar * sz + 1);
+ }
+ else
+ {
+ b.resize(Sys_eqs::nvar * tot + 1, Sys_eqs::nvar * sz);
+ }
}
// Calculate the starting point
@@ -211,7 +211,7 @@ class DCPSE_scheme: public MatrixAssembler
*
*/
variable_b(particles_type & parts)
- :parts(parts)
+ :parts(parts)
{}
/*! \brief Get the b term on a grid point
@@ -240,15 +240,15 @@ class DCPSE_scheme: public MatrixAssembler
// A and B must have the same rows
if (row != row_b)
{
- std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " the term B and the Matrix A for Ax=B must contain the same number of rows\n";
- return;
+ std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " the term B and the Matrix A for Ax=B must contain the same number of rows\n";
+ return;
+ }
+ if (row_b != p_map.size_local() * Sys_eqs::nvar) {
+ std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " your system is underdetermined you set "
+ << row_b << " conditions " << " but i am expecting " << p_map.size_local() * Sys_eqs::nvar
+ << std::endl;
+ return;
}
- if (row_b != p_map.size_local() * Sys_eqs::nvar) {
- std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " your system is underdetermined you set "
- << row_b << " conditions " << " but i am expecting " << p_map.size_local() * Sys_eqs::nvar
- << std::endl;
- return;
- }
// Indicate all the non zero rows
openfpm::vector<unsigned char> nz_rows;
@@ -256,11 +256,11 @@ class DCPSE_scheme: public MatrixAssembler
for (size_t i = 0 ; i < trpl.size() ; i++)
{
- if (trpl.get(i).row() - s_pnt*Sys_eqs::nvar >= nz_rows.size())
- {
- std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " It seems that you are setting colums that does not exist \n";
- }
- if (trpl.get(i).value() != 0)
+ if (trpl.get(i).row() - s_pnt*Sys_eqs::nvar >= nz_rows.size())
+ {
+ std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " It seems that you are setting colums that does not exist \n";
+ }
+ if (trpl.get(i).value() != 0)
{nz_rows.get(trpl.get(i).row() - s_pnt*Sys_eqs::nvar) = true;}
}
@@ -275,7 +275,7 @@ class DCPSE_scheme: public MatrixAssembler
for (size_t i = 0 ; i < trpl.size() ; i++)
{
- if (trpl.get(i).value() != 0)
+ if (trpl.get(i).value() != 0)
{nz_cols.get(trpl.get(i).col()) = true;}
}
@@ -322,18 +322,18 @@ class DCPSE_scheme: public MatrixAssembler
template<typename solType, typename exp1, typename ... othersExp>
void copy_nested(solType & x, unsigned int & comp, exp1 exp, othersExp ... exps)
{
- copy_impl(x,exp,comp);
- comp++;
+ copy_impl(x,exp,comp);
+ comp++;
- copy_nested(x,comp,exps ...);
+ copy_nested(x,comp,exps ...);
}
template<typename solType, typename exp1>
void copy_nested(solType & x, unsigned int & comp, exp1 exp)
{
- copy_impl(x,exp,comp);
- comp++;
+ copy_impl(x,exp,comp);
+ comp++;
}
public:
@@ -378,11 +378,11 @@ public:
template<typename ... expr_type>
void solve(expr_type ... exps)
{
- if (sizeof...(exps) != Sys_eqs::nvar)
- {std::cerr << __FILE__ << ":" << __LINE__ << " Error the number of properties you gave does not match the solution in\
+ if (sizeof...(exps) != Sys_eqs::nvar)
+ {std::cerr << __FILE__ << ":" << __LINE__ << " Error the number of properties you gave does not match the solution in\
dimensionality, I am expecting " << Sys_eqs::nvar <<
- " properties " << std::endl;};
- typename Sys_eqs::solver_type solver;
+ " properties " << std::endl;};
+ typename Sys_eqs::solver_type solver;
// umfpack_solver<double> solver;
auto x = solver.solve(getA(opt),getB(opt));
@@ -423,10 +423,10 @@ public:
template<typename SolverType, typename ... expr_type>
void try_solve_with_solver(SolverType & solver, expr_type ... exps)
{
- if (sizeof...(exps) != Sys_eqs::nvar)
- {std::cerr << __FILE__ << ":" << __LINE__ << " Error the number of properties you gave does not match the solution in\
+ if (sizeof...(exps) != Sys_eqs::nvar)
+ {std::cerr << __FILE__ << ":" << __LINE__ << " Error the number of properties you gave does not match the solution in\
dimensionality, I am expecting " << Sys_eqs::nvar <<
- " properties " << std::endl;};
+ " properties " << std::endl;};
auto x = solver.try_solve(getA(opt),getB(opt));
@@ -445,7 +445,7 @@ public:
*
*/
DCPSE_scheme(particles_type & part, options_solver opt = options_solver::STANDARD)
- :parts(part),p_map(part.getDecomposition(),0),row(0),row_b(0),opt(opt)
+ :parts(part),p_map(part.getDecomposition(),0),row(0),row_b(0),opt(opt)
{
p_map.resize(part.size_local());
@@ -470,8 +470,8 @@ public:
*/
template<typename T, typename index_type, unsigned int prp_id>
void impose(const T & op , openfpm::vector<index_type> & subset,
- const prop_id<prp_id> & num,
- eq_id id = eq_id())
+ const prop_id<prp_id> & num,
+ eq_id id = eq_id())
{
auto itd = subset.template getIteratorElements<0>();
@@ -497,8 +497,8 @@ public:
*/
template<typename T, typename index_type, typename RHS_type, typename sfinae = typename std::enable_if< !std::is_fundamental<RHS_type>::type::value >::type >
void impose(const T & op , openfpm::vector<index_type> & subset,
- const RHS_type & rhs,
- eq_id id = eq_id())
+ const RHS_type & rhs,
+ eq_id id = eq_id())
{
auto itd = subset.template getIteratorElements<0>();
@@ -521,9 +521,9 @@ public:
*
*/
template<typename T, typename index_type> void impose(const T & op ,
- openfpm::vector<index_type> & subset,
- const typename Sys_eqs::stype num,
- eq_id id = eq_id())
+ openfpm::vector<index_type> & subset,
+ const typename Sys_eqs::stype num,
+ eq_id id = eq_id())
{
auto itd = subset.template getIteratorElements<0>();
@@ -669,15 +669,17 @@ public:
// get the particle
auto key = it.get();
-/* if (key == 298 && create_vcluster().rank() == 1)
+/*
+ if (key == 298 && create_vcluster().rank() == 1)
{
int debug = 0;
debug++;
- }*/
+ }
+*/
// Calculate the non-zero colums
typename Sys_eqs::stype coeff = 1.0;
- op.template value_nz<Sys_eqs>(p_map,key,cols,coeff,id);
+ op.template value_nz<Sys_eqs>(p_map,key,cols,coeff,0);
// indicate if the diagonal has been set
bool is_diag = false;
@@ -693,11 +695,11 @@ public:
if (trpl.last().row() == trpl.last().col())
is_diag = true;
- if (trpl.last().col() == 1323)
+/* if (trpl.last().col() == 1323)
{
- int debug = 0;
- debug++;
- }
+ int debug = 0;
+ debug++;
+ }*/
// std::cout << "(" << trpl.last().row() << "," << trpl.last().col() << "," << trpl.last().value() << ")" << "\n";
}
@@ -737,4 +739,4 @@ public:
-#endif //OPENFPM_PDATA_DCPSE_SOLVER_HPP
+#endif //OPENFPM_PDATA_DCPSE_SOLVER_HPP
\ No newline at end of file
diff --git a/src/DCPSE_op/DCPSE_copy/Dcpse.hpp b/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
index e4523a8e1617e16c9ba5a1b2efd2fcffcbdb948e..8eb2768ea60b978820e34c682cf0a3a4285a6713 100644
--- a/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
+++ b/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
@@ -103,7 +103,7 @@ public:
particles.template getProp<prp>(xqK) += computeKernel(normalizedArg, a);
DrawKernelKounter++;
}
- //std::cout<<"Number of Neighbours: "<<DrawKernelKounter<<std::endl;
+// std::cout<<"Number of Neighbours: "<<DrawKernelKounter<<std::endl;
}
template<unsigned int prp>
diff --git a/src/DCPSE_op/DCPSE_op.hpp b/src/DCPSE_op/DCPSE_op.hpp
index a37c443ed9923e2f237cad87647e52a3452b5b7a..c443dd0a76b12820877a6902c88c9a1e40e530d2 100644
--- a/src/DCPSE_op/DCPSE_op.hpp
+++ b/src/DCPSE_op/DCPSE_op.hpp
@@ -673,6 +673,7 @@ public:
p.get(0) = 1;
dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor, opt);
+
}
template<typename operand_type>
@@ -683,6 +684,14 @@ public:
return vector_dist_expression_op<operand_type,dcpse_type,VECT_DCPSE>(arg,*(dcpse_type *)dcpse);
}
+
+ template<unsigned int prp, typename particles_type>
+ void DrawKernel(particles_type &particles,int k)
+ {
+ auto dcpse2 = (Dcpse<particles_type::dims,particles_type> *)dcpse;
+ dcpse2->template DrawKernel<prp>(particles,k);
+
+ }
};
@@ -701,6 +710,12 @@ public:
p.get(1) = 1;
dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor, opt);
+
+ Dcpse<particles_type::dims,particles_type> * dcpse_ptr = (Dcpse<particles_type::dims,particles_type> *)dcpse;
+
+ new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor, opt);
+ dcpse_ptr++;
+
}
template<typename operand_type>
@@ -711,6 +726,17 @@ public:
return vector_dist_expression_op<operand_type,dcpse_type,VECT_DCPSE>(arg,*(dcpse_type *)dcpse);
}
+
+ template<unsigned int prp, typename particles_type>
+ void DrawKernel(particles_type &particles,int k)
+ {
+ auto dcpse2 = (Dcpse<particles_type::dims,particles_type> *)dcpse;
+ dcpse2->template DrawKernel<prp>(particles,k);
+
+ }
+
+
+
};
class Derivative_z
{
@@ -776,6 +802,18 @@ public:
return vector_dist_expression_op<operand_type,dcpse_type,VECT_DCPSE_V>(arg,*(dcpse_type(*)[operand_type::vtype::dims])dcpse);
}
+
+ template<unsigned int prp, typename particles_type>
+ void DrawKernel(particles_type &particles,int k)
+ {
+ Dcpse<particles_type::dims,particles_type> * dcpse_ptr = (Dcpse<particles_type::dims,particles_type> *)dcpse;
+
+ for (int i = 0 ; i < particles_type::dims ; i++)
+ {
+ dcpse_ptr[i].template DrawKernel<prp>(particles,i,k);
+ }
+
+ }
};
class Curl2D
@@ -988,6 +1026,12 @@ public:
p.get(1) = 1;
dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,dcpse_oversampling_factor, opt);
+
+ Dcpse<particles_type::dims,particles_type> * dcpse_ptr = (Dcpse<particles_type::dims,particles_type> *)dcpse;
+
+ new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor, opt);
+ dcpse_ptr++;
+
}
template<typename operand_type>
@@ -998,6 +1042,15 @@ public:
return vector_dist_expression_op<operand_type,dcpse_type,VECT_DCPSE>(arg,*(dcpse_type *)dcpse);
}
+
+ template<unsigned int prp, typename particles_type>
+ void DrawKernel(particles_type &particles,int k)
+ {
+ Dcpse<particles_type::dims,particles_type> * dcpse_ptr = (Dcpse<particles_type::dims,particles_type> *)dcpse;
+
+ dcpse_ptr[0].template DrawKernel<prp>(particles,k);
+
+ }
};
diff --git a/src/DCPSE_op/DCPSE_op_test.cpp b/src/DCPSE_op/DCPSE_op_test.cpp
index fc760366c43abab1f60c1e7a482b70cada5d73be..3948ac13d348d5842bfc362e94d0afcbeaf8ea02 100644
--- a/src/DCPSE_op/DCPSE_op_test.cpp
+++ b/src/DCPSE_op/DCPSE_op_test.cpp
@@ -31,6 +31,7 @@ const bool equations3d1E::boundary[] = {NON_PERIODIC, NON_PERIODIC};
const bool equations3d3E::boundary[] = {NON_PERIODIC, NON_PERIODIC};
const bool equations2d1E::boundary[] = {NON_PERIODIC, NON_PERIODIC};
const bool equations2d2E::boundary[] = {NON_PERIODIC, NON_PERIODIC};
+const bool equations2d3E::boundary[] = {NON_PERIODIC, NON_PERIODIC};
const bool equations2d1pE::boundary[] = {PERIODIC, NON_PERIODIC};
const bool equations2d2pE::boundary[] = {PERIODIC, NON_PERIODIC};
const bool equations2d3pE::boundary[] = {PERIODIC, NON_PERIODIC};
@@ -1539,8 +1540,14 @@ 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.1);
- double rCut = 3.1*spacing;
+ double rCut =spacing*(3.1);
+ double ord = 3;
+ double sampling = 2.3;
+ double rCut2 = 3.1*spacing;
+ double ord2 = 2;
+ double sampling2 = 1.9;
+
+ Ghost<2, double> ghost(rCut);
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);
@@ -1637,7 +1644,6 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
if (up.isInside(xp) == true) {
up_p.add();
up_p.last().get<0>() = p.getKey();
- Particles. template getProp<1>(p) = 12.0/spacing;// -2.0*Particles.getProp<0>(p)/(spacing*spacing) - 12.0/spacing;
}
else if (down.isInside(xp) == true) {
dw_p.add();
@@ -1671,188 +1677,48 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
}
++it2;
}
-/* //test for copy
- for(int j=0;j<up_p1.size();j++)
- {
- auto p1=up_p1.get<0>(j);
- Point<2, double> xp = Particles.getPos(p1);
- Particles.getProp<3>(p1) = xp[0];
- }
- for(int j=0;j<l_p1.size();j++)
- {
- auto p1=l_p1.get<0>(j);
- Point<2, double> xp = Particles.getPos(p1);
- Particles.getProp<3>(p1) = xp[1];
- }
- for(int j=0;j<r_p1.size();j++)
- {
- auto p1=r_p1.get<0>(j);
- Point<2, double> xp = Particles.getPos(p1);
- Particles.getProp<3>(p1) = xp[1];
- }
- for(int j=0;j<dw_p1.size();j++)
- {
- auto p1=dw_p1.get<0>(j);
- Point<2, double> xp = Particles.getPos(p1);
- Particles.getProp<3>(p1)= xp[0];
- }
- //Copy Module
- for(int j=0;j<up_p1.size();j++)
- { auto p=up_p.get<0>(j);
- auto p1=up_p1.get<0>(j);
- Particles.getProp<3>(p) = Particles.getProp<3>(p1);
- if(j==0)
- { auto p=l_p.get<0>(l_p.size()-1);
- auto p2=l_p.get<0>(l_p.size()-2);
- Particles.getProp<3>(p) = Particles.getProp<3>(p1);
- Particles.getProp<3>(p2) = Particles.getProp<3>(p1);
- }
- if(j==up_p1.size()-1)
- { auto p=r_p.get<0>(r_p.size()-1);
- auto p2=r_p.get<0>(r_p.size()-2);
- Particles.getProp<3>(p) = Particles.getProp<3>(p1);
- Particles.getProp<3>(p2) = Particles.getProp<3>(p1);
- }
- }
-
- for(int j=0;j<l_p1.size();j++)
- {
- auto p=l_p.get<0>(j+2);
- auto p1=l_p1.get<0>(j);
- Particles.getProp<3>(p) = Particles.getProp<3>(p1);
- }
- for(int j=0;j<r_p1.size();j++)
- {
- auto p=r_p.get<0>(j+2);
- auto p1=r_p1.get<0>(j);
- Particles.getProp<3>(p) = Particles.getProp<3>(p1);
- }
- for(int j=0;j<dw_p1.size();j++)
- { auto p=dw_p.get<0>(j);
- auto p1=dw_p1.get<0>(j);
- Particles.getProp<3>(p)= Particles.getProp<3>(p1);
- if(j==0)
- { auto p=l_p.get<0>(0);
- auto p2=l_p.get<0>(1);
- Particles.getProp<3>(p) = Particles.getProp<3>(p1);
- Particles.getProp<3>(p2) = Particles.getProp<3>(p1);
- }
- if(j==dw_p1.size()-1)
- { auto p=r_p.get<0>(0);
- auto p2=r_p.get<0>(1);
- Particles.getProp<3>(p) = Particles.getProp<3>(p1);
- Particles.getProp<3>(p2) = Particles.getProp<3>(p1);
- }
- }*/
-
-/*
- auto Sf = getV<0>(Particles);
- auto W = getV<1>(Particles);
- auto dW = getV<2>(Particles);
- auto RHS = getV<3>(Particles);
- auto Wnew = getV<4>(Particles);
- auto V = getV<5>(Particles);
-*/
for(int j=0;j<up_p.size();j++) {
auto p = up_p.get<0>(j);
Particles.getProp<1>(p) = -12.0/spacing;
}
-// Particles.write_frame("Re1000-1e-3-Lid_sf",0);
-// Derivative_x Dx(Particles, 2, rCut,1);
-// Derivative_y Dy(Particles, 2, rCut,1);
-// Gradient Grad(Particles, 2, rCut,1);
-// Laplacian Lap(Particles, 2, rCut,1.9,support_options::RADIUS);
- Laplacian Lap(Particles, 2, rCut,1.9);
-// Curl2D Curl(Particles, 2, rCut, 1);
+ Laplacian Lap(Particles, ord2, rCut2,sampling,support_options::RADIUS);
+ //Gradient Grad(Particles, ord, rCut,sampling,support_options::RADIUS);
+ Derivative_xy Dxy(Particles,ord, 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);
+
auto its = Particles.getDomainIterator();
int ctr=0;
while (its.isNext()) {
auto p = its.get();
- Lap.DrawKernel<3>(Particles, p.getKey());
+ 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());
+ //Grad.DrawKernel<4>(Particles, p.getKey());
+
Particles.write_frame("LapKer",ctr);
for(int j=0;j<BULK.size();j++) {
auto p1 = BULK.get<0>(j);
+ Particles.getProp<0>(p1) = 0;
+ Particles.getProp<1>(p1) = 0;
+ Particles.getProp<2>(p1) = 0;
Particles.getProp<3>(p1) = 0;
+ Particles.getProp<4>(p1) = 0;
+ Particles.getProp<5>(p1)[0] = 0;
+ Particles.getProp<5>(p1)[1] = 0;
+
}
++its;
ctr++;
}
-/* double dt=0.003;
- double nu=0.01;
- int n=5;
- std::cout<<"Init Done"<<std::endl;
- for(int i =0; i<=n ;i++)
- { dW=Dx(Sf)*Dy(W)-Dy(Sf)*Dx(W)+nu*Lap(W);
- //Lap.DrawKernel<3>(Particles,837);
- Wnew=W+dt*dW;
- W=Wnew;
- //Copy Module
- for(int j=0;j<up_p1.size();j++)
- { auto p=up_p.get<0>(j);
- auto p1=up_p1.get<0>(j);
- Particles.getProp<1>(p) = -2.0*Particles.getProp<0>(p1)/(spacing*spacing)-12.0/spacing;
- if(j==0)
- { auto p=l_p.get<0>(l_p.size()-1);
- auto p2=l_p.get<0>(l_p.size()-2);
- //Particles.getProp<1>(p) = 0;//-2.0*Particles.getProp<0>(p1)/(spacing*spacing);
- Particles.getProp<1>(p2) = 0;//-2.0*Particles.getProp<0>(p1)/(spacing*spacing);
- }
- if(j==up_p1.size()-1)
- { auto p=r_p.get<0>(r_p.size()-1);
- auto p2=r_p.get<0>(r_p.size()-2);
- //Particles.getProp<1>(p) = 0;//-2.0*Particles.getProp<0>(p1)/(spacing*spacing);
- Particles.getProp<1>(p2) = 0;//-2.0*Particles.getProp<0>(p1)/(spacing*spacing);
- }
- }
- for(int j=0;j<l_p1.size();j++)
- {
- auto p=l_p.get<0>(j+2);
- auto p1=l_p1.get<0>(j);
- Particles.getProp<1>(p) = -2.0*Particles.getProp<0>(p1)/(spacing*spacing);
- }
- for(int j=0;j<r_p1.size();j++)
- {
- auto p=r_p.get<0>(j+2);
- auto p1=r_p1.get<0>(j);
- Particles.getProp<1>(p) = -2.0*Particles.getProp<0>(p1)/(spacing*spacing);
- }
- for(int j=0;j<dw_p1.size();j++)
- { auto p=dw_p.get<0>(j);
- auto p1=dw_p1.get<0>(j);
- Particles.getProp<1>(p) = -2.0*Particles.getProp<0>(p1)/(spacing*spacing);
- if(j==0)
- { auto p=l_p.get<0>(0);
- auto p2=l_p.get<0>(1);
- //Particles.getProp<1>(p) = 0;//-2.0*Particles.getProp<0>(p1)/(spacing*spacing);
- Particles.getProp<1>(p2) = 0;//-2.0*Particles.getProp<0>(p1)/(spacing*spacing);
- }
- if(j==dw_p1.size()-1)
- { auto p=r_p.get<0>(0);
- auto p2=r_p.get<0>(1);
- // Particles.getProp<1>(p) = 0;//-2.0*Particles.getProp<0>(p1)/(spacing*spacing);
- Particles.getProp<1>(p2) = 0;//-2.0*Particles.getProp<0>(p1)/(spacing*spacing);
- }
- }
- std::cout<<"W Done"<<std::endl;
-
- DCPSE_scheme<equations,decltype(Particles)> Solver(2*rCut, Particles);
- auto Sf_Poisson = -Lap(Sf);
- Solver.impose(Sf_Poisson, bulk, prop_id<1>());
- Solver.impose(Sf, up_p, 0);
- Solver.impose(Sf, dw_p,0);
- Solver.impose(Sf, l_p,0);
- Solver.impose(Sf, r_p, 0);
- Solver.solve(Sf);
- V=Curl(Sf);
- std::cout<<"Poisson Solved"<<std::endl;
- Particles.write_frame("Re1000-1e-3-Lid_sf",i);
- //if (i%10==0)
- std::cout<<i<<std::endl;
- }*/
+
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
BOOST_AUTO_TEST_CASE(dcpse_poisson_Robin_anal) {
diff --git a/src/DCPSE_op/DCPSE_op_test3d.cpp b/src/DCPSE_op/DCPSE_op_test3d.cpp
index 33985a2382a0ba845f04b6960cdd654aa5775c98..c22c283327f4b399c29e141fd72744aed3fa02fb 100644
--- a/src/DCPSE_op/DCPSE_op_test3d.cpp
+++ b/src/DCPSE_op/DCPSE_op_test3d.cpp
@@ -63,9 +63,11 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests3)
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);
@@ -174,6 +176,13 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests3)
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;
@@ -207,8 +216,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests3)
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);
- petsc_solver<double> solverPetsc;
- solverPetsc.setRestart(500);
+
Solver.solve_with_solver(solverPetsc,V_star[0],V_star[1],V_star[2]);
std::cout << "Stokes Solved" << std::endl;
Particles.ghost_get<2>();
@@ -225,8 +233,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests3)
SolverH.impose(H, l_p,0);
SolverH.impose(H, f_p,0);
SolverH.impose(H, b_p,0);
- petsc_solver<double> solverPetsc2;
- solverPetsc2.setRestart(500);
+
SolverH.solve_with_solver(solverPetsc2,H);
Particles.ghost_get<5>();
std::cout << "Helmholtz Solved" << std::endl;
diff --git a/src/DCPSE_op/EqnsStruct.hpp b/src/DCPSE_op/EqnsStruct.hpp
index 1f7132d55b082c1b8fae1cb95bc7614304ce8002..aaf2156540c566a8de53aad1d14720b555483374 100644
--- a/src/DCPSE_op/EqnsStruct.hpp
+++ b/src/DCPSE_op/EqnsStruct.hpp
@@ -228,6 +228,30 @@ struct equations2d2E {
typedef umfpack_solver<double> solver_type;
};
+struct equations2d3E {
+ //! dimensionaly of the equation ( 3D problem ...)
+ static const unsigned int dims = 2;
+ //! number of fields in the system
+ static const unsigned int nvar = 3;
+
+ //! 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;
+
+ typedef umfpack_solver<double> solver_type;
+};
+
struct equations2d1pE {
//! dimensionaly of the equation ( 3D problem ...)
diff --git a/src/DCPSE_op/Fast_test.cpp b/src/DCPSE_op/Fast_test.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b07bf218781ec5c008c5364061ff267d4c450cea
--- /dev/null
+++ b/src/DCPSE_op/Fast_test.cpp
@@ -0,0 +1,205 @@
+//
+// Created by Abhinav Singh on 25.04.20.
+//
+#include "config.h"
+
+#define BOOST_TEST_DYN_LINK
+
+#include "util/util_debug.hpp"
+#include <boost/test/unit_test.hpp>
+#include <iostream>
+#include "DCPSE_op.hpp"
+#include "DCPSE_Solver.hpp"
+#include "Operators/Vector/vector_dist_operators.hpp"
+#include "Vector/vector_dist_subset.hpp"
+#include "EqnsStruct.hpp"
+BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
+ BOOST_AUTO_TEST_CASE(dcpse_Kernels) {
+ 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);
+ double rCut =spacing*(2.8);
+ double ord = 2;
+ double sampling = 1.6;
+ double rCut2 = 3.1*spacing;
+ double ord2 = 2;
+ double sampling2 = 1.9;
+
+ Ghost<2, double> ghost(rCut2);
+ 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);
+
+ auto it = Particles.getGridIterator(sz);
+ while (it.isNext()) {
+ Particles.add();
+ auto key = it.get();
+ double x = key.get(0) * it.getSpacing(0);
+ Particles.getLastPos()[0] = x;
+ double y = key.get(1) * it.getSpacing(1);
+ Particles.getLastPos()[1] = y;
+ ++it;
+ }
+
+ Particles.map();
+ Particles.ghost_get<0>();
+
+ openfpm::vector<aggregate<int>> bulk;
+ openfpm::vector<aggregate<int>> up_p;
+ openfpm::vector<aggregate<int>> dw_p;
+ openfpm::vector<aggregate<int>> l_p;
+ openfpm::vector<aggregate<int>> r_p;
+
+ openfpm::vector<aggregate<int>> up_p1;
+ openfpm::vector<aggregate<int>> dw_p1;
+ openfpm::vector<aggregate<int>> l_p1;
+ openfpm::vector<aggregate<int>> r_p1;
+
+ openfpm::vector<aggregate<int>> BULK;
+
+
+
+ // Here fill up the boxes for particle detection.
+
+ Box<2, double> up({box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getHigh(1) + spacing / 2.0});
+
+ Box<2, double> up_d({box.getLow(0) + spacing / 2.0, box.getHigh(1) - 3 * spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
+ Box<2, double> down({box.getLow(0) + spacing / 2.0, box.getLow(1) - spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0});
+
+ Box<2, double> down_u({box.getLow(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getLow(1) + 3 * spacing / 2.0});
+
+ Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
+ {box.getLow(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
+
+ Box<2, double> left_r({box.getLow(0) + spacing / 2.0, box.getLow(1) + 3 * spacing / 2.0},
+ {box.getLow(0) + 3 * spacing / 2.0, box.getHigh(1) - 3 * spacing / 2.0});
+
+ Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
+
+ Box<2, double> right_l({box.getHigh(0) - 3 * spacing / 2.0, box.getLow(1) + 3 * spacing / 2.0},
+ {box.getHigh(0) - spacing / 2.0, box.getHigh(1) - 3 * spacing / 2.0});
+ Box<2, double> Bulk_box({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) +spacing / 2.0});
+
+
+ openfpm::vector<Box<2, double>> boxes;
+ boxes.add(up);
+ boxes.add(up_d);
+ boxes.add(down);
+ boxes.add(down_u);
+ boxes.add(left);
+ boxes.add(left_r);
+ boxes.add(right);
+ boxes.add(right_l);
+ boxes.add(Bulk_box);
+ VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
+ vtk_box.add(boxes);
+ vtk_box.write("box_sf.vtk");
+ // Particles.write_frame("Re1000-1e-3-Lid_sf",0);
+
+ auto it2 = Particles.getDomainIterator();
+
+ while (it2.isNext()) {
+ auto p = it2.get();
+ Point<2, double> xp = Particles.getPos(p);
+ Particles.getProp<0>(p) =0;//-M_PI*M_PI*sin(M_PI*xp[0])*sin(M_PI*xp[1]);
+ Particles.getProp<1>(p) =0;//sin(M_PI*xp[0])*sin(M_PI*xp[1]);
+ Particles.getProp<2>(p) =0.0;
+ Particles.getProp<3>(p) =0.0;
+ Particles.getProp<4>(p) =0.0;
+ Particles.getProp<5>(p)[0] =0.0;
+ Particles.getProp<5>(p)[1] =0.0;
+
+ BULK.add();
+ BULK.last().get<0>() = p.getKey();
+
+ if (up.isInside(xp) == true) {
+ up_p.add();
+ up_p.last().get<0>() = p.getKey();
+ }
+ else if (down.isInside(xp) == true) {
+ dw_p.add();
+ dw_p.last().get<0>() = p.getKey();
+ }
+ else if (left.isInside(xp) == true) {
+ l_p.add();
+ l_p.last().get<0>() = p.getKey();
+ } else if (right.isInside(xp) == true) {
+ r_p.add();
+ r_p.last().get<0>() = p.getKey();
+ }
+ else if (Bulk_box.isInside(xp) == true)
+ {
+ if (up_d.isInside(xp) == true) {
+ up_p1.add();
+ up_p1.last().get<0>() = p.getKey();
+ }
+ else if (down_u.isInside(xp) == true) {
+ dw_p1.add();
+ dw_p1.last().get<0>() = p.getKey();
+ }else if (left_r.isInside(xp) == true) {
+ l_p1.add();
+ l_p1.last().get<0>() = p.getKey();
+ } else if (right_l.isInside(xp) == true) {
+ r_p1.add();
+ r_p1.last().get<0>() = p.getKey();
+ }
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ }
+ ++it2;
+ }
+
+ for(int j=0;j<up_p.size();j++) {
+ auto p = up_p.get<0>(j);
+ Particles.getProp<1>(p) = -12.0/spacing;
+ }
+
+ Laplacian Lap(Particles, ord2, rCut2,sampling,support_options::RADIUS);
+ //Gradient Grad(Particles, ord, rCut,sampling,support_options::RADIUS);
+ //Derivative_xy Dxy(Particles,ord, 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);
+
+
+ auto its = Particles.getDomainIterator();
+ int ctr=0;
+ while (its.isNext()) {
+ auto p = its.get();
+ 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());
+ //Grad.DrawKernel<4>(Particles, p.getKey());
+
+ Particles.write_frame("LapKer",ctr);
+ for(int j=0;j<BULK.size();j++) {
+ auto p1 = BULK.get<0>(j);
+ Particles.getProp<0>(p1) = 0;
+ Particles.getProp<1>(p1) = 0;
+ Particles.getProp<2>(p1) = 0;
+ Particles.getProp<3>(p1) = 0;
+ Particles.getProp<4>(p1) = 0;
+ Particles.getProp<5>(p1)[0] = 0;
+ Particles.getProp<5>(p1)[1] = 0;
+
+ }
+ ++its;
+ ctr++;
+ }
+
+
+ }
+
+
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/src/Operators/Vector/vector_dist_operators.hpp b/src/Operators/Vector/vector_dist_operators.hpp
index 203656ea17ad5b6274ca005d05bc921293e44425..a3ba492b194b9dd50146924a1d2c688cff467d1c 100644
--- a/src/Operators/Vector/vector_dist_operators.hpp
+++ b/src/Operators/Vector/vector_dist_operators.hpp
@@ -922,6 +922,10 @@ class vector_dist_expression_op<exp1,boost::mpl::int_<n>,VECT_COMP>
int comp[n];
int var_id = 0;
+ void setVarId(int var_id)
+ {
+ this->var_id = var_id;
+ }
typedef vector_dist_expression_op<exp1,boost::mpl::int_<n>,VECT_COMP> myself;