main.cpp 40.7 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
 *
incardon's avatar
incardon committed
15
 * \htmlonly
incardon's avatar
incardon committed
16
 * <a href="#" onclick="hide_show('vector-video-3')" >Simulation video 1</a><br>
incardon's avatar
incardon committed
17
 * <div style="display:none" id="vector-video-3">
incardon's avatar
incardon committed
18
 * <video id="vid3" width="1200" height="576" controls> <source src="http://openfpm.mpi-cbg.de/web/images/examples/7_SPH_dlb/sph_speed.mp4" type="video/mp4"></video>
incardon's avatar
incardon committed
19
 * </div>
incardon's avatar
incardon committed
20
 * <a href="#" onclick="hide_show('vector-video-4')" >Simulation video 2</a><br>
incardon's avatar
incardon committed
21
 * <div style="display:none" id="vector-video-4">
incardon's avatar
incardon committed
22
23
24
25
26
27
28
29
30
 * <video id="vid4" width="1200" height="576" controls> <source src="http://openfpm.mpi-cbg.de/web/images/examples/7_SPH_dlb/sph_speed2.mp4" type="video/mp4"></video>
 * </div>
 * <a href="#" onclick="hide_show('vector-video-15')" >Simulation dynamic load balancing video 1</a><br>
 * <div style="display:none" id="vector-video-15">
 * <video id="vid15" width="1200" height="576" controls> <source src="http://openfpm.mpi-cbg.de/web/images/examples/7_SPH_dlb/sph_dlb.mp4" type="video/mp4"></video>
 * </div>
 * <a href="#" onclick="hide_show('vector-video-16')" >Simulation dynamic load balancing video 2</a><br>
 * <div style="display:none" id="vector-video-16">
 * <video id="vid16" width="1200" height="576" controls> <source src="http://openfpm.mpi-cbg.de/web/images/examples/7_SPH_dlb/sph_dlb2.mp4" type="video/mp4"></video>
incardon's avatar
incardon committed
31
32
33
34
 * </div>
 * \endhtmlonly
 *
 * \htmlonly
incardon's avatar
incardon committed
35
36
37
38
 * <img src="http://ppmcore.mpi-cbg.de/web/images/examples/7_SPH_dlb/dam_break_all.jpg"/>
 * \endhtmlonly
 *
 * ## Inclusion ## {#e7_sph_inclusion}
incardon's avatar
incardon committed
39
40
 *
 * In order to use distributed vectors in our code we have to include the file Vector/vector_dist.hpp
incardon's avatar
incardon committed
41
42
 * we also include DrawParticles that has nice utilities to draw particles in parallel accordingly
 * to simple shapes
incardon's avatar
incardon committed
43
44
45
46
47
 *
 * \snippet Vector/7_SPH_dlb/main.cpp inclusion
 *
 */

incardon's avatar
incardon committed
48
49
50
//#define SE_CLASS1
//#define STOP_ON_ERROR

incardon's avatar
incardon committed
51
52
53
//! \cond [inclusion] \endcond
#include "Vector/vector_dist.hpp"
#include <math.h>
incardon's avatar
incardon committed
54
#include "Draw/DrawParticles.hpp"
incardon's avatar
incardon committed
55
56
//! \cond [inclusion] \endcond

incardon's avatar
incardon committed
57
58
59
/*!
 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balacing
 *
incardon's avatar
incardon committed
60
 * ## SPH simulation {#e7_sph_parameters}
incardon's avatar
incardon committed
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
 *
 * 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
incardon's avatar
incardon committed
85
86
 * can be found regarding SPH in general in the original Monghagan SPH paper.
 * In this example we use the sligtly modified version
incardon's avatar
incardon committed
87
88
 * 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.
incardon's avatar
incardon committed
89
90
91
92
93
 *
 * ### Parameters {#e7_sph_parameters}
 *
 * Based on the equation
 * reported before several constants must be defined.
incardon's avatar
incardon committed
94
95
96
97
98
99
 *
 * \snippet Vector/7_SPH_dlb/main.cpp sim parameters
 *
 */

/*! \cond [sim parameters] \endcond */
incardon's avatar
incardon committed
100

incardon's avatar
incardon committed
101
// A constant to indicate boundary particles
incardon's avatar
incardon committed
102
103
#define BOUNDARY 0

incardon's avatar
incardon committed
104
105
// A constant to indicate fluid particles
#define FLUID 1
incardon's avatar
incardon committed
106

incardon's avatar
incardon committed
107
// initial spacing between particles dp in the formulas
incardon's avatar
incardon committed
108
const double dp = 0.0085;
incardon's avatar
incardon committed
109
// Maximum height of the fluid water
incardon's avatar
incardon committed
110
// is going to be calculated and filled later on
incardon's avatar
incardon committed
111
double h_swl = 0.0;
incardon's avatar
incardon committed
112

