From d4c1d1b712ce5d640780ec4a02113c6aba25d76c Mon Sep 17 00:00:00 2001
From: absingh <absingh@mpi-cbg.de>
Date: Mon, 6 Apr 2020 15:03:03 +0200
Subject: [PATCH] update
---
src/DCPSE_op/DCPSE_copy/DcpseRhs.hpp | 10 +--
src/DCPSE_op/DCPSE_copy/MonomialBasis.hpp | 4 +-
src/DCPSE_op/DCPSE_op_test.cpp | 102 +++++++++++++++++++++-
src/DCPSE_op/DCPSE_op_test2.cpp | 36 ++++++--
4 files changed, 135 insertions(+), 17 deletions(-)
diff --git a/src/DCPSE_op/DCPSE_copy/DcpseRhs.hpp b/src/DCPSE_op/DCPSE_copy/DcpseRhs.hpp
index 6e80bd6e..b3f1a3ca 100644
--- a/src/DCPSE_op/DCPSE_copy/DcpseRhs.hpp
+++ b/src/DCPSE_op/DCPSE_copy/DcpseRhs.hpp
@@ -50,11 +50,11 @@ MatrixType &DcpseRhs<dim>::getVector(MatrixType &b)
const Monomial<dim> dm = derivatives.getElement(i);
b(i, 0) = sign * dm.evaluate(Point<dim, T>(0));
}
- //Choosing a(0,0) as a free parameter can let us set b(0,0) for numerical robustness
- //if (b(0,0) == 0.0 && sign == 1)
- //{
- // b(0,0) = 100;
- //}
+ //Choosing a(0,0) for even order as a free parameter can let us set b(0,0) for numerical robustness
+/* if (b(0,0) == 0.0 && sign == 1)
+ {
+ b(0,0) = 100;
+ }*/
diff --git a/src/DCPSE_op/DCPSE_copy/MonomialBasis.hpp b/src/DCPSE_op/DCPSE_copy/MonomialBasis.hpp
index 8f6cf6c6..c9390a9d 100644
--- a/src/DCPSE_op/DCPSE_copy/MonomialBasis.hpp
+++ b/src/DCPSE_op/DCPSE_copy/MonomialBasis.hpp
@@ -124,9 +124,9 @@ void MonomialBasis<dim>::generateBasis(std::vector<unsigned int> m, unsigned int
grid_key_dx_iterator_sub_bc<dim> it(grid, start, stop, bc);
// Finally compute alpha_min
- //unsigned char alphaMin = static_cast<unsigned char>(!(mSum % 2)); // if mSum is even, alpha_min must be 1
+ unsigned char alphaMin = static_cast<unsigned char>(!(mSum % 2)); // if mSum is even, alpha_min must be 1
//std::cout<<"AlphaMin: "<<alphaMin<<std::endl;
- unsigned char alphaMin = 0; // we want to always have 1 in the basis
+ //unsigned char alphaMin = 0; // we want to always have 1 in the basis
while (it.isNext())
{
diff --git a/src/DCPSE_op/DCPSE_op_test.cpp b/src/DCPSE_op/DCPSE_op_test.cpp
index abaf96b3..f91d5e22 100644
--- a/src/DCPSE_op/DCPSE_op_test.cpp
+++ b/src/DCPSE_op/DCPSE_op_test.cpp
@@ -874,7 +874,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] = {150,150};
+ 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);
@@ -1316,7 +1316,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
//https://fenicsproject.org/docs/dolfin/1.4.0/python/demo/documented/neumann-poisson/python/documentation.html
// int rank;
// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
- const size_t sz[2] = {150,150};
+ const size_t sz[2] = {81,81};
Box<2, double> box({0, 0}, {1.0, 1.0});
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing = box.getHigh(0) / (sz[0] - 1);
@@ -2080,7 +2080,105 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
//auto flux = Dx(v) + v;
}
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ BOOST_AUTO_TEST_CASE(dcpse_op_test_lap) {
+// int rank;
+// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+ size_t edgeSemiSize = 160;
+ const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
+ Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
+ size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
+ double spacing[2];
+ spacing[0] = 2 * M_PI / (sz[0] - 1);
+ spacing[1] = 2 * M_PI / (sz[1] - 1);
+ Ghost<2, double> ghost(spacing[0] * 3);
+ double rCut = 2.0 * spacing[0];
+ BOOST_TEST_MESSAGE("Init vector_dist...");
+ double sigma2 = spacing[0] * spacing[1] / (2 * 4);
+
+ vector_dist<2, double, aggregate<double, double, double, double, VectorS<2, double>>> domain(0, box,
+ bc,
+ ghost);
+
+ //Init_DCPSE(domain)
+ BOOST_TEST_MESSAGE("Init domain...");
+// std::random_device rd{};
+// std::mt19937 rng{rd()};
+ std::mt19937 rng{6666666};
+
+ std::normal_distribution<> gaussian{0, sigma2};
+
+ auto it = domain.getGridIterator(sz);
+ size_t pointId = 0;
+ size_t counter = 0;
+ double minNormOne = 999;
+ while (it.isNext()) {
+ domain.add();
+ auto key = it.get();
+ mem_id k0 = key.get(0);
+ double x = k0 * spacing[0];
+ domain.getLastPos()[0] = x;//+ gaussian(rng);
+ mem_id k1 = key.get(1);
+ double y = k1 * spacing[1];
+ domain.getLastPos()[1] = y;//+gaussian(rng);
+ // Here fill the function value
+ domain.template getLastProp<0>() = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
+ // Here fill the validation value for Lap
+ domain.template getLastProp<1>() = - sin(domain.getLastPos()[0]) - sin(domain.getLastPos()[1]);
+ domain.template getLastProp<4>()[0] = -sin(domain.getLastPos()[0]);
+ domain.template getLastProp<4>()[1] = -sin(domain.getLastPos()[1]);
+
+// domain.template getLastProp<2>() = 2 * x;
+// domain.template getLastProp<2>() = 1;
+
+ ++counter;
+ ++it;
+ }
+ BOOST_TEST_MESSAGE("Sync domain across processors...");
+
+ domain.map();
+ domain.ghost_get<0>();
+
+ //Derivative_x Dx(domain, 2, rCut);
+ //Derivative_y Dy(domain, 2, rCut);
+ //Gradient Grad(domain, 2, rCut);
+ Laplacian Lap(domain, 2, rCut,1.9);
+ auto v = getV<1>(domain);
+ auto P = getV<0>(domain);
+ auto vv = getV<2>(domain);
+ auto errv = getV<3>(domain);
+
+
+
+ vv = Lap(P);
+ auto it2 = domain.getDomainIterator();
+
+ double worst = 0.0;
+
+ while (it2.isNext()) {
+ auto p = it2.get();
+
+ if (fabs(domain.getProp<1>(p) - domain.getProp<2>(p)) > worst) {
+ worst = fabs(domain.getProp<1>(p) - domain.getProp<2>(p));
+ }
+
+ ++it2;
+ }
+
+ std::cout << "Maximum Error: " << worst << std::endl;
+ errv=v-vv;
+ domain.deleteGhost();
+ domain.write("v");
+
+ //std::cout << demangle(typeid(decltype(v)).name()) << "\n";
+
+ //Debug<decltype(expr)> a;
+
+ //typedef decltype(expr)::blabla blabla;
+
+ //auto err = Dx + Dx;
+ }
BOOST_AUTO_TEST_CASE(dcpse_slice) {
const size_t sz[2] = {31,31};
diff --git a/src/DCPSE_op/DCPSE_op_test2.cpp b/src/DCPSE_op/DCPSE_op_test2.cpp
index 2d5f101c..9ff87917 100644
--- a/src/DCPSE_op/DCPSE_op_test2.cpp
+++ b/src/DCPSE_op/DCPSE_op_test2.cpp
@@ -76,7 +76,7 @@ const bool equations1d::boundary[] = {NON_PERIODIC, NON_PERIODIC};
BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
BOOST_AUTO_TEST_CASE(dcpse_Lid_Stokes) {
- const size_t sz[2] = {31,31};
+ const size_t sz[2] = {81,81};
Box<2, double> box({0, 0}, {1,1});
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing = box.getHigh(0) / (sz[0] - 1);
@@ -109,7 +109,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
auto V = getV<1>(Particles);
auto V_star = getV<2>(Particles);
auto RHS = getV<3>(Particles);
- auto V_BC = getV<4>(Particles);
+ auto V_t = getV<4>(Particles);
auto H = getV<5>(Particles);
auto temp=getV<6>(Particles);
@@ -172,7 +172,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
}
++it2;
}
- V_BC=V;
+ V_t=V;
eq_id vx,vy;
@@ -186,6 +186,26 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
Laplacian Lap(Particles, 2, rCut, 1.9);
Advection Adv(Particles, 2, rCut, 1.9);
Divergence Div(Particles, 2, rCut, 1.9);
+
+
+ //starting the simulation at a nice *continuous* place
+ V_t=1e-3*(1e-2*Lap(V)-Adv(V,V));
+ RHS=Div(V_t);
+ DCPSE_scheme<equations1d,decltype(Particles)> Solver(2 * rCut, Particles,options_solver::LAGRANGE_MULTIPLIER);
+ auto Pressure_Poisson = Lap(P);
+ auto D_y=Dy(P);
+ auto D_x=Dx(P);
+ Solver.impose(Pressure_Poisson,bulk,prop_id<3>());
+ Solver.impose(D_y, up_p,0);
+ Solver.impose(D_x, r_p, 0);
+ Solver.impose(-D_y, dw_p,0);
+ Solver.impose(-D_x, l_p,0);
+ Solver.solve(P);
+ std::cout << "Poisson Solved" << std::endl;
+ V_star = V + (V_t - 1e-3*Grad(P));
+ V = V_star;
+
+
int n=10;
double nu=1e-2;
Particles.write_frame("Stokes",0);
@@ -239,7 +259,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
BOOST_AUTO_TEST_CASE(dcpse_Lid_normal) {
- const size_t sz[2] = {31,31};
+ const size_t sz[2] = {81,81};
Box<2, double> box({0, 0}, {1,1});
size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
double spacing = box.getHigh(0) / (sz[0] - 1);
@@ -342,11 +362,11 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
Laplacian Lap(Particles, 2, rCut, 1.9);
Advection Adv(Particles, 2, rCut, 1.9);
Divergence Div(Particles, 2, rCut, 1.9);
- double dt=3e-3;
+ double dt=5e-4;
int n=50;
double nu=1e-2;
dV=dt*(nu*Lap(V)-Adv(V,V));
- DCPSE_scheme<equations2,decltype(Particles)> Solver(2 * rCut, Particles,options_solver::LAGRANGE_MULTIPLIER);
+ DCPSE_scheme<equations1d,decltype(Particles)> Solver(2 * rCut, Particles,options_solver::LAGRANGE_MULTIPLIER);
auto Pressure_Poisson = Lap(P);
auto D_y=Dy(P);
auto D_x=Dx(P);
@@ -364,7 +384,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
for(int i=1; i<=n ;i++)
{ dV=dt*(nu*Lap(V)-Adv(V,V));
RHS=1/dt*Div(dV);
- DCPSE_scheme<equations2,decltype(Particles)> Solver(2 * rCut, Particles,options_solver::LAGRANGE_MULTIPLIER);
+ DCPSE_scheme<equations1d,decltype(Particles)> Solver(2 * rCut, Particles,options_solver::LAGRANGE_MULTIPLIER);
auto Pressure_Poisson = Lap(P);
auto D_y=Dy(P);
auto D_x=Dx(P);
@@ -401,7 +421,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
Particles.write_frame("Re100-3e-3-Lid",i);
std::cout<<i<<std::endl;
if (i==100)
- dt=5e-4;
+ dt=1e-4;
}
}
--
GitLab