From 2f7ef00207e993b0ed71bcd7eafec8fbf4867906 Mon Sep 17 00:00:00 2001
From: absingh <absingh@mpi-cbg.de>
Date: Wed, 5 Feb 2020 11:53:43 +0100
Subject: [PATCH] Added DCPSE Scheme Solver
---
src/DCPSE_op/DCPSE_Solver.hpp | 100 ++++++++++-
src/DCPSE_op/DCPSE_copy/Dcpse.hpp | 132 ++++++++++++++
src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp | 2 +-
src/DCPSE_op/DCPSE_op.hpp | 33 +++-
src/DCPSE_op/DCPSE_op_test.cpp | 199 ++++++++++++++++-----
src/FiniteDifference/FDScheme.hpp | 2 +-
6 files changed, 408 insertions(+), 60 deletions(-)
diff --git a/src/DCPSE_op/DCPSE_Solver.hpp b/src/DCPSE_op/DCPSE_Solver.hpp
index bcd06035..82ecae7c 100644
--- a/src/DCPSE_op/DCPSE_Solver.hpp
+++ b/src/DCPSE_op/DCPSE_Solver.hpp
@@ -56,6 +56,10 @@ class DCPSE_scheme: public MatrixAssembler
//! row on b
size_t row_b;
+ //! Total number of points
+ size_t tot;
+
+
/*! \brief Construct the gmap structure
*
*/
@@ -75,7 +79,7 @@ class DCPSE_scheme: public MatrixAssembler
for (size_t i = 0 ; i < v_cl.getProcessUnitID() ; i++)
s_pnt += pnt.get(i);
- size_t tot = sz;
+ tot = sz;
v_cl.sum(tot);
v_cl.execute();
@@ -134,7 +138,7 @@ class DCPSE_scheme: public MatrixAssembler
* \return the scalar
*
*/
- typename Sys_eqs::stype get(size_t & key)
+ typename Sys_eqs::stype get(size_t key)
{
return scal;
}
@@ -167,12 +171,60 @@ class DCPSE_scheme: public MatrixAssembler
* \return the scalar
*
*/
- typename Sys_eqs::stype get(size_t & key)
+ inline typename Sys_eqs::stype get(size_t key)
{
return parts.template getProp<prp_id>(key);
}
};
+ /*! \brief Check if the Matrix is consistent
+ *
+ */
+ void consistency()
+ {
+ openfpm::vector<triplet> & trpl = A.getMatrixTriplets();
+
+ // 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";
+
+ // Indicate all the non zero rows
+ openfpm::vector<unsigned char> nz_rows;
+ nz_rows.resize(row_b);
+
+ for (size_t i = 0 ; i < trpl.size() ; i++)
+ nz_rows.get(trpl.get(i).row() - s_pnt*Sys_eqs::nvar) = true;
+
+ // Indicate all the non zero colums
+ // This check can be done only on single processor
+
+ Vcluster<> & v_cl = create_vcluster();
+ if (v_cl.getProcessingUnits() == 1)
+ {
+ openfpm::vector<unsigned> nz_cols;
+ nz_cols.resize(row_b);
+
+ for (size_t i = 0 ; i < trpl.size() ; i++)
+ nz_cols.get(trpl.get(i).col()) = true;
+
+ // all the rows must have a non zero element
+ for (size_t i = 0 ; i < nz_rows.size() ; i++)
+ {
+ if (nz_rows.get(i) == false)
+ {
+ std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix row " << i << " is not filled " << " equation: " << "\n";
+ }
+ }
+
+ // all the colums must have a non zero element
+ for (size_t i = 0 ; i < nz_cols.size() ; i++)
+ {
+ if (nz_cols.get(i) == false)
+ std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " Ill posed matrix colum " << i << " is not filled\n";
+ }
+ }
+ }
+
public:
/* void * dcpse_solver;
@@ -208,7 +260,7 @@ public:
void solve(expr_type exp)
{
umfpack_solver<double> solver;
- auto x = solver.solve(A,b);
+ auto x = solver.solve(getA(),getB());
auto parts = exp.getVector();
@@ -235,7 +287,7 @@ public:
*
*/
DCPSE_scheme(const typename Sys_eqs::stype r_cut, particles_type & part)
- :parts(part),p_map(0,part.getDecomposition().getDomain(),part.getDecomposition().periodicity(),Ghost<particles_type::dims,typename particles_type::stype>(r_cut))
+ :parts(part),p_map(0,part.getDecomposition().getDomain(),part.getDecomposition().periodicity(),Ghost<particles_type::dims,typename particles_type::stype>(r_cut)),row(0),row_b(0)
{
p_map.resize(part.size_local());
@@ -262,7 +314,7 @@ public:
const prop_id<prp_id> & num,
long int id = 0)
{
- auto itd = subset.getIterator();
+ auto itd = subset.template getIteratorElements<0>();
variable_b<prp_id> vb(parts);
@@ -289,13 +341,45 @@ public:
const typename Sys_eqs::stype num,
long int id = 0)
{
- auto itd = subset.getIterator();
+ auto itd = subset.template getIteratorElements<0>();
constant_b b(num);
impose_git(op,b,id,itd);
}
+ /*! \brief produce the Matrix
+ *
+ * \return the Sparse matrix produced
+ *
+ */
+ typename Sys_eqs::SparseMatrix_type & getA()
+ {
+#ifdef SE_CLASS1
+ consistency();
+#endif
+ A.resize(tot*Sys_eqs::nvar,tot*Sys_eqs::nvar,
+ p_map.size_local()*Sys_eqs::nvar,
+ p_map.size_local()*Sys_eqs::nvar);
+
+ return A;
+
+ }
+
+ /*! \brief produce the B vector
+ *
+ * \return the vector produced
+ *
+ */
+ typename Sys_eqs::Vector_type & getB()
+ {
+#ifdef SE_CLASS1
+ consistency();
+#endif
+
+ return b;
+ }
+
/*! \brief Impose an operator
*
* This function impose an operator on a particular grid region to produce the system
@@ -320,7 +404,7 @@ public:
auto it = it_d;
- std::unordered_map<long int,float> cols;
+ std::unordered_map<long int,typename particles_type::stype> cols;
// iterate all particles points
while (it.isNext())
diff --git a/src/DCPSE_op/DCPSE_copy/Dcpse.hpp b/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
index e5fb5d35..45802c18 100644
--- a/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
+++ b/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
@@ -60,6 +60,117 @@ public:
}
}
+ template<unsigned int prp>
+ void DrawKernel(vector_type &particles, int k)
+ {
+ EMatrix<T, Eigen::Dynamic, 1> &a = localCoefficients[k];
+ Support<dim, T, part_type> support = localSupports[k];
+ auto eps = localEps[k];
+
+ size_t xpK = k;
+ Point<dim, T> xp = support.getReferencePoint();
+ for (auto &xqK : support.getKeys())
+ {
+ Point<dim, T> xq = particles.getPos(xqK);
+ Point<dim, T> normalizedArg = (xp - xq) / eps;
+
+ particles.template getProp<prp>(xqK) += computeKernel(normalizedArg, a);
+ }
+ }
+
+ template<unsigned int prp>
+ void DrawKernel(vector_type &particles, int k, int i)
+ {
+ EMatrix<T, Eigen::Dynamic, 1> &a = localCoefficients[k];
+ Support<dim, T, part_type> support = localSupports[k];
+ auto eps = localEps[k];
+
+ size_t xpK = k;
+ Point<dim, T> xp = support.getReferencePoint();
+ for (auto &xqK : support.getKeys())
+ {
+ Point<dim, T> xq = particles.getPos(xqK);
+ Point<dim, T> normalizedArg = (xp - xq) / eps;
+
+ particles.template getProp<prp>(xqK)[i] += computeKernel(normalizedArg, a);
+ }
+ }
+
+ void checkMomenta(vector_type &particles)
+ {
+ openfpm::vector<aggregate<double,double>> momenta;
+ openfpm::vector<aggregate<double,double>> momenta_accu;
+
+ momenta.resize(monomialBasis.size());
+ momenta_accu.resize(monomialBasis.size());
+
+ for (int i = 0 ; i < momenta.size() ; i++)
+ {
+ momenta.template get<0>(i) = 3000000000.0;
+ momenta.template get<1>(i) = -3000000000.0;
+ }
+
+ auto it = particles.getDomainIterator();
+ auto coefficientsIt = localCoefficients.begin();
+ auto supportsIt = localSupports.begin();
+ auto epsIt = localEps.begin();
+ while (it.isNext())
+ {
+ double eps = *epsIt;
+
+ for (int i = 0 ; i < momenta.size() ; i++)
+ {
+ momenta_accu.template get<0>(i) = 0.0;
+ }
+
+ Support<dim, T, part_type> support = *supportsIt;
+ size_t xpK = support.getReferencePointKey();
+ Point<dim, T> xp = support.getReferencePoint();
+ for (auto &xqK : support.getKeys())
+ {
+ Point<dim, T> xq = particles.getPos(xqK);
+ Point<dim, T> normalizedArg = (xp - xq) / eps;
+ EMatrix<T, Eigen::Dynamic, 1> &a = *coefficientsIt;
+
+ int counter = 0;
+ for (const Monomial<dim> &m : monomialBasis.getElements())
+ {
+ T mbValue = m.evaluate(normalizedArg);
+
+
+ momenta_accu.template get<0>(counter) += mbValue * computeKernel(normalizedArg, a);
+
+ ++counter;
+ }
+
+ }
+
+ for (int i = 0 ; i < momenta.size() ; i++)
+ {
+ if (momenta_accu.template get<0>(i) < momenta.template get<0>(i))
+ {
+ momenta.template get<0>(i) = momenta_accu.template get<0>(i);
+ }
+
+ if (momenta_accu.template get<1>(i) > momenta.template get<1>(i))
+ {
+ momenta.template get<1>(i) = momenta_accu.template get<0>(i);
+ }
+ }
+
+ //
+ ++it;
+ ++coefficientsIt;
+ ++supportsIt;
+ ++epsIt;
+ }
+
+ for (int i = 0 ; i < momenta.size() ; i++)
+ {
+ std::cout << "MOMENTA: " << monomialBasis.getElements()[i] << "Min: " << momenta.template get<0>(i) << " " << "Max: " << momenta.template get<1>(i) << std::endl;
+ }
+ }
+
/**
* Computes the value of the differential operator on all the particles,
* using the f values stored at the fValuePos position in the aggregate
@@ -170,6 +281,24 @@ public:
return localSupports[key.getKey()].getKeys()[j];
}
+
+ T getSign()
+ {
+ T sign = 1.0;
+ if (differentialOrder % 2 == 0) {
+ sign = -1;
+ }
+
+ return sign;
+ }
+
+ T getEpsilonPrefactor(const vect_dist_key_dx &key)
+ {
+ double eps = localEps[key.getKey()];
+
+ return pow(eps, differentialOrder);
+ }
+
/**
* Computes the value of the differential operator on all the particles,
* using the f values stored at the fValuePos position in the aggregate
@@ -178,6 +307,9 @@ public:
* @tparam DfValuePos Position in the aggregate of the Df values to store.
* @param particles The set of particles to iterate over.
*/
+
+
+
template<typename op_type>
auto computeDifferentialOperator(const vect_dist_key_dx &key,
op_type &o1) -> decltype(is_scalar<std::is_fundamental<decltype(o1.value(
diff --git a/src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp b/src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp
index 385dc52d..01bb533d 100644
--- a/src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp
+++ b/src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp
@@ -152,7 +152,7 @@ std::vector<size_t> SupportBuilder<dim, T, Prop>::getPointsInSetOfCells(std::set
reord pr;
- pr.dist = xp.distance(xp);
+ pr.dist = xp.distance(xq);
pr.offset = el;
rp.add(pr);
diff --git a/src/DCPSE_op/DCPSE_op.hpp b/src/DCPSE_op/DCPSE_op.hpp
index 0a9770aa..4a6ef592 100644
--- a/src/DCPSE_op/DCPSE_op.hpp
+++ b/src/DCPSE_op/DCPSE_op.hpp
@@ -126,7 +126,8 @@ public:
o1.init();
}
- template<typename r_type=VectorS<dims,stype> > inline r_type value(const vect_dist_key_dx & key) const
+ template<typename r_type= typename std::remove_reference<decltype(o1.value(vect_dist_key_dx(0)))>::type>
+ inline r_type value(const vect_dist_key_dx & key) const
{
//typedef typename std::remove_reference<decltype(o1.value(key))>::type::blabla blabla;
@@ -153,7 +154,10 @@ public:
auto coeff_dc = dcp[i].getCoeffNN(key,j);
auto k = dcp[i].getIndexNN(key,j);
- cols[p_map. template getProp<0>(k)] += coeff_dc * coeff;
+
+ cols[p_map. template getProp<0>(k)] += coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
+
+ cols[p_map. template getProp<0>(key)] += dcp[i].getSign() * coeff_dc * coeff / dcp[i].getEpsilonPrefactor(key);
}
}
}
@@ -406,9 +410,34 @@ public:
return vector_dist_expression_op<operand_type1,std::pair<operand_type2,dcpse_type>,VECT_DCPSE_V_DOT>(arg,arg2,*(dcpse_type(*)[operand_type2::vtype::dims])dcpse);
}
+
+ template<typename particles_type>
+ void checkMomenta(particles_type &particles)
+ {
+ 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].checkMomenta(particles);
+ }
+
+ }
+ 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 Derivative_xy
{
diff --git a/src/DCPSE_op/DCPSE_op_test.cpp b/src/DCPSE_op/DCPSE_op_test.cpp
index 2edb5eb7..c047210c 100644
--- a/src/DCPSE_op/DCPSE_op_test.cpp
+++ b/src/DCPSE_op/DCPSE_op_test.cpp
@@ -51,7 +51,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
BOOST_AUTO_TEST_CASE(dcpse_op_solver) {
// int rank;
// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
- const size_t sz[2] = {100, 100};
+ const size_t sz[2] = {250, 250};
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);
@@ -59,7 +59,8 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
double rCut = 2.0 * spacing;
BOOST_TEST_MESSAGE("Init vector_dist...");
- vector_dist<2, double, aggregate<double,double,double>> domain(0, box, bc, ghost);
+ vector_dist<2, double, aggregate<double,double,double,double>> domain(0, box, bc, ghost);
+
//Init_DCPSE(domain)
BOOST_TEST_MESSAGE("Init domain...");
@@ -97,6 +98,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
auto v = getV<0>(domain);
auto sol = getV<2>(domain);
+ auto anasol = getV<3>(domain);
// Here fill me
@@ -123,14 +125,13 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
vtk_box.add(boxes);
vtk_box.write("vtk_box.vtk");
- domain.write("particles");
auto it2 = domain.getDomainIterator();
while (it2.isNext()) {
auto p = it2.get();
Point<2, double> xp = domain.getPos(p);
-
+ domain.getProp<3>(p)=1+xp[0]*xp[0]+2*xp[1]*xp[1];
if (up.isInside(xp) == true) {
up_p.add();
up_p.last().get<0>() = p.getKey();
@@ -146,7 +147,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
} else if (right.isInside(xp) == true) {
r_p.add();
r_p.last().get<0>() = p.getKey();
- domain.getProp<1>(p) = 1 + xp.get(1)*xp.get(1);
+ domain.getProp<1>(p) = 2 + 2*xp.get(1)*xp.get(1);
} else {
bulk.add();
bulk.last().get<0>() = p.getKey();
@@ -164,7 +165,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
- Solver.impose(eq1, bulk, 0);
+ Solver.impose(eq1, bulk, 6);
Solver.impose(v, up_p, prop_id<1>());
Solver.impose(v, dw_p, prop_id<1>());
Solver.impose(v, l_p, prop_id<1>());
@@ -172,53 +173,27 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Solver.solve(sol);
-
-/* P = Solver.solve(b);
-
- //auto v = petsc_solver(Matrix.getA(),Matrix.getB())
-
-
-
- auto v = getV<1>(domain);
- auto P = getV<0>(domain);
- auto dv = getV<3>(domain);
-
- vv=Sol_Lap(RHS_Vec);
- v2=Sol_Dx(RHS_Vec);
- //vv=Lap(P);
- //dv=Lap(v);//+Dy(P);
- dv=Adv(v,v);//+Dy(P);
- auto it2 = domain.getDomainIterator();
-
double worst1 = 0.0;
- while (it2.isNext())
- {
- auto p = it2.get();
-
- std::cout << "VALS: " << domain.getProp<3>(p)[0] << " " << domain.getProp<4>(p)[0] << std::endl;
- std::cout << "VALS: " << domain.getProp<3>(p)[1] << " " << domain.getProp<4>(p)[1] << std::endl;
+ it2 = domain.getDomainIterator();
- if (fabs(domain.getProp<3>(p)[0] - domain.getProp<4>(p)[0]) > worst1)
- {
- worst1 = fabs(domain.getProp<3>(p)[0] - domain.getProp<4>(p)[0]);
+ while (it2.isNext()) {
+ auto p = it2.get();
+ if (fabs(domain.getProp<3>(p) - domain.getProp<2>(p)) >= worst1) {
+ worst1 = fabs(domain.getProp<3>(p) - domain.getProp<2>(p));
}
+ domain.getProp<3>(p) = fabs(domain.getProp<3>(p) - domain.getProp<2>(p));
+
++it2;
}
std::cout << "Maximum Error: " << worst1 << std::endl;
- 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;*/
+ domain.write("particles");
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -337,7 +312,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
BOOST_TEST_MESSAGE("Init vector_dist...");
double sigma2 = spacing[0] * spacing[1] / (2 * 4);
- vector_dist<2, double, aggregate<double, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>>> domain(
+ vector_dist<2, double, aggregate<double, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>,double>> domain(
0, box, bc, ghost);
//Init_DCPSE(domain)
@@ -388,11 +363,11 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
domain.map();
domain.ghost_get<0>();
- Derivative_x Dx(domain, 2, rCut);
- Derivative_y Dy(domain, 2, rCut);
- Gradient Grad(domain, 2, rCut);
- Laplacian Lap(domain, 2, rCut, 3);
- Advection Adv(domain, 3, rCut, 3);
+ //Derivative_x Dx(domain, 2, rCut);
+ //Derivative_y Dy(domain, 2, rCut);
+ //Gradient Grad(domain, 2, rCut);
+ //Laplacian Lap(domain, 2, rCut, 3);
+ Advection Adv(domain, 3, rCut, 2);
auto v = getV<1>(domain);
auto P = getV<0>(domain);
auto dv = getV<3>(domain);
@@ -411,11 +386,14 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
while (it2.isNext()) {
auto p = it2.get();
- std::cout << "VALS: " << domain.getProp<3>(p)[0] << " " << domain.getProp<4>(p)[0] << std::endl;
- std::cout << "VALS: " << domain.getProp<3>(p)[1] << " " << domain.getProp<4>(p)[1] << std::endl;
+ //std::cout << "VALS: " << domain.getProp<3>(p)[0] << " " << domain.getProp<4>(p)[0] << std::endl;
+ //std::cout << "VALS: " << domain.getProp<3>(p)[1] << " " << domain.getProp<4>(p)[1] << std::endl;
+
+ domain.getProp<0>(p)=std::sqrt((domain.getProp<3>(p)[0] - domain.getProp<4>(p)[0])*(domain.getProp<3>(p)[0] - domain.getProp<4>(p)[0])+(domain.getProp<3>(p)[1] - domain.getProp<4>(p)[1])*(domain.getProp<3>(p)[1] - domain.getProp<4>(p)[1]));
if (fabs(domain.getProp<3>(p)[0] - domain.getProp<4>(p)[0]) > worst1) {
worst1 = fabs(domain.getProp<3>(p)[0] - domain.getProp<4>(p)[0]);
+
}
++it2;
@@ -423,6 +401,9 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
std::cout << "Maximum Error: " << worst1 << std::endl;
+ //Adv.checkMomenta(domain);
+ Adv.DrawKernel<2>(domain,0);
+
domain.deleteGhost();
domain.write("v");
@@ -534,6 +515,128 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
//auto err = Dx + Dx;
}
+ BOOST_AUTO_TEST_CASE(dcpse_test_diffusion) {
+ const size_t sz[2] = {100, 100};
+ 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);
+ Ghost<2, double> ghost(spacing * 3);
+ double rCut = 2.0 * spacing;
+ BOOST_TEST_MESSAGE("Init vector_dist...");
+
+ vector_dist<2, double, aggregate<double,double,double[2]>> domain(0, box, bc, ghost);
+
+ //Init_DCPSE(domain)
+ BOOST_TEST_MESSAGE("Init domain...");
+
+ auto it = domain.getGridIterator(sz);
+ while (it.isNext()) {
+ domain.add();
+
+ auto key = it.get();
+ double x = key.get(0) * it.getSpacing(0);
+ domain.getLastPos()[0] = x;
+ double y = key.get(1) * it.getSpacing(1);
+ domain.getLastPos()[1] = y;
+
+ ++it;
+ }
+ BOOST_TEST_MESSAGE("Sync domain across processors...");
+
+ domain.map();
+ domain.ghost_get<0>();
+
+ //Derivative_x Dx(domain, 2, rCut);
+ //Derivative_y Dy(domain, 2, rCut);
+ //Gradient Grad(domain, 2, rCut);
+ Laplacian Lap(domain, 2, rCut, 3);
+ //Advection Adv(domain, 3, rCut, 3);
+ //Solver Sol_Lap(Lap),Sol_Dx(Dx);
+ DCPSE_scheme<equations,decltype(domain)> Solver(2 * rCut, domain);
+
+ 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;
+
+ auto v = getV<0>(domain);
+ auto dc = getV<1>(domain);
+ // Here fill me
+
+ Box<2, double> up({box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
+
+ Box<2, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0});
+
+ Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
+ Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
+ {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
+
+ openfpm::vector<Box<2, double>> boxes;
+ boxes.add(up);
+ boxes.add(down);
+ boxes.add(left);
+ boxes.add(right);
+
+ // Create a writer and write
+ //VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
+ //vtk_box.add(boxes);
+ //vtk_box.write("vtk_box.vtk");
+
+ domain.write("Diffusion");
+
+ auto it2 = domain.getDomainIterator();
+
+ while (it2.isNext()) {
+ auto p = it2.get();
+ Point<2, double> xp = domain.getPos(p);
+
+ if (up.isInside(xp) == true) {
+ up_p.add();
+ up_p.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = 0;
+ } else if (down.isInside(xp) == true) {
+ dw_p.add();
+ dw_p.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = 0;
+ } else if (left.isInside(xp) == true) {
+ l_p.add();
+ l_p.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = 0;
+ } else if (right.isInside(xp) == true) {
+ r_p.add();
+ r_p.last().get<0>() = p.getKey();
+ domain.getProp<1>(p) = 0;
+ } else {
+ bulk.add();
+ bulk.last().get<0>() = p.getKey();
+ if(xp[0]==0.5 && xp[1]==0.5)
+ domain.getProp<1>(p) = 0;
+
+ }
+
+ ++it2;
+ }
+
+
+ double t=0;
+ while (t<2)
+ {
+ dc=Lap(v);
+ t=t+0.1;
+ v=dc;
+ domain.deleteGhost();
+ domain.write("Diffusion");
+ }
+
+ //auto flux = Dx(v) + v;
+
+ }
+
BOOST_AUTO_TEST_SUITE_END()
diff --git a/src/FiniteDifference/FDScheme.hpp b/src/FiniteDifference/FDScheme.hpp
index d86a4acf..38f67526 100644
--- a/src/FiniteDifference/FDScheme.hpp
+++ b/src/FiniteDifference/FDScheme.hpp
@@ -793,7 +793,7 @@ public:
const Box<Sys_eqs::dims,typename Sys_eqs::stype> & domain,
const typename Sys_eqs::b_grid & b_g)
:pd(pd),gs(b_g.getGridInfoVoid()),g_map(b_g,stencil,pd),row(0),row_b(0)
- {
+ {
Initialize(domain);
}
--
GitLab