diff --git a/mpi_include b/mpi_include
index 16981421116cb2ce7f9dcf3dba32df467827856c..4439b1ca3ff6f7e0872a08610a3d8f99ee1db3a9 100644
--- a/mpi_include
+++ b/mpi_include
@@ -1 +1 @@
--I
\ No newline at end of file
+-I/home/i-bird/MPI/include
\ No newline at end of file
diff --git a/mpi_libs b/mpi_libs
index 0519ecba6ea913e21689ec692e81e9e4973fbf73..3a37a5cd2100d56757bcde4ef6ae508b3f7c0c68 100644
--- a/mpi_libs
+++ b/mpi_libs
@@ -1 +1 @@
-
\ No newline at end of file
+-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
diff --git a/src/DCPSE_op/DCPSE_Solver.hpp b/src/DCPSE_op/DCPSE_Solver.hpp
index 301907a4089bfc35128efda86e11de025ebff6c7..ce939d080aa971b5d0229bb5f374b265892413d2 100644
--- a/src/DCPSE_op/DCPSE_Solver.hpp
+++ b/src/DCPSE_op/DCPSE_Solver.hpp
@@ -376,8 +376,8 @@ public:
* \param b_g object grid that will store the solution
*
*/
- DCPSE_scheme(const typename Sys_eqs::stype r_cut, particles_type & part, options_solver opt = options_solver::STANDARD)
- :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),opt(opt)
+ 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)
{
p_map.resize(part.size_local());
diff --git a/src/DCPSE_op/DCPSE_copy/Dcpse.hpp b/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
index 76d0c8c606198c9367b961d7b1507ebe3747c4bf..e4523a8e1617e16c9ba5a1b2efd2fcffcbdb948e 100644
--- a/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
+++ b/src/DCPSE_op/DCPSE_copy/Dcpse.hpp
@@ -59,7 +59,7 @@ private:
vector_type & particles;
- typename vector_type::stype rcut_verlet = -1.0;
+ support_options opt;
public:
@@ -70,12 +70,12 @@ public:
unsigned int convergenceOrder,
T rCut,
T supportSizeFactor = 1,
- T rcut_verlet = -1.0) :
- particles(particles),
+ support_options opt = support_options::N_PARTICLES)
+ :particles(particles),
differentialSignature(differentialSignature),
differentialOrder(Monomial<dim>(differentialSignature).order()),
monomialBasis(differentialSignature.asArray(), convergenceOrder),
- rcut_verlet(rcut_verlet)
+ opt(opt)
{
if (supportSizeFactor < 1) {
initializeAdaptive(particles, convergenceOrder, rCut);
@@ -411,7 +411,7 @@ private:
const T condVTOL = 1e2;
// Get the points in the support of the DCPSE kernel and store the support for reuse
- Support<dim, T, part_type> support = supportBuilder.getSupport(it, requiredSupportSize,rcut_verlet);
+ Support<dim, T, part_type> support = supportBuilder.getSupport(it, requiredSupportSize,opt);
EMatrix<T, Eigen::Dynamic, Eigen::Dynamic> V(support.size(), monomialBasis.size());
// Vandermonde matrix computation
@@ -469,7 +469,7 @@ private:
auto it = particles.getDomainIterator();
while (it.isNext()) {
// Get the points in the support of the DCPSE kernel and store the support for reuse
- Support<dim, T, part_type> support = supportBuilder.getSupport(it, requiredSupportSize,rcut_verlet);
+ Support<dim, T, part_type> support = supportBuilder.getSupport(it, requiredSupportSize,opt);
EMatrix<T, Eigen::Dynamic, Eigen::Dynamic> V(support.size(), monomialBasis.size());
/* Some Debug code
if (it.get().getKey() == 5564)
diff --git a/src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp b/src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp
index 3fdcc49689d338e13c6cee4cefa2b2382731624c..8dfd39a4836c0bee4c6c6503e92819a922595dfd 100644
--- a/src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp
+++ b/src/DCPSE_op/DCPSE_copy/SupportBuilder.hpp
@@ -12,6 +12,12 @@
#include <Vector/vector_dist.hpp>
#include "Support.hpp"
+enum support_options
+{
+ N_PARTICLES,
+ RADIUS
+};
+
template<typename vector_type>
class SupportBuilder
{
@@ -19,13 +25,14 @@ private:
vector_type &domain;
CellList<vector_type::dims, typename vector_type::stype, Mem_fast<HeapMemory, local_index_>> cellList;
const Point<vector_type::dims, unsigned int> differentialSignature;
+ typename vector_type::stype rCut;
public:
SupportBuilder(vector_type &domain, Point<vector_type::dims, unsigned int> differentialSignature, typename vector_type::stype rCut);
SupportBuilder(vector_type &domain, unsigned int differentialSignature[vector_type::dims], typename vector_type::stype rCut);
- Support<vector_type::dims, typename vector_type::stype, typename vector_type::value_type> getSupport(vector_dist_iterator itPoint, unsigned int requiredSize, typename vector_type::stype rcut_verlet);
+ Support<vector_type::dims, typename vector_type::stype, typename vector_type::value_type> getSupport(vector_dist_iterator itPoint, unsigned int requiredSize, support_options opt);
private:
size_t getCellLinId(const grid_key_dx<vector_type::dims> &cellKey);
@@ -36,7 +43,7 @@ private:
void enlargeSetOfCellsUntilSize(std::set<grid_key_dx<vector_type::dims>> &set, unsigned int requiredSize);
- std::vector<size_t> getPointsInSetOfCells(std::set<grid_key_dx<vector_type::dims>> set, vect_dist_key_dx & p, size_t requiredSupportSize, typename vector_type::stype rcut_verlet);
+ std::vector<size_t> getPointsInSetOfCells(std::set<grid_key_dx<vector_type::dims>> set, vect_dist_key_dx & p, size_t requiredSupportSize, support_options opt);
bool isCellKeyInBounds(grid_key_dx<vector_type::dims> key);
};
@@ -45,14 +52,17 @@ private:
template<typename vector_type>
SupportBuilder<vector_type>::SupportBuilder(vector_type &domain, const Point<vector_type::dims, unsigned int> differentialSignature,
- typename vector_type::stype rCut) : domain(domain), differentialSignature(differentialSignature)
+ typename vector_type::stype rCut)
+:domain(domain),
+ differentialSignature(differentialSignature),
+ rCut(rCut)
{
cellList = domain.template getCellList<CellList<vector_type::dims, typename vector_type::stype, Mem_fast<HeapMemory, local_index_>>>(rCut);
}
template<typename vector_type>
Support<vector_type::dims, typename vector_type::stype, typename vector_type::value_type> SupportBuilder<vector_type>
-::getSupport(vector_dist_iterator itPoint, unsigned int requiredSize, typename vector_type::stype rcut_verlet)
+::getSupport(vector_dist_iterator itPoint, unsigned int requiredSize, support_options opt)
{
// Get spatial position from point iterator
vect_dist_key_dx p = itPoint.get();
@@ -67,7 +77,7 @@ Support<vector_type::dims, typename vector_type::stype, typename vector_type::va
enlargeSetOfCellsUntilSize(supportCells, requiredSize + 1); // NOTE: this +1 is because we then remove the point itself
// Now return all the points from the support into a vector
- std::vector<size_t> supportKeys = getPointsInSetOfCells(supportCells,p,requiredSize,rcut_verlet);
+ std::vector<size_t> supportKeys = getPointsInSetOfCells(supportCells,p,requiredSize,opt);
std::remove(supportKeys.begin(), supportKeys.end(), p.getKey());
return Support<vector_type::dims, typename vector_type::stype, typename vector_type::value_type>(domain, p.getKey(), supportKeys);
}
@@ -127,7 +137,7 @@ template<typename vector_type>
std::vector<size_t> SupportBuilder<vector_type>::getPointsInSetOfCells(std::set<grid_key_dx<vector_type::dims>> set,
vect_dist_key_dx & p,
size_t requiredSupportSize,
- typename vector_type::stype rcut_verlet)
+ support_options opt)
{
struct reord
{
@@ -163,19 +173,13 @@ std::vector<size_t> SupportBuilder<vector_type>::getPointsInSetOfCells(std::set<
}
}
- if (p.getKey() == 174)
- {
- int debug = 0;
- debug++;
- }
-
rp.sort();
- if (rcut_verlet >= 0)
+ if (opt == support_options::RADIUS)
{
for (int i = 0 ; i < rp.size() ; i++)
{
- if (rp.get(i).dist <= rcut_verlet)
+ 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 e412df448db71a55170dc084dd2dca258860df0a..766497d78bd7d1a450115a066cc0e49804600327 100644
--- a/src/DCPSE_op/DCPSE_op.hpp
+++ b/src/DCPSE_op/DCPSE_op.hpp
@@ -663,13 +663,13 @@ class Derivative_x
public:
template<typename particles_type>
- Derivative_x(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, double rcut_verlet = -1.0)
+ Derivative_x(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, support_options opt = support_options::N_PARTICLES)
{
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(0) = 1;
- dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor, rcut_verlet);
+ dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor, opt);
}
template<typename operand_type>
@@ -691,13 +691,13 @@ class Derivative_y
public:
template<typename particles_type>
- Derivative_y(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, double rcut_verlet = -1.0)
+ Derivative_y(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, support_options opt = support_options::N_PARTICLES)
{
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(1) = 1;
- dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor, rcut_verlet);
+ dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor, opt);
}
template<typename operand_type>
@@ -717,13 +717,13 @@ class Derivative_z
public:
template<typename particles_type>
- Derivative_z(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, double rcut_verlet = -1.0)
+ Derivative_z(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, support_options opt = support_options::N_PARTICLES)
{
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(2) = 1;
- dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor, rcut_verlet);
+ dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor, opt);
}
template<typename operand_type>
@@ -748,7 +748,7 @@ class Gradient
public:
template<typename particles_type>
- Gradient(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, double rcut_verlet = -1.0)
+ Gradient(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, support_options opt = support_options::N_PARTICLES)
{
typedef Dcpse<particles_type::dims,particles_type> DCPSE_type;
@@ -761,7 +761,7 @@ public:
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(i) = 1;
- new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor, rcut_verlet);
+ new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor, opt);
dcpse_ptr++;
}
}
@@ -782,7 +782,7 @@ class Curl2D
public:
template<typename particles_type>
- Curl2D(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, double rcut_verlet = -1.0)
+ Curl2D(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, support_options opt = support_options::N_PARTICLES)
{
typedef Dcpse<particles_type::dims,particles_type> DCPSE_type;
@@ -792,11 +792,11 @@ public:
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(1) = 1;
- new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor, rcut_verlet);
+ new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor, opt);
dcpse_ptr++;
p.zero();
p.get(0) = 1;
- new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor, rcut_verlet);
+ new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor, opt);
dcpse_ptr++;
}
@@ -818,7 +818,7 @@ class Laplacian
public:
template<typename particles_type>
- Laplacian(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, double rcut_verlet = -1.0)
+ Laplacian(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, support_options opt = support_options::N_PARTICLES)
{
typedef Dcpse<particles_type::dims,particles_type> DCPSE_type;
@@ -831,7 +831,7 @@ public:
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(i) = 2;
- new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor,rcut_verlet);
+ new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor,opt);
dcpse_ptr++;
}
}
@@ -878,7 +878,7 @@ class Divergence
public:
template<typename particles_type>
- Divergence(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, double rcut_verlet = -1.0)
+ Divergence(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, support_options opt = support_options::N_PARTICLES)
{
typedef Dcpse<particles_type::dims,particles_type> DCPSE_type;
@@ -891,7 +891,7 @@ public:
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(i) = 1;
- new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor, rcut_verlet);
+ new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor, opt);
dcpse_ptr++;
}
}
@@ -914,7 +914,7 @@ class Advection
public:
template<typename particles_type>
- Advection(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, double rcut_verlet = -1.0)
+ Advection(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, support_options opt = support_options::N_PARTICLES)
{
typedef Dcpse<particles_type::dims,particles_type> DCPSE_type;
@@ -927,7 +927,7 @@ public:
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(i) = 1;
- new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor, rcut_verlet);
+ new (dcpse_ptr) Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,oversampling_factor, opt);
dcpse_ptr++;
}
@@ -977,14 +977,14 @@ class Derivative_xy
public:
template<typename particles_type>
- Derivative_xy(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, double rcut_verlet = -1.0)
+ Derivative_xy(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, support_options opt = support_options::N_PARTICLES)
{
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(0) = 1;
p.get(1) = 1;
- dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,dcpse_oversampling_factor, rcut_verlet);
+ dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,dcpse_oversampling_factor, opt);
}
template<typename operand_type>
@@ -1006,14 +1006,14 @@ class Derivative_xx
public:
template<typename particles_type>
- Derivative_xx(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, double rcut_verlet = -1.0)
+ Derivative_xx(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, support_options opt = support_options::N_PARTICLES)
{
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(0) = 2;
p.get(1) = 0;
- dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,dcpse_oversampling_factor, rcut_verlet);
+ dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,dcpse_oversampling_factor, opt);
}
template<typename operand_type>
@@ -1035,14 +1035,14 @@ class Derivative_yy
public:
template<typename particles_type>
- Derivative_yy(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, double rcut_verlet = -1.0)
+ Derivative_yy(particles_type & parts, unsigned int ord ,typename particles_type::stype rCut,double oversampling_factor=dcpse_oversampling_factor, support_options opt = support_options::N_PARTICLES)
{
Point<particles_type::dims,unsigned int> p;
p.zero();
p.get(0) = 0;
p.get(1) = 2;
- dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,dcpse_oversampling_factor, rcut_verlet);
+ dcpse = new Dcpse<particles_type::dims,particles_type>(parts,p, ord, rCut,dcpse_oversampling_factor, opt);
}
template<typename operand_type>
diff --git a/src/DCPSE_op/DCPSE_op_test.cpp b/src/DCPSE_op/DCPSE_op_test.cpp
index 1e15b935451d75007785ab319b771e7c0f3a1a85..fa9d04aef010c5b5b0641966c9c757952b9dfdce 100644
--- a/src/DCPSE_op/DCPSE_op_test.cpp
+++ b/src/DCPSE_op/DCPSE_op_test.cpp
@@ -64,52 +64,6 @@ struct equations2d1 {
typedef Vector<double> Vector_type;
};
-
-//! Specify the general characteristic of system to solve
-struct equationsp {
- //! dimensionaly of the equation ( 3D problem ...)
- static const unsigned int dims = 2;
- //! number of fields in the system
- static const unsigned int nvar = 1;
-
- //! 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;
-};
-
-struct equations2dp {
- //! dimensionaly of the equation ( 3D problem ...)
- static const unsigned int dims = 2;
- //! number of fields in the system
- static const unsigned int nvar = 2;
-
- //! 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;
-};
-
//Change accordingly as per test
const bool equations2d1::boundary[] = {NON_PERIODIC, NON_PERIODIC};
const bool equations::boundary[] = {NON_PERIODIC, NON_PERIODIC};
@@ -194,15 +148,15 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
double lambda = 0.1;
double delmu = -1;
- Derivative_x Dx(Particles, 2, rCut,1.9,3.1*spacing);
- Derivative_y Dy(Particles, 2, rCut,1.9,3.1*spacing);
- Derivative_xy Dxy(Particles, 2, rCut,1.9,3.1*spacing);
- Derivative_xx Dxx(Particles, 2, rCut,1.9,3.1*spacing);
- Derivative_yy Dyy(Particles, 2, rCut,1.9,3.1*spacing);
- Gradient Grad(Particles, 2, rCut,1.9,3.1*spacing);
- Laplacian Lap(Particles, 2, rCut,1.9,3.1*spacing);
- Advection Adv(Particles, 2, rCut,1.9,3.1*spacing);
- Divergence Div(Particles, 2, rCut,1.9,3.1*spacing);
+ Derivative_x Dx(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS);
+ Derivative_y Dy(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS);
+ Derivative_xy Dxy(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS);
+ Derivative_xx Dxx(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS);
+ Derivative_yy Dyy(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS);
+ Gradient Grad(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS);
+ Laplacian Lap(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS);
+ Advection Adv(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS);
+ Divergence Div(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS);
// Here fill up the boxes for particle boundary detection.
@@ -309,7 +263,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
for(int i=1; i<=n ;i++)
{ RHS[x]=Grad(P)+dV[x];
RHS[y]=Grad(P)+dV[y];
- DCPSE_scheme<equations2d1,decltype(Particles)> Solver(2 * rCut, Particles);
+ DCPSE_scheme<equations2d1,decltype(Particles)> Solver( Particles);
auto Stokes1 = nu*Lap(V[x]) + 0.5*nu*(Dx(f1)+Dx(f2)+Dx(f3)+Dy(f4)+Dy(f5)+Dy(f6));
auto Stokes2 = nu*Lap(V[y]) + 0.5*nu*(Dx(f1)+Dx(f2)+Dx(f3)+Dy(f4)+Dy(f5)+Dy(f6));
Solver.impose(Stokes1,bulk,RHS[0],vx);
@@ -394,7 +348,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
std::cout<<"dV Done"<<std::endl;
RHS = Div(dV);
std::cout<<"RHS Done"<<std::endl;
- DCPSE_scheme<equations,decltype(Particles)> Solver(2 * rCut, Particles);
+ DCPSE_scheme<equations,decltype(Particles)> Solver( Particles);
auto Pressure_Poisson = Lap(H);
auto D_x = Dx(H);
auto D_y = Dy(H);
@@ -617,7 +571,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Particles.template getProp<9>(bulk.template get<0>(i)) = Particles_subset.getProp<0>(i);
}
//subset_create<0,1,2,4>(Particles,Particles_subset,bulk);
- DCPSE_scheme<equations,decltype(Particles_subset)> Solver(2*rCut, Particles_subset);
+ DCPSE_scheme<equations,decltype(Particles_subset)> Solver(Particles_subset);
auto Pressure_Poisson = Lap_sub(P_bulk);
Solver.impose(Pressure_Poisson, bulk_1,prop_id<0>());
Solver.impose(P_bulk, ref_bulk, 1);
@@ -661,7 +615,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
}
//RHS=1.0/dt*Div(dV);
std::cout<<"RHS Done"<<std::endl;
- DCPSE_scheme<equations,decltype(Particles)> Solver(2*rCut, Particles);
+ DCPSE_scheme<equations,decltype(Particles)> Solver(Particles);
auto Pressure_Poisson = Lap(P);
auto D_x = Dx(P);
auto D_y = Dy(P);
@@ -949,7 +903,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
// 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,3.1*spacing);
+ Laplacian Lap(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS);
// Curl2D Curl(Particles, 2, rCut, 1);
auto its = Particles.getDomainIterator();
@@ -1074,10 +1028,10 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
domain.map();
domain.ghost_get<0>();
- Derivative_x Dx(domain, 2, rCut,1.9,3.1*spacing);
- Derivative_y Dy(domain, 2, rCut,1.9,3.1*spacing);
+ Derivative_x Dx(domain, 2, 3.1*spacing,1.9,support_options::RADIUS);
+ Derivative_y Dy(domain, 2, 3.1*spacing,1.9,support_options::RADIUS);
//Gradient Grad(domain, 2, rCut);
- Laplacian Lap(domain, 2, rCut, 1.9,3.1*spacing);
+ Laplacian Lap(domain, 2, 3.1*spacing,1.9,support_options::RADIUS);
//Advection Adv(domain, 3, rCut, 3);
//Solver Sol_Lap(Lap),Sol_Dx(Dx);
@@ -1159,7 +1113,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
}
++it2;
}
- DCPSE_scheme<equations,decltype(domain)> Solver(2 * rCut, domain);
+ DCPSE_scheme<equations,decltype(domain)> Solver( domain);
auto Poisson = Lap(v);
auto D_x = Dx(v);
auto D_y = Dy(v);
@@ -1307,7 +1261,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
}
++it2;
}
- DCPSE_scheme<equations,decltype(domain)> Solver(2 * rCut, domain);
+ DCPSE_scheme<equations,decltype(domain)> Solver( domain);
auto Poisson = Lap(v);
auto D_x = Dx(v);
auto D_y = Dy(v);
@@ -1344,8 +1298,8 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Box<2, double> box({0, 0}, {1.0, 1.0});
size_t bc[2] = {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 = 2.0 * spacing;
BOOST_TEST_MESSAGE("Init vector_dist...");
vector_dist<2, double, aggregate<double,double,double,double,double,VectorS<2, double>>> domain(0, box, bc, ghost);
@@ -1371,9 +1325,9 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
domain.ghost_get<0>();
- Laplacian Lap(domain, 2, rCut, 1.9,3.1*spacing);
+ Laplacian Lap(domain, 2, 3.1*spacing, 1.9, support_options::RADIUS);
- DCPSE_scheme<equationsp,decltype(domain)> Solver(2 * rCut, domain);
+ DCPSE_scheme<equations,decltype(domain)> Solver( domain);
openfpm::vector<aggregate<int>> bulk;
openfpm::vector<aggregate<int>> up_p;
@@ -1443,6 +1397,8 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Solver.impose(v, up_p, 0);
Solver.impose(v, dw_p, 0);
Solver.solve(v);
+
+ domain.ghost_get<0>();
anasol=-Lap(v);
double worst1 = 0.0;
@@ -1455,7 +1411,6 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
}
std::cout << "Maximum Auto Error: " << worst1 << std::endl;
- domain.ghost_get<0>();
domain.write("Poisson_Periodic");
}
@@ -1502,8 +1457,8 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
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);
- DCPSE_scheme<equations,decltype(domain)> Solver(2 * rCut, domain);
+ //DCPSE_scheme<equations,decltype(domain)> Solver( domain);
+ DCPSE_scheme<equations,decltype(domain)> Solver( domain);
openfpm::vector<aggregate<int>> bulk;
openfpm::vector<aggregate<int>> up_p;
@@ -1646,8 +1601,8 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Laplacian Lap(domain, 3, rCut, 2);
//Advection Adv(domain, 3, rCut, 3);
//Solver Sol_Lap(Lap),Sol_Dx(Dx);
- //DCPSE_scheme<equations,decltype(domain)> Solver(2 * rCut, domain);
- DCPSE_scheme<equations,decltype(domain)> Solver(2 * rCut, domain,options_solver::LAGRANGE_MULTIPLIER);
+ //DCPSE_scheme<equations,decltype(domain)> Solver( domain);
+ DCPSE_scheme<equations,decltype(domain)> Solver( domain,options_solver::LAGRANGE_MULTIPLIER);
openfpm::vector<aggregate<int>> bulk;
openfpm::vector<aggregate<int>> up_p;
@@ -1786,7 +1741,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
Laplacian Lap(domain, 2, rCut, 2);
//Advection Adv(domain, 3, rCut, 3);
//Solver Sol_Lap(Lap),Sol_Dx(Dx);
- DCPSE_scheme<equations,decltype(domain)> Solver(2 * rCut, domain);
+ DCPSE_scheme<equations,decltype(domain)> Solver( domain);
openfpm::vector<aggregate<int>> bulk;
openfpm::vector<aggregate<int>> up_p;
@@ -2289,7 +2244,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
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);
+ DCPSE_scheme<equations,decltype(domain)> Solver( domain);
openfpm::vector<aggregate<int>> bulk;
openfpm::vector<aggregate<int>> up_p;
@@ -2593,10 +2548,10 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
domain.map();
domain.ghost_get<0>();
- Derivative_x Dx(domain, 2, rCut,1.9,3.1*spacing);
- Derivative_y Dy(domain, 2, rCut,1.9,3.1*spacing);
+ Derivative_x Dx(domain, 2, 3.1*spacing,1.9,support_options::RADIUS);
+ Derivative_y Dy(domain, 2, 3.1*spacing,1.9,support_options::RADIUS);
//Gradient Grad(domain, 2, rCut);
- Laplacian Lap(domain, 2, rCut, 1.9,3.1*spacing);
+ Laplacian Lap(domain, 2, 3.1*spacing,1.9,support_options::RADIUS);
//Advection Adv(domain, 3, rCut, 3);
//Solver Sol_Lap(Lap),Sol_Dx(Dx);
@@ -2700,7 +2655,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
vx.setId(0);
vy.setId(1);
- DCPSE_scheme<equations2d1,decltype(domain)> Solver(2 * rCut, domain);
+ DCPSE_scheme<equations2d1,decltype(domain)> Solver( domain);
auto Poisson0 = Lap(v[0]);
auto Poisson1 = Lap(v[1]);
//auto D_x = Dx(v[1]);
diff --git a/src/DCPSE_op/DCPSE_op_test2.cpp b/src/DCPSE_op/DCPSE_op_test2.cpp
index a1114a624a124faf075c84852c556d42594a805b..c1bdb9c02d16f68f0dc4eaf291b617ab5922f630 100644
--- a/src/DCPSE_op/DCPSE_op_test2.cpp
+++ b/src/DCPSE_op/DCPSE_op_test2.cpp
@@ -179,18 +179,18 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
vx.setId(0);
vy.setId(1);
- Derivative_x Dx(Particles, 2, rCut,1.9,3.1*spacing );
- Derivative_y Dy(Particles, 2, rCut,1.9,3.1*spacing);
- Gradient Grad(Particles, 2, rCut,1.9,3.1*spacing );
- Laplacian Lap(Particles, 2, rCut, 1.9,3.1*spacing);
- Advection Adv(Particles, 2, rCut, 1.9,3.1*spacing);
- Divergence Div(Particles, 2, rCut, 1.9,3.1*spacing);
+ Derivative_x Dx(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS );
+ Derivative_y Dy(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS);
+ Gradient Grad(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS );
+ Laplacian Lap(Particles, 2, 3.1*spacing, 1.9,support_options::RADIUS);
+ Advection Adv(Particles, 2, 3.1*spacing, 1.9,support_options::RADIUS);
+ Divergence Div(Particles, 2, 3.1*spacing, 1.9,support_options::RADIUS);
/* 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);
+ DCPSE_scheme<equations1d,decltype(Particles)> Solver( Particles,options_solver::LAGRANGE_MULTIPLIER);
auto Pressure_Poisson = Lap(P);
auto D_y=Dy(P);
auto D_x=Dx(P);
@@ -210,7 +210,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
Particles.write_frame("Stokes",0);
for(int i=1; i<=n ;i++)
{ RHS=-Grad(P);
- DCPSE_scheme<equations2,decltype(Particles)> Solver(2 * rCut, Particles);
+ DCPSE_scheme<equations2,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]);
Solver.impose(Stokes1,bulk,RHS[0],vx);
@@ -226,7 +226,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
Solver.solve(V_star[0],V_star[1]);
//std::cout << "Stokes Solved" << std::endl;
RHS=Div(V_star);
- DCPSE_scheme<equations1d,decltype(Particles)> SolverH(2 * rCut, Particles,options_solver::LAGRANGE_MULTIPLIER);
+ DCPSE_scheme<equations1d,decltype(Particles)> SolverH( Particles,options_solver::LAGRANGE_MULTIPLIER);
auto Helmholtz = Lap(H);
auto D_y=Dy(H);
auto D_x=Dx(H);
@@ -379,18 +379,18 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
vx.setId(0);
vy.setId(1);
- Derivative_x Dx(Particles, 2, rCut,1.9,3.1*spacing );
- Derivative_y Dy(Particles, 2, rCut,1.9,3.1*spacing);
- Gradient Grad(Particles, 2, rCut,1.9,3.1*spacing );
- Laplacian Lap(Particles, 2, rCut, 1.9,3.1*spacing);
- Advection Adv(Particles, 2, rCut, 1.9,3.1*spacing);
- Divergence Div(Particles, 2, rCut, 1.9,3.1*spacing);
+ Derivative_x Dx(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS );
+ Derivative_y Dy(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS);
+ Gradient Grad(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS );
+ Laplacian Lap(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS);
+ Advection Adv(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS);
+ Divergence Div(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS);
/* 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);
+ DCPSE_scheme<equations1d,decltype(Particles)> Solver( Particles,options_solver::LAGRANGE_MULTIPLIER);
auto Pressure_Poisson = Lap(P);
auto D_y=Dy(P);
auto D_x=Dx(P);
@@ -410,7 +410,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
Particles.write_frame("Stokes",0);
for(int i=1; i<=n ;i++)
{ RHS=-Grad(P);
- DCPSE_scheme<equations2,decltype(Particles)> Solver(2 * rCut, Particles);
+ DCPSE_scheme<equations2,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]);
Solver.impose(Stokes1,bulk,RHS[0],vx);
@@ -426,7 +426,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
Solver.solve(V_star[0],V_star[1]);
//std::cout << "Stokes Solved" << std::endl;
RHS=Div(V_star);
- DCPSE_scheme<equations1d,decltype(Particles)> SolverH(2 * rCut, Particles,options_solver::LAGRANGE_MULTIPLIER);
+ DCPSE_scheme<equations1d,decltype(Particles)> SolverH( Particles,options_solver::LAGRANGE_MULTIPLIER);
auto Helmholtz = Lap(H);
auto D_y=Dy(H);
auto D_x=Dx(H);
@@ -583,17 +583,17 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
++it2;
}
- Derivative_x Dx(Particles, 2, rCut,1.9,3.1*spacing);
- Derivative_y Dy(Particles, 2, rCut,1.9,3.1*spacing);
- Gradient Grad(Particles, 2, rCut,1.9,3.1*spacing);
- Laplacian Lap(Particles, 2, rCut, 1.9,3.1*spacing);
- Advection Adv(Particles, 2, rCut, 1.9,3.1*spacing);
- Divergence Div(Particles, 2, rCut, 1.9,3.1*spacing);
+ Derivative_x Dx(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS);
+ Derivative_y Dy(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS);
+ Gradient Grad(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS);
+ Laplacian Lap(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS);
+ Advection Adv(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS);
+ Divergence Div(Particles, 2, 3.1*spacing,1.9,support_options::RADIUS);
double dt=1e-3;
int n=50;
double nu=1e-2;
dV=dt*(nu*Lap(V)-Adv(V,V));
- DCPSE_scheme<equations1d,decltype(Particles)> Solver(2 * rCut, Particles,options_solver::LAGRANGE_MULTIPLIER);
+ DCPSE_scheme<equations1d,decltype(Particles)> Solver( Particles,options_solver::LAGRANGE_MULTIPLIER);
auto Pressure_Poisson = Lap(P);
auto D_y=Dy(P);
auto D_x=Dx(P);
@@ -611,7 +611,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<equations1d,decltype(Particles)> Solver(2 * rCut, Particles,options_solver::LAGRANGE_MULTIPLIER);
+ DCPSE_scheme<equations1d,decltype(Particles)> Solver( Particles,options_solver::LAGRANGE_MULTIPLIER);
auto Pressure_Poisson = Lap(P);
auto D_y=Dy(P);
auto D_x=Dx(P);
@@ -1958,7 +1958,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests2)
vector_dist_expression_op<void, void, VECT_COPY_1_TO_N> vcp_dl(LL_p, LL_p_ref.get<0>(0));
vector_dist_expression_op<void, void, VECT_COPY_1_TO_N> vcp_dr(LR_p, LR_p_ref.get<0>(0));
- DCPSE_scheme<equations2, decltype(Particles)> Solver(2 * rCut, Particles);
+ DCPSE_scheme<equations2, decltype(Particles)> Solver( Particles);
auto Pressure_Poisson = Lap(P);
Solver.impose(Pressure_Poisson, bulk, prop_id<3>());
Solver.impose(vcp_up, up_p, 0);
diff --git a/src/DCPSE_op/DCPSE_op_test3d.cpp b/src/DCPSE_op/DCPSE_op_test3d.cpp
index 6d633eaac1007f7f0db54af7119d62e1f3f09f87..483a79773520402e3b8f15e016b37f99bc925331 100644
--- a/src/DCPSE_op/DCPSE_op_test3d.cpp
+++ b/src/DCPSE_op/DCPSE_op_test3d.cpp
@@ -75,7 +75,7 @@ const bool equations3d2::boundary[] = {NON_PERIODIC, NON_PERIODIC};
BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests3)
- BOOST_AUTO_TEST_CASE(dcpse_sphere) {
+/* BOOST_AUTO_TEST_CASE(dcpse_sphere) {
const size_t sz[3] = {31,31,31};
Sphere<3, double> sphere(0, 5);
double spacing = 0.1;
@@ -108,7 +108,7 @@ BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests3)
Particles.ghost_get<0>();
Particles.write_frame("Sphere",i);\
- }
+ }*/
BOOST_AUTO_TEST_SUITE_END()