main.cpp 45.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
 *
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
 * </div>
incardon's avatar
incardon committed
32 33 34 35 36 37 38 39 40 41 42 43
 * <a href="#" onclick="hide_show('vector-video-17')" >Simulation countour prospective 1</a><br>
 * <div style="display:none" id="vector-video-17">
 * <video id="vid17" width="1200" height="576" controls> <source src="http://openfpm.mpi-cbg.de/web/images/examples/7_SPH_dlb/sph_zoom.mp4" type="video/mp4"></video>
 * </div>
 * <a href="#" onclick="hide_show('vector-video-18')" >Simulation countour prospective 2</a><br>
 * <div style="display:none" id="vector-video-18">
 * <video id="vid18" width="1200" height="576" controls> <source src="http://openfpm.mpi-cbg.de/web/images/examples/7_SPH_dlb/sph_back.mp4" type="video/mp4"></video>
 * </div>
 * <a href="#" onclick="hide_show('vector-video-19')" >Simulation countour prospective 3</a><br>
 * <div style="display:none" id="vector-video-19">
 * <video id="vid19" width="1200" height="576" controls> <source src="http://openfpm.mpi-cbg.de/web/images/examples/7_SPH_dlb/sph_all.mp4" type="video/mp4"></video>
 * </div>
incardon's avatar
incardon committed
44 45 46
 * \endhtmlonly
 *
 * \htmlonly
incardon's avatar
incardon committed
47 48 49 50
 * <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
51 52
 *
 * In order to use distributed vectors in our code we have to include the file Vector/vector_dist.hpp
incardon's avatar
incardon committed
53 54
 * we also include DrawParticles that has nice utilities to draw particles in parallel accordingly
 * to simple shapes
incardon's avatar
incardon committed
55 56 57 58 59
 *
 * \snippet Vector/7_SPH_dlb/main.cpp inclusion
 *
 */

incardon's avatar
incardon committed
60 61 62
//#define SE_CLASS1
//#define STOP_ON_ERROR

incardon's avatar
incardon committed
63 64 65
//! \cond [inclusion] \endcond
#include "Vector/vector_dist.hpp"
#include <math.h>
incardon's avatar
incardon committed
66
#include "Draw/DrawParticles.hpp"
incardon's avatar
incardon committed
67 68
//! \cond [inclusion] \endcond

incardon's avatar
incardon committed
69 70 71
/*!
 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balacing
 *
incardon's avatar
incardon committed
72
 * ## SPH simulation {#e7_sph_parameters}
incardon's avatar
incardon committed
73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96
 *
 * 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
97 98
 * 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
99 100
 * 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
101 102 103 104 105
 *
 * ### Parameters {#e7_sph_parameters}
 *
 * Based on the equation
 * reported before several constants must be defined.
incardon's avatar
incardon committed
106 107 108 109 110 111
 *
 * \snippet Vector/7_SPH_dlb/main.cpp sim parameters
 *
 */

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

incardon's avatar
incardon committed
113
// A constant to indicate boundary particles
incardon's avatar
incardon committed
114 115
#define BOUNDARY 0

incardon's avatar
incardon committed
116 117
// A constant to indicate fluid particles
#define FLUID 1
incardon's avatar
incardon committed
118

incardon's avatar
incardon committed
119
// initial spacing between particles dp in the formulas
incardon's avatar
incardon committed
120
const double dp = 0.0085;
incardon's avatar
incardon committed
121
// Maximum height of the fluid water
incardon's avatar
incardon committed
122
// is going to be calculated and filled later on
incardon's avatar
incardon committed
123
double h_swl = 0.0;
incardon's avatar
incardon committed
124

incardon's avatar
incardon committed
125
// c_s in the formulas (constant used to calculate the sound speed)
incardon's avatar
incardon committed
126
const double coeff_sound = 20.0;
incardon's avatar
incardon committed
127 128

// gamma in the formulas
incardon's avatar
incardon committed
129
const double gamma_ = 7.0;
incardon's avatar
incardon committed
130 131

