From 528ec2063faf7fe233ad56b8b9f5532a21637188 Mon Sep 17 00:00:00 2001
From: Pietro Incardona <incardon@mpi-cbg.de>
Date: Fri, 27 Jan 2017 22:22:05 +0100
Subject: [PATCH] examples SPH update
---
example/Vector/0_simple/main.cpp | 1 +
.../4_multiphase_celllist_verlet/main.cpp | 4 +-
example/Vector/7_SPH_dlb/main.cpp | 778 ++++++++++++------
openfpm_numerics | 2 +-
openfpm_pdata.doc | 2 +-
src/Decomposition/CartDecomposition.hpp | 7 +-
.../Distribution/ParMetisDistribution.hpp | 29 +-
src/Graph/CartesianGraphFactory.hpp | 23 +-
src/Graph/ids.hpp | 17 +
src/Makefile.am | 2 +-
src/Vector/vector_dist.hpp | 20 +-
src/Vector/vector_dist_MP_unit_tests.hpp | 16 +-
src/Vector/vector_dist_comm.hpp | 23 +-
.../vector_dist_multiphase_functions.hpp | 30 +-
src/Vector/vector_dist_unit_test.hpp | 2 +-
15 files changed, 667 insertions(+), 289 deletions(-)
diff --git a/example/Vector/0_simple/main.cpp b/example/Vector/0_simple/main.cpp
index a2c9b7b2..b17023e1 100644
--- a/example/Vector/0_simple/main.cpp
+++ b/example/Vector/0_simple/main.cpp
@@ -11,6 +11,7 @@
* \subpage Vector_5_md_vl_sym
* \subpage Vector_5_md_vl_sym_crs
* \subpage Vector_6_complex_usage
+ * \subpage Vector_7_sph_dlb
*
*/
diff --git a/example/Vector/4_multiphase_celllist_verlet/main.cpp b/example/Vector/4_multiphase_celllist_verlet/main.cpp
index e426b889..65b92e93 100644
--- a/example/Vector/4_multiphase_celllist_verlet/main.cpp
+++ b/example/Vector/4_multiphase_celllist_verlet/main.cpp
@@ -317,7 +317,7 @@ int main(int argc, char* argv[])
CL_phase1 = phases.get(1).getCellListSym(r_cut);
// This function create a Verlet-list between phases 0 and 1
- NN_ver01 = createVerletSym(phases.get(0),CL_phase1,r_cut);
+ NN_ver01 = createVerletSym(phases.get(0),phases.get(1),CL_phase1,r_cut);
// Get an iterator over the real and ghost particles
it = phases.get(0).getDomainAndGhostIterator();
@@ -419,7 +419,7 @@ int main(int argc, char* argv[])
CL_all = createCellListSymM<2>(phases,r_cut);
// Type of the multiphase Verlet-list
- typedef decltype(createVerletSymM<2>(phases.get(0),CL_all,r_cut)) verlet_type;
+ typedef decltype(createVerletSymM<2>(phases.get(0),phases,CL_all,r_cut)) verlet_type;
// for each phase we create one Verlet-list that contain the neighborhood
// from all the phases
diff --git a/example/Vector/7_SPH_dlb/main.cpp b/example/Vector/7_SPH_dlb/main.cpp
index 9a6dfc1b..2c6f10ce 100644
--- a/example/Vector/7_SPH_dlb/main.cpp
+++ b/example/Vector/7_SPH_dlb/main.cpp
@@ -1,5 +1,5 @@
/*!
- * \page Vector_7_sph_dlb Vector 7 SPH Dam break simulation with Dynamic load balacing
+ * \page Vector_7_sph_dlb Vector 7 SPH Dam break simulation with Dynamic load balacing
*
*
* [TOC]
@@ -8,11 +8,15 @@
* # SPH with Dynamic load Balancing # {#SPH_dlb}
*
*
- * This example show the classical SPH Dam break simulation with Dynamic load balancing
+ * This example show the classical SPH Dam break simulation with Load Balancing and Dynamic load balancing. With
+ * Load balancing and Dynamic load balancing we indicate the possibility of the system to re-adapt the domain
+ * decomposition to keep all the processor load and reduce idle time.
*
* ## inclusion ## {#e0_v_inclusion}
*
* In order to use distributed vectors in our code we have to include the file Vector/vector_dist.hpp
+ * we also include DrawParticles that has nice utilities to draw particles in parallel accordingly
+ * to simple shapes
*
* \snippet Vector/7_SPH_dlb/main.cpp inclusion
*
@@ -21,56 +25,211 @@
//! \cond [inclusion] \endcond
#include "Vector/vector_dist.hpp"
#include <math.h>
+#include "Draw/DrawParticles.hpp"
//! \cond [inclusion] \endcond
-#include "Draw/DrawParticles.hpp"
+/*!
+ * \page Vector_7_sph_dlb Vector 7 SPH Dam break simulation with Dynamic load balacing
+ *
+ * ## Parameters {#e7_sph_parameters}
+ *
+ * The SPH formulation used in this example code follow these equations
+ *
+ * \f$\frac{dv_a}{dt} = - \sum_{b = NN(a) } m_b \left(\frac{P_a + P_b}{\rho_a \rho_b} + \Pi_{ab} \right) \nabla_{a} W_{ab} + g \tag{1} \f$
+ *
+ * \f$\frac{d\rho_a}{dt} = \sum_{b = NN(a) } m_b v_{ab} \cdot \nabla_{a} W_{ab} \tag{2} \f$
+ *
+ * \f$ P_a = b \left[ \left( \frac{\rho_a}{\rho_{0}} \right)^{\gamma} - 1 \right] \tag{3} \f$
+ *
+ * with
+ *
+ * \f$ \Pi_{ab} = \begin{cases} - \frac {\alpha \bar{c_{ab}} \mu_{ab} }{\bar{\rho_{ab}} } & v_{ab} \cdot r_{ab} > 0 \\ 0 & v_{ab} \cdot r_{ab} < 0 \end{cases} \tag{4}\f$
+ *
+ * and the constants defined as
+ *
+ * \f$ b = \frac{c_{s}^{2} \rho_0}{\gamma} \tag{5} \f$
+ *
+ * \f$ c_s = \sqrt{g \cdot h_{swl}} \tag{6} \f$
+ *
+ * While the particle kernel support is given by
+ *
+ * \f$ H = \sqrt{3 \cdot dp} \tag{7} \f$
+ *
+ * Explain the equations is out of the context of this tutorial. An introduction
+ * can be found in the original Monghagan SPH paper. In this example we use the version
+ * used by Dual-SPH (http://www.dual.sphysics.org/). A summary of the equation and constants can be founded in
+ * their User Manual and the XML user Manual.
+ * In the following we define all the constants required by the simulation
+ *
+ * \snippet Vector/7_SPH_dlb/main.cpp sim parameters
+ *
+ */
+
+/*! \cond [sim parameters] \endcond */
+// A constant to indicate boundary particles
#define BOUNDARY 0
-#define FLUID 1
-double lower_z = 2.0;
+// A constant to indicate fluid particles
+#define FLUID 1
+// initial spacing between particles dp in the formulas
const double dp = 0.0085;
+// Maximum height of the fluid water
+// is coing to be calculated and filled later on
double h_swl = 0.0;
+
+// in the formulas indicated with c_s (constant used to calculate the sound speed)
const double coeff_sound = 20.0;
+
+// gamma in the formulas
const double gamma_ = 7.0;
-// sqrt(3.0*dp*dp)
+
+// sqrt(3.0*dp*dp) support of the kernel
const double H = 0.0147224318643;
+
+// Eta in the formulas
const double Eta2 = 0.01 * H*H;
+
+
const double visco = 0.1;
double cbar = 0.0;
+
+// Mass of the fluid particles
const double MassFluid = 0.000614125;
+
+// Mass of the boundary particles
const double MassBound = 0.000614125;
-const double t_end = 1.0;
+
+// End simulation time
+const double t_end = 1.5;
+
+// Gravity acceleration
const double gravity = 9.81;
+
+// Reference densitu 1000Kg/m^3
const double rho_zero = 1000.0;
+
+// Filled later require h_swl, it is b in the formulas
double B = 0.0;
+
+// Constant used to define time integration
const double CFLnumber = 0.2;
+
+// Minimum T
const double DtMin = 0.00001;
+// Minimum Rho allowed
+const double RhoMin = 700.0;
+
+// Maximum Rho allowed
+const double RhoMax = 1300.0;
+
// Filled in initialization
double max_fluid_height = 0.0;
+// Properties
+
+// FLUID or BOUNDARY
const size_t type = 0;
+
+// Density
const int rho = 1;
+
+// Density at step n-1
const int rho_prev = 2;
+
+// Pressure
const int Pressure = 3;
+
+// Delta rho calculated in the force calculation
const int drho = 4;
+
+// calculated force
const int force = 5;
+
+// velocity
const int velocity = 6;
+
+// velocity at previous step
const int velocity_prev = 7;
-typedef vector_dist<3,double,aggregate<size_t,double,double,double,double,double[3],double[3],double[3]>> particles;
+// Type of the vector containing particles
+typedef vector_dist<3,double,aggregate<size_t,double, double, double, double, double[3], double[3], double[3]>> particles;
+// | | | | | | | |
+// | | | | | | | |
+// type density density Pressure delta force velocity velocity
+// at n-1 density at n - 1
-inline void EqState(particles & vd)
+
+/*! \cond [sim parameters] \endcond */
+
+/*! \brief Linear model
+ *
+ * The linear model count each particle as weight one
+ *
+ */
+struct ModelCustom
+{
+ size_t factor = 1;
+
+ template<typename Decomposition, typename vector> inline void addComputation(Decomposition & dec, const vector & vd, size_t v, size_t p)
+ {
+ if (vd.template getProp<type>(p) == FLUID)
+ {
+ dec.addComputationCost(v,3);
+ }
+ else
+ {
+ dec.addComputationCost(v,2);
+ }
+
+ }
+
+ template<typename Decomposition> inline void applyModel(Decomposition & dec, size_t v)
+ {
+ dec.setSubSubDomainComputationCost(v, dec.getSubSubDomainComputationCost(v) * dec.getSubSubDomainComputationCost(v));
+ }
+};
+
+/*! \brief Linear model
+ *
+ * The linear model count each particle as weight one
+ *
+ */
+struct ModelCustom1
{
-// double min = 100000000.0;
-// double max = 0.0;
-// double accum = 0.0;
-// size_t n = 0;
+ size_t factor = 1;
+
+ template<typename Decomposition, typename vector> inline void addComputation(Decomposition & dec, const vector & vd, size_t v, size_t p)
+ {
+
+ dec.addComputationCost(v,100);
+ }
-// Vcluster & v_cl = create_vcluster();
+ template<typename Decomposition> inline void applyModel(Decomposition & dec, size_t v)
+ {
+ }
+};
+
+/*!
+ * \page Vector_7_sph_dlb Vector 7 SPH Dam break simulation with Dynamic load balacing
+ *
+ * ## Equation of state and SPH Kernels {#e7_sph_equation_state}
+ *
+ * This function implement the formula 3 in the set of equations. It calculate the
+ * pressure of each particle based on the local density of each particle.
+ *
+ * \snippet Vector/7_SPH_dlb/main.cpp eq_state_and_ker
+ *
+ */
+
+/*! \cond [eq_state_and_ker] \endcond */
+
+
+inline void EqState(particles & vd)
+{
auto it = vd.getDomainIterator();
while (it.isNext())
@@ -82,40 +241,27 @@ inline void EqState(particles & vd)
vd.template getProp<Pressure>(a) = B*( rho_frac*rho_frac*rho_frac*rho_frac*rho_frac*rho_frac*rho_frac - 1.0);
- /// DEBUG
-
-/* if (vd.template getProp<Pressure>(a) < min)
- min = vd.template getProp<Pressure>(a);
-
- if (vd.template getProp<Pressure>(a) > max)
- max = vd.template getProp<Pressure>(a);
-
- if (vd.template getProp<Pressure>(a) > 2849.0)
- std::cout << "Particle: " << Point<3,double>(vd.getPos(a)).toString() << std::endl;
-
- accum += vd.template getProp<Pressure>(a);
- n++;*/
-
++it;
}
+}
-/* v_cl.max(max);
- v_cl.min(min);
- v_cl.sum(accum);
- v_cl.sum(n);
+/*! \cond [eq_state_and_ker] \endcond */
- v_cl.execute();
+/*!
+ * \page Vector_7_sph_dlb Vector 7 SPH Dam break simulation with Dynamic load balancing
+ *
+ * This function define the Cubic kernel or \f$ W_{ab} \f$ in the set of equations. The cubic kernel is
+ * defined as
+ *
+ * \f$ \begin{cases} 1.0 - \frac{3}{2} q^2 + \frac{3}{4} q^3 & 0 < q < 1 \\ (2 - q)^3 & 1 < q < 2 \\ 0 & q > 2 \end{cases} \f$
+ *
+ * \snippet Vector/7_SPH_dlb/main.cpp kernel_sph
+ *
+ */
- std::cout << "Max: " << max << " min: " << min << " accum: " << accum/n << " " << B << " n: " << n << std::endl;*/
-}
+/*! \cond [kernel_sph] \endcond */
-const double c1 = -3.0/M_PI/H/H/H/H;
-const double d1 = 9.0/4.0/M_PI/H/H/H/H;
-const double c2 = -3.0/4.0/M_PI/H/H/H/H;
const double a2 = 1.0/M_PI/H/H/H;
-const double a2_4 = 0.25*a2;
-// Filled later
-double W_dap = 0.0;
inline double Wab(double r)
{
@@ -129,6 +275,29 @@ inline double Wab(double r)
return 0.0;
}
+/*! \cond [kernel_sph] \endcond */
+
+/*!
+ * \page Vector_7_sph_dlb Vector 7 SPH Dam break simulation with Dynamic load balancing
+ *
+ * This function define the derivative of the Cubic kernel function \f$ W_{ab} \f$ in the set of equations.
+ *
+ * \f$ \nabla W_{ab} = \beta (x,y,z) \f$
+ *
+ * \f$ \beta = \begin{cases} (c_1 q + d_1 q^2) & 0 < q < 1 \\ c_2 (2 - q)^2 & 1 < q < 2 \end{cases} \f$
+ *
+ * \snippet Vector/7_SPH_dlb/main.cpp kernel_sph_der
+ *
+ */
+
+/*! \cond [kernel_sph_der] \endcond */
+
+const double c1 = -3.0/M_PI/H/H/H/H;
+const double d1 = 9.0/4.0/M_PI/H/H/H/H;
+const double c2 = -3.0/4.0/M_PI/H/H/H/H;
+const double a2_4 = 0.25*a2;
+// Filled later
+double W_dap = 0.0;
inline void DWab(Point<3,double> & dx, Point<3,double> & DW, double r, bool print)
{
@@ -160,8 +329,26 @@ inline void DWab(Point<3,double> & dx, Point<3,double> & DW, double r, bool prin
}
}
+/*! \cond [kernel_sph_der] \endcond */
+
+/*!
+ * \page Vector_7_sph_dlb Vector 7 SPH Dam break simulation with Dynamic load balancing
+ *
+ * This function define the Tensile term. An explanation of the Tensile term is out of the
+ * context of this tutorial, but in brief is an additional repulsive term that avoid the particles
+ * to get enough near. Can be considered at small scale like a repulsive force that avoid
+ * particles to get too close like the Lennard-Jhonned potential at atomistic level. A good
+ * reference is the Monaghan paper "SPH without a Tensile Instability"
+ *
+ * \snippet Vector/7_SPH_dlb/main.cpp tensile_term
+ *
+ *
+ */
+
+/*! \cond [tensile_term] \endcond */
+
// Tensile correction
-inline double Tensile(double r, double rhoa, double rhob, double prs1, double prs2, bool print)
+inline double Tensile(double r, double rhoa, double rhob, double prs1, double prs2)
{
const double qq=r/H;
//-Cubic Spline kernel
@@ -187,13 +374,26 @@ inline double Tensile(double r, double rhoa, double rhob, double prs1, double pr
const double tensilp1=(prs1/(rhoa*rhoa))*(prs1>0? 0.01: -0.2);
const double tensilp2=(prs2/(rhob*rhob))*(prs2>0? 0.01: -0.2);
-// if (print == true)
-// std::cout << "fab " << fab << " tensilp1: " << tensilp1 << " tensilp2: " << tensilp2 << std::endl;
-
return (fab*(tensilp1+tensilp2));
}
-inline double Pi(const Point<3,double> & dr, double rr2, Point<3,double> & dv, double rhoa, double rhob, double massb, double & visc, bool print)
+/*! \cond [tensile_term] \endcond */
+
+
+/*!
+ *
+ * \page Vector_7_sph_dlb Vector 7 SPH Dam break simulation with Dynamic load balancing
+ *
+ * This function is the implementation of the viscous term \f$ \Pi_{ab} \f$
+ *
+ * \snippet Vector/7_SPH_dlb/main.cpp viscous_term
+ *
+ *
+ */
+
+/*! \cond [viscous_term] \endcond */
+
+inline double Pi(const Point<3,double> & dr, double rr2, Point<3,double> & dv, double rhoa, double rhob, double massb, double & visc)
{
const double dot = dr.get(0)*dv.get(0) + dr.get(1)*dv.get(1) + dr.get(2)*dv.get(2);
const double dot_rr2 = dot/(rr2+Eta2);
@@ -205,169 +405,189 @@ inline double Pi(const Point<3,double> & dr, double rr2, Point<3,double> & dv, d
const float robar=(rhoa+rhob)*0.5f;
const float pi_visc=(-visco*cbar*amubar/robar);
- if (print == true)
- std::cout << " visco: " << visco << " " << cbar << " " << amubar << " " << robar << " ";
-
return pi_visc;
}
else
return 0.0;
}
+/*! \cond [viscous_term] \endcond */
+
+/*!
+ *
+ * \page Vector_7_sph_dlb Vector 7 SPH Dam break simulation with Dynamic load balancing
+ *
+ * ## Force calculation {#e7_force_calc}
+ *
+ * Calculate forces. It calculate equation 1 and 2 in the set of formulas
+ *
+ * \snippet Vector/7_SPH_dlb/main.cpp calc_forces
+ *
+ *
+ */
+
+/*! \cond [calc_forces] \endcond */
+
template<typename CellList> inline double calc_forces(particles & vd, CellList & NN, double & max_visc)
{
auto part = vd.getDomainIterator();
-
- Point<3,double> ForceMax({ 0.0, 0.0, 0.0});
- Point<3,double> ForceMin({10000000.0,1000000.0,1000000,0});
-
double visc = 0;
-// std::cout << "c1: " << c1 << " c2: " << c2 << " d1: " << d1 << std::endl;
+ // Update the cell-list
vd.updateCellList(NN);
+ // For each particle ...
while (part.isNext())
{
+ // ... a
auto a = part.get();
// Get the position xp of the particle
Point<3,double> xa = vd.getPos(a);
- if (vd.getProp<type>(a) != FLUID)
- {
- ++part;
- continue;
- }
-
+ // Take the mass of the particle dependently if it is FLUID or BOUNDARY
double massa = (vd.getProp<type>(a) == FLUID)?MassFluid:MassBound;
+
+ // Get the density of the of the particle a
double rhoa = vd.getProp<rho>(a);
+
+ // Get the pressure of the particle a
double Pa = vd.getProp<Pressure>(a);
+
+ // Get the Velocity of the particle a
Point<3,double> va = vd.getProp<velocity>(a);
- // Reset the force counter
+ // Reset the force counter (- gravity on zeta direction)
vd.template getProp<force>(a)[0] = 0.0;
vd.template getProp<force>(a)[1] = 0.0;
vd.template getProp<force>(a)[2] = -gravity;
vd.template getProp<drho>(a) = 0.0;
- size_t cnt = 0;
+ // We threat FLUID particle differently from BOUNDARY PARTICLES ...
+ if (vd.getProp<type>(a) != FLUID)
+ {
+ // If it is a boundary particle calculate the delta rho based on equation 2
+ // This require to run across the neighborhoods particles of a
+ auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(a)));
-// std::cout << "---------------------" << std::endl;
-// std::cout << vd.getPos(a)[0] << " " << vd.getPos(a)[1] << " " << vd.getPos(a)[2] << std::endl;
+ // For each neighborhood particle
+ while (Np.isNext() == true)
+ {
+ // ... q
+ auto b = Np.get();
- // Get an iterator over the neighborhood particles of p
- auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(a)));
+ // Get the position xp of the particle
+ Point<3,double> xb = vd.getPos(b);
- while (Np.isNext() == true)
- {
- // ... q
- auto b = Np.get();
+ // if (p == q) skip this particle
+ if (a.getKey() == b) {++Np; continue;};
- // Get the position xp of the particle
- Point<3,double> xb = vd.getPos(b);
+ // get the mass of the particle
+ double massb = (vd.getProp<type>(b) == FLUID)?MassFluid:MassBound;
- // if (p == q) skip this particle
- if (a.getKey() == b) {++Np; continue;};
+ // Get the velocity of the particle b
+ Point<3,double> vb = vd.getProp<velocity>(b);
- double massb = (vd.getProp<type>(b) == FLUID)?MassFluid:MassBound;
- Point<3,double> vb = vd.getProp<velocity>(b);
- double Pb = vd.getProp<Pressure>(b);
- double rhob = vd.getProp<rho>(b);
+ // Get the pressure and density of particle b
+ double Pb = vd.getProp<Pressure>(b);
+ double rhob = vd.getProp<rho>(b);
- // Get the distance between p and q
- Point<3,double> dr = xa - xb;
- // take the norm of this vector
- double r2 = norm2(dr);
+ // Get the distance between p and q
+ Point<3,double> dr = xa - xb;
+ // take the norm of this vector
+ double r2 = norm2(dr);
- if (r2 < 4.0*H*H)
- {
- double r = sqrt(r2);
+ // If the particles interact ...
+ if (r2 < 4.0*H*H)
+ {
+ // ... calculate delta rho
+ double r = sqrt(r2);
- Point<3,double> v_rel = va - vb;
+ Point<3,double> dv = va - vb;
- Point<3,double> DW;
- DWab(dr,DW,r,false);
+ Point<3,double> DW;
+ DWab(dr,DW,r,false);
- double factor = - massb*((vd.getProp<Pressure>(a) + vd.getProp<Pressure>(b)) / (rhoa * rhob) + Tensile(r,rhoa,rhob,Pa,Pb,false) + Pi(dr,r2,v_rel,rhoa,rhob,massb,visc,false));
-/*
- if (create_vcluster().getProcessUnitID() == 0)
- {
- std::cout << "PARTICLE: " << dr.toString() << std::endl;
- std::cout << "Pressure: " << Pa << " " << Pb << std::endl;
- std::cout << "Density: " << rhoa << " " << rhob << std::endl;
- std::cout << "FACTOR: " << factor << std::endl;
- std::cout << "Tensile: " << Tensile(r,rhoa,rhob,Pa,Pb,false) << std::endl;
- std::cout << "DW: " << DW.get(0) << " " << DW.get(1) << " " << DW.get(2) << std::endl;
- }*/
-
- vd.getProp<force>(a)[0] += factor * DW.get(0);
- vd.getProp<force>(a)[1] += factor * DW.get(1);
- vd.getProp<force>(a)[2] += factor * DW.get(2);
-
- if (xa.get(0) > 0.0085 && xa.get(0) < 0.0105 &&
- xa.get(1) > 0.0085 && xa.get(1) < 0.0105 &&
- xa.get(2) > 0.0085 && xa.get(2) < 0.0105)
- {
- std::cout << "POSITION: " << xb.toString() << std::endl;
- std::cout << "FORCE: " << factor*DW.get(0) << " " << factor*DW.get(1) << " " << factor*DW.get(2) << std::endl;
- std::cout << "FORCE TERM: " << Tensile(r,rhoa,rhob,Pa,Pb,true) << " " << Pi(dr,r2,v_rel,rhoa,rhob,massb,visc,true) << std::endl;
- cnt++;
- }
+ const double dot = dr.get(0)*dv.get(0) + dr.get(1)*dv.get(1) + dr.get(2)*dv.get(2);
+ const double dot_rr2 = dot/(r2+Eta2);
+ max_visc=std::max(dot_rr2,max_visc);
- vd.getProp<drho>(a) += massb*(v_rel.get(0)*DW.get(0)+v_rel.get(1)*DW.get(1)+v_rel.get(2)*DW.get(2));
+ vd.getProp<drho>(a) += massb*(dv.get(0)*DW.get(0)+dv.get(1)*DW.get(1)+dv.get(2)*DW.get(2));
+ }
-/* if (vd.getProp<type>(a) == FLUID)
- {
- std::cout << "DELTA RHO: " << massb*(v_rel.get(0)*DW.get(0)+v_rel.get(1)*DW.get(1)+v_rel.get(2)*DW.get(2)) << " " << v_rel.get(0) << " " << v_rel.get(1) << " " << v_rel.get(2) << " VISC: " << Pi(dr,r2,v_rel,rhoa,rhob,massb,visc,true) << std::endl;
- }*/
+ ++Np;
}
-
- ++Np;
}
-
- if (xa.get(0) > 0.0085 && xa.get(0) < 0.0105 &&
- xa.get(1) > 0.0085 && xa.get(1) < 0.0105 &&
- xa.get(2) > 0.0085 && xa.get(2) < 0.0105)
+ else
{
- std::cout << "FORCE FINAL: " << vd.getProp<force>(a)[0] << " " << vd.getProp<force>(a)[1] << " " << vd.getProp<force>(a)[2] << " " << cnt << std::endl;
- }
+ // If it is a fluid particle calculate based on equation 1 and 2
-/* if (vd.getProp<type>(a) == FLUID)
- {
- std::cout << "DELTA DENSITY: " << vd.getProp<drho>(a) << std::endl;
- std::cout << "+++++++++++++++++++++++++++++++++++" << std::endl;
- }*/
+ // Get an iterator over the neighborhood particles of p
+ auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(a)));
- ++part;
+ // For each neighborhood particle
+ while (Np.isNext() == true)
+ {
+ // ... q
+ auto b = Np.get();
- if (Point<3,double>(vd.getProp<force>(a)).norm() > ForceMax.norm() )
- {
- ForceMax = Point<3,double>(vd.getProp<force>(a));
-// std::cout << "ForceMax: " << ForceMax.toString() << " " << a.getKey() << std::endl;
+ // Get the position xp of the particle
+ Point<3,double> xb = vd.getPos(b);
- Point<3,double> p({0.01,0.0,0.0});
- Point<3,double> DW;
+ // if (p == q) skip this particle
+ if (a.getKey() == b) {++Np; continue;};
- DWab(p,DW,0.01,false);
+ double massb = (vd.getProp<type>(b) == FLUID)?MassFluid:MassBound;
+ Point<3,double> vb = vd.getProp<velocity>(b);
+ double Pb = vd.getProp<Pressure>(b);
+ double rhob = vd.getProp<rho>(b);
-// std::cout << DW.get(0) << " " << DW.get(1) << " " << DW.get(2) << std::endl;
- }
+ // Get the distance between p and q
+ Point<3,double> dr = xa - xb;
+ // take the norm of this vector
+ double r2 = norm2(dr);
- if (Point<3,double>(vd.getProp<force>(a)).norm() < ForceMin.norm() )
- {
- ForceMin = Point<3,double>(vd.getProp<force>(a));
+ // if they interact
+ if (r2 < 4.0*H*H)
+ {
+ double r = sqrt(r2);
+
+ Point<3,double> v_rel = va - vb;
+
+ Point<3,double> DW;
+ DWab(dr,DW,r,false);
+
+ double factor = - massb*((vd.getProp<Pressure>(a) + vd.getProp<Pressure>(b)) / (rhoa * rhob) + Tensile(r,rhoa,rhob,Pa,Pb) + Pi(dr,r2,v_rel,rhoa,rhob,massb,visc));
+
+ vd.getProp<force>(a)[0] += factor * DW.get(0);
+ vd.getProp<force>(a)[1] += factor * DW.get(1);
+ vd.getProp<force>(a)[2] += factor * DW.get(2);
+
+ vd.getProp<drho>(a) += massb*(v_rel.get(0)*DW.get(0)+v_rel.get(1)*DW.get(1)+v_rel.get(2)*DW.get(2));
+ }
+
+ ++Np;
+ }
}
+
+ ++part;
}
+}
- // Get the maximum viscosity term across processors
- Vcluster & v_cl = create_vcluster();
- v_cl.max(visc);
- v_cl.execute();
- max_visc = visc;
+/*! \cond [calc_forces] \endcond */
-// std::cout << "---------------------------------" << std::endl;
-}
+/*!
+ *
+ * \page Vector_7_sph_dlb Vector 7 SPH Dam break simulation with Dynamic load balancing
+ *
+ * This function calculate the Maximum acceleration and velocity across the particles.
+ *
+ * \snippet Vector/7_SPH_dlb/main.cpp max_acc_vel
+ *
+ *
+ */
+/*! \cond [max_acc_vel] \endcond */
void max_acceleration_and_velocity(particles & vd, double & max_acc, double & max_vel)
{
@@ -396,7 +616,30 @@ void max_acceleration_and_velocity(particles & vd, double & max_acc, double & ma
max_vel = sqrt(max_vel);
}
-double dt_old = 0.0;
+/*! \cond [max_acc_vel] \endcond */
+
+/*!
+ *
+ * \page Vector_7_sph_dlb Vector 7 SPH Dam break simulation with Dynamic load balancing
+ *
+ * In this example we are using Dynamic time-stepping. The Dynamic time stepping is
+ * calculated with the Courant-Friedrich-Lewy condition. See Monaghan 1992 "Smoothed Particle Hydrodynamic"
+ *
+ * \f$ \delta t = CFL \cdot min(t_f,t_{cv}) \f$
+ *
+ * where
+ *
+ * \f$ \delta t_f = min \sqrt{h/f_a}\f$
+ *
+ * \f$ \delta t_{cv} = min \frac{h}{c_s + max \left| \frac{hv_{ab} \cdot r_{ab}}{r_{ab}^2} \right|} \f$
+ *
+ *
+ * \snippet Vector/7_SPH_dlb/main.cpp dyn_stepping
+ *
+ *
+ */
+
+/*! \cond [dyn_stepping] \endcond */
double calc_deltaT(particles & vd, double ViscDtMax)
{
@@ -410,29 +653,51 @@ double calc_deltaT(particles & vd, double ViscDtMax)
//-dt2 combines the Courant and the viscous time-step controls.
const double dt_cv = H/(std::max(cbar,Maxvel*10.) + H*ViscDtMax);
-// std::cout << "dt_f " << dt_f << " dt_cv: " << dt_cv << " Maxvel: " << Maxvel << std::endl;
-
//-dt new value of time step.
double dt=double(CFLnumber)*std::min(dt_f,dt_cv);
if(dt<double(DtMin))
dt=double(DtMin);
- if (dt_old != dt)
- {
- std::cout << "Dt changed to " << dt << " MaxVel: " << Maxvel << " Maxacc: " << Maxacc << " lower_z: " << lower_z << std::endl;
-
- dt_old = dt;
- }
-
-// std::cout << "Returned dt: " << dt << std::endl;
return dt;
}
+/*! \cond [dyn_stepping] \endcond */
+
+/*!
+ *
+ * \page Vector_7_sph_dlb Vector 7 SPH Dam break simulation with Dynamic load balancing
+ *
+ * This function perform verlet integration accordingly to the Verlet time stepping scheme
+ *
+ * \f$ v_a^{n+1} = v_a^{n-1} + 2 \delta t F_a^{n} \f$
+ *
+ * \f$ r_a^{n+1} = \delta t V_a^n + 0.5 \delta t^2 F_a^n \f$
+ *
+ * \f$ \rho_a^{n+1} = \rho_a^{n-1} + 2 \delta t D_a^n \f$
+ *
+ * Every N Verlet steps the euler stepping scheme is choosen to avoid instabilities
+ *
+ * \f$ v_a^{n+1} = v_a^{n} + \delta t F_a^n \f$
+ *
+ * \f$ r_a^{n+1} = r_a^{n} + \delta t V_a^n + 0.5 delta t^2 F_a^n \f$
+ *
+ * \f$ \rho_a^n + \delta t D_a^n \f$
+ *
+ * More the integration this function also check that no particles go outside the simulation
+ * domain or their density go dangerously out of range
+ *
+ * \snippet Vector/7_SPH_dlb/main.cpp verlet_int
+ *
+ *
+ */
+
+/*! \cond [verlet_int] \endcond */
+
openfpm::vector<size_t> to_remove;
-bool particle_out = true;
+size_t cnt = 0;
-void verlet_int(particles & vd, double dt)
+void verlet_int(particles & vd, double dt, bool VerletStep)
{
to_remove.clear();
@@ -448,7 +713,26 @@ void verlet_int(particles & vd, double dt)
if (vd.template getProp<type>(a) == BOUNDARY)
{
+ double rhop = vd.template getProp<rho>(a);
+
// Update only the density
+ if (VerletStep == true)
+ {
+ vd.template getProp<velocity>(a)[0] = 0.0;
+ vd.template getProp<velocity>(a)[1] = 0.0;
+ vd.template getProp<velocity>(a)[2] = 0.0;
+ vd.template getProp<rho>(a) = vd.template getProp<rho_prev>(a) + dt2*vd.template getProp<drho>(a);
+ }
+ else
+ {
+ vd.template getProp<velocity>(a)[0] = 0.0;
+ vd.template getProp<velocity>(a)[1] = 0.0;
+ vd.template getProp<velocity>(a)[2] = 0.0;
+ vd.template getProp<rho>(a) = vd.template getProp<rho>(a) + dt*vd.template getProp<drho>(a);
+ }
+
+ vd.template getProp<rho_prev>(a) = rhop;
+
++part;
continue;
}
@@ -458,75 +742,62 @@ void verlet_int(particles & vd, double dt)
double dy = vd.template getProp<velocity>(a)[1]*dt + vd.template getProp<force>(a)[1]*dt205;
double dz = vd.template getProp<velocity>(a)[2]*dt + vd.template getProp<force>(a)[2]*dt205;
-// bool outrhop=(rhopnew<RhopOutMin||rhopnew>RhopOutMax);
-
-// if (vd.getPos(a)[0] > 0.0085 && vd.getPos(a)[0] < 0.0105 &&
-// vd.getPos(a)[1] > 0.0085 && vd.getPos(a)[1] < 0.0105 &&
-// vd.getPos(a)[2] > 0.0085 && vd.getPos(a)[2] < 0.0105)
-// {
-// std::cout << "DeltaX: " << dx << " " << dy << " " << dz << std::endl;
-/* std::cout << "FORCE: " << vd.template getProp<force>(a)[0] << " "
- << vd.template getProp<force>(a)[1] << " "
- << vd.template getProp<force>(a)[2] <<
- " DENSITY: " << vd.template getProp<rho_prev>(a) + dt2*vd.template getProp<drho>(a) <<
- " PRESSURE: " << vd.template getProp<Pressure>(a) <<
- " POSITION Z: " << vd.getPos(a)[2]
- << std::endl;*/
- // std::cout << "FORCE TERM: " << Tensile(r,rhoa,rhob,Pa,Pb,true) << " " << massb*Pi(dr,r2,v_rel,rhoa,rhob,massa,massb) << std::endl;
-// }
-
vd.getPos(a)[0] += dx;
vd.getPos(a)[1] += dy;
vd.getPos(a)[2] += dz;
- if (vd.getPos(a)[2] < lower_z)
- lower_z = vd.getPos(a)[2];
-
- if (vd.getPos(a)[0] < 0.000263878 || vd.getPos(a)[1] < 0.000263878 || vd.getPos(a)[2] < 0.000263878 ||
- vd.getPos(a)[0] > 0.000263878+1.59947 || vd.getPos(a)[1] > 0.000263878+0.672972 || vd.getPos(a)[2] > 0.000263878+0.903944)
- {
-// std::cout << "Particle out" << std::endl;
- to_remove.add(a.getKey());
-
- particle_out = true;
- }
-
double velX = vd.template getProp<velocity>(a)[0];
double velY = vd.template getProp<velocity>(a)[1];
double velZ = vd.template getProp<velocity>(a)[2];
double rhop = vd.template getProp<rho>(a);
- vd.template getProp<velocity>(a)[0] = vd.template getProp<velocity_prev>(a)[0] + vd.template getProp<force>(a)[0]*dt2;
- vd.template getProp<velocity>(a)[1] = vd.template getProp<velocity_prev>(a)[1] + vd.template getProp<force>(a)[1]*dt2;
- vd.template getProp<velocity>(a)[2] = vd.template getProp<velocity_prev>(a)[2] + vd.template getProp<force>(a)[2]*dt2;
- vd.template getProp<rho>(a) = vd.template getProp<rho_prev>(a) + dt2*vd.template getProp<drho>(a);
+ if (VerletStep == true)
+ {
+ vd.template getProp<velocity>(a)[0] = vd.template getProp<velocity_prev>(a)[0] + vd.template getProp<force>(a)[0]*dt2;
+ vd.template getProp<velocity>(a)[1] = vd.template getProp<velocity_prev>(a)[1] + vd.template getProp<force>(a)[1]*dt2;
+ vd.template getProp<velocity>(a)[2] = vd.template getProp<velocity_prev>(a)[2] + vd.template getProp<force>(a)[2]*dt2;
+ vd.template getProp<rho>(a) = vd.template getProp<rho_prev>(a) + dt2*vd.template getProp<drho>(a);
+ }
+ else
+ {
+ vd.template getProp<velocity>(a)[0] = vd.template getProp<velocity>(a)[0] + vd.template getProp<force>(a)[0]*dt;
+ vd.template getProp<velocity>(a)[1] = vd.template getProp<velocity>(a)[1] + vd.template getProp<force>(a)[1]*dt;
+ vd.template getProp<velocity>(a)[2] = vd.template getProp<velocity>(a)[2] + vd.template getProp<force>(a)[2]*dt;
+ vd.template getProp<rho>(a) = vd.template getProp<rho>(a) + dt*vd.template getProp<drho>(a);
+ }
+
+ // Check if there are particles to remove
+
+ if (vd.getPos(a)[0] < 0.000263878 || vd.getPos(a)[1] < 0.000263878 || vd.getPos(a)[2] < 0.000263878 ||
+ vd.getPos(a)[0] > 0.000263878+1.59947 || vd.getPos(a)[1] > 0.000263878+0.672972 || vd.getPos(a)[2] > 0.000263878+0.903944 ||
+ vd.template getProp<rho>(a) < RhoMin || vd.template getProp<rho>(a) > RhoMax)
+ {
+ std::cout << "Particle_out" << std::endl;
+ to_remove.add(a.getKey());
+ }
vd.template getProp<velocity_prev>(a)[0] = velX;
vd.template getProp<velocity_prev>(a)[1] = velY;
vd.template getProp<velocity_prev>(a)[2] = velZ;
vd.template getProp<rho_prev>(a) = rhop;
-// std::cout << "VELOCITY: " << vd.template getProp<velocity>(a)[0] << " " << vd.template getProp<velocity>(a)[1] << " " << vd.template getProp<velocity>(a)[2] << std::endl;
-// std::cout << "++++++++++++++++++++++++++++++++++++++++++" << std::endl;
-
-/* if (particle_out == true)
- {
- std::cout << "PARTICLE DENSITY: " << vd.template getProp<rho>(a) << " PARTICLE PRESSURE: " << vd.template getProp<Pressure>(a) << " Delta rho: " << vd.template getProp<drho>(a) << " VELOCITY: " << vd.template getProp<velocity>(a)[0] << " " << vd.template getProp<velocity>(a)[1] << " " << vd.template getProp<velocity>(a)[2] << std::endl;
- }*/
-
++part;
}
vd.remove(to_remove,0);
+
+ cnt++;
}
+/*! \cond [verlet_int] \endcond */
+
int main(int argc, char* argv[])
{
-
/*!
- * \page Vector_7_sph_dlb Vector 7 SPH Dam break simulation with Dynamic load balacing
*
- * ## Initialization ##
+ * \page Vector_7_sph_dlb Vector 7 SPH Dam break simulation with Dynamic load balancing
+ *
+ * ## Main function ##
*
* Here we Initialize the library, we create a Box that define our domain, boundary conditions, ghost
*
@@ -542,8 +813,8 @@ int main(int argc, char* argv[])
openfpm_init(&argc,&argv);
// Here we define our domain a 2D box with internals from 0 to 1.0 for x and y
- Box<3,double> domain({-0.05,-0.05,-0.05},{2.0070,1.0040,1.0040});
- size_t sz[3] = {243,125,125};
+ Box<3,double> domain({-0.05,-0.05,-0.05},{1.7010,0.7065,0.5025});
+ size_t sz[3] = {207,90,66};
// Fill W_dap
W_dap = 1.0/Wab(H/1.5);
@@ -557,7 +828,7 @@ int main(int argc, char* argv[])
//! \cond [Initialization and parameters] \endcond
/*!
- * \page Vector_7_SPH_dlb Vector 7 SPH Dam break simulation with Dynamic load balacing
+ * \page Vector_7_SPH_dlb Vector 7 SPH Dam break simulation with Dynamic load balancing
*
* ## %Vector create ##
*
@@ -573,15 +844,14 @@ int main(int argc, char* argv[])
//! \cond [vector inst] \endcond
- particles vd(0,domain,bc,g);
+ particles vd(0,domain,bc,g,DEC_GRAN(4096));
//! \cond [vector inst] \endcond
// the scalar is the element at position 0 in the aggregate
const int type = 0;
-// Box<3,double> fluid_box({dp/2.0,dp/2.0,dp/2.0},{0.4+dp/2.0,0.67-dp/2.0,0.3+dp/2.0});
- Box<3,double> fluid_box({0.2,0.2,0.298},{0.2+2*dp,0.2+2*dp,0.298+2*dp});
+ Box<3,double> fluid_box({dp/2.0,dp/2.0,dp/2.0},{0.4+dp/2.0,0.67-dp/2.0,0.3+dp/2.0});
// first we create Fluid particles
// Fluid particles are created
@@ -592,8 +862,6 @@ int main(int argc, char* argv[])
B = (coeff_sound)*(coeff_sound)*gravity*h_swl*rho_zero / gamma_;
cbar = coeff_sound * sqrt(gravity * h_swl);
-// std::cout << "MAX FLUID: " << max_fluid_height << std::endl;
-
while (fluid_it.isNext())
{
vd.add();
@@ -613,8 +881,6 @@ int main(int argc, char* argv[])
vd.template getLastProp<Pressure>() = rho_zero * gravity * (max_fluid_height - fluid_it.get().get(2));
- std::cout << "B: " << B << std::endl;
-
vd.template getLastProp<rho>() = pow(vd.template getLastProp<Pressure>() / B + 1, 1.0/gamma_) * rho_zero;
vd.template getLastProp<rho_prev>() = vd.template getLastProp<rho>();
vd.template getLastProp<velocity>()[0] = 0.0;
@@ -690,7 +956,31 @@ int main(int argc, char* argv[])
}
vd.map();
- vd.ghost_get<rho,Pressure>();
+ vd.getDecomposition().write("Decomposition_before_load_bal");
+
+ // Now that we fill the vector with particles
+ ModelCustom md;
+
+ vd.addComputationCosts(md);
+ vd.getDecomposition().getDistribution().write("BEFORE_DECOMPOSE");
+ vd.getDecomposition().decompose();
+ vd.map();
+
+ vd.addComputationCosts(md);
+ vd.getDecomposition().getDistribution().write("AFTER_DECOMPOSE1");
+
+ vd.getDecomposition().rebalance(1);
+
+ vd.map();
+ vd.getDecomposition().getDistribution().write("AFTER_DECOMPOSE2");
+
+ std::cout << "N particles: " << vd.size_local() << " " << create_vcluster().getProcessUnitID() << " " << "Get processor Load " << vd.getDecomposition().getDistribution().getProcessorLoad() << std::endl;
+
+ vd.write("Geometry");
+ vd.getDecomposition().write("Decomposition_after_load_bal");
+ vd.getDecomposition().getDistribution().write("Distribution_load_bal");
+
+ vd.ghost_get<type,rho,Pressure,velocity>();
auto NN = vd.getCellList(2*H);
@@ -699,27 +989,44 @@ int main(int argc, char* argv[])
size_t write = 0;
size_t it = 0;
+ size_t it_reb = 0;
double t = 0.0;
while (t <= t_end)
{
- const size_t type = 0;
- const int rho = 1;
- const int Pressure = 2;
- const int drho = 3;
- const int force = 4;
- const int velocity = 5;
- const int velocity_prev = 6;
+ timer it_time;
+
+ ////// Do rebalancing every 200 timesteps
+ it_reb++;
+ if (it_reb == 10)
+ {
+ vd.map();
+
+ it_reb = 0;
+ ModelCustom md;
+ vd.addComputationCosts(md);
+ vd.getDecomposition().rebalance(1);
+
+ std::cout << "REBALANCED " << std::endl;
+ }
vd.map();
- vd.ghost_get<type,rho,Pressure,velocity,velocity_prev>();
+ vd.ghost_get<type,rho,Pressure,velocity>();
// Calculate pressure from the density
EqState(vd);
double max_visc = 0.0;
+ it_time.start();
+
// Calc forces
calc_forces(vd,NN,max_visc);
+ it_time.stop();
+
+ // Get the maximum viscosity term across processors
+ Vcluster & v_cl = create_vcluster();
+ v_cl.max(max_visc);
+ v_cl.execute();
// Calculate delta t integration
double dt = calc_deltaT(vd,max_visc);
@@ -727,7 +1034,14 @@ int main(int argc, char* argv[])
// std::cout << "Calculate deltaT: " << dt << " " << DtMin << std::endl;
// VerletStep
- verlet_int(vd,dt);
+ it++;
+ if (it < 40)
+ verlet_int(vd,dt,true);
+ else
+ {
+ verlet_int(vd,dt,false);
+ it = 0;
+ }
t += dt;
@@ -736,6 +1050,12 @@ int main(int argc, char* argv[])
vd.write("Geometry",write);
write++;
+
+ std::cout << "TIME: " << t << " write " << it_time.getwct() << " " << v_cl.getProcessUnitID() << " " << cnt << std::endl;
+ }
+ else
+ {
+ std::cout << "TIME: " << t << " " << it_time.getwct() << " " << v_cl.getProcessUnitID() << " " << cnt << std::endl;
}
}
@@ -746,7 +1066,7 @@ int main(int argc, char* argv[])
//! \cond [finalize] \endcond
/*!
- * \page Vector_7_SPH_dlb Vector 7 SPH Dam break simulation with Dynamic load balacing
+ * \page Vector_7_SPH_dlb Vector 7 SPH Dam break simulation with Dynamic load balancing
*
* ## Full code ## {#code_e0_sim}
*
diff --git a/openfpm_numerics b/openfpm_numerics
index 2a805e31..9e49bd95 160000
--- a/openfpm_numerics
+++ b/openfpm_numerics
@@ -1 +1 @@
-Subproject commit 2a805e31cd0bc1e5fe89c8f6a03ba22b1691ad2d
+Subproject commit 9e49bd9594b1b1c0d9cdd19e609aaba7a163d907
diff --git a/openfpm_pdata.doc b/openfpm_pdata.doc
index 21a5e06b..335cc0cb 100644
--- a/openfpm_pdata.doc
+++ b/openfpm_pdata.doc
@@ -38,7 +38,7 @@ PROJECT_NAME = "OpenFPM_pdata"
# could be handy for archiving the generated documentation or if some version
# control system is used.
-PROJECT_NUMBER = 0.7.0
+PROJECT_NUMBER = 0.8.0
# Using the PROJECT_BRIEF tag one can provide an optional one line description
# for a project that appears at the top of each page and should give viewer a
diff --git a/src/Decomposition/CartDecomposition.hpp b/src/Decomposition/CartDecomposition.hpp
index d9e5e91c..be5dbdb9 100755
--- a/src/Decomposition/CartDecomposition.hpp
+++ b/src/Decomposition/CartDecomposition.hpp
@@ -351,11 +351,11 @@ public:
for (size_t i = 0; i < dist.getNSubSubDomains(); i++)
{
- dist.setMigrationCost(i, norm * migration * dist.getSubSubDomainComputationCost(i));
+ dist.setMigrationCost(i, norm * migration /* * dist.getSubSubDomainComputationCost(i) */);
for (size_t s = 0; s < dist.getNSubSubDomainNeighbors(i); s++)
{
- dist.setCommunicationCost(i, s, 1 * dist.getSubSubDomainComputationCost(i) * ts);
+ dist.setCommunicationCost(i, s, 1 * /* dist.getSubSubDomainComputationCost(i) * */ ts);
}
prev += dist.getNSubSubDomainNeighbors(i);
}
@@ -994,6 +994,7 @@ public:
sub_domains.clear();
box_nn_processor.clear();
fine_s.clear();
+ loc_box.clear();
nn_prcs<dim, T>::reset();
ie_ghost<dim, T>::reset();
ie_loc_ghost<dim, T>::reset();
@@ -1122,7 +1123,7 @@ public:
*/
inline size_t getSubSubDomainComputationCost(size_t id)
{
- return dist.getComputationCost(id);
+ return dist.getSubSubDomainComputationCost(id);
}
/*! \brief Operator to access the size of the sub-graph
diff --git a/src/Decomposition/Distribution/ParMetisDistribution.hpp b/src/Decomposition/Distribution/ParMetisDistribution.hpp
index 3bddfc98..6790e49f 100644
--- a/src/Decomposition/Distribution/ParMetisDistribution.hpp
+++ b/src/Decomposition/Distribution/ParMetisDistribution.hpp
@@ -118,28 +118,33 @@ class ParMetisDistribution
for (size_t i = 0; i <= Np; i++)
vtxdist.get(i) = n_vtxdist.get(i);
+ openfpm::vector<size_t> cnt;
+ cnt.resize(Np);
+
// Renumber the main graph and re-create the map
- for (size_t p = 0; p < (size_t)Np; p++)
+/* for (size_t p = 0; p < (size_t)Np; p++)
{
size_t i = 0;
for (rid j = vtxdist.get(p); j < vtxdist.get(p + 1); ++j, i++)
{
setMapId(j, v_per_proc.get(p).get(i));
- gp.vertex(v_per_proc.get(p).get(i).id).template get<nm_v::id>() = j.id;
+ gp.vertex(v_per_proc.get(p).get(i).id).template get<nm_v::id>() = cnt.get(p) + vtxdist.get(p).id;
+ cnt.get(p)++;
}
- }
- }
-
- void createMapsFromGlobalGraph(openfpm::vector<size_t> & vtxdist)
- {
-/* openfpm::vector<size_t> cnt_np;
+ }*/
- for (size_t i = 0 ; i < gp.getNVertex() ; i++)
+ for (size_t i = 0 ; i < gp.getNVertex(); ++i)
{
- cnt_np(gp.template vertex<nm_v::proc_id>)++;
+ size_t pid = gp.template vertex_p<nm_v::proc_id>(i);
- gp.setMapId()
- }*/
+ rid j = rid(vtxdist.get(pid).id + cnt.get(pid));
+ gid gi = gid(i);
+
+ gp.template vertex_p<nm_v::id>(i) = j.id;
+ cnt.get(pid)++;
+
+ setMapId(j,gi);
+ }
}
/*! \brief operator to access the vertex by mapped position
diff --git a/src/Graph/CartesianGraphFactory.hpp b/src/Graph/CartesianGraphFactory.hpp
index a3d8d275..d9638e03 100755
--- a/src/Graph/CartesianGraphFactory.hpp
+++ b/src/Graph/CartesianGraphFactory.hpp
@@ -76,6 +76,9 @@ struct fill_id<dim, G_v, NO_VERTEX_ID>
template<unsigned int dim, int lin_id, typename dT, typename G_v, typename v, int impl>
class fill_prop
{
+ //! Domain
+ const Box<dim,dT> & domain;
+
//! Reference to an array containing the spacing
const dT (&szd)[dim];
@@ -91,8 +94,8 @@ class fill_prop
public:
//! Fill the object from where to take the properties
- fill_prop(G_v & g_v, const dT (&szd)[dim], grid_key_dx<dim> & gk, const grid_sm<dim, void> & gs) :
- szd(szd), gk(gk), g_v(g_v), gs(gs)
+ fill_prop(G_v & g_v, const dT (&szd)[dim], grid_key_dx<dim> & gk, const grid_sm<dim, void> & gs, const Box<dim,dT> & domain)
+ :domain(domain), szd(szd), gk(gk), g_v(g_v), gs(gs)
{
}
@@ -102,7 +105,7 @@ public:
{
typedef typename boost::fusion::result_of::at<v, boost::mpl::int_<T::value>>::type t_val;
- g_v.template get<t_val::value>() = gk.get(T::value) * szd[T::value];
+ g_v.template get<t_val::value>() = gk.get(T::value) * szd[T::value] + domain.getLow(T::value);
fill_id<dim, G_v, lin_id>::fill(g_v, gk, gs);
}
};
@@ -140,7 +143,7 @@ class fill_prop<dim, lin_id, dT, G_v, v, 0>
public:
//! Fill the object from where to take the properties
- fill_prop(G_v & g_v, const dT (&szd)[dim], grid_key_dx<dim> & gk, const grid_sm<dim, void> & gs)
+ fill_prop(G_v & g_v, const dT (&szd)[dim], grid_key_dx<dim> & gk, const grid_sm<dim, void> & gs, const Box<dim,dT> & domain)
{
}
@@ -180,6 +183,8 @@ public:
template<unsigned int dim, int lin_id, typename dT, typename G_v, typename v>
class fill_prop<dim, lin_id, dT, G_v, v, 2>
{
+ //! Domain
+ const Box<dim,dT> & domain;
//! Reference to an array containing the spacing
const dT (&szd)[dim];
@@ -196,8 +201,8 @@ class fill_prop<dim, lin_id, dT, G_v, v, 2>
public:
//! Fill the object from where to take the properties
- fill_prop(G_v & g_v, const dT (&szd)[dim], grid_key_dx<dim> & gk, const grid_sm<dim, void> & gs) :
- szd(szd), gk(gk), g_v(g_v), gs(gs)
+ fill_prop(G_v & g_v, const dT (&szd)[dim], grid_key_dx<dim> & gk, const grid_sm<dim, void> & gs, const Box<dim,dT> & domain)
+ :domain(domain), szd(szd), gk(gk), g_v(g_v), gs(gs)
{
}
@@ -212,7 +217,7 @@ public:
g_v.template get<t_val::value>()[i] = 0.0;
for (size_t i = 0 ; i < dim ; i++)
- g_v.template get<t_val::value>()[i] = gk.get(i) * static_cast<float>(szd[i]);
+ g_v.template get<t_val::value>()[i] = gk.get(i) * static_cast<float>(szd[i]) + domain.getLow(i);
fill_id<dim, G_v, lin_id>::fill(g_v, gk, gs);
}
@@ -330,7 +335,7 @@ public:
// vertex spatial properties functor
- fill_prop<dim, lin_id, T, decltype(gp.vertex(g.LinId(key))), typename to_boost_vmpl<pos...>::type, fill_prop_by_type<sizeof...(pos), p, Graph, pos...>::value> flp(obj, szd, key, g);
+ fill_prop<dim, lin_id, T, decltype(gp.vertex(g.LinId(key))), typename to_boost_vmpl<pos...>::type, fill_prop_by_type<sizeof...(pos), p, Graph, pos...>::value> flp(obj, szd, key, g, dom);
// fill properties
@@ -451,7 +456,7 @@ public:
// vertex spatial properties functor
- fill_prop<dim, lin_id, T, decltype(gp.vertex(g.LinId(key))), typename to_boost_vmpl<pos...>::type, fill_prop_by_type<sizeof...(pos), p, Graph, pos...>::value> flp(obj, szd, key, g);
+ fill_prop<dim, lin_id, T, decltype(gp.vertex(g.LinId(key))), typename to_boost_vmpl<pos...>::type, fill_prop_by_type<sizeof...(pos), p, Graph, pos...>::value> flp(obj, szd, key, g, dom);
// fill properties
diff --git a/src/Graph/ids.hpp b/src/Graph/ids.hpp
index 63c3dc79..ea4826a1 100644
--- a/src/Graph/ids.hpp
+++ b/src/Graph/ids.hpp
@@ -17,8 +17,17 @@
*/
struct rid
{
+ //! id
idx_t id;
+ //! Cosntructor
+ rid(size_t id)
+ :id(id)
+ {}
+
+ //! Constructor
+ rid() {}
+
inline bool operator<=(const rid & r) const
{
return id <= r.id;
@@ -79,6 +88,14 @@ struct rid
struct gid
{
size_t id;
+
+ //! Constructor
+ gid(){};
+
+ //! Constructor
+ gid(size_t id)
+ :id(id)
+ {}
};
/*! Here we define different the remapped-id
diff --git a/src/Makefile.am b/src/Makefile.am
index 23d87159..ca89186e 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -14,7 +14,7 @@ nobase_include_HEADERS = Decomposition/CartDecomposition.hpp Decomposition/CartD
example.mk \
Decomposition/Distribution/metis_util.hpp Decomposition/Distribution/SpaceDistribution.hpp Decomposition/Distribution/parmetis_dist_util.hpp Decomposition/Distribution/parmetis_util.hpp Decomposition/Distribution/MetisDistribution.hpp Decomposition/Distribution/ParMetisDistribution.hpp Decomposition/Distribution/DistParMetisDistribution.hpp dec_optimizer.hpp SubdomainGraphNodes.hpp \
Graph/ids.hpp Graph/dist_map_graph.hpp Graph/DistGraphFactory.hpp \
- DLB/DLB.hpp
+ DLB/DLB.hpp DLB/LB_Model.hpp
lib_LIBRARIES = libofpm_pdata.a
libofpm_pdata_a_SOURCES = lib/pdata.cpp
diff --git a/src/Vector/vector_dist.hpp b/src/Vector/vector_dist.hpp
index 031a0d43..0034c7a5 100644
--- a/src/Vector/vector_dist.hpp
+++ b/src/Vector/vector_dist.hpp
@@ -30,10 +30,13 @@
#include "data_type/aggregate.hpp"
#include "NN/VerletList/VerletList.hpp"
#include "vector_dist_comm.hpp"
+#include "DLB/LB_Model.hpp"
#define NO_ID false
#define ID true
+#define DEC_GRAN(gr) ((size_t)gr << 32)
+
// Perform a ghost get or a ghost put
#define GET 1
#define PUT 2
@@ -236,6 +239,9 @@ public:
check_new(this,8,VECTOR_DIST_EVENT,4);
#endif
+ if (opt >> 32 != 0)
+ this->setDecompositionGranularity(opt >> 32);
+
check_parameters(box);
init_structures(np);
@@ -1022,21 +1028,19 @@ public:
g_m--;
}
- /*! \brief Add the computation cost on the decomposition comming from the particles
+ /*! \brief Add the computation cost on the decomposition coming from the particles
*
*/
- inline void addComputationCosts()
+ template <typename Model=ModelLin>inline void addComputationCosts(Model md=Model())
{
- CellDecomposer_sm<dim, St> cdsm;
+ CellDecomposer_sm<dim, St, shift<dim,St>> cdsm;
Decomposition & dec = getDecomposition();
cdsm.setDimensions(dec.getDomain(), dec.getGrid().getSize(), 0);
for (size_t i = 0; i < getDecomposition().getNSubSubDomains(); i++)
- {
dec.setSubSubDomainComputationCost(i, 1);
- }
auto it = getDomainIterator();
@@ -1044,11 +1048,15 @@ public:
{
size_t v = cdsm.getCell(this->getPos(it.get()));
- dec.addComputationCost(v, 1);
+ md.addComputation(dec,*this,v,it.get().getKey());
++it;
}
+ // Go throught all the sub-sub-domains and apply the model
+
+ for (size_t i = 0 ; i < dec.getDistribution().getNSubSubDomains(); i++)
+ md.applyModel(dec,i);
}
/*! \brief Output particle position and properties
diff --git a/src/Vector/vector_dist_MP_unit_tests.hpp b/src/Vector/vector_dist_MP_unit_tests.hpp
index 44797b43..4f797a77 100644
--- a/src/Vector/vector_dist_MP_unit_tests.hpp
+++ b/src/Vector/vector_dist_MP_unit_tests.hpp
@@ -85,7 +85,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_multiphase_cell_list_test )
auto CL_phase1 = phases.get(1).getCellList(r_cut);
// This function create a Verlet-list between phases 0 and 1
- auto NN_ver01 = createVerlet(phases.get(0),CL_phase1,r_cut);
+ auto NN_ver01 = createVerlet(phases.get(0),phases.get(1),CL_phase1,r_cut);
// Check NNver0_1
@@ -128,7 +128,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_multiphase_cell_list_test )
auto CL_all = createCellListM<2>(phases,r_cut);
// This create a Verlet-list between phase 0 and all the other phases
- auto NNver0_all = createVerletM<2>(phases.get(0),CL_all,r_cut);
+ auto NNver0_all = createVerletM<2>(phases.get(0),phases,CL_all,r_cut);
it = phases.get(0).getDomainIterator();
@@ -235,7 +235,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_multiphase_cell_list_sym_test )
auto CL_phase1 = phases.get(1).getCellListSym(r_cut);
// This function create a Verlet-list between phases 0 and 1
- auto NN_ver01 = createVerletSym(phases.get(0),CL_phase1,r_cut);
+ auto NN_ver01 = createVerletSym(phases.get(0),phases.get(1),CL_phase1,r_cut);
// Check NNver0_1
@@ -305,15 +305,15 @@ BOOST_AUTO_TEST_CASE( vector_dist_multiphase_cell_list_sym_test )
// This function create an "Empty" Multiphase Cell List
auto CL_all = createCellListSymM<2>(phases,r_cut);
- typedef decltype(createVerletSymM<2>(phases.get(0),CL_all,r_cut)) verlet_type;
+ typedef decltype(createVerletSymM<2>(phases.get(0),phases,CL_all,r_cut)) verlet_type;
verlet_type NNver_all[4];
// This create a Verlet-list between phase all phases to all the other phases
- NNver_all[0] = createVerletSymM<2>(phases.get(0),CL_all,r_cut);
- NNver_all[1] = createVerletSymM<2>(phases.get(1),CL_all,r_cut);
- NNver_all[2] = createVerletSymM<2>(phases.get(2),CL_all,r_cut);
- NNver_all[3] = createVerletSymM<2>(phases.get(3),CL_all,r_cut);
+ NNver_all[0] = createVerletSymM<2>(phases.get(0),phases,CL_all,r_cut);
+ NNver_all[1] = createVerletSymM<2>(phases.get(1),phases,CL_all,r_cut);
+ NNver_all[2] = createVerletSymM<2>(phases.get(2),phases,CL_all,r_cut);
+ NNver_all[3] = createVerletSymM<2>(phases.get(3),phases,CL_all,r_cut);
// all phases to all phases
diff --git a/src/Vector/vector_dist_comm.hpp b/src/Vector/vector_dist_comm.hpp
index 355035dd..d79031be 100644
--- a/src/Vector/vector_dist_comm.hpp
+++ b/src/Vector/vector_dist_comm.hpp
@@ -8,8 +8,6 @@
#ifndef SRC_VECTOR_VECTOR_DIST_COMM_HPP_
#define SRC_VECTOR_VECTOR_DIST_COMM_HPP_
-#define V_SUB_UNIT_FACTOR 64
-
#define SKIP_LABELLING 512
#define KEEP_PROPERTIES 512
@@ -47,6 +45,9 @@ inline static size_t compute_options(size_t opt)
template<unsigned int dim, typename St, typename prop, typename Decomposition = CartDecomposition<dim,St>, typename Memory = HeapMemory>
class vector_dist_comm
{
+ //! Number of units for each sub-domain
+ size_t v_sub_unit_factor = 64;
+
//! definition of the send vector for position
typedef openfpm::vector<Point<dim, St>, Memory> send_pos_vector;
@@ -784,14 +785,24 @@ public:
{
}
- /*! \brief Get the number of minimum sub-domain
+ /*! \brief Get the number of minimum sub-domain per processor
*
* \return minimum number
*
*/
- static size_t getDefaultNsubsub()
+ size_t getDecompositionGranularity()
+ {
+ return v_sub_unit_factor;
+ }
+
+ /*! \brief Set the minimum number of sub-domain per processor
+ *
+ * \param n_sub
+ *
+ */
+ void setDecompositionGranularity(size_t n_sub)
{
- return V_SUB_UNIT_FACTOR;
+ this->v_sub_unit_factor = n_sub;
}
/*! \brief Initialize the decomposition
@@ -826,7 +837,7 @@ public:
// Get the number of processor and calculate the number of sub-domain
// for decomposition
size_t n_proc = v_cl.getProcessingUnits();
- size_t n_sub = n_proc * getDefaultNsubsub();
+ size_t n_sub = n_proc * getDecompositionGranularity();
// Calculate the maximum number (before merging) of sub-domain on
// each dimension
diff --git a/src/Vector/vector_dist_multiphase_functions.hpp b/src/Vector/vector_dist_multiphase_functions.hpp
index 4a3c6e32..543eb68e 100644
--- a/src/Vector/vector_dist_multiphase_functions.hpp
+++ b/src/Vector/vector_dist_multiphase_functions.hpp
@@ -11,20 +11,25 @@
#include "NN/CellList/CellListM.hpp"
#include "NN/VerletList/VerletListM.hpp"
-template<typename Vector,typename CL, typename T> VerletList<Vector::dims,typename Vector::stype,FAST,shift<Vector::dims,typename Vector::stype>> createVerlet(Vector & v, CL & cl, T r_cut)
+template<typename Vector, typename CL, typename T> VerletList<Vector::dims,typename Vector::stype,FAST,shift<Vector::dims,typename Vector::stype>> createVerlet(Vector & v, Vector & v1, CL & cl, T r_cut)
{
VerletList<Vector::dims,typename Vector::stype,FAST,shift<Vector::dims,typename Vector::stype>> ver;
- ver.Initialize(cl,r_cut,v.getPosVector(),v.size_local());
+ ver.Initialize(cl,r_cut,v.getPosVector(),v1.getPosVector(),v.size_local());
return ver;
}
-template<unsigned int sh_byte, typename Vector,typename CL, typename T> VerletListM<Vector::dims,typename Vector::stype,sh_byte,shift<Vector::dims,typename Vector::stype>> createVerletM(Vector & v, CL & cl, T r_cut)
+template<unsigned int sh_byte, typename Vector , typename Vector1,typename CL, typename T> VerletListM<Vector::dims,typename Vector::stype,sh_byte,CL,shift<Vector::dims,typename Vector::stype>> createVerletM(Vector & v, Vector1 & phases, CL & cl, T r_cut)
{
- VerletListM<Vector::dims,typename Vector::stype,sh_byte,shift<Vector::dims,typename Vector::stype>> ver;
+ VerletListM<Vector::dims,typename Vector::stype,sh_byte,CL,shift<Vector::dims,typename Vector::stype>> ver;
- ver.Initialize(cl,r_cut,v.getPosVector(),v.size_local());
+ openfpm::vector<pos_v<Vector::dims,typename Vector::stype>> v_phases;
+
+ for (size_t i = 0 ; i < phases.size() ; i++)
+ v_phases.add(pos_v<Vector::dims,typename Vector::stype>(phases.get(i).getPosVector()));
+
+ ver.Initialize(cl,r_cut,v.getPosVector(),v_phases,v.size_local());
return ver;
}
@@ -65,20 +70,25 @@ template<unsigned int nbit, typename Vector, typename T> CellListM<Vector::dims,
/////// Symmetric version
-template<typename Vector,typename CL, typename T> VerletList<Vector::dims,typename Vector::stype,FAST,shift<Vector::dims,typename Vector::stype>> createVerletSym(Vector & v, CL & cl, T r_cut)
+template<typename Vector,typename CL, typename T> VerletList<Vector::dims,typename Vector::stype,FAST,shift<Vector::dims,typename Vector::stype>> createVerletSym(Vector & v, Vector & v1, CL & cl, T r_cut)
{
VerletList<Vector::dims,typename Vector::stype,FAST,shift<Vector::dims,typename Vector::stype>> ver;
- ver.Initialize(cl,r_cut,v.getPosVector(),v.size_local());
+ ver.Initialize(cl,r_cut,v.getPosVector(),v1.getPosVector(),v.size_local());
return ver;
}
-template<unsigned int sh_byte, typename Vector,typename CL, typename T> VerletListM<Vector::dims,typename Vector::stype,sh_byte,shift<Vector::dims,typename Vector::stype>> createVerletSymM(Vector & v, CL & cl, T r_cut)
+template<unsigned int sh_byte, typename Vector, typename Vector1 ,typename CL, typename T> VerletListM<Vector::dims,typename Vector::stype,sh_byte,CL,shift<Vector::dims,typename Vector::stype>> createVerletSymM(Vector & v, Vector1 & phases, CL & cl, T r_cut)
{
- VerletListM<Vector::dims,typename Vector::stype,sh_byte,shift<Vector::dims,typename Vector::stype>> ver;
+ VerletListM<Vector::dims,typename Vector::stype,sh_byte,CL,shift<Vector::dims,typename Vector::stype>> ver;
+
+ openfpm::vector<pos_v<Vector::dims,typename Vector::stype>> v_phases;
+
+ for (size_t i = 0 ; i < phases.size() ; i++)
+ v_phases.add(pos_v<Vector::dims,typename Vector::stype>(phases.get(i).getPosVector()));
- ver.Initialize(cl,r_cut,v.getPosVector(),v.size_local(),VL_SYMMETRIC);
+ ver.Initialize(cl,r_cut,v.getPosVector(),v_phases,v.size_local(),VL_SYMMETRIC);
return ver;
}
diff --git a/src/Vector/vector_dist_unit_test.hpp b/src/Vector/vector_dist_unit_test.hpp
index aefd55ab..9102beb4 100644
--- a/src/Vector/vector_dist_unit_test.hpp
+++ b/src/Vector/vector_dist_unit_test.hpp
@@ -105,7 +105,7 @@ void Test2D_ghost(Box<2,float> & box)
typedef Point_test<float> p;
// Get the default minimum number of sub-sub-domain per processor (granularity of the decomposition)
- size_t n_sub = vector_dist<2,float, Point_test<float>, CartDecomposition<2,float> >::getDefaultNsubsub() * v_cl.getProcessingUnits();
+ size_t n_sub = 64 * v_cl.getProcessingUnits();
// Convert the request of having a minimum n_sub number of sub-sub domain into grid decompsition of the space
size_t sz = CartDecomposition<2,float>::getDefaultGrid(n_sub);
--
GitLab