main.cpp 28.5 KB
Newer Older
incardon's avatar
incardon committed
1
/*!
incardon's avatar
incardon committed
2
 * \page Vector_7_sph_dlb Vector 7 SPH Dam break simulation with Dynamic load balacing
incardon's avatar
incardon committed
3
4
5
6
7
8
9
10
 *
 *
 * [TOC]
 *
 *
 * # SPH with Dynamic load Balancing # {#SPH_dlb}
 *
 *
incardon's avatar
incardon committed
11
12
13
 * 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.
incardon's avatar
incardon committed
14
15
16
17
 *
 * ## inclusion ## {#e0_v_inclusion}
 *
 * In order to use distributed vectors in our code we have to include the file Vector/vector_dist.hpp
incardon's avatar
incardon committed
18
19
 * we also include DrawParticles that has nice utilities to draw particles in parallel accordingly
 * to simple shapes
incardon's avatar
incardon committed
20
21
22
23
24
25
26
27
 *
 * \snippet Vector/7_SPH_dlb/main.cpp inclusion
 *
 */

//! \cond [inclusion] \endcond
#include "Vector/vector_dist.hpp"
#include <math.h>
incardon's avatar
incardon committed
28
#include "Draw/DrawParticles.hpp"
incardon's avatar
incardon committed
29
30
//! \cond [inclusion] \endcond

incardon's avatar
incardon committed
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
/*!
 * \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 */
incardon's avatar
incardon committed
69

incardon's avatar
incardon committed
70
// A constant to indicate boundary particles
incardon's avatar
incardon committed
71
72
#define BOUNDARY 0

incardon's avatar
incardon committed
73
74
// A constant to indicate fluid particles
#define FLUID 1
incardon's avatar
incardon committed
75

incardon's avatar
incardon committed
76
// initial spacing between particles dp in the formulas
incardon's avatar
incardon committed
77
const double dp = 0.0085;
incardon's avatar
incardon committed
78
79
// Maximum height of the fluid water
// is coing to be calculated and filled later on
incardon's avatar
incardon committed
80
double h_swl = 0.0;
incardon's avatar
incardon committed
81
82

// in the formulas indicated with c_s (constant used to calculate the sound speed)
incardon's avatar
incardon committed
83
const double coeff_sound = 20.0;
incardon's avatar
incardon committed
84
85

// gamma in the formulas
incardon's avatar
incardon committed
86
const double gamma_ = 7.0;
incardon's avatar
incardon committed
87
88

// sqrt(3.0*dp*dp) support of the kernel
incardon's avatar
incardon committed
89
const double H = 0.0147224318643;
incardon's avatar
incardon committed
90
91

// Eta in the formulas
incardon's avatar
incardon committed
92
const double Eta2 = 0.01 * H*H;
incardon's avatar
incardon committed
93
94


incardon's avatar
incardon committed
95
96
const double visco = 0.1;
double cbar = 0.0;
incardon's avatar
incardon committed
97
98

// Mass of the fluid particles
incardon's avatar
incardon committed
99
const double MassFluid = 0.000614125;
incardon's avatar
incardon committed
100
101

// Mass of the boundary particles
incardon's avatar
incardon committed
102
const double MassBound = 0.000614125;
incardon's avatar
incardon committed
103
104
105
106
107

// End simulation time
const double t_end = 1.5;

// Gravity acceleration
incardon's avatar
incardon committed
108
const double gravity = 9.81;
incardon's avatar
incardon committed
109
110

// Reference densitu 1000Kg/m^3
incardon's avatar
incardon committed
111
const double rho_zero = 1000.0;
incardon's avatar
incardon committed
112
113

// Filled later require h_swl, it is b in the formulas
incardon's avatar
incardon committed
114
double B = 0.0;
incardon's avatar
incardon committed
115
116

// Constant used to define time integration
incardon's avatar
incardon committed
117
const double CFLnumber = 0.2;
incardon's avatar
incardon committed
118
119

// Minimum T
incardon's avatar
incardon committed
120
121
const double DtMin = 0.00001;

incardon's avatar
incardon committed
122
123
124
125
126
127
// Minimum Rho allowed
const double RhoMin = 700.0;

// Maximum Rho allowed
const double RhoMax = 1300.0;

incardon's avatar
incardon committed
128
129
130
// Filled in initialization
double max_fluid_height = 0.0;