// sqrt(3.0*dp*dp) support of the kernel
incardon's avatar
incardon committed
132
const double H = 0.0147224318643;
incardon's avatar
incardon committed
133 134

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

incardon's avatar
incardon committed
137
// alpha in the formula
incardon's avatar
incardon committed
138
const double visco = 0.1;
incardon's avatar
incardon committed
139 140

// cbar in the formula (calculated later)
incardon's avatar
incardon committed
141
double cbar = 0.0;
incardon's avatar
incardon committed
142 143

// Mass of the fluid particles
incardon's avatar
incardon committed
144
const double MassFluid = 0.000614125;
incardon's avatar
incardon committed
145 146

// Mass of the boundary particles
incardon's avatar
incardon committed
147
const double MassBound = 0.000614125;
incardon's avatar
incardon committed
148 149

// End simulation time
150 151 152
#ifdef TEST_RUN
const double t_end = 0.001;
#else
incardon's avatar
incardon committed
153
const double t_end = 1.5;
154
#endif
incardon's avatar
incardon committed
155 156

// Gravity acceleration
incardon's avatar
incardon committed
157
const double gravity = 9.81;
incardon's avatar
incardon committed
158 159

// Reference densitu 1000Kg/m^3
incardon's avatar
incardon committed
160
const double rho_zero = 1000.0;
incardon's avatar
incardon committed
161 162

// Filled later require h_swl, it is b in the formulas
incardon's avatar
incardon committed
163
double B = 0.0;
incardon's avatar
incardon committed
164 165

// Constant used to define time integration
incardon's avatar
incardon committed
166
const double CFLnumber = 0.2;
incardon's avatar
incardon committed
167 168

// Minimum T
incardon's avatar
incardon committed
169 170
const double DtMin = 0.00001;

incardon's avatar
incardon committed
171 172 173 174 175 176
// Minimum Rho allowed
const double RhoMin = 700.0;

// Maximum Rho allowed
const double RhoMax = 1300.0;

incardon's avatar
incardon committed
177 178 179
// Filled in initialization
double max_fluid_height = 0.0;

incardon's avatar
incardon committed
180 181 182
// Properties

// FLUID or BOUNDARY
incardon's avatar
incardon committed
183
const size_t type = 0;
incardon's avatar
incardon committed
184 185

// Density
incardon's avatar
incardon committed
186
const int rho = 1;
incardon's avatar
incardon committed
187 188

// Density at step n-1
incardon's avatar
incardon committed
189
const int rho_prev = 2;
incardon's avatar
incardon committed
190 191

// Pressure
incardon's avatar
incardon committed
192
const int Pressure = 3;
incardon's avatar
incardon committed
193 194

// Delta rho calculated in the force calculation
incardon's avatar
incardon committed
195
const int drho = 4;
incardon's avatar
incardon committed
196 197

// calculated force
incardon's avatar
incardon committed
198
const int force = 5;
incardon's avatar
incardon committed
199 200

// velocity
incardon's avatar
incardon committed
201
const int velocity = 6;
incardon's avatar
incardon committed
202 203

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

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

incardon's avatar
incardon committed
208 209
/*! \cond [vector_dist_def] \endcond */

incardon's avatar
incardon committed
210 211 212 213 214 215
// 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
216

incardon's avatar
incardon committed
217
/*! \cond [vector_dist_def] \endcond */
incardon's avatar
incardon committed
218

incardon's avatar
incardon committed
219
/*! \cond [model custom] \endcond */
incardon's avatar
incardon committed
220 221 222

struct ModelCustom
{
incardon's avatar
incardon committed
223
	template<typename Decomposition, typename vector> inline void addComputation(Decomposition & dec,
224
			                                                                     vector & vd,
incardon's avatar
incardon committed
225 226
																				 size_t v,
																				 size_t p)
incardon's avatar
incardon committed
227 228
	{
		if (vd.template getProp<type>(p) == FLUID)
incardon's avatar
incardon committed
229
			dec.addComputationCost(v,4);
incardon's avatar
incardon committed
230
		else
incardon's avatar
incardon committed
231
			dec.addComputationCost(v,3);
incardon's avatar
incardon committed
232 233 234 235 236 237
	}

	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
238 239 240 241 242