incardon's avatar
incardon committed
113
// c_s in the formulas (constant used to calculate the sound speed)
incardon's avatar
incardon committed
114
const double coeff_sound = 20.0;
incardon's avatar
incardon committed
115
116

// gamma in the formulas
incardon's avatar
incardon committed
117
const double gamma_ = 7.0;
incardon's avatar
incardon committed
118
119

// sqrt(3.0*dp*dp) support of the kernel
incardon's avatar
incardon committed
120
const double H = 0.0147224318643;
incardon's avatar
incardon committed
121
122

// Eta in the formulas
incardon's avatar
incardon committed
123
const double Eta2 = 0.01 * H*H;
incardon's avatar
incardon committed
124

incardon's avatar
incardon committed
125
// alpha in the formula
incardon's avatar
incardon committed
126
const double visco = 0.1;
incardon's avatar
incardon committed
127
128

// cbar in the formula (calculated later)
incardon's avatar
incardon committed
129
double cbar = 0.0;
incardon's avatar
incardon committed
130
131

// Mass of the fluid particles
incardon's avatar
incardon committed
132
const double MassFluid = 0.000614125;
incardon's avatar
incardon committed
133
134

// Mass of the boundary particles
incardon's avatar
incardon committed
135
const double MassBound = 0.000614125;
incardon's avatar
incardon committed
136
137
138
139
140

// End simulation time
const double t_end = 1.5;

// Gravity acceleration
incardon's avatar
incardon committed
141
const double gravity = 9.81;
incardon's avatar
incardon committed
142
143

// Reference densitu 1000Kg/m^3
incardon's avatar
incardon committed
144
const double rho_zero = 1000.0;
incardon's avatar
incardon committed
145
146

// Filled later require h_swl, it is b in the formulas
incardon's avatar
incardon committed
147
double B = 0.0;
incardon's avatar
incardon committed
148
149

// Constant used to define time integration
incardon's avatar
incardon committed
150
const double CFLnumber = 0.2;
incardon's avatar
incardon committed
151
152

// Minimum T
incardon's avatar
incardon committed
153
154
const double DtMin = 0.00001;

incardon's avatar
incardon committed
155
156
157
158
159
160
// Minimum Rho allowed
const double RhoMin = 700.0;

// Maximum Rho allowed
const double RhoMax = 1300.0;

incardon's avatar
incardon committed
161
162
163
// Filled in initialization
double max_fluid_height = 0.0;

incardon's avatar
incardon committed
164
165
166
// Properties

// FLUID or BOUNDARY
incardon's avatar
incardon committed
167
const size_t type = 0;
incardon's avatar
incardon committed
168
169

// Density
incardon's avatar
incardon committed
170
const int rho = 1;
incardon's avatar
incardon committed
171
172

// Density at step n-1
incardon's avatar
incardon committed
173
const int rho_prev = 2;
incardon's avatar
incardon committed
174
175

// Pressure
incardon's avatar
incardon committed
176
const int Pressure = 3;
incardon's avatar
incardon committed
177
178

// Delta rho calculated in the force calculation
incardon's avatar
incardon committed
179
const int drho = 4;
incardon's avatar
incardon committed
180
181

// calculated force
incardon's avatar
incardon committed
182
const int force = 5;
incardon's avatar
incardon committed
183
184

// velocity
incardon's avatar
incardon committed
185
const int velocity = 6;
incardon's avatar
incardon committed
186
187

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

incardon's avatar
incardon committed
190
191
/*! \cond [sim parameters] \endcond */

incardon's avatar
incardon committed
192
193
/*! \cond [vector_dist_def] \endcond */

incardon's avatar
incardon committed
194
195
196
197
198
199
// 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
200

incardon's avatar
incardon committed
201
/*! \cond [vector_dist_def] \endcond */
incardon's avatar
incardon committed
202

incardon's avatar
incardon committed
203
/*! \cond [model custom] \endcond */
incardon's avatar
incardon committed
204
205
206

struct ModelCustom
{
incardon's avatar
incardon committed
207
	template<typename Decomposition, typename vector> inline void addComputation(Decomposition & dec,
incardon's avatar
incardon committed
208
			                                                                     vector & vd,
incardon's avatar
incardon committed
209
210
																				 size_t v,
																				 size_t p)
incardon's avatar
incardon committed
211
212
	{
		if (vd.template getProp<type>(p) == FLUID)
incardon's avatar
incardon committed
213
			dec.addComputationCost(v,4);
incardon's avatar
incardon committed
214
		else
incardon's avatar
incardon committed
215
			dec.addComputationCost(v,3);
incardon's avatar
incardon committed
216
217
218
219
220
221
	}