incardon's avatar
incardon committed
131
132
133
// Properties

// FLUID or BOUNDARY
incardon's avatar
incardon committed
134
const size_t type = 0;
incardon's avatar
incardon committed
135
136

// Density
incardon's avatar
incardon committed
137
const int rho = 1;
incardon's avatar
incardon committed
138
139

// Density at step n-1
incardon's avatar
incardon committed
140
const int rho_prev = 2;
incardon's avatar
incardon committed
141
142

// Pressure
incardon's avatar
incardon committed
143
const int Pressure = 3;
incardon's avatar
incardon committed
144
145

// Delta rho calculated in the force calculation
incardon's avatar
incardon committed
146
const int drho = 4;
incardon's avatar
incardon committed
147
148

// calculated force
incardon's avatar
incardon committed
149
const int force = 5;
incardon's avatar
incardon committed
150
151

// velocity
incardon's avatar
incardon committed
152
const int velocity = 6;
incardon's avatar
incardon committed
153
154

// velocity at previous step
incardon's avatar
incardon committed
155
156
const int velocity_prev = 7;

incardon's avatar
incardon committed
157
158
159
160
161
162
// 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
incardon's avatar
incardon committed
163

incardon's avatar
incardon committed
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200

/*! \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
incardon's avatar
incardon committed
201
{
incardon's avatar
incardon committed
202
203
204
205
206
207
208
	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);
	}
incardon's avatar
incardon committed
209

incardon's avatar
incardon committed
210
211
	template<typename Decomposition> inline void applyModel(Decomposition & dec, size_t v)
	{
incardon's avatar
incardon committed
212

incardon's avatar
incardon committed
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
	}
};

/*!
 * \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)
{
incardon's avatar
incardon committed
233
234
235
236
237
238
239
240
241
242
243
244
245
	auto it = vd.getDomainIterator();

	while (it.isNext())
	{
		auto a = it.get();

		double rho_a = vd.template getProp<rho>(a);
		double rho_frac = rho_a / rho_zero;

		vd.template getProp<Pressure>(a) = B*( rho_frac*rho_frac*rho_frac*rho_frac*rho_frac*rho_frac*rho_frac - 1.0);

		++it;
	}
incardon's avatar
incardon committed
246
}
incardon's avatar
incardon committed
247

incardon's avatar
incardon committed
248
/*! \cond [eq_state_and_ker] \endcond */
incardon's avatar
incardon committed
249

incardon's avatar
incardon committed
250
251
252
253
254
255
256
257
258
259
260
/*!
 * \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
 *
 */
incardon's avatar
incardon committed
261

incardon's avatar
incardon committed
262
/*! \cond [kernel_sph] \endcond */
incardon's avatar
incardon committed
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277

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;
}

incardon's avatar
incardon committed
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
/*! \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;
incardon's avatar
incardon committed
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331

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;
	}
}

incardon's avatar
incardon committed
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
/*! \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 */

incardon's avatar
incardon committed
350
// Tensile correction
incardon's avatar
incardon committed
351
inline double Tensile(double r, double rhoa, double rhob, double prs1, double prs2)
incardon's avatar
incardon committed
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
{
	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));
}

incardon's avatar
incardon committed
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
/*! \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)
incardon's avatar
incardon committed
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
{
	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;
}

incardon's avatar
incardon committed
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
/*! \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 */