	double distributionTol()
	{
		return 1.01;
	}
incardon's avatar
incardon committed
243 244
};

incardon's avatar
incardon committed
245
/*! \cond [model custom] \endcond */
incardon's avatar
incardon committed
246 247 248 249

/*!
 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balacing
 *
incardon's avatar
incardon committed
250
 * ### Equation of state {#e7_sph_equation_state}
incardon's avatar
incardon committed
251 252 253 254 255 256 257 258 259 260 261 262 263
 *
 * 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
264 265 266 267 268 269 270 271 272 273 274 275 276
	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
277
}
incardon's avatar
incardon committed
278

incardon's avatar
incardon committed
279
/*! \cond [eq_state_and_ker] \endcond */
incardon's avatar
incardon committed
280

incardon's avatar
incardon committed
281 282 283
/*!
 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
 *
incardon's avatar
incardon committed
284 285 286
 * ### 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
287 288 289 290 291 292 293
 * 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
294

incardon's avatar
incardon committed
295
/*! \cond [kernel_sph] \endcond */
incardon's avatar
incardon committed
296 297 298 299 300 301 302 303 304 305 306 307 308 309 310

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
311 312 313 314 315
/*! \cond [kernel_sph] \endcond */

/*!
 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
 *
incardon's avatar
incardon committed
316
 * This function define the gradient of the Cubic kernel function \f$ W_{ab} \f$.
incardon's avatar
incardon committed
317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333
 *
 * \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
334 335 336 337 338

inline void DWab(Point<3,double> & dx, Point<3,double> & DW, double r, bool print)
{
	const double qq=r/H;

339 340 341
    double qq2 = qq * qq;
    double fac1 = (c1*qq + d1*qq2)/r;
    double b1 = (qq < 1.0)?1.0f:0.0f;
incardon's avatar
incardon committed
342

343 344 345
    double wqq = (2.0 - qq);
    double fac2 = c2 * wqq * wqq / r;
    double b2 = (qq >= 1.0 && qq < 2.0)?1.0f:0.0f;
incardon's avatar
incardon committed
346

347 348 349 350 351
    double factor = (b1*fac1 + b2*fac2);

    DW.get(0) = factor * dx.get(0);
    DW.get(1) = factor * dx.get(1);
    DW.get(2) = factor * dx.get(2);
incardon's avatar
incardon committed
352 353
}

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


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

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

incardon's avatar
incardon committed
465
	// For each particle ...
incardon's avatar
incardon committed
466 467
	while (part.isNext())
	{
incardon's avatar
incardon committed
468
		// ... a
incardon's avatar
incardon committed
469 470 471 472 473
		auto a = part.get();

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

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

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

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

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

incardon's avatar
incardon committed
486
		// Reset the force counter (- gravity on zeta direction)
incardon's avatar
incardon committed
487 488 489 490 491
		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
492 493 494 495 496 497
		// 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
498

incardon's avatar
incardon committed
499 500 501 502 503
			// For each neighborhood particle
			while (Np.isNext() == true)
			{
				// ... q
				auto b = Np.get();
incardon's avatar
incardon committed
504

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

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

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

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

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

incardon's avatar
incardon committed
521 522 523 524
				// 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
525

incardon's avatar
incardon committed
526 527 528 529 530
				// If the particles interact ...
				if (r2 < 4.0*H*H)
				{
					// ... calculate delta rho
					double r = sqrt(r2);
incardon's avatar
incardon committed
531

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

incardon's avatar
incardon committed
534 535
					Point<3,double> DW;
					DWab(dr,DW,r,false);
incardon's avatar
incardon committed
536

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

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

incardon's avatar
incardon committed
544
				++Np;
incardon's avatar
incardon committed
545 546
			}
		}
incardon's avatar
incardon committed
547
		else
incardon's avatar
incardon committed
548
		{
incardon's avatar
incardon committed
549
			// If it is a fluid particle calculate based on equation 1 and 2
incardon's avatar
incardon committed
550

incardon's avatar
incardon committed
551 552
			// 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
553

incardon's avatar
incardon committed
554 555 556 557 558
			// For each neighborhood particle
			while (Np.isNext() == true)
			{
				// ... q
				auto b = Np.get();
incardon's avatar
incardon committed
559

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

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

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

incardon's avatar
incardon committed
571 572 573 574
				// 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
575

incardon's avatar
incardon committed
576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596
				// 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
597
		}
incardon's avatar
incardon committed
598 599

		++part;
incardon's avatar
incardon committed
600
	}
incardon's avatar
incardon committed
601
}
incardon's avatar
incardon committed
602

incardon's avatar
incardon committed
603
/*! \cond [calc_forces] \endcond */
incardon's avatar
incardon committed
604

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

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

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
646 647 648 649 650

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

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

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
698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715
/*! \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
716
 * \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
717
 *
incardon's avatar
incardon committed
718
 * \f$ \rho_a^{n+1} = \rho_a^n + \delta t D_a^n \f$
incardon's avatar
incardon committed
719
 *
incardon's avatar
incardon committed
720
 * This function also check that no particles go outside the simulation
incardon's avatar
incardon committed
721 722
 * 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
723 724 725 726 727 728 729 730
 *
 * \snippet Vector/7_SPH_dlb/main.cpp verlet_int
 *
 *
 */

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