	template<typename Decomposition> inline void applyModel(Decomposition & dec, size_t v)
	{
		dec.setSubSubDomainComputationCost(v, dec.getSubSubDomainComputationCost(v) * dec.getSubSubDomainComputationCost(v));
	}
incardon's avatar
incardon committed
222
223
224
225
226

	double distributionTol()
	{
		return 1.01;
	}
incardon's avatar
incardon committed
227
228
};

incardon's avatar
incardon committed
229
/*! \cond [model custom] \endcond */
incardon's avatar
incardon committed
230
231
232
233

/*!
 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balacing
 *
incardon's avatar
incardon committed
234
 * ### Equation of state {#e7_sph_equation_state}
incardon's avatar
incardon committed
235
236
237
238
239
240
241
242
243
244
245
246
247
 *
 * 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
248
249
250
251
252
253
254
255
256
257
258
259
260
	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
261
}
incardon's avatar
incardon committed
262

incardon's avatar
incardon committed
263
/*! \cond [eq_state_and_ker] \endcond */
incardon's avatar
incardon committed
264

incardon's avatar
incardon committed
265
266
267
/*!
 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
 *
incardon's avatar
incardon committed
268
269
270
 * ### Cubic SPH kernel and derivatives {#e7_sph_kernel}
 *
 * This function define the Cubic kernel or \f$ W_{ab} \f$. The cubic kernel is
incardon's avatar
incardon committed
271
272
273
274
275
276
277
 * 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
278

incardon's avatar
incardon committed
279
/*! \cond [kernel_sph] \endcond */
incardon's avatar
incardon committed
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294

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
295
296
297
298
299
/*! \cond [kernel_sph] \endcond */

/*!
 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
 *
incardon's avatar
incardon committed
300
 * This function define the gradient of the Cubic kernel function \f$ W_{ab} \f$.
incardon's avatar
incardon committed
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
 *
 * \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
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348

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
349
350
351
352
353
/*! \cond [kernel_sph_der] \endcond */

/*!
 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
 *
incardon's avatar
incardon committed
354
355
 * ### Tensile correction {#e7_sph_tensile}
 *
incardon's avatar
incardon committed
356
357
 * 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
incardon's avatar
incardon committed
358
 * to get too near. Can be considered at small scale like a repulsive force that avoid
incardon's avatar
incardon committed
359
360
361
362
363
364
365
366
367
368
 * 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
369
// Tensile correction
incardon's avatar
incardon committed
370
inline double Tensile(double r, double rhoa, double rhob, double prs1, double prs2)
incardon's avatar
incardon committed
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
{
	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
399
400
401
402
403
404
405
/*! \cond [tensile_term] \endcond */


/*!
 *
 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
 *
incardon's avatar
incardon committed
406
407
408
 * ### Viscous term {#e7_sph_viscous}
 *
 * This function implement the viscous term \f$ \Pi_{ab} \f$
incardon's avatar
incardon committed
409
410
411
412
413
414
415
416
417
 *
 * \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
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
{
	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
435
436
437
438
439
440
/*! \cond [viscous_term] \endcond */

/*!
 *
 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
 *
incardon's avatar
incardon committed
441
 * ### Force calculation {#e7_force_calc}
incardon's avatar
incardon committed
442
443
444
445
446
447
448
449
450
451
 *
 * 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
452
453
454
455
456
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
457
	// Update the cell-list
incardon's avatar
incardon committed
458
459
	vd.updateCellList(NN);

incardon's avatar
incardon committed
460
	// For each particle ...
incardon's avatar
incardon committed
461
462
	while (part.isNext())
	{
incardon's avatar
incardon committed
463
		// ... a
incardon's avatar
incardon committed
464
465
466
467
468
		auto a = part.get();

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

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

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

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

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

incardon's avatar
incardon committed
481
		// Reset the force counter (- gravity on zeta direction)
incardon's avatar
incardon committed
482
483
484
485
486
		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
487
488
489
490
491
492
		// 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
493

incardon's avatar
incardon committed
494
495
496
497
498
			// For each neighborhood particle
			while (Np.isNext() == true)
			{
				// ... q
				auto b = Np.get();
incardon's avatar
incardon committed
499

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

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

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

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

incardon's avatar
incardon committed
512
513
514
				// 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
515

incardon's avatar
incardon committed
516
517
518
519
				// 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
520

incardon's avatar
incardon committed
521
522
523
524
525
				// If the particles interact ...
				if (r2 < 4.0*H*H)
				{
					// ... calculate delta rho
					double r = sqrt(r2);
incardon's avatar
incardon committed
526

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

incardon's avatar
incardon committed
529
530
					Point<3,double> DW;
					DWab(dr,DW,r,false);
incardon's avatar
incardon committed
531

incardon's avatar
incardon committed
532
533
534
					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
535

incardon's avatar
incardon committed
536
537
					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
538

incardon's avatar
incardon committed
539
				++Np;
incardon's avatar
incardon committed
540
541
			}
		}
incardon's avatar
incardon committed
542
		else
incardon's avatar
incardon committed
543
		{
incardon's avatar
incardon committed
544
			// If it is a fluid particle calculate based on equation 1 and 2
incardon's avatar
incardon committed
545

incardon's avatar
incardon committed
546
547
			// 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
548

incardon's avatar
incardon committed
549
550
551
552
553
			// For each neighborhood particle
			while (Np.isNext() == true)
			{
				// ... q
				auto b = Np.get();
incardon's avatar
incardon committed
554

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

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

incardon's avatar
incardon committed
561
562
563
564
				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
565

incardon's avatar
incardon committed
566
567
568
569
				// 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
570

incardon's avatar
incardon committed
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
				// 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
592
		}
incardon's avatar
incardon committed
593
594

		++part;
incardon's avatar
incardon committed
595
	}
incardon's avatar
incardon committed
596
}
incardon's avatar
incardon committed
597

incardon's avatar
incardon committed
598
/*! \cond [calc_forces] \endcond */
incardon's avatar
incardon committed
599