incardon's avatar
incardon committed
431
432
433
434
435
template<typename CellList> inline double calc_forces(particles & vd, CellList & NN, double & max_visc)
{
	auto part = vd.getDomainIterator();
	double visc = 0;

incardon's avatar
incardon committed
436
	// Update the cell-list
incardon's avatar
incardon committed
437
438
	vd.updateCellList(NN);

incardon's avatar
incardon committed
439
	// For each particle ...
incardon's avatar
incardon committed
440
441
	while (part.isNext())
	{
incardon's avatar
incardon committed
442
		// ... a
incardon's avatar
incardon committed
443
444
445
446
447
		auto a = part.get();

		// Get the position xp of the particle
		Point<3,double> xa = vd.getPos(a);

incardon's avatar
incardon committed
448
		// Take the mass of the particle dependently if it is FLUID or BOUNDARY
incardon's avatar
incardon committed
449
		double massa = (vd.getProp<type>(a) == FLUID)?MassFluid:MassBound;
incardon's avatar
incardon committed
450
451

		// Get the density of the of the particle a
incardon's avatar
incardon committed
452
		double rhoa = vd.getProp<rho>(a);
incardon's avatar
incardon committed
453
454

		// Get the pressure of the particle a
incardon's avatar
incardon committed
455
		double Pa = vd.getProp<Pressure>(a);
incardon's avatar
incardon committed
456
457

		// Get the Velocity of the particle a
incardon's avatar
incardon committed
458
459
		Point<3,double> va = vd.getProp<velocity>(a);

incardon's avatar
incardon committed
460
		// Reset the force counter (- gravity on zeta direction)
incardon's avatar
incardon committed
461
462
463
464
465
		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;

incardon's avatar
incardon committed
466
467
468
469
470
471
		// 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)));
incardon's avatar
incardon committed
472

incardon's avatar
incardon committed
473
474
475
476
477
			// For each neighborhood particle
			while (Np.isNext() == true)
			{
				// ... q
				auto b = Np.get();
incardon's avatar
incardon committed
478

incardon's avatar
incardon committed
479
480
				// Get the position xp of the particle
				Point<3,double> xb = vd.getPos(b);
incardon's avatar
incardon committed
481

incardon's avatar
incardon committed
482
483
				// if (p == q) skip this particle
				if (a.getKey() == b)	{++Np; continue;};
incardon's avatar
incardon committed
484

incardon's avatar
incardon committed
485
486
				// get the mass of the particle
				double massb = (vd.getProp<type>(b) == FLUID)?MassFluid:MassBound;
incardon's avatar
incardon committed
487

incardon's avatar
incardon committed
488
489
				// Get the velocity of the particle b
				Point<3,double> vb = vd.getProp<velocity>(b);
incardon's avatar
incardon committed
490

incardon's avatar
incardon committed
491
492
493
				// Get the pressure and density of particle b
				double Pb = vd.getProp<Pressure>(b);
				double rhob = vd.getProp<rho>(b);
incardon's avatar
incardon committed
494

incardon's avatar
incardon committed
495
496
497
498
				// Get the distance between p and q
				Point<3,double> dr = xa - xb;
				// take the norm of this vector
				double r2 = norm2(dr);
incardon's avatar
incardon committed
499

incardon's avatar
incardon committed
500
501
502
503
504
				// If the particles interact ...
				if (r2 < 4.0*H*H)
				{
					// ... calculate delta rho
					double r = sqrt(r2);
incardon's avatar
incardon committed
505

incardon's avatar
incardon committed
506
					Point<3,double> dv = va - vb;
incardon's avatar
incardon committed
507

incardon's avatar
incardon committed
508
509
					Point<3,double> DW;
					DWab(dr,DW,r,false);
incardon's avatar
incardon committed
510

incardon's avatar
incardon committed
511
512
513
					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);
incardon's avatar
incardon committed
514

incardon's avatar
incardon committed
515
516
					vd.getProp<drho>(a) += massb*(dv.get(0)*DW.get(0)+dv.get(1)*DW.get(1)+dv.get(2)*DW.get(2));
				}
incardon's avatar
incardon committed
517

incardon's avatar
incardon committed
518
				++Np;
incardon's avatar
incardon committed
519
520
			}
		}
incardon's avatar
incardon committed
521
		else
incardon's avatar
incardon committed
522
		{
incardon's avatar
incardon committed
523
			// If it is a fluid particle calculate based on equation 1 and 2
incardon's avatar
incardon committed
524

incardon's avatar
incardon committed
525
526
			// Get an iterator over the neighborhood particles of p
			auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(a)));
incardon's avatar
incardon committed
527

incardon's avatar
incardon committed
528
529
530
531
532
			// For each neighborhood particle
			while (Np.isNext() == true)
			{
				// ... q
				auto b = Np.get();
incardon's avatar
incardon committed
533

incardon's avatar
incardon committed
534
535
				// Get the position xp of the particle
				Point<3,double> xb = vd.getPos(b);
incardon's avatar
incardon committed
536

incardon's avatar
incardon committed
537
538
				// if (p == q) skip this particle
				if (a.getKey() == b)	{++Np; continue;};
incardon's avatar
incardon committed
539

incardon's avatar
incardon committed
540
541
542
543
				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);
