/*! * \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 * * \endhtmlonly * * */ //#define SE_CLASS1 //#define STOP_ON_ERROR #include "Vector/vector_dist.hpp" #include #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(); }