64 #include "Vector/vector_dist.hpp"
66 #include "Draw/DrawParticles.hpp"
120 const double dp = 0.0085;
126 const double coeff_sound = 20.0;
129 const double gamma_ = 7.0;
132 const double H = 0.0147224318643;
135 const double Eta2 = 0.01 * H*H;
138 const double visco = 0.1;
144 const double MassFluid = 0.000614125;
147 const double MassBound = 0.000614125;
151 const double t_end = 0.001;
153 const double t_end = 1.5;
157 const double gravity = 9.81;
160 const double rho_zero = 1000.0;
166 const double CFLnumber = 0.2;
169 const double DtMin = 0.00001;
172 const double RhoMin = 700.0;
175 const double RhoMax = 1300.0;
178 double max_fluid_height = 0.0;
183 const size_t type = 0;
189 const int rho_prev = 2;
192 const int Pressure = 3;
201 const int velocity = 6;
204 const int velocity_prev = 7;
223 template<
typename Decomposition,
typename vector>
inline void addComputation(
Decomposition & dec,
228 if (vd.template getProp<type>(p) == FLUID)
229 dec.addComputationCost(v,4);
231 dec.addComputationCost(v,3);
234 template<
typename Decomposition>
inline void applyModel(
Decomposition & dec,
size_t v)
236 dec.setSubSubDomainComputationCost(v, dec.getSubSubDomainComputationCost(v) * dec.getSubSubDomainComputationCost(v));
239 double distributionTol()
270 double rho_a = vd.template getProp<rho>(a);
271 double rho_frac = rho_a / rho_zero;
273 vd.template getProp<Pressure>(a) = B*( rho_frac*rho_frac*rho_frac*rho_frac*rho_frac*rho_frac*rho_frac - 1.0);
297 const double a2 = 1.0/M_PI/H/H/H;
299 inline double Wab(
double r)
304 return (1.0 - 3.0/2.0*r*r + 3.0/4.0*r*r*r)*a2;
306 return (1.0/4.0*(2.0 - r*r)*(2.0 - r*r)*(2.0 - r*r))*a2;
328 const double c1 = -3.0/M_PI/H/H/H/H;
329 const double d1 = 9.0/4.0/M_PI/H/H/H/H;
330 const double c2 = -3.0/4.0/M_PI/H/H/H/H;
331 const double a2_4 = 0.25*a2;
339 double qq2 = qq * qq;
340 double fac1 = (c1*qq + d1*qq2)/r;
341 double b1 = (qq < 1.0)?1.0f:0.0f;
343 double wqq = (2.0 - qq);
344 double fac2 = c2 * wqq * wqq / r;
345 double b2 = (qq >= 1.0 && qq < 2.0)?1.0f:0.0f;
347 double factor = (b1*fac1 + b2*fac2);
349 DW.
get(0) = factor * dx.
get(0);
350 DW.
get(1) = factor * dx.
get(1);
351 DW.
get(2) = factor * dx.
get(2);
375 inline double Tensile(
double r,
double rhoa,
double rhob,
double prs1,
double prs2)
383 double wqq2=wqq1*wqq1;
385 wab=a2_4*(wqq2*wqq1);
392 wab=a2*(1.0f-1.5f*wqq2+0.75f*wqq3);
396 double fab=wab*W_dap;
398 const double tensilp1=(prs1/(rhoa*rhoa))*(prs1>0? 0.01: -0.2);
399 const double tensilp2=(prs2/(rhob*rhob))*(prs2>0? 0.01: -0.2);
401 return (fab*(tensilp1+tensilp2));
425 const double dot_rr2 = dot/(rr2+Eta2);
426 visc=std::max(dot_rr2,visc);
430 const float amubar=H*dot_rr2;
431 const float robar=(rhoa+rhob)*0.5f;
432 const float pi_visc=(-visco*cbar*amubar/robar);
457 template<
typename CellList>
inline double calc_forces(
particles & vd,
CellList & NN,
double & max_visc)
466 while (part.isNext())
475 double massa = (vd.
getProp<type>(a) == FLUID)?MassFluid:MassBound;
478 double rhoa = vd.
getProp<rho>(a);
481 double Pa = vd.
getProp<Pressure>(a);
487 vd.template getProp<force>(a)[0] = 0.0;
488 vd.template getProp<force>(a)[1] = 0.0;
489 vd.template getProp<force>(a)[2] = -gravity;
490 vd.template getProp<drho>(a) = 0.0;
493 if (vd.
getProp<type>(a) != FLUID)
497 auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.
getPos(a)));
500 while (Np.isNext() ==
true)
509 if (a.getKey() == b) {++Np;
continue;};
512 double massb = (vd.
getProp<type>(b) == FLUID)?MassFluid:MassBound;
518 double Pb = vd.
getProp<Pressure>(b);
519 double rhob = vd.
getProp<rho>(b);
524 double r2 = norm2(dr);
537 const double dot = dr.get(0)*dv.
get(0) + dr.get(1)*dv.
get(1) + dr.get(2)*dv.
get(2);
538 const double dot_rr2 = dot/(r2+Eta2);
539 max_visc=std::max(dot_rr2,max_visc);
552 auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.
getPos(a)));
555 while (Np.isNext() ==
true)
564 if (a.getKey() == b) {++Np;
continue;};
566 double massb = (vd.
getProp<type>(b) == FLUID)?MassFluid:MassBound;
568 double Pb = vd.
getProp<Pressure>(b);
569 double rhob = vd.
getProp<rho>(b);
574 double r2 = norm2(dr);
586 double factor = - massb*((vd.
getProp<Pressure>(a) + vd.
getProp<Pressure>(b)) / (rhoa * rhob) + Tensile(r,rhoa,rhob,Pa,Pb) + Pi(dr,r2,v_rel,rhoa,rhob,massb,visc));
588 vd.
getProp<force>(a)[0] += factor * DW.
get(0);
589 vd.
getProp<force>(a)[1] += factor * DW.
get(1);
590 vd.
getProp<force>(a)[2] += factor * DW.
get(2);
621 void max_acceleration_and_velocity(
particles & vd,
double & max_acc,
double & max_vel)
626 while (part.isNext())
631 double acc2 = norm2(acc);
634 double vel2 = norm2(vel);
644 max_acc = sqrt(max_acc);
645 max_vel = sqrt(max_vel);
647 Vcluster & v_cl = create_vcluster();
678 double calc_deltaT(
particles & vd,
double ViscDtMax)
682 max_acceleration_and_velocity(vd,Maxacc,Maxvel);
685 const double dt_f = (Maxacc)?sqrt(H/Maxacc):std::numeric_limits<int>::max();
688 const double dt_cv = H/(std::max(cbar,Maxvel*10.) + H*ViscDtMax);
691 double dt=double(CFLnumber)*std::min(dt_f,dt_cv);
735 void verlet_int(
particles & vd,
double dt)
743 double dt205 = dt*dt*0.5;
747 while (part.isNext())
753 if (vd.template getProp<type>(a) == BOUNDARY)
756 double rhop = vd.template getProp<rho>(a);
759 vd.template getProp<velocity>(a)[0] = 0.0;
760 vd.template getProp<velocity>(a)[1] = 0.0;
761 vd.template getProp<velocity>(a)[2] = 0.0;
762 vd.template getProp<rho>(a) = vd.template getProp<rho_prev>(a) + dt2*vd.template getProp<drho>(a);
764 vd.template getProp<rho_prev>(a) = rhop;
771 double dx = vd.template getProp<velocity>(a)[0]*dt + vd.template getProp<force>(a)[0]*dt205;
772 double dy = vd.template getProp<velocity>(a)[1]*dt + vd.template getProp<force>(a)[1]*dt205;
773 double dz = vd.template getProp<velocity>(a)[2]*dt + vd.template getProp<force>(a)[2]*dt205;
779 double velX = vd.template getProp<velocity>(a)[0];
780 double velY = vd.template getProp<velocity>(a)[1];
781 double velZ = vd.template getProp<velocity>(a)[2];
782 double rhop = vd.template getProp<rho>(a);
784 vd.template getProp<velocity>(a)[0] = vd.template getProp<velocity_prev>(a)[0] + vd.template getProp<force>(a)[0]*dt2;
785 vd.template getProp<velocity>(a)[1] = vd.template getProp<velocity_prev>(a)[1] + vd.template getProp<force>(a)[1]*dt2;
786 vd.template getProp<velocity>(a)[2] = vd.template getProp<velocity_prev>(a)[2] + vd.template getProp<force>(a)[2]*dt2;
787 vd.template getProp<rho>(a) = vd.template getProp<rho_prev>(a) + dt2*vd.template getProp<drho>(a);
790 if (vd.
getPos(a)[0] < 0.000263878 || vd.
getPos(a)[1] < 0.000263878 || vd.
getPos(a)[2] < 0.000263878 ||
791 vd.
getPos(a)[0] > 0.000263878+1.59947 || vd.
getPos(a)[1] > 0.000263878+0.672972 || vd.
getPos(a)[2] > 0.000263878+0.903944 ||
792 vd.template getProp<rho>(a) < RhoMin || vd.template getProp<rho>(a) > RhoMax)
794 to_remove.add(a.getKey());
797 vd.template getProp<velocity_prev>(a)[0] = velX;
798 vd.template getProp<velocity_prev>(a)[1] = velY;
799 vd.template getProp<velocity_prev>(a)[2] = velZ;
800 vd.template getProp<rho_prev>(a) = rhop;
812 void euler_int(
particles & vd,
double dt)
820 double dt205 = dt*dt*0.5;
824 while (part.isNext())
830 if (vd.template getProp<type>(a) == BOUNDARY)
833 double rhop = vd.template getProp<rho>(a);
836 vd.template getProp<velocity>(a)[0] = 0.0;
837 vd.template getProp<velocity>(a)[1] = 0.0;
838 vd.template getProp<velocity>(a)[2] = 0.0;
839 vd.template getProp<rho>(a) = vd.template getProp<rho>(a) + dt*vd.template getProp<drho>(a);
841 vd.template getProp<rho_prev>(a) = rhop;
848 double dx = vd.template getProp<velocity>(a)[0]*dt + vd.template getProp<force>(a)[0]*dt205;
849 double dy = vd.template getProp<velocity>(a)[1]*dt + vd.template getProp<force>(a)[1]*dt205;
850 double dz = vd.template getProp<velocity>(a)[2]*dt + vd.template getProp<force>(a)[2]*dt205;
856 double velX = vd.template getProp<velocity>(a)[0];
857 double velY = vd.template getProp<velocity>(a)[1];
858 double velZ = vd.template getProp<velocity>(a)[2];
859 double rhop = vd.template getProp<rho>(a);
861 vd.template getProp<velocity>(a)[0] = vd.template getProp<velocity>(a)[0] + vd.template getProp<force>(a)[0]*dt;
862 vd.template getProp<velocity>(a)[1] = vd.template getProp<velocity>(a)[1] + vd.template getProp<force>(a)[1]*dt;
863 vd.template getProp<velocity>(a)[2] = vd.template getProp<velocity>(a)[2] + vd.template getProp<force>(a)[2]*dt;
864 vd.template getProp<rho>(a) = vd.template getProp<rho>(a) + dt*vd.template getProp<drho>(a);
867 if (vd.
getPos(a)[0] < 0.000263878 || vd.
getPos(a)[1] < 0.000263878 || vd.
getPos(a)[2] < 0.000263878 ||
868 vd.
getPos(a)[0] > 0.000263878+1.59947 || vd.
getPos(a)[1] > 0.000263878+0.672972 || vd.
getPos(a)[2] > 0.000263878+0.903944 ||
869 vd.template getProp<rho>(a) < RhoMin || vd.template getProp<rho>(a) > RhoMax)
871 to_remove.add(a.getKey());
874 vd.template getProp<velocity_prev>(a)[0] = velX;
875 vd.template getProp<velocity_prev>(a)[1] = velY;
876 vd.template getProp<velocity_prev>(a)[2] = velZ;
877 vd.template getProp<rho_prev>(a) = rhop;
916 template<
typename Vector,
typename CellList>
917 inline void sensor_pressure(
Vector & vd,
922 Vcluster & v_cl = create_vcluster();
926 for (
size_t i = 0 ; i < probes.size() ; i++)
928 float press_tmp = 0.0f;
932 if (vd.getDecomposition().isLocal(probes.get(i)) ==
true)
938 auto itg = NN.template getNNIterator<NO_CHECK>(NN.getCell(probes.get(i)));
944 if (vd.template getProp<type>(q) != FLUID)
955 double r = sqrt(norm2(xp - xq));
957 double ker = Wab(r) * (MassFluid / rho_zero);
964 press_tmp += vd.template getProp<Pressure>(q) * ker;
975 press_tmp = 1.0 / tot_ker * press_tmp;
986 press_t.last().add(press_tmp);
992 int main(
int argc,
char* argv[])
1012 openfpm_init(&argc,&argv);
1018 probes.add({0.8779,0.3,0.02});
1019 probes.add({0.754,0.31,0.02});
1022 Box<3,double> domain({-0.05,-0.05,-0.05},{1.7010,0.7065,0.5025});
1023 size_t sz[3] = {207,90,66};
1026 W_dap = 1.0/Wab(H/1.5);
1029 size_t bc[3]={NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
1073 particles vd(0,domain,bc,g,DEC_GRAN(512));
1115 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});
1121 max_fluid_height = fluid_it.getBoxMargins().
getHigh(2);
1122 h_swl = fluid_it.getBoxMargins().
getHigh(2) - fluid_it.getBoxMargins().
getLow(2);
1123 B = (coeff_sound)*(coeff_sound)*gravity*h_swl*rho_zero / gamma_;
1124 cbar = coeff_sound * sqrt(gravity * h_swl);
1127 while (fluid_it.isNext())
1133 vd.getLastPos()[0] = fluid_it.get().get(0);
1134 vd.getLastPos()[1] = fluid_it.get().get(1);
1135 vd.getLastPos()[2] = fluid_it.get().get(2);
1138 vd.template getLastProp<type>() = FLUID;
1147 vd.template getLastProp<Pressure>() = rho_zero * gravity * (max_fluid_height - fluid_it.get().get(2));
1149 vd.template getLastProp<rho>() = pow(vd.template getLastProp<Pressure>() / B + 1, 1.0/gamma_) * rho_zero;
1150 vd.template getLastProp<rho_prev>() = vd.template getLastProp<rho>();
1151 vd.template getLastProp<velocity>()[0] = 0.0;
1152 vd.template getLastProp<velocity>()[1] = 0.0;
1153 vd.template getLastProp<velocity>()[2] = 0.0;
1155 vd.template getLastProp<velocity_prev>()[0] = 0.0;
1156 vd.template getLastProp<velocity_prev>()[1] = 0.0;
1157 vd.template getLastProp<velocity_prev>()[2] = 0.0;
1191 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});
1192 Box<3,double> recipient2({dp,dp,dp},{1.6-dp/2.0,0.67-dp/2.0,0.4+dp/2.0});
1194 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});
1195 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});
1196 Box<3,double> obstacle3({0.9+dp,0.24,0.0},{1.02,0.36,0.45});
1199 holes.add(recipient2);
1200 holes.add(obstacle1);
1203 while (bound_box.isNext())
1207 vd.getLastPos()[0] = bound_box.get().get(0);
1208 vd.getLastPos()[1] = bound_box.get().get(1);
1209 vd.getLastPos()[2] = bound_box.get().get(2);
1211 vd.template getLastProp<type>() = BOUNDARY;
1212 vd.template getLastProp<rho>() = rho_zero;
1213 vd.template getLastProp<rho_prev>() = rho_zero;
1214 vd.template getLastProp<velocity>()[0] = 0.0;
1215 vd.template getLastProp<velocity>()[1] = 0.0;
1216 vd.template getLastProp<velocity>()[2] = 0.0;
1218 vd.template getLastProp<velocity_prev>()[0] = 0.0;
1219 vd.template getLastProp<velocity_prev>()[1] = 0.0;
1220 vd.template getLastProp<velocity_prev>()[2] = 0.0;
1247 while (obstacle_box.isNext())
1251 vd.getLastPos()[0] = obstacle_box.get().get(0);
1252 vd.getLastPos()[1] = obstacle_box.get().get(1);
1253 vd.getLastPos()[2] = obstacle_box.get().get(2);
1255 vd.template getLastProp<type>() = BOUNDARY;
1256 vd.template getLastProp<rho>() = rho_zero;
1257 vd.template getLastProp<rho_prev>() = rho_zero;
1258 vd.template getLastProp<velocity>()[0] = 0.0;
1259 vd.template getLastProp<velocity>()[1] = 0.0;
1260 vd.template getLastProp<velocity>()[2] = 0.0;
1262 vd.template getLastProp<velocity_prev>()[0] = 0.0;
1263 vd.template getLastProp<velocity_prev>()[1] = 0.0;
1264 vd.template getLastProp<velocity_prev>()[2] = 0.0;
1359 vd.addComputationCosts(md);
1360 vd.getDecomposition().decompose();
1365 vd.ghost_get<type,rho,Pressure,velocity>();
1367 auto NN = vd.getCellList(2*H);
1393 Vcluster & v_cl = create_vcluster();
1404 vd.addComputationCosts(md);
1405 vd.getDecomposition().decompose();
1408 std::cout <<
"REBALANCED " << std::endl;
1416 double max_visc = 0.0;
1418 vd.ghost_get<type,rho,Pressure,velocity>();
1421 calc_forces(vd,NN,max_visc);
1428 double dt = calc_deltaT(vd,max_visc);
1445 sensor_pressure(vd,NN,press_t,probes);
1447 vd.write(
"Geometry",write);
1451 std::cout <<
"TIME: " << t <<
" write " << it_time.
getwct() <<
" " << v_cl.
getProcessUnitID() <<
" " << cnt << std::endl;
1456 std::cout <<
"TIME: " << t <<
" " << it_time.
getwct() <<
" " << v_cl.
getProcessUnitID() <<
" " << 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
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.
Sparse Matrix implementation stub object when OpenFPM is compiled with no linear algebra support...
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
This class implement the point shape in an N-dimensional space.
void max(T &num)
Get the maximum number across all processors (or reduction with infinity norm)
double getwct()
Return the elapsed real time.
void updateCellList(CellL &cell_list, bool no_se3=false)
Update a cell list using the stored particles.
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.
const T & get(size_t i) const
Get coordinate.
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.
Model for Dynamic load balancing.
Class for FAST cell list implementation.
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.