incardon's avatar
incardon committed
544

incardon's avatar
incardon committed
545
546
547
548
				// Get the distance between p and q
				Point<3,double> dr = xa - xb;
				// take the norm of this vector
				double r2 = norm2(dr);
incardon's avatar
incardon committed
549

incardon's avatar
incardon committed
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
				// 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);

					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));

					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);

					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));
				}

				++Np;
			}
incardon's avatar
incardon committed
571
		}
incardon's avatar
incardon committed
572
573

		++part;
incardon's avatar
incardon committed
574
	}
incardon's avatar
incardon committed
575
}
incardon's avatar
incardon committed
576

incardon's avatar
incardon committed
577
/*! \cond [calc_forces] \endcond */
incardon's avatar
incardon committed
578

incardon's avatar
incardon committed
579
580
581
582
583
584
585
586
587
588
/*!
 *
 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
 *
 * This function calculate the Maximum acceleration and velocity across the particles.
 *
 * \snippet Vector/7_SPH_dlb/main.cpp max_acc_vel
 *
 *
 */
incardon's avatar
incardon committed
589

incardon's avatar
incardon committed
590
/*! \cond [max_acc_vel] \endcond */
incardon's avatar
incardon committed
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618

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<force>(a));
		double acc2 = norm2(acc);

		Point<3,double> vel(vd.getProp<velocity>(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);
}

incardon's avatar
incardon committed
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
/*! \cond [max_acc_vel] \endcond */

/*!
 *
 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
 *
 * In this example we are using Dynamic time-stepping. The Dynamic time stepping is
 * calculated with the Courant-Friedrich-Lewy condition. See Monaghan 1992 "Smoothed Particle Hydrodynamic"
 *
 * \f$ \delta t = CFL \cdot min(t_f,t_{cv}) \f$
 *
 * where
 *
 * \f$ \delta t_f = min \sqrt{h/f_a}\f$
 *
 * \f$  \delta t_{cv} = min \frac{h}{c_s + max \left| \frac{hv_{ab} \cdot r_{ab}}{r_{ab}^2} \right|} \f$
 *
 *
 * \snippet Vector/7_SPH_dlb/main.cpp dyn_stepping
 *
 *
 */

/*! \cond [dyn_stepping] \endcond */
incardon's avatar
incardon committed
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663

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<int>::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<double(DtMin))
		dt=double(DtMin);

	return dt;
}

incardon's avatar
incardon committed
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
/*! \cond [dyn_stepping] \endcond */

/*!
 *
 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
 *
 * This function perform verlet integration accordingly to the Verlet time stepping scheme
 *
 * \f$ v_a^{n+1} = v_a^{n-1} + 2 \delta t F_a^{n} \f$
 *
 * \f$ r_a^{n+1} = \delta t V_a^n + 0.5 \delta t^2 F_a^n \f$
 *
 * \f$ \rho_a^{n+1} = \rho_a^{n-1} + 2 \delta t D_a^n \f$
 *
 * Every N Verlet steps the euler stepping scheme is choosen to avoid instabilities
 *
 * \f$ v_a^{n+1} = v_a^{n} + \delta t F_a^n \f$
 *
 * \f$ r_a^{n+1} = r_a^{n} + \delta t V_a^n + 0.5 delta t^2 F_a^n \f$
 *
 * \f$ \rho_a^n + \delta t D_a^n \f$
 *
 * More the integration this function also check that no particles go outside the simulation
 * domain or their density go dangerously out of range
 *
 * \snippet Vector/7_SPH_dlb/main.cpp verlet_int
 *
 *
 */

/*! \cond [verlet_int] \endcond */

incardon's avatar
incardon committed
696
697
openfpm::vector<size_t> to_remove;

incardon's avatar
incardon committed
698
size_t cnt = 0;
incardon's avatar
incardon committed
699

