Commit 528ec206 authored by incardon's avatar incardon

examples SPH update

parent 065f9a27
......@@ -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
*
*/
......
......@@ -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
......
/*!
* \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