incardon's avatar
incardon committed
600
601
602
603
/*!
 *
 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
 *
incardon's avatar
incardon committed
604
605
 * ### Integration and dynamic time integration {#e7_delta_time_t}
 *
incardon's avatar
incardon committed
606
 * This function calculate the Maximum acceleration and velocity across the particles.
incardon's avatar
incardon committed
607
 * It is used to calculate a dynamic time-stepping.
incardon's avatar
incardon committed
608
609
610
611
612
 *
 * \snippet Vector/7_SPH_dlb/main.cpp max_acc_vel
 *
 *
 */
incardon's avatar
incardon committed
613

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

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
641
642
643
644
645

	Vcluster & v_cl = create_vcluster();
	v_cl.max(max_acc);
	v_cl.max(max_vel);
	v_cl.execute();
incardon's avatar
incardon committed
646
647
}

incardon's avatar
incardon committed
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
/*! \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
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692

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
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
/*! \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$
 *
incardon's avatar
incardon committed
711
 * \f$ r_a^{n+1} = r_a^{n} + \delta t V_a^n + \frac{1}{2} \delta t^2 F_a^n \f$
incardon's avatar
incardon committed
712
 *
incardon's avatar
incardon committed
713
 * \f$ \rho_a^{n+1} = \rho_a^n + \delta t D_a^n \f$
incardon's avatar
incardon committed
714
 *
incardon's avatar
incardon committed
715
 * This function also check that no particles go outside the simulation
incardon's avatar
incardon committed
716
717
 * domain or their density go dangerously out of range. If a particle go out of range is removed
 * from the simulation
incardon's avatar
incardon committed
718
719
720
721
722
723
724
725
 *
 * \snippet Vector/7_SPH_dlb/main.cpp verlet_int
 *
 *
 */

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

incardon's avatar
incardon committed
726
727
openfpm::vector<size_t> to_remove;

incardon's avatar
incardon committed
728
size_t cnt = 0;
incardon's avatar
incardon committed
729

incardon's avatar
incardon committed
730
void verlet_int(particles & vd, double dt)
incardon's avatar
incardon committed
731
{
incardon's avatar
incardon committed
732
	// list of the particle to remove
incardon's avatar
incardon committed
733
734
	to_remove.clear();

incardon's avatar
incardon committed
735
	// particle iterator
incardon's avatar
incardon committed
736
737
738
739
740
	auto part = vd.getDomainIterator();

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

incardon's avatar
incardon committed
741
	// For each particle ...
incardon's avatar
incardon committed
742
743
	while (part.isNext())
	{
incardon's avatar
incardon committed
744
		// ... a
incardon's avatar
incardon committed
745
746
		auto a = part.get();

incardon's avatar
incardon committed
747
		// if the particle is boundary
incardon's avatar
incardon committed
748
749
		if (vd.template getProp<type>(a) == BOUNDARY)
		{
incardon's avatar
incardon committed
750
			// Update rho
incardon's avatar
incardon committed
751
752
			double rhop = vd.template getProp<rho>(a);

incardon's avatar
incardon committed
753
			// Update only the density
incardon's avatar
incardon committed
754
755
756
757
	    	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);
incardon's avatar
incardon committed
758
759
760

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

incardon's avatar
incardon committed
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
			++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
779
780
781
782
783
784
785
786
787
    	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);

	    // 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<rho>(a) < RhoMin || vd.template getProp<rho>(a) > RhoMax)
incardon's avatar
incardon committed
788
	    {
incardon's avatar
incardon committed
789
	                   to_remove.add(a.getKey());
incardon's avatar
incardon committed
790
	    }
incardon's avatar
incardon committed
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833

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

	// remove the particles
	vd.remove(to_remove,0);

	// increment the iteration counter
	cnt++;
}