incardon's avatar
incardon committed
700
void verlet_int(particles & vd, double dt, bool VerletStep)
incardon's avatar
incardon committed
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
{
	to_remove.clear();

	// Calculate the maximum acceleration
	auto part = vd.getDomainIterator();

	double dt205 = dt*dt*0.5;
	double dt2 = dt*2.0;

	while (part.isNext())
	{
		auto a = part.get();

		if (vd.template getProp<type>(a) == BOUNDARY)
		{
incardon's avatar
incardon committed
716
717
			double rhop = vd.template getProp<rho>(a);

incardon's avatar
incardon committed
718
			// Update only the density
incardon's avatar
incardon committed
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
		    if (VerletStep == true)
		    {
		    	vd.template getProp<velocity>(a)[0] = 0.0;
		    	vd.template getProp<velocity>(a)[1] = 0.0;
		    	vd.template getProp<velocity>(a)[2] = 0.0;
		    	vd.template getProp<rho>(a) = vd.template getProp<rho_prev>(a) + dt2*vd.template getProp<drho>(a);
		    }
		    else
		    {
		    	vd.template getProp<velocity>(a)[0] = 0.0;
		    	vd.template getProp<velocity>(a)[1] = 0.0;
		    	vd.template getProp<velocity>(a)[2] = 0.0;
		    	vd.template getProp<rho>(a) = vd.template getProp<rho>(a) + dt*vd.template getProp<drho>(a);
		    }

		    vd.template getProp<rho_prev>(a) = rhop;

incardon's avatar
incardon committed
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
			++part;
			continue;
		}

		//-Calculate displacement and update position / Calcula desplazamiento y actualiza posicion.
		double dx = vd.template getProp<velocity>(a)[0]*dt + vd.template getProp<force>(a)[0]*dt205;
	    double dy = vd.template getProp<velocity>(a)[1]*dt + vd.template getProp<force>(a)[1]*dt205;
	    double dz = vd.template getProp<velocity>(a)[2]*dt + vd.template getProp<force>(a)[2]*dt205;

	    vd.getPos(a)[0] += dx;
	    vd.getPos(a)[1] += dy;
	    vd.getPos(a)[2] += dz;

	    double velX = vd.template getProp<velocity>(a)[0];
	    double velY = vd.template getProp<velocity>(a)[1];
	    double velZ = vd.template getProp<velocity>(a)[2];
	    double rhop = vd.template getProp<rho>(a);

incardon's avatar
incardon committed
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
	    if (VerletStep == true)
	    {
	    	vd.template getProp<velocity>(a)[0] = vd.template getProp<velocity_prev>(a)[0] + vd.template getProp<force>(a)[0]*dt2;
	    	vd.template getProp<velocity>(a)[1] = vd.template getProp<velocity_prev>(a)[1] + vd.template getProp<force>(a)[1]*dt2;
	    	vd.template getProp<velocity>(a)[2] = vd.template getProp<velocity_prev>(a)[2] + vd.template getProp<force>(a)[2]*dt2;
	    	vd.template getProp<rho>(a) = vd.template getProp<rho_prev>(a) + dt2*vd.template getProp<drho>(a);
	    }
	    else
	    {
	    	vd.template getProp<velocity>(a)[0] = vd.template getProp<velocity>(a)[0] + vd.template getProp<force>(a)[0]*dt;
	    	vd.template getProp<velocity>(a)[1] = vd.template getProp<velocity>(a)[1] + vd.template getProp<force>(a)[1]*dt;
	    	vd.template getProp<velocity>(a)[2] = vd.template getProp<velocity>(a)[2] + vd.template getProp<force>(a)[2]*dt;
	    	vd.template getProp<rho>(a) = vd.template getProp<rho>(a) + dt*vd.template getProp<drho>(a);
	    }

	    // Check if there are particles to remove

	    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<rho>(a) < RhoMin || vd.template getProp<rho>(a) > RhoMax)
	    {
	    	std::cout << "Particle_out" << std::endl;
	                   to_remove.add(a.getKey());
	    }
incardon's avatar
incardon committed
778
779
780
781
782
783
784
785
786
787

	    vd.template getProp<velocity_prev>(a)[0] = velX;
	    vd.template getProp<velocity_prev>(a)[1] = velY;
	    vd.template getProp<velocity_prev>(a)[2] = velZ;
	    vd.template getProp<rho_prev>(a) = rhop;

		++part;
	}

	vd.remove(to_remove,0);
incardon's avatar
incardon committed
788
789

	cnt++;
incardon's avatar
incardon committed
790
791
}

incardon's avatar
incardon committed
792
793
/*! \cond [verlet_int] \endcond */