incardon's avatar
incardon committed
731 732
openfpm::vector<size_t> to_remove;

incardon's avatar
incardon committed
733
size_t cnt = 0;
incardon's avatar
incardon committed
734

incardon's avatar
incardon committed
735
void verlet_int(particles & vd, double dt)
incardon's avatar
incardon committed
736
{
incardon's avatar
incardon committed
737
	// list of the particle to remove
incardon's avatar
incardon committed
738 739
	to_remove.clear();

incardon's avatar
incardon committed
740
	// particle iterator
incardon's avatar
incardon committed
741 742 743 744 745
	auto part = vd.getDomainIterator();

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

incardon's avatar
incardon committed
746
	// For each particle ...
incardon's avatar
incardon committed
747 748
	while (part.isNext())
	{
incardon's avatar
incardon committed
749
		// ... a
incardon's avatar
incardon committed
750 751
		auto a = part.get();

incardon's avatar
incardon committed
752
		// if the particle is boundary
incardon's avatar
incardon committed
753 754
		if (vd.template getProp<type>(a) == BOUNDARY)
		{
incardon's avatar
incardon committed
755
			// Update rho
incardon's avatar
incardon committed
756 757
			double rhop = vd.template getProp<rho>(a);

incardon's avatar
incardon committed
758
			// Update only the density
incardon's avatar
incardon committed
759 760 761 762
	    	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
763 764 765

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

incardon's avatar
incardon committed
766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783
			++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
784 785 786 787 788 789 790 791 792
    	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
793
	    {
incardon's avatar
incardon committed
794
	                   to_remove.add(a.getKey());
incardon's avatar
incardon committed
795
	    }
incardon's avatar
incardon committed
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 834 835 836 837 838

	    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
839 840
	    	vd.template getProp<rho>(a) = vd.template getProp<rho>(a) + dt*vd.template getProp<drho>(a);

incardon's avatar
incardon committed
841 842 843 844 845
		    vd.template getProp<rho_prev>(a) = rhop;

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

incardon's avatar
incardon committed
847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866
		//-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
867 868 869 870 871 872
	    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
873 874 875 876 877 878 879 880 881

	    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
882
	// remove the particles
incardon's avatar
incardon committed
883
	vd.remove(to_remove,0);
incardon's avatar
incardon committed
884

incardon's avatar
incardon committed
885
	// increment the iteration counter
incardon's avatar
incardon committed
886
	cnt++;
incardon's avatar
incardon committed
887 888
}

incardon's avatar
incardon committed
889 890
/*! \cond [verlet_int] \endcond */

incardon's avatar
incardon committed
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 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 984 985 986 987 988 989 990 991
/*!
 *
 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
 *
 * ### Probes/sensors {#e7_sph_prob_sens}
 *
 * This function show how to create a pressure sensor/probe on a set of specified points. To do this
 * from the cell-list we just get an iterator across the neighborhood points of the sensors and we
 * calculate the pressure profile. On the other hand because the sensor is in the processor domain
 * of only one processor, only one processor must do this calculation. We will use the function isLocal
 * to determine which processor contain the probe and only such processor will do the calculation.
 *
 * \warning This type of calculation is suitable if the number of probes is small (like 10) and pressure is not
 * calculated every time step. In case the number of
 * probes is comparable to the number of particles or the pressure is calculated every time-step than we suggest
 *  to create a set of "probe" particles
 *
 *
 * \snippet Vector/7_SPH_dlb/main.cpp sens_press
 *
 *
 */

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

template<typename Vector, typename CellList>
inline void sensor_pressure(Vector & vd,
                            CellList & NN,
                            openfpm::vector<openfpm::vector<double>> & press_t,
                            openfpm::vector<Point<3,double>> & probes)
{
    Vcluster & v_cl = create_vcluster();

    press_t.add();

    for (size_t i = 0 ; i < probes.size() ; i++)
    {
        float press_tmp = 0.0f;
        float tot_ker = 0.0;

        // if the probe is inside the processor domain
		if (vd.getDecomposition().isLocal(probes.get(i)) == true)
		{
			// Get the position of the probe i
			Point<3,double> xp = probes.get(i);

			// get the iterator over the neighbohood particles of the probes position
			auto itg = NN.template getNNIterator<NO_CHECK>(NN.getCell(probes.get(i)));
			while (itg.isNext())
			{
				auto q = itg.get();

				// Only the fluid particles are importants
				if (vd.template getProp<type>(q) != FLUID)
				{
					++itg;
					continue;
				}

				// Get the position of the neighborhood particle q
				Point<3,double> xq = vd.template getPos(q);

				// Calculate the contribution of the particle to the pressure
				// of the probe
				double r = sqrt(norm2(xp - xq));

				double ker = Wab(r) * (MassFluid / rho_zero);

				// Also keep track of the calculation of the summed
				// kernel
				tot_ker += ker;

				// Add the total pressure contribution
				press_tmp += vd.template getProp<Pressure>(q) * ker;

				// next neighborhood particle
				++itg;
			}

			// We calculate the pressure normalizing the
			// sum over all kernels
			if (tot_ker == 0.0)
				press_tmp = 0.0;
			else
				press_tmp = 1.0 / tot_ker * press_tmp;

		}

		// This is not necessary in principle, but if you
		// want to make all processor aware of the history of the calculated
		// pressure we have to execute this
		v_cl.sum(press_tmp);
		v_cl.execute();

		// We add the calculated pressure into the history
		press_t.last().add(press_tmp);
	}
}

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

incardon's avatar
incardon committed
992 993 994 995
int main(int argc, char* argv[])
{
	/*!
	 *
incardon's avatar
incardon committed
996 997
	 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
	 *
incardon's avatar
incardon committed
998
	 * ## Main function {#e7_sph_main}
incardon's avatar
incardon committed
999
	 *
incardon's avatar
incardon committed
1000 1001
	 * Here we Initialize the library, we create a Box that define our domain, boundary conditions and ghost. We also create
	 * a vector that contain two probes to measure pressure
incardon's avatar
incardon committed
1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013
	 *
	 * \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);

incardon's avatar
incardon committed
1014 1015 1016 1017 1018 1019 1020
	// It contain for each time-step the value detected by the probes
	openfpm::vector<openfpm::vector<double>> press_t;
	openfpm::vector<Point<3,double>> probes;

	probes.add({0.8779,0.3,0.02});
	probes.add({0.754,0.31,0.02});

incardon's avatar
incardon committed
1021
	// Here we define our domain a 2D box with internals from 0 to 1.0 for x and y
incardon's avatar
incardon committed
1022 1023
	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
1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036

	// 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
1037
	 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
incardon's avatar
incardon committed
1038
	 *
incardon's avatar
incardon committed
1039
	 * ### Vector create {#e7_sph_vcreate}
incardon's avatar
incardon committed
1040
	 *
incardon's avatar
incardon committed
1041
	 * Here we define a distributed vector in 3D, we use the particles type that we defined previously.
incardon's avatar
incardon committed
1042 1043 1044
	 * Each particle contain the following properties
	 * * **type** Type of the particle
	 * * **rho** Density of the particle
incardon's avatar
incardon committed
1045 1046 1047 1048 1049 1050 1051 1052
 	 * * **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
1053 1054
	 * In this case the vector contain 0 particles initially
	 *
incardon's avatar
incardon committed
1055 1056 1057 1058 1059 1060 1061 1062 1063 1064
	 * \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
1065
	 *
incardon's avatar
incardon committed
1066 1067
	 * \snippet Vector/7_SPH_dlb/main.cpp vector inst
	 * \snippet Vector/7_SPH_dlb/main.cpp vector_dist_def
incardon's avatar
incardon committed
1068 1069 1070 1071 1072
	 *
	 */

	//! \cond [vector inst] \endcond

incardon's avatar
incardon committed
1073
	particles vd(0,domain,bc,g,DEC_GRAN(512));
incardon's avatar
incardon committed
1074 1075 1076

	//! \cond [vector inst] \endcond

incardon's avatar
incardon committed
1077 1078 1079
	/*!
	 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
	 *
incardon's avatar
incardon committed
1080
	 * ### Draw particles and initialization ## {#e7_sph_draw_part_init}
incardon's avatar
incardon committed
1081
	 *
incardon's avatar
incardon committed
1082 1083 1084 1085
	 * 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
1086
	 *
incardon's avatar
incardon committed
1087
	 *  ### Draw Fluid ### {#e7_sph_draw_part_fluid}
incardon's avatar
incardon committed
1088
	 *
incardon's avatar
incardon committed
1089 1090 1091 1092
	 * 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
1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112
	 *
	 *  \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
1113 1114
	// 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
1115
	Box<3,double> fluid_box({dp/2.0,dp/2.0,dp/2.0},{0.4+dp/2.0,0.67-dp/2.0,0.3+dp/2.0});
incardon's avatar
incardon committed
1116

incardon's avatar
incardon committed
1117
	// return an iterator to the fluid particles to add to vd
incardon's avatar
incardon committed
1118
	auto fluid_it = DrawParticles::DrawBox(vd,sz,domain,fluid_box);
incardon's avatar
incardon committed
1119 1120

	// here we fill some of the constants needed by the simulation
incardon's avatar
incardon committed
1121 1122 1123 1124 1125
	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
1126
	// for each particle inside the fluid box ...
incardon's avatar
incardon committed
1127 1128
	while (fluid_it.isNext())
	{
incardon's avatar
incardon committed
1129
		// ... add a particle ...
incardon's avatar
incardon committed
1130 1131
		vd.add();

incardon's avatar
incardon committed
1132
		// ... and set it position ...
incardon's avatar
incardon committed
1133 1134 1135 1136
		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
1137
		// and its type.
incardon's avatar
incardon committed
1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158
		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
1159
		// next fluid particle
incardon's avatar
incardon committed
1160 1161 1162
		++fluid_it;
	}

incardon's avatar
incardon committed
1163 1164 1165 1166 1167 1168 1169
	//! \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
1170 1171 1172
	 * 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
1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189
	 *
	 * \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
1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224
	// 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
1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243
	//! \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
1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269

	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
1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281

	//! \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
1282
	 * processor that has particles in white has few particles and all of them are non fluid.
incardon's avatar
incardon committed
1283 1284 1285 1286 1287 1288
	 * 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
1289 1290 1291 1292 1293 1294 1295
	 * 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
1296 1297 1298
	 *
	 * \snippet Vector/7_SPH_dlb/main.cpp model custom
	 *
incardon's avatar
incardon committed
1299 1300 1301 1302
	 *  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
1303
	 *
incardon's avatar
incardon committed
1304
	 * \f$ w_v =  4 N_{fluid} + 3 N_{boundary} \f$
incardon's avatar
incardon committed
1305
	 *
incardon's avatar
incardon committed
1306 1307 1308
	 * 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
1309
	 * A second cycle is performed in order to calculate a complex function of this number (for example squaring).
incardon's avatar
incardon committed
1310 1311 1312 1313 1314 1315 1316 1317 1318 1319 1320
	 *
	 * 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
1321 1322 1323 1324
	 * 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
1325
	 *       more about that please follow the video tutorials
incardon's avatar
incardon committed
1326
	 *
incardon's avatar
incardon committed
1327
	 * \htmlonly
incardon's avatar
incardon committed
1328
	 * <a href="#" onclick="hide_show('vector-video-6')" >Dynamic load balancing the theory part1</a><br>
incardon's avatar
incardon committed
1329 1330 1331
	 * <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>
incardon's avatar
incardon committed
1332
	 * <a href="#" onclick="hide_show('vector-video-7')" >Dynamic load balancing the theory part2</a><br>
incardon's avatar
incardon committed
1333 1334 1335
	 * <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>
incardon's avatar
incardon committed
1336
	 * <a href="#" onclick="hide_show('vector-video-8')" >Dynamic load balancing practice part1</a><br>
incardon's avatar
incardon committed
1337
	 * <div style="display:none" id="vector-video-8">
incardon's avatar
incardon committed
1338 1339 1340 1341 1342
	 * <video id="vid8" width="1200" height="576" controls> <source src="http://openfpm.mpi-cbg.de/upload/video/dlb-3.mp4" type="video/mp4"></video>
	 * </div>
	 * <a href="#" onclick="hide_show('vector-video-9')" >Dynamic load balancing practice part2</a><br>
	 * <div style="display:none" id="vector-video-9">
	 * <video id="vid9" width="1200" height="576" controls> <source src="http://openfpm.mpi-cbg.de/upload/video/dlb-4.mp4" type="video/mp4"></video>
incardon's avatar
incardon committed
1343 1344 1345
	 * </div>
	 * \endhtmlonly
	 *
incardon's avatar
incardon committed
1346 1347 1348 1349 1350 1351 1352 1353 1354
	 * \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
1355 1356 1357 1358 1359 1360 1361 1362

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

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

incardon's avatar
incardon committed
1363
	//! \cond [load balancing] \endcond
incardon's avatar
incardon committed
1364 1365

	vd.ghost_get<type,rho,Pressure,velocity>();
incardon's avatar
incardon committed
1366 1367 1368 1369 1370

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

	// Evolve

incardon's avatar
incardon committed
1371 1372 1373 1374 1375 1376 1377
	/*!
	 * \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
incardon's avatar
incardon committed
1378 1379
	 * and velocity. After 200 time-step we do a re-balancing. We save the configuration
	 * and we calculate the pressure on the probe position every 0.01 seconds
incardon's avatar
incardon committed
1380 1381 1382 1383 1384 1385
	 *
	 * \snippet Vector/7_SPH_dlb/main.cpp main loop
	 *
	 */

	//! \cond [main loop] \endcond
incardon's avatar
incardon committed
1386 1387 1388

	size_t write = 0;
	size_t it = 0;
incardon's avatar
incardon committed
1389
	size_t it_reb = 0;
incardon's avatar
incardon committed
1390 1391 1392
	double t = 0.0;
	while (t <= t_end)
	{
incardon's avatar
incardon committed
1393
		Vcluster & v_cl = create_vcluster();
incardon's avatar
incardon committed
1394 1395 1396 1397
		timer it_time;

		////// Do rebalancing every 200 timesteps
		it_reb++;
incardon's avatar
incardon committed
1398
		if (it_reb == 200)
incardon's avatar
incardon committed
1399 1400 1401
		{
			vd.map();

incardon's avatar
incardon committed
1402
			it_reb = 0;
incardon's avatar
incardon committed
1403 1404
			ModelCustom md;
			vd.addComputationCosts(md);
incardon's avatar
incardon committed
1405
			vd.getDecomposition().decompose();
incardon's avatar
incardon committed
1406

incardon's avatar
incardon committed
1407 1408
			if (v_cl.getProcessUnitID() == 0)
				std::cout << "REBALANCED " << std::endl;
incardon's avatar
incardon committed
1409
		}
incardon's avatar
incardon committed
1410 1411

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

incardon's avatar
incardon committed
1413 1414 1415 1416 1417
		// Calculate pressure from the density
		EqState(vd);

		double max_visc = 0.0;

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

incardon's avatar
incardon committed
1420 1421
		// Calc forces
		calc_forces(vd,NN,max_visc);
incardon's avatar
incardon committed
1422 1423 1424 1425

		// Get the maximum viscosity term across processors
		v_cl.max(max_visc);
		v_cl.execute();
incardon's avatar
incardon committed
1426 1427 1428 1429

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

incardon's avatar
incardon committed
1430
		// VerletStep or euler step
incardon's avatar
incardon committed
1431 1432
		it++;
		if (it < 40)
incardon's avatar
incardon committed
1433
			verlet_int(vd,dt);
incardon's avatar
incardon committed
1434 1435
		else
		{
incardon's avatar
incardon committed
1436
			euler_int(vd,dt);
incardon's avatar
incardon committed
1437 1438
			it = 0;
		}
incardon's avatar
incardon committed
1439 1440 1441 1442 1443

		t += dt;

		if (write < t*100)
		{
incardon's avatar
incardon committed
1444 1445
			// calculate the pressure at the sensor points
			sensor_pressure(vd,NN,press_t,probes);
incardon's avatar
incardon committed
1446

incardon's avatar
incardon committed
1447
			vd.write_frame("Geometry",write);
incardon's avatar
incardon committed
1448
			write++;
incardon's avatar
incardon committed
1449

incardon's avatar
incardon committed
1450 1451
			if (v_cl.getProcessUnitID() == 0)
				std::cout << "TIME: " << t << "  write " << it_time.getwct() << "   " << v_cl.getProcessUnitID() << "   " << cnt << std::endl;
incardon's avatar
incardon committed
1452 1453 1454
		}
		else
		{
incardon's avatar
incardon committed
1455 1456
			if (v_cl.getProcessUnitID() == 0)
				std::cout << "TIME: " << t << "  " << it_time.getwct() << "   " << v_cl.getProcessUnitID() << "   " << cnt << std::endl;
incardon's avatar
incardon committed
1457 1458 1459
		}
	}

incardon's avatar
incardon committed
1460 1461 1462 1463 1464 1465 1466 1467 1468 1469 1470 1471 1472 1473 1474
	//! \cond [main loop] \endcond

	/*!
	 *
	 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
	 *
	 * ## Finalize ## {#finalize_e0_sim}
	 *
	 *
	 *  At the very end of the program we have always de-initialize the library
	 *
	 * \snippet Vector/7_SPH_dlb/main.cpp finalize
	 *
	 */

incardon's avatar
incardon committed
1475 1476 1477 1478 1479 1480 1481
	//! \cond [finalize] \endcond

	openfpm_finalize();

	//! \cond [finalize] \endcond

	/*!
incardon's avatar
incardon committed
1482
	 * \page Vector_7_sph_dlb Vector 7 SPH Dam break  simulation with Dynamic load balancing
incardon's avatar
incardon committed
1483
	 *
incardon's avatar
incardon committed
1484
	 * ## Full code ## {#code_e7_sph_dlb}
incardon's avatar
incardon committed
1485
	 *
incardon's avatar
incardon committed
1486
	 * \include Vector/7_SPH_dlb/main.cpp
incardon's avatar
incardon committed
1487 1488 1489 1490
	 *
	 */
}