56 #include "Vector/vector_dist.hpp"
58 #include "Draw/DrawParticles.hpp"
218 const double dp = 0.0085;
224 const double coeff_sound = 20.0;
227 const double gamma_ = 7.0;
230 const double H = 0.0147224318643;
232 const double FourH2 = 4.0 * H*H;
235 const double Eta2 = 0.01 * H*H;
238 const double visco = 0.1;
244 const double MassFluid = 0.000614125;
247 const double MassBound = 0.000614125;
251 const double t_end = 1.5;
253 const double t_end = 0.001;
257 const double gravity = 9.81;
260 const double rho_zero = 1000.0;
266 const double CFLnumber = 0.2;
269 const double DtMin = 0.00001;
272 const double RhoMin = 700.0;
275 const double RhoMax = 1300.0;
278 double max_fluid_height = 0.0;
283 const size_t type = 0;
286 const size_t nn_num = 1;
292 const int rho_prev = 3;
295 const int Pressure = 4;
304 const int velocity = 7;
307 const int velocity_prev = 8;
314 typedef vector_dist<3,double,aggregate<int, int,double, double, double, double, double[3], double[3], double[3]> >
particles;
325 template<
typename Decomposition,
typename vector>
inline void addComputation(
Decomposition & dec, vector & vd,
size_t v,
size_t p)
327 if (vd.template getProp<type>(p) == FLUID )
328 dec.addComputationCost(v,4);
330 dec.addComputationCost(v,3);
333 template<
typename Decomposition>
inline void applyModel(
Decomposition & dec,
size_t v)
335 dec.setSubSubDomainComputationCost(v, dec.getSubSubDomainComputationCost(v) * dec.getSubSubDomainComputationCost(v));
338 float distributionTol()
347 template<
typename Decomposition,
typename vector>
inline void addComputation(
Decomposition & dec, vector & vd,
size_t v,
size_t p)
349 dec.addComputationCost(v,vd.template getProp<nn_num>(p) + 4);
352 template<
typename Decomposition>
inline void applyModel(
Decomposition & dec,
size_t v)
356 float distributionTol()
370 double rho_a = vd.template getProp<rho>(a);
371 double rho_frac = rho_a / rho_zero;
373 vd.template getProp<Pressure>(a) = B*( rho_frac*rho_frac*rho_frac*rho_frac*rho_frac*rho_frac*rho_frac - 1.0);
379 const double a2 = 1.0/M_PI/H/H/H;
381 inline double Wab(
double r)
386 return (1.0 - 3.0/2.0*r*r + 3.0/4.0*r*r*r)*a2;
388 return (1.0/4.0*(2.0 - r*r)*(2.0 - r*r)*(2.0 - r*r))*a2;
393 const double c1 = -3.0/M_PI/H/H/H/H;
394 const double d1 = 9.0/4.0/M_PI/H/H/H/H;
395 const double c2 = -3.0/4.0/M_PI/H/H/H/H;
396 const double a2_4 = 0.25*a2;
404 double qq2 = qq * qq;
405 double fac1 = (c1*qq + d1*qq2)/r;
406 double b1 = (qq < 1.0)?1.0f:0.0f;
408 double wqq = (2.0 - qq);
409 double fac2 = c2 * wqq * wqq / r;
410 double b2 = (qq >= 1.0 && qq < 2.0)?1.0f:0.0f;
412 double factor = (b1*fac1 + b2*fac2);
414 DW.
get(0) = factor * dx.
get(0);
415 DW.
get(1) = factor * dx.
get(1);
416 DW.
get(2) = factor * dx.
get(2);
419 inline double Tensile(
double r,
double rhoa,
double rhob,
double prs1,
double prs2)
427 double wqq2=wqq1*wqq1;
429 wab=a2_4*(wqq2*wqq1);
436 wab=a2*(1.0f-1.5f*wqq2+0.75f*wqq3);
440 double fab=wab*W_dap;
442 const double tensilp1=(prs1/(rhoa*rhoa))*(prs1>0? 0.01: -0.2);
443 const double tensilp2=(prs2/(rhob*rhob))*(prs2>0? 0.01: -0.2);
445 return (fab*(tensilp1+tensilp2));
452 const double dot_rr2 = dot/(rr2+Eta2);
453 visc=std::max(dot_rr2,visc);
457 const float amubar=H*dot_rr2;
458 const float robar=(rhoa+rhob)*0.5f;
459 const float pi_visc=(-visco*cbar*amubar/robar);
468 template<
typename VerletList>
inline double calc_forces(
particles & vd,
VerletList & NN,
double & max_visc)
480 vd.template getProp<force>(p)[0] = 0.0;
481 vd.template getProp<force>(p)[1] = 0.0;
482 vd.template getProp<force>(p)[2] = -gravity;
483 vd.template getProp<drho>(p) = 0.0;
494 while (itg2.isNext())
500 vd.template getProp<force>(p)[0] = 0.0;
501 vd.template getProp<force>(p)[1] = 0.0;
502 vd.template getProp<force>(p)[2] = 0.0;
503 vd.template getProp<drho>(p) = 0.0;
516 while (part.isNext())
525 size_t typea = vd.
getProp<type>(a);
528 double massa = (vd.
getProp<type>(a) == FLUID)?MassFluid:MassBound;
531 double rhoa = vd.
getProp<rho>(a);
534 double Pa = vd.
getProp<Pressure>(a);
540 auto Np = NN.template getNNIterator<NO_CHECK>(a);
545 while (Np.isNext() ==
true)
556 float r2 = norm2(dr);
559 if (r2 < FourH2 && r2 > 1e-18)
563 unsigned int typeb = vd.
getProp<type>(b);
565 double massb = (typeb == FLUID)?MassFluid:MassBound;
567 double Pb = vd.
getProp<Pressure>(b);
568 double rhob = vd.
getProp<rho>(b);
577 double factor = - ((Pa + Pb) / (rhoa * rhob) + Tensile(r,rhoa,rhob,Pa,Pb) + Pi(dr,r2,v_rel,rhoa,rhob,massb,visc));
580 factor = (typea == BOUNDARY && typeb == BOUNDARY)?0.0f:factor;
582 vd.
getProp<force>(a)[0] += massb * factor * DW.
get(0);
583 vd.
getProp<force>(a)[1] += massb * factor * DW.
get(1);
584 vd.
getProp<force>(a)[2] += massb * factor * DW.
get(2);
586 vd.
getProp<force>(b)[0] -= massa * factor * DW.
get(0);
587 vd.
getProp<force>(b)[1] -= massa * factor * DW.
get(1);
588 vd.
getProp<force>(b)[2] -= massa * factor * DW.
get(2);
590 double scal = (v_rel.get(0)*DW.
get(0)+v_rel.get(1)*DW.
get(1)+v_rel.get(2)*DW.
get(2));
593 scal = (typea == BOUNDARY && typeb == BOUNDARY)?0.0f:scal;
595 vd.
getProp<drho>(a) += massb*scal;
596 vd.
getProp<drho>(b) += massa*scal;
613 vd.template ghost_put<add_,drho,force>();
618 void max_acceleration_and_velocity(
particles & vd,
double & max_acc,
double & max_vel)
623 while (part.isNext())
628 double acc2 = norm2(acc);
631 double vel2 = norm2(vel);
641 max_acc = sqrt(max_acc);
642 max_vel = sqrt(max_vel);
644 Vcluster & v_cl = create_vcluster();
650 double calc_deltaT(
particles & vd,
double ViscDtMax)
654 max_acceleration_and_velocity(vd,Maxacc,Maxvel);
657 const double dt_f = (Maxacc)?sqrt(H/Maxacc):std::numeric_limits<int>::max();
660 const double dt_cv = H/(std::max(cbar,Maxvel*10.) + H*ViscDtMax);
663 double dt=double(CFLnumber)*std::min(dt_f,dt_cv);
676 void verlet_int(
particles & vd,
double dt,
double & max_disp)
685 double dt205 = dt*dt*0.5;
693 while (part.isNext())
699 if (vd.template getProp<type>(a) == BOUNDARY)
702 double rhop = vd.template getProp<rho>(a);
705 vd.template getProp<velocity>(a)[0] = 0.0;
706 vd.template getProp<velocity>(a)[1] = 0.0;
707 vd.template getProp<velocity>(a)[2] = 0.0;
708 double rhonew = vd.template getProp<rho_prev>(a) + dt2*vd.template getProp<drho>(a);
709 vd.template getProp<rho>(a) = (rhonew < rho_zero)?rho_zero:rhonew;
711 vd.template getProp<rho_prev>(a) = rhop;
718 double dx = vd.template getProp<velocity>(a)[0]*dt + vd.template getProp<force>(a)[0]*dt205;
719 double dy = vd.template getProp<velocity>(a)[1]*dt + vd.template getProp<force>(a)[1]*dt205;
720 double dz = vd.template getProp<velocity>(a)[2]*dt + vd.template getProp<force>(a)[2]*dt205;
727 double d2 = dx*dx + dy*dy + dz*dz;
729 max_disp = (max_disp > d2)?max_disp:d2;
732 double velX = vd.template getProp<velocity>(a)[0];
733 double velY = vd.template getProp<velocity>(a)[1];
734 double velZ = vd.template getProp<velocity>(a)[2];
735 double rhop = vd.template getProp<rho>(a);
737 vd.template getProp<velocity>(a)[0] = vd.template getProp<velocity_prev>(a)[0] + vd.template getProp<force>(a)[0]*dt2;
738 vd.template getProp<velocity>(a)[1] = vd.template getProp<velocity_prev>(a)[1] + vd.template getProp<force>(a)[1]*dt2;
739 vd.template getProp<velocity>(a)[2] = vd.template getProp<velocity_prev>(a)[2] + vd.template getProp<force>(a)[2]*dt2;
740 vd.template getProp<rho>(a) = vd.template getProp<rho_prev>(a) + dt2*vd.template getProp<drho>(a);
743 if (vd.
getPos(a)[0] < 0.000263878 || vd.
getPos(a)[1] < 0.000263878 || vd.
getPos(a)[2] < 0.000263878 ||
744 vd.
getPos(a)[0] > 0.000263878+1.59947 || vd.
getPos(a)[1] > 0.000263878+0.672972 || vd.
getPos(a)[2] > 0.000263878+0.903944 ||
745 vd.template getProp<rho>(a) < RhoMin || vd.template getProp<rho>(a) > RhoMax)
747 to_remove.add(a.getKey());
755 vd.template getProp<velocity_prev>(a)[0] = velX;
756 vd.template getProp<velocity_prev>(a)[1] = velY;
757 vd.template getProp<velocity_prev>(a)[2] = velZ;
758 vd.template getProp<rho_prev>(a) = rhop;
765 Vcluster & v_cl = create_vcluster();
769 max_disp = sqrt(max_disp);
778 void euler_int(
particles & vd,
double dt,
double & max_disp)
787 double dt205 = dt*dt*0.5;
793 while (part.isNext())
799 if (vd.template getProp<type>(a) == BOUNDARY)
802 double rhop = vd.template getProp<rho>(a);
805 vd.template getProp<velocity>(a)[0] = 0.0;
806 vd.template getProp<velocity>(a)[1] = 0.0;
807 vd.template getProp<velocity>(a)[2] = 0.0;
808 double rhonew = vd.template getProp<rho>(a) + dt*vd.template getProp<drho>(a);
809 vd.template getProp<rho>(a) = (rhonew < rho_zero)?rho_zero:rhonew;
811 vd.template getProp<rho_prev>(a) = rhop;
818 double dx = vd.template getProp<velocity>(a)[0]*dt + vd.template getProp<force>(a)[0]*dt205;
819 double dy = vd.template getProp<velocity>(a)[1]*dt + vd.template getProp<force>(a)[1]*dt205;
820 double dz = vd.template getProp<velocity>(a)[2]*dt + vd.template getProp<force>(a)[2]*dt205;
826 double d2 = dx*dx + dy*dy + dz*dz;
828 max_disp = (max_disp > d2)?max_disp:d2;
830 double velX = vd.template getProp<velocity>(a)[0];
831 double velY = vd.template getProp<velocity>(a)[1];
832 double velZ = vd.template getProp<velocity>(a)[2];
833 double rhop = vd.template getProp<rho>(a);
835 vd.template getProp<velocity>(a)[0] = vd.template getProp<velocity>(a)[0] + vd.template getProp<force>(a)[0]*dt;
836 vd.template getProp<velocity>(a)[1] = vd.template getProp<velocity>(a)[1] + vd.template getProp<force>(a)[1]*dt;
837 vd.template getProp<velocity>(a)[2] = vd.template getProp<velocity>(a)[2] + vd.template getProp<force>(a)[2]*dt;
838 vd.template getProp<rho>(a) = vd.template getProp<rho>(a) + dt*vd.template getProp<drho>(a);
841 if (vd.
getPos(a)[0] < 0.000263878 || vd.
getPos(a)[1] < 0.000263878 || vd.
getPos(a)[2] < 0.000263878 ||
842 vd.
getPos(a)[0] > 0.000263878+1.59947 || vd.
getPos(a)[1] > 0.000263878+0.672972 || vd.
getPos(a)[2] > 0.000263878+0.903944 ||
843 vd.template getProp<rho>(a) < RhoMin || vd.template getProp<rho>(a) > RhoMax)
845 to_remove.add(a.getKey());
848 vd.template getProp<velocity_prev>(a)[0] = velX;
849 vd.template getProp<velocity_prev>(a)[1] = velY;
850 vd.template getProp<velocity_prev>(a)[2] = velZ;
851 vd.template getProp<rho_prev>(a) = rhop;
856 Vcluster & v_cl = create_vcluster();
860 max_disp = sqrt(max_disp);
866 int main(
int argc,
char* argv[])
869 openfpm_init(&argc,&argv);
872 Box<3,double> domain({-0.05,-0.05,-0.05},{1.7010,0.7065,0.5025});
873 size_t sz[3] = {207,90,66};
876 W_dap = 1.0/Wab(H/1.5);
879 size_t bc[3]={NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
883 double skin = 0.25 * 2*H;
884 double r_gskin = 2*H + skin;
900 particles vd(0,domain,bc,g,BIND_DEC_TO_GHOST);
906 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});
912 max_fluid_height = fluid_it.getBoxMargins().
getHigh(2);
913 h_swl = fluid_it.getBoxMargins().
getHigh(2) - fluid_it.getBoxMargins().
getLow(2);
914 B = (coeff_sound)*(coeff_sound)*gravity*h_swl*rho_zero / gamma_;
915 cbar = coeff_sound * sqrt(gravity * h_swl);
918 while (fluid_it.isNext())
929 vd.template getLastProp<type>() = FLUID;
938 vd.template getLastProp<Pressure>() = rho_zero * gravity * (max_fluid_height - fluid_it.get().get(2));
940 vd.template getLastProp<rho>() = pow(vd.template getLastProp<Pressure>() / B + 1, 1.0/gamma_) * rho_zero;
941 vd.template getLastProp<rho_prev>() = vd.template getLastProp<rho>();
942 vd.template getLastProp<velocity>()[0] = 0.0;
943 vd.template getLastProp<velocity>()[1] = 0.0;
944 vd.template getLastProp<velocity>()[2] = 0.0;
946 vd.template getLastProp<velocity_prev>()[0] = 0.0;
947 vd.template getLastProp<velocity_prev>()[1] = 0.0;
948 vd.template getLastProp<velocity_prev>()[2] = 0.0;
955 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});
956 Box<3,double> recipient2({dp,dp,dp},{1.6-dp/2.0,0.67-dp/2.0,0.4+dp/2.0});
958 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});
959 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});
963 holes.add(recipient2);
964 holes.add(obstacle1);
967 while (bound_box.isNext())
975 vd.template getLastProp<type>() = BOUNDARY;
976 vd.template getLastProp<rho>() = rho_zero;
977 vd.template getLastProp<rho_prev>() = rho_zero;
978 vd.template getLastProp<velocity>()[0] = 0.0;
979 vd.template getLastProp<velocity>()[1] = 0.0;
980 vd.template getLastProp<velocity>()[2] = 0.0;
982 vd.template getLastProp<velocity_prev>()[0] = 0.0;
983 vd.template getLastProp<velocity_prev>()[1] = 0.0;
984 vd.template getLastProp<velocity_prev>()[2] = 0.0;
991 while (obstacle_box.isNext())
995 vd.
getLastPos()[0] = obstacle_box.get().get(0);
996 vd.
getLastPos()[1] = obstacle_box.get().get(1);
997 vd.
getLastPos()[2] = obstacle_box.get().get(2);
999 vd.template getLastProp<type>() = BOUNDARY;
1000 vd.template getLastProp<rho>() = rho_zero;
1001 vd.template getLastProp<rho_prev>() = rho_zero;
1002 vd.template getLastProp<velocity>()[0] = 0.0;
1003 vd.template getLastProp<velocity>()[1] = 0.0;
1004 vd.template getLastProp<velocity>()[2] = 0.0;
1006 vd.template getLastProp<velocity_prev>()[0] = 0.0;
1007 vd.template getLastProp<velocity_prev>()[1] = 0.0;
1008 vd.template getLastProp<velocity_prev>()[2] = 0.0;
1022 auto & v_cl = create_vcluster();
1025 std::cout <<
"SUM: " << tot_part << std::endl;
1034 vd.
ghost_get<type,rho,Pressure,velocity>();
1044 double tot_disp = 0.0;
1048 Vcluster & v_cl = create_vcluster();
1054 if (2*tot_disp >= skin)
1071 std::cout <<
"REBALANCED " << std::endl;
1078 vd.
ghost_get<type,rho,Pressure,velocity>();
1087 std::cout <<
"RECONSTRUCT Verlet " << std::endl;
1094 vd.
ghost_get<type,rho,Pressure,velocity>(SKIP_LABELLING);
1099 double max_visc = 0.0;
1102 calc_forces(vd,NN,max_visc);
1109 double dt = calc_deltaT(vd,max_visc);
1116 verlet_int(vd,dt,max_disp);
1119 euler_int(vd,dt,max_disp);
1123 tot_disp += max_disp;
1132 vd.
write_frame(
"Geometry",write,VTK_WRITER | FORMAT_BINARY);
1134 vd.
ghost_get<type,rho,Pressure,velocity>(SKIP_LABELLING);
1138 std::cout <<
"TIME: " << t <<
" write " << it_time.
getwct() <<
" " << v_cl.
getProcessUnitID() <<
" TOT disp: " << tot_disp <<
" " << cnt << std::endl;
1143 std::cout <<
"TIME: " << t <<
" " << it_time.
getwct() <<
" " << v_cl.
getProcessUnitID() <<
" TOT disp: " << tot_disp <<
" " << cnt << std::endl;
void sum(T &num)
Sum the numbers across all processors and get the result.
auto getProp(vect_dist_key_dx vec_key) -> decltype(v_prp.template get< id >(vec_key.getKey()))
Get the property of an element.
T getLow(int i) const
get the i-coordinate of the low bound interval of the box
openfpm::vector_key_iterator_seq< typename vrl::Mem_type_type::loc_index > getParticleIteratorCRS(vrl &NN)
Get a special particle iterator able to iterate across particles using symmetric crossing scheme...
size_t getProcessUnitID()
Get the process unit id.
static PointIteratorSkin< dim, T, Decomposition > DrawSkin(vector_dist< dim, T, aggr, layout, layout_base, Decomposition > &vd, size_t(&sz)[dim], Box< dim, T > &domain, Box< dim, T > &sub_A, Box< dim, T > &sub_B)
Draw particles in a box B excluding the area of a second box A (B - A)
void execute()
Execute all the requests.
Class for Verlet list implementation.
bool write_frame(std::string out, size_t iteration, int opt=VTK_WRITER)
Output particle position and properties.
Second model for dynamic load balancing.
auto getPos(vect_dist_key_dx vec_key) -> decltype(v_pos.template get< 0 >(vec_key.getKey()))
Get the position of an element.
T getHigh(int i) const
get the high interval of the box
Decomposition & getDecomposition()
Get the decomposition.
VerletL getVerletCrs(St r_cut)
for each particle get the symmetric verlet list
This class implement the point shape in an N-dimensional space.
vector_dist_iterator getGhostIterator() const
Get the iterator across the position of the ghost particles.
void max(T &num)
Get the maximum number across all processors (or reduction with infinity norm)
void deleteGhost()
Delete the particles on the ghost.
double getwct()
Return the elapsed real time.
Implementation of VCluster class.
This class define the domain decomposition interface.
void remove(openfpm::vector< size_t > &keys, size_t start=0)
Remove a set of elements from the distributed vector.
size_t size_local() const
return the local size of the vector
void updateVerlet(VerletList< dim, St, Mem_type, shift< dim, St > > &ver, St r_cut, size_t opt=VL_NON_SYMMETRIC)
for each particle get the verlet list
const T & get(size_t i) const
Get coordinate.
void addComputationCosts(const self &vd, Model md=Model())
Add the computation cost on the decomposition coming from the particles.
void map(size_t opt=NONE)
It move all the particles that does not belong to the local processor to the respective processor...
void ghost_get(size_t opt=WITH_POSITION)
It synchronize the properties and position of the ghost particles.
void add()
Add local particle.
This class represent an N-dimensional box.
vector_dist_iterator getDomainIterator() const
Get an iterator that traverse the particles in the domain.
vect_dist_key_dx get()
Get the actual key.
auto getLastPos() -> decltype(v_pos.template get< 0 >(0))
Get the position of the last element.
Model for Dynamic load balancing.
Class for cpu time benchmarking.
static PointIterator< dim, T, Decomposition > DrawBox(vector_dist< dim, T, aggr, layout, layout_base, Decomposition > &vd, size_t(&sz)[dim], Box< dim, T > &domain, Box< dim, T > &sub)
Draw particles in a box.