incardon's avatar
incardon committed
794
795
796
797
int main(int argc, char* argv[])
{
	/*!
	 *
incardon's avatar
incardon committed
798
799
800
	 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
	 *
	 * ## Main function ##
incardon's avatar
incardon committed
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
	 *
	 * Here we Initialize the library, we create a Box that define our domain, boundary conditions, ghost
	 *
	 * \see \ref e0_s_init
	 *
	 * \snippet Vector/7_SPH_dlb/main.cpp Initialization and parameters
	 *
	 */

	//! \cond [Initialization and parameters] \endcond

    // 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
incardon's avatar
incardon committed
816
817
	Box<3,double> domain({-0.05,-0.05,-0.05},{1.7010,0.7065,0.5025});
	size_t sz[3] = {207,90,66};
incardon's avatar
incardon committed
818
819
820
821
822
823
824
825
826
827
828
829
830

	// 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};

	// extended boundary around the domain, and the processor domain
	Ghost<3,double> g(2*H);
	
	//! \cond [Initialization and parameters] \endcond

	/*!
incardon's avatar
incardon committed
831
	 * \page Vector_7_SPH_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
incardon's avatar
incardon committed
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
	 *
	 * ## %Vector create ##
	 *
	 * Here we define a distributed vector in 3D, containing 3 properties, a
	 * scalar double, a vector double[3], and a tensor or rank 2 double[3][3].
	 * In this case the vector contain 0 particles initially
	 *
	 * \see \ref vector_inst
	 *
	 * \snippet Vector/1_celllist/main.cpp vector inst
	 *
	 */

	//! \cond [vector inst] \endcond

incardon's avatar
incardon committed
847
	particles vd(0,domain,bc,g,DEC_GRAN(4096));
incardon's avatar
incardon committed
848
849
850
851
852
853

	//! \cond [vector inst] \endcond

	// the scalar is the element at position 0 in the aggregate
	const int type = 0;

incardon's avatar
incardon committed
854
	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});
incardon's avatar
incardon committed
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958

	// first we create Fluid particles
	// Fluid particles are created

	auto fluid_it = DrawParticles::DrawBox(vd,sz,domain,fluid_box);
	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);

	while (fluid_it.isNext())
	{
		vd.add();

		vd.getLastPos()[0] = fluid_it.get().get(0);
		vd.getLastPos()[1] = fluid_it.get().get(1);
		vd.getLastPos()[2] = fluid_it.get().get(2);

		vd.template getLastProp<type>() = 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<Pressure>() = rho_zero * gravity *  (max_fluid_height - fluid_it.get().get(2));

		vd.template getLastProp<rho>() = pow(vd.template getLastProp<Pressure>() / B + 1, 1.0/gamma_) * rho_zero;
		vd.template getLastProp<rho_prev>() = vd.template getLastProp<rho>();
		vd.template getLastProp<velocity>()[0] = 0.0;
		vd.template getLastProp<velocity>()[1] = 0.0;
		vd.template getLastProp<velocity>()[2] = 0.0;

		vd.template getLastProp<velocity_prev>()[0] = 0.0;
		vd.template getLastProp<velocity_prev>()[1] = 0.0;
		vd.template getLastProp<velocity_prev>()[2] = 0.0;

		++fluid_it;
	}

	// 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<Box<3,double>> 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<type>() = BOUNDARY;
		vd.template getLastProp<rho>() = rho_zero;
		vd.template getLastProp<rho_prev>() = rho_zero;
		vd.template getLastProp<velocity>()[0] = 0.0;
		vd.template getLastProp<velocity>()[1] = 0.0;
		vd.template getLastProp<velocity>()[2] = 0.0;

		vd.template getLastProp<velocity_prev>()[0] = 0.0;
		vd.template getLastProp<velocity_prev>()[1] = 0.0;
		vd.template getLastProp<velocity_prev>()[2] = 0.0;

		++bound_box;
	}

	// Obstacle

	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<type>() = BOUNDARY;
		vd.template getLastProp<rho>() = rho_zero;
		vd.template getLastProp<rho_prev>() = rho_zero;
		vd.template getLastProp<velocity>()[0] = 0.0;
		vd.template getLastProp<velocity>()[1] = 0.0;
		vd.template getLastProp<velocity>()[2] = 0.0;

		vd.template getLastProp<velocity_prev>()[0] = 0.0;
		vd.template getLastProp<velocity_prev>()[1] = 0.0;
		vd.template getLastProp<velocity_prev>()[2] = 0.0;

		++obstacle_box;
	}

	vd.map();