void euler_int(particles & vd, double dt)
{
	// 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;

	// For each particle ...
	while (part.isNext())
	{
		// ... a
		auto a = part.get();

		// if the particle is boundary
		if (vd.template getProp<type>(a) == BOUNDARY)
		{
			// Update rho
			double rhop = vd.template getProp<rho>(a);

			// Update only the density
	    	vd.template getProp<velocity>(a)[0] = 0.0;
	    	vd.template getProp<velocity>(a)[1] = 0.0;
	    	vd.template getProp<velocity>(a)[2] = 0.0;
incardon's avatar
incardon committed
834
835
	    	vd.template getProp<rho>(a) = vd.template getProp<rho>(a) + dt*vd.template getProp<drho>(a);

incardon's avatar
incardon committed
836
837
838
839
840
		    vd.template getProp<rho_prev>(a) = rhop;

			++part;
			continue;
		}
incardon's avatar
incardon committed
841

incardon's avatar
incardon committed
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
		//-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);

    	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 the particle go out of range in space and in density
incardon's avatar
incardon committed
862
863
864
865
866
867
	    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)
	    {
	                   to_remove.add(a.getKey());
	    }
incardon's avatar
incardon committed
868
869
870
871
872
873
874
875
876

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

incardon's avatar
incardon committed
877
	// remove the particles
incardon's avatar
incardon committed
878
	vd.remove(to_remove,0);
incardon's avatar
incardon committed
879

incardon's avatar
incardon committed
880
	// increment the iteration counter
incardon's avatar
incardon committed
881
	cnt++;
incardon's avatar
incardon committed
882
883
}

incardon's avatar
incardon committed
884
885
/*! \cond [verlet_int] \endcond */

incardon's avatar
incardon committed
886
887
888
889
int main(int argc, char* argv[])
{
	/*!
	 *
incardon's avatar
incardon committed
890
891
	 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
	 *
incardon's avatar
incardon committed
892
	 * ## Main function {#e7_sph_main}
incardon's avatar
incardon committed
893
	 *
incardon's avatar
incardon committed
894
	 * Here we Initialize the library, we create a Box that define our domain, boundary conditions and ghost
incardon's avatar
incardon committed
895
896
897
898
899
900
901
902
903
904
905
906
907
	 *
	 * \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
908
909
	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
910
911
912
913
914
915
916
917
918
919
920
921
922

	// 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
923
	 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
incardon's avatar
incardon committed
924
	 *
incardon's avatar
incardon committed
925
	 * ### Vector create {#e7_sph_vcreate}
incardon's avatar
incardon committed
926
	 *
incardon's avatar
incardon committed
927
	 * Here we define a distributed vector in 3D, we use the particles type that we defined previously.
incardon's avatar
incardon committed
928
929
930
	 * Each particle contain the following properties
	 * * **type** Type of the particle
	 * * **rho** Density of the particle
incardon's avatar
incardon committed
931
932
933
934
935
936
937
938
 	 * * **rho_prev** Density at previous timestep
     * * **Pressure** Pressure of the particle
 	 * * **drho** Derivative of the density over time
	 * * **force** acceleration of the particles
	 * * **velocity** velocity of the particles
	 * * **velocity_prev** velocity of the particles at previous time-step
	 *
	 *
incardon's avatar
incardon committed
939
940
	 * In this case the vector contain 0 particles initially
	 *
incardon's avatar
incardon committed
941
942
943
944
945
946
947
948
949
950
	 * \see \ref e0_s_vector_inst
	 *
	 * The option DEC_GRAN(512) is related to the Load-Balancing decomposition
	 * granularity. It indicate that the space must be decomposed by at least
	 *
	 * \f$ N_{subsub} = 512 \cdot N_p \f$
	 *
	 * Where \f$ N_{subsub} \f$ is the number of sub-sub-domain in which the space
	 * must be decomposed and \f$ N_p \f$ is the number of processors. (The concept
	 * of sub-sub-domain will be explained leter)
incardon's avatar
incardon committed
951
	 *
incardon's avatar
incardon committed
952
953
	 * \snippet Vector/7_SPH_dlb/main.cpp vector inst
	 * \snippet Vector/7_SPH_dlb/main.cpp vector_dist_def
incardon's avatar
incardon committed
954
955
956
957
958
	 *
	 */

	//! \cond [vector inst] \endcond

incardon's avatar
incardon committed
959
	particles vd(0,domain,bc,g,DEC_GRAN(512));
incardon's avatar
incardon committed
960
961
962

	//! \cond [vector inst] \endcond

incardon's avatar
incardon committed
963
964
965
	/*!
	 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
	 *
incardon's avatar
incardon committed
966
	 * ### Draw particles and initialization ## {#e7_sph_draw_part_init}
incardon's avatar
incardon committed
967
	 *
incardon's avatar
incardon committed
968
969
970
971
	 * In this part we initialize the problem creating particles. In order to do it we use the class DrawParticles. Because some of
	 * the simulation constants require the maximum height \f$ h_{swl} \f$ of the fluid to be calculated
	 *  and the maximum fluid height is determined at runtime, some of the constants just after we create the
	 *  fluid particles
incardon's avatar
incardon committed
972
	 *
incardon's avatar
incardon committed
973
	 *  ### Draw Fluid ### {#e7_sph_draw_part_fluid}
incardon's avatar
incardon committed
974
	 *
incardon's avatar
incardon committed
975
976
977
978
	 * The Function DrawParticles::DrawBox return an iterator that can be used to create particle in a predefined
	 * box (smaller than the simulation domain) with a predefined spacing.
	 * We start drawing the fluid particles, the initial pressure is initialized accordingly to the
	 * Hydrostatic pressure given by:
incardon's avatar
incardon committed
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
	 *
	 *  \f$ P = \rho_{0} g (h_{max} - z) \f$
	 *
	 * Where \f$ h_{max} \f$ is the maximum height of the fluid.
	 * The density instead is given by the equation (3). Assuming \f$ \rho \f$ constant to
	 * \f$ \rho_{0} \f$ in the Hydrostatic equation is a good approximation. Velocity is
	 * initialized to zero.
	 *
	 * \see \ref e0_s_vector_inst
	 *
	 * \htmlonly
	 * <img src="http://ppmcore.mpi-cbg.de/web/images/examples/7_SPH_dlb/fluid.jpg"/>
	 * \endhtmlonly
	 *
	 * \snippet Vector/7_SPH_dlb/main.cpp draw fluid
	 *
	 */

	//! \cond [draw fluid] \endcond

incardon's avatar
incardon committed
999
1000
	// 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
incardon's avatar
incardon committed
1001
	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
1002

incardon's avatar
incardon committed
1003
	// return an iterator to the fluid particles to add to vd
incardon's avatar
incardon committed
1004
	auto fluid_it = DrawParticles::DrawBox(vd,sz,domain,fluid_box);
incardon's avatar
incardon committed
1005
1006

	// here we fill some of the constants needed by the simulation
incardon's avatar
incardon committed
1007
1008
1009
1010
1011
	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);

