/*!
* \page Vector_7_sph_dlb_opt Vector 7 SPH Dam break simulation with Dynamic load balacing (Optimized version)
*
*
* [TOC]
*
*
* # SPH with Dynamic load Balancing # {#SPH_dlb}
*
* This is just a rework of the SPH Dam break simulation optimized to get better performance we will focus on the
* optimization and differences with the previous example
*
* \see \ref Vector_7_sph_dlb
*
* \htmlonly
* Simulation video 1
*
*
* Simulation video 2
*
*
* Simulation dynamic load balancing video 1
*
*
* Simulation dynamic load balancing video 2
*
*
* \endhtmlonly
*
* \htmlonly
*
#include "Draw/DrawParticles.hpp"
/*!
* \page Vector_7_sph_dlb_opt Vector 7 SPH Dam break simulation with Dynamic load balacing (Optimized version)
*
* ## Using verlet list with skin{#e7_sph_dlb_opt}
*
* The first optimization that we operate is the usage of verlet list. The verlet are reconstructed only when
* the maximum displacement is bigger than the half skin. Because we have to calculate
* the maximum displacement the verlet and euler integration has been modified to do this.
* The function accept a reference to max_disp that is filled with the maximum displacement calculated
* from these functions.
*
* \snippet Vector/7_SPH_dlb_opt/main.cpp verlet_new_arg
*
* \snippet Vector/7_SPH_dlb_opt/main.cpp euler_new_arg
*
*
* The variable is reset inside verlet and euler time integration function
*
* \snippet Vector/7_SPH_dlb_opt/main.cpp reset_max_disp
*
* while iteration across particle the maximum displacement is saved inside the variable
* max_disp
*
* \snippet Vector/7_SPH_dlb_opt/main.cpp calc_max_disp
*
* We also have to be careful that if we are removing particles we have to reconstruct the verlet list,
* so we set it to a really big number
*
* \snippet Vector/7_SPH_dlb_opt/main.cpp big_number_set
*
* Because the maximum displacement has to be calculated across processors, we use the
* function max in Vcluster to calculate the maximum displacement across processors.
*
* \snippet Vector/7_SPH_dlb_opt/main.cpp max_across_proc
*
* We also calculate the skin part and ghost plus skin. Consider also that the ghost
* must be extended to ghost + skin so r_gskin
*
* \snippet Vector/7_SPH_dlb_opt/main.cpp skin_calc
*
* As we explained before, we update the verlet only if particles move more than the skin.
* In case they move more than the skin we do first a map to redistribute the particles and in the
* meanwhile we check if it is a good moment to rebalance. We decided to combine these two steps
* because in case we rebalance we have anyway to reconstruct the Verler-list. Than we calculate
* the pressure for all the particles, refresh the ghost, update the Verlet-list and reset the
* total displacement. In case the the total displacement does not overshoot the skin we just
* calculate the pressure for all the particles and refresh the ghost. We must use the option
* **SKIP_LABELLING**
*
* \snippet Vector/7_SPH_dlb_opt/main.cpp update_verlet
*
* We pass the max_displacement variable to verlet_int and euler_int function. We also add
* the maximum displacement per iteration to the total maximum displacement
*
* \snippet Vector/7_SPH_dlb_opt/main.cpp pass_ver_eu
*
*
*
* ## Symmetric interaction (Crossing scheme){#e7_sph_dlb_sym}
*
* Symmetric interaction give the possibility to reduce the computation by half and speed-up
* your simulation. To do this we have to do some changes into the function calc forces.
* Symmetric interaction require to write on the ghost area. So at the beginning of the function
* we reset the ghost part. In the meanwhile because we have the external force gravity that
* operate per particles, we set this force at the beginning.
*
* \warning The requirement to set per particle external forces outside the particle loop
* come from the symmetric scheme. Suppose to have in pseudocode this
* \code{.unparsed}
1 for each particles p
2 reset the force for p
3 for each neighborhood particle q of p
4 calculate the force p-q
5 add the contribution to p
6 add the contribution to q
* \endcode
* suppose we are on particle p=0 and calculate the force with q=10 we add the
* contribution to p and q. Unfortunately accordingly to this cycle when we reach
* particle q = 10 we reset what we previously calculated. So we have to write
* \code{.unparsed}
*
1 for each particles p
2 reset the force for p
3 for each particles p
4 for each neighborhood particle q of p
5 calculate the force p-q
6 add the contribution to p
7 add the contribution to q
\endcode
*
* With this code we set the per particle external force to gravity and reset the derivative of the
* density for the domain particles
*
* \snippet Vector/7_SPH_dlb_opt/main.cpp reset_particles
*
* With this code we reset the force and derivative of the density of the particles on the ghost part
*
* \snippet Vector/7_SPH_dlb_opt/main.cpp reset_particles2
*
* Small changes must be done to iterate over the neighborhood particles
*
* \snippet Vector/7_SPH_dlb_opt/main.cpp nn_part
*
* skip the self interaction
*
* \snippet Vector/7_SPH_dlb_opt/main.cpp skip_self
*
* This is instead an important change (and honestly it took quite some hour of debuging to
* discover the problem). In case we are on boundary particle (p = boundary particle) and
* calculating an interaction with a particle q = fluid particle we have to remeber that we have
* also to calculate the force for q (not only drho)
*
* \snippet Vector/7_SPH_dlb_opt/main.cpp symmetry
*
* for a fluid particle instead we calculate p-q interaction and we add the contribution
* to p and q. Because we do not integrate over the boundary particles we can also avoid to
* check that q is a boundary particle
*
* \snippet Vector/7_SPH_dlb_opt/main.cpp symmetry2
*
* After the calculation cycle we have to merge the forces and delta density calculated on
* the ghost with the real particles.
*
* \snippet Vector/7_SPH_dlb_opt/main.cpp ghost_put
*
* It is important when we construct our vector of particles to pass the option **BIND_DEC_TO_GHOST**.
* To use symmetric calculation in parallel environment the decomposition must be consistent with the
* cell decomposition of the space.
*
* \snippet Vector/7_SPH_dlb_opt/main.cpp important_option
*
* To construct a Verlet-list using the CRS scheme we use the following function
*
* \snippet Vector/7_SPH_dlb_opt/main.cpp get_verlet_crs
*
* while to update the verlet list we use the following
*
* \snippet Vector/7_SPH_dlb_opt/main.cpp update_verlet_crs
*
* ## Using re-decompose instead of decompose {#e7_sph_dlb_opt_red}
*
* Using redecompose instead of decompose produce less jumping decomposition during the simulation
*
* \htmlonly
* Video 1
*
*
*
*
* Video 2
*
*
*
*
* \endhtmlonly
*
*/
// A constant to indicate boundary particles
#define BOUNDARY 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 going to be calculated and filled later on
double h_swl = 0.0;
// c_s in the formulas (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) support of the kernel
const double H = 0.0147224318643;
// Eta in the formulas
const double Eta2 = 0.01 * H*H;
// alpha in the formula
const double visco = 0.1;
// cbar in the formula (calculated later)
double cbar = 0.0;
// Mass of the fluid particles
const double MassFluid = 0.000614125;
// Mass of the boundary particles
const double MassBound = 0.000614125;
// End simulation time
#ifndef TEST_RUN
const double t_end = 1.5;
#else
const double t_end = 0.001;
#endif
// 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;
/*! \cond [sim parameters] \endcond */
/*! \cond [vector_dist_def] \endcond */
// Type of the vector containing particles
typedef vector_dist<3,double,aggregate > particles;
// | | | | | | | |
// | | | | | | | |
// type density density Pressure delta force velocity velocity
// at n-1 density at n - 1
struct ModelCustom
{
template inline void addComputation(Decomposition & dec, vector & vd, size_t v, size_t p)
{
if (vd.template getProp(p) == FLUID)
dec.addComputationCost(v,4);
else
dec.addComputationCost(v,3);
}
template inline void applyModel(Decomposition & dec, size_t v)
{
dec.setSubSubDomainComputationCost(v, dec.getSubSubDomainComputationCost(v) * dec.getSubSubDomainComputationCost(v));
}
double distributionTol()
{
return 1.01;
}
};
inline void EqState(particles & vd)
{
auto it = vd.getDomainIterator();
while (it.isNext())
{
auto a = it.get();
double rho_a = vd.template getProp(a);
double rho_frac = rho_a / rho_zero;
vd.template getProp(a) = B*( rho_frac*rho_frac*rho_frac*rho_frac*rho_frac*rho_frac*rho_frac - 1.0);
++it;
}
}
const double a2 = 1.0/M_PI/H/H/H;
inline double Wab(double r)
{
r /= H;
if (r < 1.0)
return (1.0 - 3.0/2.0*r*r + 3.0/4.0*r*r*r)*a2;
else if (r < 2.0)
return (1.0/4.0*(2.0 - r*r)*(2.0 - r*r)*(2.0 - r*r))*a2;
else
return 0.0;
}
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)
{
const double qq=r/H;
if (qq < 1.0)
{
double qq2 = qq * qq;
double fac = (c1*qq + d1*qq2)/r;
DW.get(0) = fac*dx.get(0);
DW.get(1) = fac*dx.get(1);
DW.get(2) = fac*dx.get(2);
}
else if (qq < 2.0)
{
double wqq = (2.0 - qq);
double fac = c2 * wqq * wqq / r;
DW.get(0) = fac * dx.get(0);
DW.get(1) = fac * dx.get(1);
DW.get(2) = fac * dx.get(2);
}
else
{
DW.get(0) = 0.0;
DW.get(1) = 0.0;
DW.get(2) = 0.0;
}
}
inline double Tensile(double r, double rhoa, double rhob, double prs1, double prs2)
{
const double qq=r/H;
//-Cubic Spline kernel
double wab;
if(r>H)
{
double wqq1=2.0f-qq;
double wqq2=wqq1*wqq1;
wab=a2_4*(wqq2*wqq1);
}
else
{
double wqq2=qq*qq;
double wqq3=wqq2*qq;
wab=a2*(1.0f-1.5f*wqq2+0.75f*wqq3);
}
//-Tensile correction.
double fab=wab*W_dap;
fab*=fab; fab*=fab; //fab=fab^4
const double tensilp1=(prs1/(rhoa*rhoa))*(prs1>0? 0.01: -0.2);
const double tensilp2=(prs2/(rhob*rhob))*(prs2>0? 0.01: -0.2);
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)
{
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);
visc=std::max(dot_rr2,visc);
if(dot < 0)
{
const float amubar=H*dot_rr2;
const float robar=(rhoa+rhob)*0.5f;
const float pi_visc=(-visco*cbar*amubar/robar);
return pi_visc;
}
else
return 0.0;
}
template inline double calc_forces(particles & vd, VerletList & NN, double & max_visc)
{
/*! \cond [reset_particles] \endcond */
// Reset the ghost
auto itg = vd.getDomainIterator();
while (itg.isNext())
{
auto p = itg.get();
// Reset force
// Reset the force counter (- gravity on zeta direction)
vd.template getProp(p)[0] = 0.0;
vd.template getProp(p)[1] = 0.0;
vd.template getProp(p)[2] = -gravity;
vd.template getProp(p) = 0.0;
++itg;
}
/*! \cond [reset_particles] \endcond */
/*! \cond [reset_particles2] \endcond */
auto itg2 = vd.getGhostIterator();
while (itg2.isNext())
{
auto p = itg2.get();
// Reset force
// Reset the force counter (- gravity on zeta direction)
vd.template getProp(p)[0] = 0.0;
vd.template getProp(p)[1] = 0.0;
vd.template getProp(p)[2] = 0.0;
vd.template getProp(p) = 0.0;
++itg2;
}
/*! \cond [reset_particles2] \endcond */
// Get an iterator over particles
auto part = vd.getParticleIteratorCRS(NN.getInternalCellList());
double visc = 0;
// 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);
// Take the mass of the particle dependently if it is FLUID or BOUNDARY
double massa = (vd.getProp(a) == FLUID)?MassFluid:MassBound;
// Get the density of the of the particle a
double rhoa = vd.getProp(a);
// Get the pressure of the particle a
double Pa = vd.getProp(a);
// Get the Velocity of the particle a
Point<3,double> va = vd.getProp(a);
// We threat FLUID particle differently from BOUNDARY PARTICLES ...
if (vd.getProp(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
//! \cond [nn_part] \endcond
auto Np = NN.template getNNIterator(a);
//! \cond [nn_part] \endcond
// For each neighborhood particle
while (Np.isNext() == true)
{
// ... q
auto b = Np.get();
// Get the position xp of the particle
Point<3,double> xb = vd.getPos(b);
//! \cond [skip_self] \endcond
// if (p == q) skip this particle
if (a == b) {++Np; continue;};
//! \cond [skip_self] \endcond
// get the mass of the particle
double massb = (vd.getProp(b) == FLUID)?MassFluid:MassBound;
// Get the velocity of the particle b
Point<3,double> vb = vd.getProp(b);
// Get the pressure and density of particle b
double Pb = vd.getProp(b);
double rhob = vd.getProp(b);
// Get the distance between p and q
Point<3,double> dr = xa - xb;
// take the norm of this vector
double r2 = norm2(dr);
// If the particles interact ...
if (r2 < 4.0*H*H)
{
// ... calculate delta rho
double r = sqrt(r2);
Point<3,double> dv = va - vb;
Point<3,double> DW;
DWab(dr,DW,r,false);
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);
double scal = (dv.get(0)*DW.get(0)+dv.get(1)*DW.get(1)+dv.get(2)*DW.get(2));
vd.getProp(a) += massb*scal;
vd.getProp(b) += massa*scal;
//! \cond [symmetry] \endcond
// If it is a fluid we have to calculate force
if (vd.getProp(b) == FLUID)
{
Point<3,double> v_rel = va - vb;
double factor = - ((vd.getProp(a) + vd.getProp(b)) / (rhoa * rhob) + Tensile(r,rhoa,rhob,Pa,Pb) + Pi(dr,r2,v_rel,rhoa,rhob,massb,visc));
vd.getProp(b)[0] -= massa * factor * DW.get(0);
vd.getProp(b)[1] -= massa * factor * DW.get(1);
vd.getProp(b)[2] -= massa * factor * DW.get(2);
}
//! \cond [symmetry] \endcond
}
++Np;
}
}
else
{
// If it is a fluid particle calculate based on equation 1 and 2
// Get an iterator over the neighborhood particles of p
auto Np = NN.template getNNIterator(a);
// For each neighborhood particle
while (Np.isNext() == true)
{
// ... q
auto b = Np.get();
// Get the position xp of the particle
Point<3,double> xb = vd.getPos(b);
// if (p == q) skip this particle
if (a == b) {++Np; continue;};
double massb = (vd.getProp(b) == FLUID)?MassFluid:MassBound;
Point<3,double> vb = vd.getProp(b);
double Pb = vd.getProp(b);
double rhob = vd.getProp(b);
// Get the distance between p and q
Point<3,double> dr = xa - xb;
// take the norm of this vector
double r2 = norm2(dr);
// 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);
//! \cond [symmetry2] \endcond
double factor = - ((vd.getProp(a) + vd.getProp(b)) / (rhoa * rhob) + Tensile(r,rhoa,rhob,Pa,Pb) + Pi(dr,r2,v_rel,rhoa,rhob,massb,visc));
vd.getProp(a)[0] += massb * factor * DW.get(0);
vd.getProp(a)[1] += massb * factor * DW.get(1);
vd.getProp(a)[2] += massb * factor * DW.get(2);
vd.getProp(b)[0] -= massa * factor * DW.get(0);
vd.getProp(b)[1] -= massa * factor * DW.get(1);
vd.getProp(b)[2] -= massa * factor * DW.get(2);
double scal = (v_rel.get(0)*DW.get(0)+v_rel.get(1)*DW.get(1)+v_rel.get(2)*DW.get(2));
vd.getProp(a) += massb*scal;
vd.getProp(b) += massa*scal;
//! \cond [symmetry2] \endcond
}
++Np;
}
}
++part;
}
//! \cond [ghost_put] \endcond
vd.template ghost_put();
//! \cond [ghost_put] \endcond
}
void max_acceleration_and_velocity(particles & vd, double & max_acc, double & max_vel)
{
// Calculate the maximum acceleration
auto part = vd.getDomainIterator();
while (part.isNext())
{
auto a = part.get();
Point<3,double> acc(vd.getProp(a));
double acc2 = norm2(acc);
Point<3,double> vel(vd.getProp(a));
double vel2 = norm2(vel);
if (vel2 >= max_vel)
max_vel = vel2;
if (acc2 >= max_acc)
max_acc = acc2;
++part;
}
max_acc = sqrt(max_acc);
max_vel = sqrt(max_vel);
Vcluster & v_cl = create_vcluster();
v_cl.max(max_acc);
v_cl.max(max_vel);
v_cl.execute();
}
double calc_deltaT(particles & vd, double ViscDtMax)
{
double Maxacc = 0.0;
double Maxvel = 0.0;
max_acceleration_and_velocity(vd,Maxacc,Maxvel);
//-dt1 depends on force per unit mass.
const double dt_f = (Maxacc)?sqrt(H/Maxacc):std::numeric_limits::max();
//-dt2 combines the Courant and the viscous time-step controls.
const double dt_cv = H/(std::max(cbar,Maxvel*10.) + H*ViscDtMax);
//-dt new value of time step.
double dt=double(CFLnumber)*std::min(dt_f,dt_cv);
if(dt to_remove;
size_t cnt = 0;
/*! \cond [verlet_new_arg] \endcond */
void verlet_int(particles & vd, double dt, double & max_disp)
/*! \cond [verlet_new_arg] \endcond */
{
// list of the particle to remove
to_remove.clear();
// particle iterator
auto part = vd.getDomainIterator();
double dt205 = dt*dt*0.5;
double dt2 = dt*2.0;
/*! \cond [reset_max_disp] \endcond */
max_disp = 0;
/*! \cond [reset_max_disp] \endcond */
// For each particle ...
while (part.isNext())
{
// ... a
auto a = part.get();
// if the particle is boundary
if (vd.template getProp(a) == BOUNDARY)
{
// Update rho
double rhop = vd.template getProp(a);
// Update only the density
vd.template getProp(a)[0] = 0.0;
vd.template getProp(a)[1] = 0.0;
vd.template getProp(a)[2] = 0.0;
vd.template getProp(a) = vd.template getProp(a) + dt2*vd.template getProp(a);
vd.template getProp(a) = rhop;
++part;
continue;
}
//-Calculate displacement and update position / Calcula desplazamiento y actualiza posicion.
double dx = vd.template getProp(a)[0]*dt + vd.template getProp(a)[0]*dt205;
double dy = vd.template getProp(a)[1]*dt + vd.template getProp(a)[1]*dt205;
double dz = vd.template getProp(a)[2]*dt + vd.template getProp(a)[2]*dt205;
vd.getPos(a)[0] += dx;
vd.getPos(a)[1] += dy;
vd.getPos(a)[2] += dz;
/*! \cond [calc_max_disp] \endcond */
double d2 = dx*dx + dy*dy + dz*dz;
max_disp = (max_disp > d2)?max_disp:d2;
/*! \cond [calc_max_disp] \endcond */
double velX = vd.template getProp(a)[0];
double velY = vd.template getProp(a)[1];
double velZ = vd.template getProp(a)[2];
double rhop = vd.template getProp(a);
vd.template getProp(a)[0] = vd.template getProp(a)[0] + vd.template getProp(a)[0]*dt2;
vd.template getProp(a)[1] = vd.template getProp(a)[1] + vd.template getProp(a)[1]*dt2;
vd.template getProp(a)[2] = vd.template getProp(a)[2] + vd.template getProp(a)[2]*dt2;
vd.template getProp(a) = vd.template getProp(a) + dt2*vd.template getProp(a);
// Check if the particle go out of range in space and in density
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(a) < RhoMin || vd.template getProp(a) > RhoMax)
{
to_remove.add(a.getKey());
/*! \cond [big_number_set] \endcond */
// if we remove something the verlet are not anymore correct and must be reconstructed
max_disp = 100.0;
/*! \cond [big_number_set] \endcond */
}
vd.template getProp(a)[0] = velX;
vd.template getProp(a)[1] = velY;
vd.template getProp(a)[2] = velZ;
vd.template getProp(a) = rhop;
++part;
}
/*! \cond [max_across_proc] \endcond */
Vcluster & v_cl = create_vcluster();
v_cl.max(max_disp);
v_cl.execute();
max_disp = sqrt(max_disp);
/*! \cond [max_across_proc] \endcond */
// remove the particles
vd.remove(to_remove,0);
// increment the iteration counter
cnt++;
}
/*! \cond [euler_new_arg] \endcond */
void euler_int(particles & vd, double dt, double & max_disp)
/*! \cond [euler_new_arg] \endcond */
{
// list of the particle to remove
to_remove.clear();
// particle iterator
auto part = vd.getDomainIterator();
double dt205 = dt*dt*0.5;
double dt2 = dt*2.0;
max_disp = 0;
// For each particle ...
while (part.isNext())
{
// ... a
auto a = part.get();
// if the particle is boundary
if (vd.template getProp(a) == BOUNDARY)
{
// Update rho
double rhop = vd.template getProp(a);
// Update only the density
vd.template getProp(a)[0] = 0.0;
vd.template getProp(a)[1] = 0.0;
vd.template getProp(a)[2] = 0.0;
vd.template getProp(a) = vd.template getProp(a) + dt*vd.template getProp(a);
vd.template getProp(a) = rhop;
++part;
continue;
}
//-Calculate displacement and update position / Calcula desplazamiento y actualiza posicion.
double dx = vd.template getProp(a)[0]*dt + vd.template getProp(a)[0]*dt205;
double dy = vd.template getProp(a)[1]*dt + vd.template getProp(a)[1]*dt205;
double dz = vd.template getProp(a)[2]*dt + vd.template getProp(a)[2]*dt205;
vd.getPos(a)[0] += dx;
vd.getPos(a)[1] += dy;
vd.getPos(a)[2] += dz;
double d2 = dx*dx + dy*dy + dz*dz;
max_disp = (max_disp > d2)?max_disp:d2;
double velX = vd.template getProp(a)[0];
double velY = vd.template getProp(a)[1];
double velZ = vd.template getProp(a)[2];
double rhop = vd.template getProp(a);
vd.template getProp(a)[0] = vd.template getProp(a)[0] + vd.template getProp(a)[0]*dt;
vd.template getProp(a)[1] = vd.template getProp(a)[1] + vd.template getProp(a)[1]*dt;
vd.template getProp(a)[2] = vd.template getProp(a)[2] + vd.template getProp(a)[2]*dt;
vd.template getProp(a) = vd.template getProp(a) + dt*vd.template getProp(a);
// Check if the particle go out of range in space and in density
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(a) < RhoMin || vd.template getProp(a) > RhoMax)
{
to_remove.add(a.getKey());
// if we remove something the verlet are not anymore correct and must be reconstructed
max_disp = 100.0;
}
vd.template getProp(a)[0] = velX;
vd.template getProp(a)[1] = velY;
vd.template getProp(a)[2] = velZ;
vd.template getProp(a) = rhop;
++part;
}
// remove the particles
vd.remove(to_remove,0);
Vcluster & v_cl = create_vcluster();
v_cl.max(max_disp);
v_cl.execute();
max_disp = sqrt(max_disp);
// increment the iteration counter
cnt++;
}
int main(int argc, char* argv[])
{
// initialize the library
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},{1.7010,0.7065,0.5025});
size_t sz[3] = {207,90,66};
// Fill W_dap
W_dap = 1.0/Wab(H/1.5);
// Here we define the boundary conditions of our problem
size_t bc[3]={NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
/*! \cond [skin_calc] \endcond */
double skin = 0.25 * 2*H;
double r_gskin = 2*H + skin;
// extended boundary around the domain, and the processor domain
// by the support of the cubic kernel
Ghost<3,double> g(r_gskin);
/*! \cond [skin_calc] \endcond */
// Eliminating the lower part of the ghost
// We are using CRS scheme
g.setLow(0,0.0);
g.setLow(1,0.0);
g.setLow(2,0.0);
/*! \cond [important_option] \endcond */
particles vd(0,domain,bc,g,BIND_DEC_TO_GHOST);
/*! \cond [important_option] \endcond */
#if 0
// You can ignore all these dp/2.0 is a trick to reach the same initialization
// of Dual-SPH that use a different criteria to draw particles
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});
// return an iterator to the fluid particles to add to vd
auto fluid_it = DrawParticles::DrawBox(vd,sz,domain,fluid_box);
// here we fill some of the constants needed by the simulation
max_fluid_height = fluid_it.getBoxMargins().getHigh(2);
h_swl = fluid_it.getBoxMargins().getHigh(2) - fluid_it.getBoxMargins().getLow(2);
B = (coeff_sound)*(coeff_sound)*gravity*h_swl*rho_zero / gamma_;
cbar = coeff_sound * sqrt(gravity * h_swl);
// for each particle inside the fluid box ...
while (fluid_it.isNext())
{
// ... add a particle ...
vd.add();
// ... and set it position ...
vd.getLastPos()[0] = fluid_it.get().get(0);
vd.getLastPos()[1] = fluid_it.get().get(1);
vd.getLastPos()[2] = fluid_it.get().get(2);
// and its type.
vd.template getLastProp() = FLUID;
// We also initialize the density of the particle and the hydro-static pressure given by
//
// rho_zero*g*h = P
//
// rho_p = (P/B + 1)^(1/Gamma) * rho_zero
//
vd.template getLastProp() = rho_zero * gravity * (max_fluid_height - fluid_it.get().get(2));
vd.template getLastProp() = pow(vd.template getLastProp() / B + 1, 1.0/gamma_) * rho_zero;
vd.template getLastProp() = vd.template getLastProp();
vd.template getLastProp()[0] = 0.0;
vd.template getLastProp()[1] = 0.0;
vd.template getLastProp()[2] = 0.0;
vd.template getLastProp()[0] = 0.0;
vd.template getLastProp()[1] = 0.0;
vd.template getLastProp()[2] = 0.0;
// next fluid particle
++fluid_it;
}
#endif
// Recipient
Box<3,double> recipient1({0.0,0.0,0.0},{1.6+dp/2.0,0.67+dp/2.0,0.4+dp/2.0});
Box<3,double> recipient2({dp,dp,dp},{1.6-dp/2.0,0.67-dp/2.0,0.4+dp/2.0});
Box<3,double> obstacle1({0.9,0.24-dp/2.0,0.0},{1.02+dp/2.0,0.36,0.45+dp/2.0});
Box<3,double> obstacle2({0.9+dp,0.24+dp/2.0,0.0},{1.02-dp/2.0,0.36-dp,0.45-dp/2.0});
Box<3,double> obstacle3({0.9+dp,0.24,0.0},{1.02,0.36,0.45});
openfpm::vector> holes;
holes.add(recipient2);
holes.add(obstacle1);
auto bound_box = DrawParticles::DrawSkin(vd,sz,domain,holes,recipient1);
while (bound_box.isNext())
{
vd.add();
vd.getLastPos()[0] = bound_box.get().get(0);
vd.getLastPos()[1] = bound_box.get().get(1);
vd.getLastPos()[2] = bound_box.get().get(2);
vd.template getLastProp() = BOUNDARY;
vd.template getLastProp() = rho_zero;
vd.template getLastProp() = rho_zero;
vd.template getLastProp()[0] = 0.0;
vd.template getLastProp()[1] = 0.0;
vd.template getLastProp()[2] = 0.0;
vd.template getLastProp()[0] = 0.0;
vd.template getLastProp()[1] = 0.0;
vd.template getLastProp()[2] = 0.0;
++bound_box;
}
auto obstacle_box = DrawParticles::DrawSkin(vd,sz,domain,obstacle2,obstacle1);
while (obstacle_box.isNext())
{
vd.add();
vd.getLastPos()[0] = obstacle_box.get().get(0);
vd.getLastPos()[1] = obstacle_box.get().get(1);
vd.getLastPos()[2] = obstacle_box.get().get(2);
vd.template getLastProp() = BOUNDARY;
vd.template getLastProp() = rho_zero;
vd.template getLastProp() = rho_zero;
vd.template getLastProp()[0] = 0.0;
vd.template getLastProp()[1] = 0.0;
vd.template getLastProp()[2] = 0.0;
vd.template getLastProp()[0] = 0.0;
vd.template getLastProp()[1] = 0.0;
vd.template getLastProp()[2] = 0.0;
++obstacle_box;
}
vd.map();
vd.write("Recipient");
openfpm_finalize();
return 0;
// Now that we fill the vector with particles
ModelCustom md;
vd.addComputationCosts(md);
vd.getDecomposition().decompose();
vd.map();
vd.ghost_get();
/*! \cond [get_verlet_crs] \endcond */
auto NN = vd.getVerletCrs(r_gskin);
/*! \cond [get_verlet_crs] \endcond */
size_t write = 0;
size_t it = 0;
size_t it_reb = 0;
double t = 0.0;
double tot_disp = 0.0;
double max_disp;
while (t <= t_end)
{
Vcluster & v_cl = create_vcluster();
timer it_time;
/*! \cond [update_verlet] \endcond */
it_reb++;
if (2*tot_disp >= skin)
{
vd.map();
if (it_reb > 200)
{
ModelCustom md;
vd.addComputationCosts(md);
vd.getDecomposition().redecompose(200);
vd.map();
it_reb = 0;
if (v_cl.getProcessUnitID() == 0)
std::cout << "REBALANCED " << std::endl;
}
// Calculate pressure from the density
EqState(vd);
vd.ghost_get();
/*! \cond [update_verlet_crs] \endcond */
vd.updateVerlet(NN,r_gskin,VL_CRS_SYMMETRIC);
/*! \cond [update_verlet_crs] \endcond */
tot_disp = 0.0;
if (v_cl.getProcessUnitID() == 0)
std::cout << "RECONSTRUCT Verlet " << std::endl;
}
else
{
// Calculate pressure from the density
EqState(vd);
vd.ghost_get(SKIP_LABELLING);
}
/*! \cond [update_verlet] \endcond */
double max_visc = 0.0;
// Calc forces
calc_forces(vd,NN,max_visc);
// Get the maximum viscosity term across processors
v_cl.max(max_visc);
v_cl.execute();
// Calculate delta t integration
double dt = calc_deltaT(vd,max_visc);
/*! \cond [pass_ver_eu] \endcond */
// VerletStep or euler step
it++;
if (it < 40)
verlet_int(vd,dt,max_disp);
else
{
euler_int(vd,dt,max_disp);
it = 0;
}
tot_disp += max_disp;
/*! \cond [pass_ver_eu] \endcond */
t += dt;
if (write < t*100)
{
vd.deleteGhost();
vd.write("Geometry",write);
vd.ghost_get(SKIP_LABELLING);
write++;
if (v_cl.getProcessUnitID() == 0)
std::cout << "TIME: " << t << " write " << it_time.getwct() << " " << v_cl.getProcessUnitID() << " TOT disp: " << tot_disp << " " << cnt << std::endl;
}
else
{
if (v_cl.getProcessUnitID() == 0)
std::cout << "TIME: " << t << " " << it_time.getwct() << " " << v_cl.getProcessUnitID() << " TOT disp: " << tot_disp << " " << cnt << std::endl;
}
}
openfpm_finalize();
}
* \endhtmlonly
*
*
*/
//#define SE_CLASS1
//#define STOP_ON_ERROR
#include "Vector/vector_dist.hpp"
#include