incardon's avatar
incardon committed
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
	vd.getDecomposition().write("Decomposition_before_load_bal");

	// Now that we fill the vector with particles
	ModelCustom md;

	vd.addComputationCosts(md);
	vd.getDecomposition().getDistribution().write("BEFORE_DECOMPOSE");
	vd.getDecomposition().decompose();
	vd.map();

	vd.addComputationCosts(md);
	vd.getDecomposition().getDistribution().write("AFTER_DECOMPOSE1");

	vd.getDecomposition().rebalance(1);

	vd.map();
	vd.getDecomposition().getDistribution().write("AFTER_DECOMPOSE2");

	std::cout << "N particles: " << vd.size_local()  << "    " << create_vcluster().getProcessUnitID() << "      " << "Get processor Load " << vd.getDecomposition().getDistribution().getProcessorLoad() << std::endl;

	vd.write("Geometry");
	vd.getDecomposition().write("Decomposition_after_load_bal");
	vd.getDecomposition().getDistribution().write("Distribution_load_bal");

	vd.ghost_get<type,rho,Pressure,velocity>();
incardon's avatar
incardon committed
984
985
986
987
988
989
990
991

	auto NN = vd.getCellList(2*H);

	// Evolve


	size_t write = 0;
	size_t it = 0;
incardon's avatar
incardon committed
992
	size_t it_reb = 0;
incardon's avatar
incardon committed
993
994
995
	double t = 0.0;
	while (t <= t_end)
	{
incardon's avatar
incardon committed
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
		timer it_time;

		////// Do rebalancing every 200 timesteps
		it_reb++;
		if (it_reb == 10)
		{
			vd.map();

			it_reb = 0;
			ModelCustom md;
			vd.addComputationCosts(md);
			vd.getDecomposition().rebalance(1);

			std::cout << "REBALANCED " << std::endl;
		}
incardon's avatar
incardon committed
1011
1012

		vd.map();
incardon's avatar
incardon committed
1013
		vd.ghost_get<type,rho,Pressure,velocity>();
incardon's avatar
incardon committed
1014
1015
1016
1017
1018
1019

		// Calculate pressure from the density
		EqState(vd);

		double max_visc = 0.0;

incardon's avatar
incardon committed
1020
1021
		it_time.start();

incardon's avatar
incardon committed
1022
1023
		// Calc forces
		calc_forces(vd,NN,max_visc);
incardon's avatar
incardon committed
1024
1025
1026
1027
1028
1029
		it_time.stop();

		// Get the maximum viscosity term across processors
		Vcluster & v_cl = create_vcluster();
		v_cl.max(max_visc);
		v_cl.execute();
incardon's avatar
incardon committed
1030
1031
1032
1033
1034
1035
1036

		// Calculate delta t integration
		double dt = calc_deltaT(vd,max_visc);

//		std::cout << "Calculate deltaT: " << dt << "   " << DtMin << std::endl;

		// VerletStep
incardon's avatar
incardon committed
1037
1038
1039
1040
1041
1042
1043
1044
		it++;
		if (it < 40)
			verlet_int(vd,dt,true);
		else
		{
			verlet_int(vd,dt,false);
			it = 0;
		}
incardon's avatar
incardon committed
1045
1046
1047
1048
1049
1050
1051
1052

		t += dt;

		if (write < t*100)
		{

			vd.write("Geometry",write);
			write++;
incardon's avatar
incardon committed
1053
1054
1055
1056
1057
1058

			std::cout << "TIME: " << t << "  write " << it_time.getwct() << "   " << v_cl.getProcessUnitID() << "   " << cnt << std::endl;
		}
		else
		{
			std::cout << "TIME: " << t << "  " << it_time.getwct() << "   " << v_cl.getProcessUnitID() << "   " << cnt << std::endl;
incardon's avatar
incardon committed
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
		}
	}

	//! \cond [finalize] \endcond

	openfpm_finalize();

	//! \cond [finalize] \endcond

	/*!
incardon's avatar
incardon committed
1069
	 * \page Vector_7_SPH_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
incardon's avatar
incardon committed
1070
1071
1072
1073
1074
1075
1076
1077
	 *
	 * ## Full code ## {#code_e0_sim}
	 *
	 * \include Vector/0_simple/main.cpp
	 *
	 */
}