incardon's avatar
incardon committed
1012
	// for each particle inside the fluid box ...
incardon's avatar
incardon committed
1013
1014
	while (fluid_it.isNext())
	{
incardon's avatar
incardon committed
1015
		// ... add a particle ...
incardon's avatar
incardon committed
1016
1017
		vd.add();

incardon's avatar
incardon committed
1018
		// ... and set it position ...
incardon's avatar
incardon committed
1019
1020
1021
1022
		vd.getLastPos()[0] = fluid_it.get().get(0);
		vd.getLastPos()[1] = fluid_it.get().get(1);
		vd.getLastPos()[2] = fluid_it.get().get(2);

incardon's avatar
incardon committed
1023
		// and its type.
incardon's avatar
incardon committed
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
		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;

incardon's avatar
incardon committed
1045
		// next fluid particle
incardon's avatar
incardon committed
1046
1047
1048
		++fluid_it;
	}

incardon's avatar
incardon committed
1049
1050
1051
1052
1053
1054
1055
	//! \cond [draw fluid] \endcond

	/*!
	 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
	 *
	 * ### Draw Recipient ###
	 *
incardon's avatar
incardon committed
1056
1057
1058
	 * Here we draw the recipient using the function DrawParticles::DrawSkin. This function can draw a set
	 * of particles inside a box A removed of a second box or an array of boxes. So all the particles in the
	 *  area included in the area A - B - C. There is no restriction that B or C must be included into A.
incardon's avatar
incardon committed
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
	 *
	 * \htmlonly
	 * <img src="http://ppmcore.mpi-cbg.de/web/images/examples/7_SPH_dlb/recipient.jpg"/>
	 * \endhtmlonly
	 *
	 * In this case A is the box defining the recipient, B is the box cutting out the internal
	 * part of the recipient, C is the hole where we will place the obstacle.
     * Because we use Dynamic boundary condition (DBC) we initialize the density
	 * to \f$ \rho_{0} \f$. It will be update over time according to equation (3) to keep
	 * the particles confined.
	 *
	 * \snippet Vector/7_SPH_dlb/main.cpp draw recipient
	 *
	 */

	//! \cond [draw recipient] \endcond

incardon's avatar
incardon committed
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
	// 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;
	}

incardon's avatar
incardon committed
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
	//! \cond [draw recipient] \endcond

	/*!
	 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
	 *
	 *  ### Draw Obstacle ###
	 *
	 * Here we draw the obstacle in the same way we draw the recipient. also for the obstacle
	 * is valid the same concept of using Dynamic boundary condition (DBC)
	 *
	 * \htmlonly
	 * <img src="http://ppmcore.mpi-cbg.de/web/images/examples/7_SPH_dlb/obstacle.jpg"/>
	 * \endhtmlonly
	 *
	 * \snippet Vector/7_SPH_dlb/main.cpp draw obstacle
	 *
	 */

	//! \cond [draw obstacle] \endcond
incardon's avatar
incardon committed
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155

	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
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167

	//! \cond [draw obstacle] \endcond

	/*!
	 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
	 *
	 * ## Load balancing and Dynamic load balancing ##
	 *
	 * ### Load Balancing ###
	 *
	 * If at this point we output the particles and we visualize where they are accordingly
	 * to their processor id we can easily see that particles are distributed unevenly. The
incardon's avatar
incardon committed
1168
	 * processor that has particles in white has few particles and all of them are non fluid.
incardon's avatar
incardon committed
1169
1170
1171
1172
1173
1174
	 * This mean that it will be almost in idle. This situation is not ideal
	 *
	 * \htmlonly
	 * <img src="http://ppmcore.mpi-cbg.de/web/images/examples/7_SPH_dlb/unbalanced_particles.jpg"/>
	 * \endhtmlonly
	 *
incardon's avatar
incardon committed
1175
1176
1177
1178
1179
1180
1181
	 * In order to reach an optimal situation we have to distribute the particles to
	 * reach a balanced situation. To do this we have to set the computation of each
	 * sub-sub-domain, redecompose the space and distribute the particles accordingly to this
	 * new configuration. To do this we need a model. A model specify how to set
	 * the computational cost for each sub-sub-domains (for example it specify if the computational cost to
	 * process a sub-sub-domain is quadratic or linear with the number of
	 * particles ...). A model look like this.
incardon's avatar
incardon committed
1182
1183
1184
	 *
	 * \snippet Vector/7_SPH_dlb/main.cpp model custom
	 *
incardon's avatar
incardon committed
1185
1186
1187
1188
	 *  Setting the the computational cost on sub-sub-domains is performed running
	 *  across the particles. For each one of them, it is calculated on which sub-sub-domain it belong.
	 *   Than the function **addComputation** is called. Inside this call we can set the weight
	 *   in the way we prefer. In this case we set the weight as:
incardon's avatar
incardon committed
1189
	 *
incardon's avatar
incardon committed
1190
	 * \f$ w_v =  4 N_{fluid} + 3 N_{boundary} \f$
incardon's avatar
incardon committed
1191
	 *
incardon's avatar
incardon committed
1192
1193
1194
	 * Where \f$ N_{fluid} \f$ Is the number of fluid particles in the sub-sub-domains and \f$ N_{boundary} \f$
	 * are the number of boundary particles. For example in our ModelCustom we square this number,
	 *  because the computation is proportional to the square of the number of particles in each sub-sub-domain.
incardon's avatar
incardon committed
1195
	 * A second cycle is performed in order to calculate a complex function of this number (for example squaring).
incardon's avatar
incardon committed
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
	 *
	 * Implicitly the communication cost is given by \f$ \frac{V_{ghost}}{V_{sub-sub}}
	 * t_s \f$, while the migration cost is given by \f$ v_{sub-sub} \f$. In general\f$ t_s \f$ is the number
	 *  of ghost get between two rebalance. In this special case where we have two type of particles,
	 * we have two different computation for each of them, this mean that fluid particles
	 * and boundary particles has different computation cost.
	 *
	 *  After filling the computational cost based on our model
	 * we can decompose the problem in computationally equal chunk for each processor.
	 * We use the function **decomposed** to redecompose the space and subsequently we use
	 *  the function map to redistribute
incardon's avatar
incardon committed
1207
1208
1209
1210
	 * the particles.
	 *
	 * \note All processors now has part of the fluid. It is good to note that the computationaly
	 *       balanced configuration does not correspond to the evenly distributed particles to know
incardon's avatar
incardon committed
1211
	 *       more about that please follow the video tutorials
incardon's avatar
incardon committed
1212
	 *
incardon's avatar
incardon committed
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
	 * \htmlonly
	 * <a href="#" onclick="hide_show('vector-video-6')" >Dynamic load balancing the theory part1</a>
	 * <div style="display:none" id="vector-video-6">
	 * <video id="vid6" width="1200" height="576" controls> <source src="http://openfpm.mpi-cbg.de/upload/video/dlb-1.mp4" type="video/mp4"></video>
	 * </div>
	 * <a href="#" onclick="hide_show('vector-video-7')" >Dynamic load balancing the theory part2</a>
	 * <div style="display:none" id="vector-video-7">
	 * <video id="vid7" width="1200" height="576" controls> <source src="http://openfpm.mpi-cbg.de/upload/video/dlb-2.mp4" type="video/mp4"></video>
	 * </div>
	 * <a href="#" onclick="hide_show('vector-video-8')" >Dynamic load balancing practice part1</a>
	 * <div style="display:none" id="vector-video-8">
	 * <video id="vid8" width="1200" height="576" controls> <source src="http://openfpm.mpi-cbg.de/upload/video/dlb-2.mp4" type="video/mp4"></video>
	 * </div>
	 * \endhtmlonly
	 *
incardon's avatar
incardon committed
1228
1229
1230
1231
1232
1233
1234
1235
1236
	 * \snippet Vector/7_SPH_dlb/main.cpp load balancing
	 *
	 * \htmlonly
	 * <img src="http://ppmcore.mpi-cbg.de/web/images/examples/7_SPH_dlb/load_balanced_particles.jpg"/>
	 * \endhtmlonly
	 *
	 */

	//! \cond [load balancing] \endcond
incardon's avatar
incardon committed
1237
1238
1239
1240
1241
1242
1243
1244

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

	vd.addComputationCosts(md);
	vd.getDecomposition().decompose();
	vd.map();

incardon's avatar
incardon committed
1245
	//! \cond [load balancing] \endcond
incardon's avatar
incardon committed
1246
1247

	vd.ghost_get<type,rho,Pressure,velocity>();
incardon's avatar
incardon committed
1248
1249
1250
1251
1252

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

	// Evolve

incardon's avatar
incardon committed
1253
1254
1255
1256
1257
1258
1259
1260
	/*!
	 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
	 *
	 * ## Main Loop ##
	 *
	 * The main loop do time integration. It calculate the pressure based on the
	 * density, than calculate the forces, than we calculate delta time, and finally update position
	 * and velocity. After 200 time-step we do a rebalancing. And we save the configuration
incardon's avatar
incardon committed
1261
	 * every 0.01 seconds
incardon's avatar
incardon committed
1262
1263
1264
1265
1266
1267
	 *
	 * \snippet Vector/7_SPH_dlb/main.cpp main loop
	 *
	 */

	//! \cond [main loop] \endcond
incardon's avatar
incardon committed
1268
1269
1270

	size_t write = 0;
	size_t it = 0;
incardon's avatar
incardon committed
1271
	size_t it_reb = 0;
incardon's avatar
incardon committed
1272
1273
1274
	double t = 0.0;
	while (t <= t_end)
	{
incardon's avatar
incardon committed
1275
		Vcluster & v_cl = create_vcluster();
incardon's avatar
incardon committed
1276
1277
1278
1279
		timer it_time;

		////// Do rebalancing every 200 timesteps
		it_reb++;
incardon's avatar
incardon committed
1280
		if (it_reb == 200)
incardon's avatar
incardon committed
1281
1282
1283
		{
			vd.map();

incardon's avatar
incardon committed
1284
			it_reb = 0;
incardon's avatar
incardon committed
1285
1286
			ModelCustom md;
			vd.addComputationCosts(md);
incardon's avatar
incardon committed
1287
			vd.getDecomposition().decompose();
incardon's avatar
incardon committed
1288

incardon's avatar
incardon committed
1289
1290
			if (v_cl.getProcessUnitID() == 0)
				std::cout << "REBALANCED " << std::endl;
incardon's avatar
incardon committed
1291
		}
incardon's avatar
incardon committed
1292
1293

		vd.map();
incardon's avatar
incardon committed
1294

incardon's avatar
incardon committed
1295
1296
1297
1298
1299
		// Calculate pressure from the density
		EqState(vd);

		double max_visc = 0.0;

incardon's avatar
incardon committed
1300
		vd.ghost_get<type,rho,Pressure,velocity>();
incardon's avatar
incardon committed
1301

incardon's avatar
incardon committed
1302
1303
		// Calc forces
		calc_forces(vd,NN,max_visc);
incardon's avatar
incardon committed
1304
1305
1306
1307

		// Get the maximum viscosity term across processors
		v_cl.max(max_visc);
		v_cl.execute();
incardon's avatar
incardon committed
1308
1309
1310
1311

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

incardon's avatar
incardon committed
1312
		// VerletStep or euler step
incardon's avatar
incardon committed
1313
1314
		it++;
		if (it < 40)
incardon's avatar
incardon committed
1315
			verlet_int(vd,dt);
incardon's avatar
incardon committed
1316
1317
		else