main_vic_petsc_opt.cpp 28 KB
Newer Older
1
/*! \page Vortex_in_cell_petsc_opt Vortex in Cell 3D (Optimization)
incardon's avatar
incardon committed
2
 *
3
 * # Vortex in Cell 3D ring with ringlets optimization # {#vic_ringlets_optimization}
4 5
 *
 * In this example we solve the Navier-Stokes equation in the vortex formulation in 3D
6 7
 * for an incompressible fluid. This example
 * has the following changes compared to \ref Vortex_in_cell_petsc
8
 *
9 10 11
 * * Constructing grid is expensive in particular with a lot of cores. For this
 *   reason we create the grid in the main function rather  than in the **comp_vel**
 *   function and **helmotz_hodge_projection**
incardon's avatar
incardon committed
12
 *
13
 *
14 15 16
 * * Constructing also FDScheme is expensive so we construct it once in the main. We set
 *   the left hand side to the poisson operator, and inside the functions **comp_vel**
 *   and **helmotz_hodge_projection** just write the right hand side with **impose_dit_b**
17
 *
18
 * \snippet Numerics/Vortex_in_cell/main_vic_petsc_opt.cpp construct grids
incardon's avatar
incardon committed
19
 *
20
 * \snippet Numerics/Vortex_in_cell/main_vic_petsc_opt.cpp create b
incardon's avatar
incardon committed
21
 *
22 23
 * Another optimization that we do is to use a Hilbert space-filling curve as sub-sub-domain
 * distribution strategy
incardon's avatar
incardon committed
24 25 26 27 28 29
 *
 */

//#define SE_CLASS1
//#define PRINT_STACKTRACE

incardon's avatar
incardon committed
30
#include "interpolation/interpolation.hpp"
incardon's avatar
incardon committed
31 32 33 34 35 36 37
#include "Grid/grid_dist_id.hpp"
#include "Vector/vector_dist.hpp"
#include "Matrix/SparseMatrix.hpp"
#include "Vector/Vector.hpp"
#include "FiniteDifference/FDScheme.hpp"
#include "Solvers/petsc_solver.hpp"
#include "interpolation/mp4_kernel.hpp"
incardon's avatar
incardon committed
38
#include "Solvers/petsc_solver_AMG_report.hpp"
39
#include "Decomposition/Distribution/SpaceDistribution.hpp"
incardon's avatar
incardon committed
40 41 42 43

constexpr int x = 0;
constexpr int y = 1;
constexpr int z = 2;
incardon's avatar
incardon committed
44
constexpr unsigned int phi = 0;
incardon's avatar
incardon committed
45

46
// The type of the grids
47 48 49 50
typedef grid_dist_id<3,float,aggregate<float[3]>,CartDecomposition<3,float,HeapMemory,SpaceDistribution<3,float>>> grid_type;

// The type of the grids
typedef grid_dist_id<3,float,aggregate<float>,CartDecomposition<3,float,HeapMemory,SpaceDistribution<3,float>>> grid_type_s;
51 52

// The type of the particles
53 54 55
typedef vector_dist<3,float,aggregate<float[3],float[3],float[3],float[3],float[3]>,memory_traits_lin<aggregate<float[3],float[3],float[3],float[3],float[3]>>::type,memory_traits_lin,CartDecomposition<3,float,HeapMemory,SpaceDistribution<3,float>>> particles_type;

typedef vector_dist<3,float,aggregate<float>,memory_traits_lin<aggregate<float>>::type,memory_traits_lin,CartDecomposition<3,float,HeapMemory,SpaceDistribution<3,float>>> particles_type_s;
incardon's avatar
incardon committed
56 57

// radius of the torus
incardon's avatar
incardon committed
58
float ringr1 = 1.0;
incardon's avatar
incardon committed
59
// radius of the core of the torus
incardon's avatar
incardon committed
60
float sigma = 1.0/3.523;
incardon's avatar
incardon committed
61 62
// Reynold number (If you want to use 7500.0 you have to use a grid 1600x400x400)
//float tgtre  = 7500.0;
incardon's avatar
incardon committed
63
float tgtre = 3000.0;
incardon's avatar
incardon committed
64

incardon's avatar
incardon committed
65
// Kinematic viscosity
incardon's avatar
incardon committed
66
float nu = 1.0/tgtre;
incardon's avatar
incardon committed
67

68
// Time step
incardon's avatar
incardon committed
69
// float dt = 0.0025 for Re 7500
incardon's avatar
incardon committed
70
float dt = 0.0125;
incardon's avatar
incardon committed
71

incardon's avatar
incardon committed
72 73 74 75 76 77
#ifdef TEST_RUN
const unsigned int nsteps = 10;
#else
const unsigned int nsteps = 10001;
#endif

78
// All the properties by index
incardon's avatar
incardon committed
79 80 81 82 83 84 85
constexpr unsigned int vorticity = 0;
constexpr unsigned int velocity = 0;
constexpr unsigned int p_vel = 1;
constexpr int rhs_part = 2;
constexpr unsigned int old_vort = 3;
constexpr unsigned int old_pos = 4;

86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
template<typename grid> void calc_and_print_max_div_and_int(grid & g_vort)
{
	g_vort.template ghost_get<vorticity>();
	auto it5 = g_vort.getDomainIterator();

	double max_vort = 0.0;

	double int_vort[3] = {0.0,0.0,0.0};

	while (it5.isNext())
	{
		auto key = it5.get();

		double tmp;

        tmp = 1.0f/g_vort.spacing(x)/2.0f*(g_vort.template get<vorticity>(key.move(x,1))[x] - g_vort.template get<vorticity>(key.move(x,-1))[x] ) +
			         	 	 1.0f/g_vort.spacing(y)/2.0f*(g_vort.template get<vorticity>(key.move(y,1))[y] - g_vort.template get<vorticity>(key.move(y,-1))[y] ) +
							 1.0f/g_vort.spacing(z)/2.0f*(g_vort.template get<vorticity>(key.move(z,1))[z] - g_vort.template get<vorticity>(key.move(z,-1))[z] );

        int_vort[x] += g_vort.template get<vorticity>(key)[x];
        int_vort[y] += g_vort.template get<vorticity>(key)[y];
        int_vort[z] += g_vort.template get<vorticity>(key)[z];

        if (tmp > max_vort)
        	max_vort = tmp;

		++it5;
	}

	Vcluster & v_cl = create_vcluster();
	v_cl.max(max_vort);
	v_cl.sum(int_vort[0]);
	v_cl.sum(int_vort[1]);
	v_cl.sum(int_vort[2]);
	v_cl.execute();

	if (v_cl.getProcessUnitID() == 0)
	{std::cout << "Max div for vorticity " << max_vort << "   Integral: " << int_vort[0] << "  " << int_vort[1] << "   " << int_vort[2] << std::endl;}
}

incardon's avatar
incardon committed
126
void init_ring(grid_type & gr, const Box<3,float> & domain)
incardon's avatar
incardon committed
127
{
128 129
	// To add some noise to the vortex ring we create two random
	// vector
incardon's avatar
incardon committed
130 131
	constexpr int nk = 32;

incardon's avatar
incardon committed
132 133
	float ak[nk];
	float bk[nk];
incardon's avatar
incardon committed
134 135 136

	for (size_t i = 0 ; i < nk ; i++)
	{
137 138
	     ak[i] = rand()/RAND_MAX;
	     bk[i] = rand()/RAND_MAX;
incardon's avatar
incardon committed
139 140
	}

141
	// We calculate the circuitation gamma
incardon's avatar
incardon committed
142 143 144
	float gamma = nu * tgtre;
	float rinv2 = 1.0f/(sigma*sigma);
	float max_vorticity = gamma*rinv2/M_PI;
incardon's avatar
incardon committed
145

146
	// We go through the grid to initialize the vortex
incardon's avatar
incardon committed
147 148 149 150 151 152 153
	auto it = gr.getDomainIterator();

	while (it.isNext())
	{
		auto key_d = it.get();
		auto key = it.getGKey(key_d);

incardon's avatar
incardon committed
154 155 156
        float tx = (key.get(x)-2)*gr.spacing(x) + domain.getLow(x);
        float ty = (key.get(y)-2)*gr.spacing(y) + domain.getLow(y);
        float tz = (key.get(z)-2)*gr.spacing(z) + domain.getLow(z);
157
        float theta1 = atan2((ty-domain.getHigh(1)/2.0f),(tz-domain.getHigh(2)/2.0f));
incardon's avatar
incardon committed
158

159

incardon's avatar
incardon committed
160
        float rad1r  = sqrt((ty-domain.getHigh(1)/2.0f)*(ty-domain.getHigh(1)/2.0f) + (tz-domain.getHigh(2)/2.0f)*(tz-domain.getHigh(2)/2.0f)) - ringr1;
incardon's avatar
incardon committed
161 162 163
        float rad1t = tx - 1.0f;
        float rad1sq = rad1r*rad1r + rad1t*rad1t;
        float radstr = -exp(-rad1sq*rinv2)*rinv2*gamma/M_PI;
incardon's avatar
incardon committed
164 165 166 167
        gr.template get<vorticity>(key_d)[x] = 0.0f;
        gr.template get<vorticity>(key_d)[y] = -radstr * cos(theta1);
        gr.template get<vorticity>(key_d)[z] = radstr * sin(theta1);

incardon's avatar
incardon committed
168 169
        // kill the axis term

incardon's avatar
incardon committed
170
        float rad1r_  = sqrt((ty-domain.getHigh(1)/2.0f)*(ty-domain.getHigh(1)/2.0f) + (tz-domain.getHigh(2)/2.0f)*(tz-domain.getHigh(2)/2.0f)) + ringr1;
incardon's avatar
incardon committed
171 172 173 174 175
        float rad1sqTILDA = rad1sq*rinv2;
        radstr = exp(-rad1sq*rinv2)*rinv2*gamma/M_PI;
        gr.template get<vorticity>(key_d)[x] = 0.0f;
        gr.template get<vorticity>(key_d)[y] = -radstr * cos(theta1);
        gr.template get<vorticity>(key_d)[z] = radstr * sin(theta1);
incardon's avatar
incardon committed
176 177 178 179 180 181 182

		++it;
	}
}

struct poisson_nn_helm
{
incardon's avatar
incardon committed
183
		//! 3D Stystem
incardon's avatar
incardon committed
184
        static const unsigned int dims = 3;
incardon's avatar
incardon committed
185
        //! We are solving for \psi that is a scalar
incardon's avatar
incardon committed
186
        static const unsigned int nvar = 1;
incardon's avatar
incardon committed
187
        //! Boundary conditions
incardon's avatar
incardon committed
188
        static const bool boundary[];
incardon's avatar
incardon committed
189
        //! type of the spatial coordinates
incardon's avatar
incardon committed
190
        typedef float stype;
incardon's avatar
incardon committed
191
        //! grid that store \psi
192
        typedef grid_type_s b_grid;
incardon's avatar
incardon committed
193
        //! Sparse matrix used to sove the linear system (we use PETSC)
incardon's avatar
incardon committed
194
        typedef SparseMatrix<double,int,PETSC_BASE> SparseMatrix_type;
incardon's avatar
incardon committed
195
        //! Vector to solve the system (PETSC)
incardon's avatar
incardon committed
196
        typedef Vector<double,PETSC_BASE> Vector_type;
incardon's avatar
incardon committed
197
        //! It is a normal grid
incardon's avatar
incardon committed
198 199 200
        static const int grid_type = NORMAL_GRID;
};

incardon's avatar
incardon committed
201
//! boundary conditions are PERIODIC
incardon's avatar
incardon committed
202 203
const bool poisson_nn_helm::boundary[] = {PERIODIC,PERIODIC,PERIODIC};

204

205
void helmotz_hodge_projection(grid_dist_id<3,float,aggregate<float>,CartDecomposition<3,float,HeapMemory,SpaceDistribution<3,float>>> & psi,
incardon's avatar
incardon committed
206 207 208 209 210 211
							  FDScheme<poisson_nn_helm> & fd,
							  grid_type & gr,
							  const Box<3,float> & domain,
							  petsc_solver<double> & solver,
							  petsc_solver<double>::return_type & x_ ,
							  bool init)
incardon's avatar
incardon committed
212 213 214 215
{
	// ghost get
	gr.template ghost_get<vorticity>();

216
	// ghost size of the psi function
incardon's avatar
incardon committed
217 218 219 220 221 222 223 224 225
    Ghost<3,long int> g(2);

	// Calculate the divergence of the vortex field
	auto it = gr.getDomainIterator();	

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

226
        psi.template get<phi>(key) = 1.0f/gr.spacing(x)/2.0f*(gr.template get<vorticity>(key.move(x,1))[x] - gr.template get<vorticity>(key.move(x,-1))[x] ) +
incardon's avatar
incardon committed
227 228 229 230 231 232
			         	 	 1.0f/gr.spacing(y)/2.0f*(gr.template get<vorticity>(key.move(y,1))[y] - gr.template get<vorticity>(key.move(y,-1))[y] ) +
							 1.0f/gr.spacing(z)/2.0f*(gr.template get<vorticity>(key.move(z,1))[z] - gr.template get<vorticity>(key.move(z,-1))[z] );

		++it;
	}

233
	calc_and_print_max_div_and_int(gr);
incardon's avatar
incardon committed
234

235
	//! \cond [create b] \endcond
236

incardon's avatar
incardon committed
237 238
	fd.new_b();
	fd.template impose_dit_b<0>(psi,psi.getDomainIterator());
incardon's avatar
incardon committed
239

240
	//! \cond [create b] \endcond
incardon's avatar
incardon committed
241

242
	timer tm_solve;
incardon's avatar
incardon committed
243 244 245 246
	if (init == true)
	{
		// Set the conjugate-gradient as method to solve the linear solver
		solver.setSolver(KSPBCGS);
247

incardon's avatar
incardon committed
248
		// Set the absolute tolerance to determine that the method is converged
249
		solver.setAbsTol(0.0001);
250

incardon's avatar
incardon committed
251 252
		// Set the maximum number of iterations
		solver.setMaxIter(500);
incardon's avatar
incardon committed
253

254 255
#ifdef USE_MULTIGRID

256 257 258 259 260 261 262 263
        solver.setPreconditioner(PCHYPRE_BOOMERAMG);
        solver.setPreconditionerAMG_nl(6);
        solver.setPreconditionerAMG_maxit(1);
        solver.setPreconditionerAMG_relax("SOR/Jacobi");
        solver.setPreconditionerAMG_cycleType("V",0,4);
        solver.setPreconditionerAMG_coarsenNodalType(0);
        solver.setPreconditionerAMG_coarsen("HMIS");

264 265 266

#endif

incardon's avatar
incardon committed
267
		tm_solve.start();
incardon's avatar
incardon committed
268

incardon's avatar
incardon committed
269 270
		// Give to the solver A and b, return x, the solution
		solver.solve(fd.getA(),x_,fd.getB());
271

incardon's avatar
incardon committed
272 273 274 275
		tm_solve.stop();
	}
	else
	{
276
		tm_solve.start();
incardon's avatar
incardon committed
277
		solver.solve(x_,fd.getB());
278
		tm_solve.stop();
incardon's avatar
incardon committed
279
	}
incardon's avatar
incardon committed
280

281 282
	// copy the solution x to the grid psi
	fd.template copy<phi>(x_,psi);
incardon's avatar
incardon committed
283

284
	psi.template ghost_get<phi>();
incardon's avatar
incardon committed
285 286 287 288 289 290 291 292 293

	// Correct the vorticity to make it divergence free

	auto it2 = gr.getDomainIterator();

	while (it2.isNext())
	{
		auto key = it2.get();

294 295 296
		gr.template get<vorticity>(key)[x] -= 1.0f/2.0f/psi.spacing(x)*(psi.template get<phi>(key.move(x,1))-psi.template get<phi>(key.move(x,-1)));
		gr.template get<vorticity>(key)[y] -= 1.0f/2.0f/psi.spacing(y)*(psi.template get<phi>(key.move(y,1))-psi.template get<phi>(key.move(y,-1)));
		gr.template get<vorticity>(key)[z] -= 1.0f/2.0f/psi.spacing(z)*(psi.template get<phi>(key.move(z,1))-psi.template get<phi>(key.move(z,-1)));
incardon's avatar
incardon committed
297 298 299 300

		++it2;
	}

301
	calc_and_print_max_div_and_int(gr);
incardon's avatar
incardon committed
302 303
}

incardon's avatar
incardon committed
304
void remesh(particles_type & vd, grid_type & gr,Box<3,float> & domain)
305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339
{
	// Remove all particles because we reinitialize in a grid like position
	vd.clear();

	// Get a grid iterator
	auto git = vd.getGridIterator(gr.getGridInfo().getSize());

	// For each grid point
	while (git.isNext())
	{
		// Get the grid point in global coordinates (key). p is in local coordinates
		auto p = git.get();
		auto key = git.get_dist();

		// Add the particle
		vd.add();

		// Assign the position to the particle in a grid like way
		vd.getLastPos()[x] = gr.spacing(x)*p.get(x) + domain.getLow(x);
		vd.getLastPos()[y] = gr.spacing(y)*p.get(y) + domain.getLow(y);
		vd.getLastPos()[z] = gr.spacing(z)*p.get(z) + domain.getLow(z);

		// Initialize the vorticity
		vd.template getLastProp<vorticity>()[x] = gr.template get<vorticity>(key)[x];
		vd.template getLastProp<vorticity>()[y] = gr.template get<vorticity>(key)[y];
		vd.template getLastProp<vorticity>()[z] = gr.template get<vorticity>(key)[z];

		// next grid point
		++git;
	}

	// redistribute the particles
	vd.map();
}

incardon's avatar
incardon committed
340

341 342
void comp_vel(grid_type_s & gr_ps,
		      grid_type & phi_v,
incardon's avatar
incardon committed
343 344 345 346 347 348
			  FDScheme<poisson_nn_helm> & fd,
		      Box<3,float> & domain,
			  grid_type & g_vort,
			  grid_type & g_vel,
			  petsc_solver<double>::return_type (& phi_s)[3],
			  petsc_solver<double> & solver)
incardon's avatar
incardon committed
349
{
350 351 352
	// We calculate and print the maximum divergence of the vorticity
	// and the integral of it
	calc_and_print_max_div_and_int(g_vort);
353

incardon's avatar
incardon committed
354 355 356
	// For each component solve a poisson
	for (size_t i = 0 ; i < 3 ; i++)
	{
357 358
		//! \cond [solve_poisson_comp] \endcond

359
		// Copy the vorticity component i in gr_ps
incardon's avatar
incardon committed
360 361 362 363 364 365 366
		auto it2 = gr_ps.getDomainIterator();

		// calculate the velocity from the curl of phi
		while (it2.isNext())
		{
			auto key = it2.get();

367
			// copy
incardon's avatar
incardon committed
368 369 370 371 372
			gr_ps.get<vorticity>(key) = g_vort.template get<vorticity>(key)[i];

			++it2;
		}

373
		// impose the poisson equation using gr_ps = vorticity for the right-hand-side (on the full domain)
incardon's avatar
incardon committed
374 375
		fd.new_b();
		fd.template impose_dit_b<phi>(gr_ps,gr_ps.getDomainIterator());
incardon's avatar
incardon committed
376

377
		solver.setAbsTol(0.01);
incardon's avatar
incardon committed
378

379
		// Get the vector that represent the right-hand-side
incardon's avatar
incardon committed
380 381
		Vector<double,PETSC_BASE> & b = fd.getB();

incardon's avatar
incardon committed
382 383 384
		timer tm_solve;
		tm_solve.start();

385 386 387
		// Give to the solver A and b in this case we are giving
		// an intitial guess phi_s calculated in the previous
		// time step
incardon's avatar
incardon committed
388
		solver.solve(phi_s[i],b);
incardon's avatar
incardon committed
389

incardon's avatar
incardon committed
390 391
		tm_solve.stop();

incardon's avatar
incardon committed
392 393
		// Calculate the residual

incardon's avatar
incardon committed
394
		solError serr;
incardon's avatar
incardon committed
395
		serr = solver.get_residual_error(phi_s[i],b);
incardon's avatar
incardon committed
396

397
		Vcluster & v_cl = create_vcluster();
incardon's avatar
incardon committed
398 399
		if (v_cl.getProcessUnitID() == 0)
		{std::cout << "Solved component " << i << "  Error: " << serr.err_inf << std::endl;}
incardon's avatar
incardon committed
400 401

		// copy the solution to grid
402
		fd.template copy<phi>(phi_s[i],gr_ps);
incardon's avatar
incardon committed
403 404 405 406 407 408 409 410

		auto it3 = gr_ps.getDomainIterator();

		// calculate the velocity from the curl of phi
		while (it3.isNext())
		{
			auto key = it3.get();

411
			phi_v.get<velocity>(key)[i] = gr_ps.get<phi>(key);
incardon's avatar
incardon committed
412 413 414 415 416

			++it3;
		}
	}

417
	phi_v.ghost_get<phi>();
incardon's avatar
incardon committed
418

incardon's avatar
incardon committed
419 420 421
	float inv_dx = 0.5f / g_vort.spacing(0);
	float inv_dy = 0.5f / g_vort.spacing(1);
	float inv_dz = 0.5f / g_vort.spacing(2);
422 423

	auto it3 = phi_v.getDomainIterator();
incardon's avatar
incardon committed
424 425 426 427 428 429

	// calculate the velocity from the curl of phi
	while (it3.isNext())
	{
		auto key = it3.get();

incardon's avatar
incardon committed
430 431 432 433 434 435
		float phi_xy = (phi_v.get<phi>(key.move(y,1))[x] - phi_v.get<phi>(key.move(y,-1))[x])*inv_dy;
		float phi_xz = (phi_v.get<phi>(key.move(z,1))[x] - phi_v.get<phi>(key.move(z,-1))[x])*inv_dz;
		float phi_yx = (phi_v.get<phi>(key.move(x,1))[y] - phi_v.get<phi>(key.move(x,-1))[y])*inv_dx;
		float phi_yz = (phi_v.get<phi>(key.move(z,1))[y] - phi_v.get<phi>(key.move(z,-1))[y])*inv_dz;
		float phi_zx = (phi_v.get<phi>(key.move(x,1))[z] - phi_v.get<phi>(key.move(x,-1))[z])*inv_dx;
		float phi_zy = (phi_v.get<phi>(key.move(y,1))[z] - phi_v.get<phi>(key.move(y,-1))[z])*inv_dy;
incardon's avatar
incardon committed
436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483

		g_vel.template get<velocity>(key)[x] = phi_zy - phi_yz;
		g_vel.template get<velocity>(key)[y] = phi_xz - phi_zx;
		g_vel.template get<velocity>(key)[z] = phi_yx - phi_xy;

		++it3;
	}
}


template<unsigned int prp> void set_zero(grid_type & gr)
{
	auto it = gr.getDomainGhostIterator();

	// calculate the velocity from the curl of phi
	while (it.isNext())
	{
		auto key = it.get();

		gr.template get<prp>(key)[0] = 0.0;
		gr.template get<prp>(key)[1] = 0.0;
		gr.template get<prp>(key)[2] = 0.0;

		++it;
	}
}


template<unsigned int prp> void set_zero(particles_type & vd)
{
	auto it = vd.getDomainIterator();

	// calculate the velocity from the curl of phi
	while (it.isNext())
	{
		auto key = it.get();

		vd.template getProp<prp>(key)[0] = 0.0;
		vd.template getProp<prp>(key)[1] = 0.0;
		vd.template getProp<prp>(key)[2] = 0.0;

		++it;
	}
}


template<typename grid> void calc_rhs(grid & g_vort, grid & g_vel, grid & g_dwp)
{
484
	// usefull constant
incardon's avatar
incardon committed
485 486
	constexpr int rhs = 0;

487 488
	// calculate several pre-factors for the stencil finite
	// difference
incardon's avatar
incardon committed
489 490 491
	float fac1 = 1.0f*nu/(g_vort.spacing(0)*g_vort.spacing(0));
	float fac2 = 1.0f*nu/(g_vort.spacing(1)*g_vort.spacing(1));
	float fac3 = 1.0f*nu/(g_vort.spacing(2)*g_vort.spacing(2));
incardon's avatar
incardon committed
492

incardon's avatar
incardon committed
493 494 495
	float fac4 = 0.5f/(g_vort.spacing(0));
	float fac5 = 0.5f/(g_vort.spacing(1));
	float fac6 = 0.5f/(g_vort.spacing(2));
incardon's avatar
incardon committed
496 497 498 499 500 501 502 503 504

	auto it = g_dwp.getDomainIterator();

	g_vort.template ghost_get<vorticity>();

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

505
		g_dwp.template get<rhs>(key)[x] = fac1*(g_vort.template get<vorticity>(key.move(x,1))[x]+g_vort.template get<vorticity>(key.move(x,-1))[x])+
incardon's avatar
incardon committed
506 507
				                          fac2*(g_vort.template get<vorticity>(key.move(y,1))[x]+g_vort.template get<vorticity>(key.move(y,-1))[x])+
										  fac3*(g_vort.template get<vorticity>(key.move(z,1))[x]+g_vort.template get<vorticity>(key.move(z,-1))[x])-
508
										  2.0f*(fac1+fac2+fac3)*g_vort.template get<vorticity>(key)[x]+
incardon's avatar
incardon committed
509 510 511 512 513 514 515
										  fac4*g_vort.template get<vorticity>(key)[x]*
										  (g_vel.template get<velocity>(key.move(x,1))[x] - g_vel.template get<velocity>(key.move(x,-1))[x])+
										  fac5*g_vort.template get<vorticity>(key)[y]*
										  (g_vel.template get<velocity>(key.move(y,1))[x] - g_vel.template get<velocity>(key.move(y,-1))[x])+
										  fac6*g_vort.template get<vorticity>(key)[z]*
										  (g_vel.template get<velocity>(key.move(z,1))[x] - g_vel.template get<velocity>(key.move(z,-1))[x]);

516
		g_dwp.template get<rhs>(key)[y] = fac1*(g_vort.template get<vorticity>(key.move(x,1))[y]+g_vort.template get<vorticity>(key.move(x,-1))[y])+
incardon's avatar
incardon committed
517 518
				                          fac2*(g_vort.template get<vorticity>(key.move(y,1))[y]+g_vort.template get<vorticity>(key.move(y,-1))[y])+
										  fac3*(g_vort.template get<vorticity>(key.move(z,1))[y]+g_vort.template get<vorticity>(key.move(z,-1))[y])-
519
										  2.0f*(fac1+fac2+fac3)*g_vort.template get<vorticity>(key)[y]+
incardon's avatar
incardon committed
520 521 522 523 524 525 526 527
										  fac4*g_vort.template get<vorticity>(key)[x]*
										  (g_vel.template get<velocity>(key.move(x,1))[y] - g_vel.template get<velocity>(key.move(x,-1))[y])+
										  fac5*g_vort.template get<vorticity>(key)[y]*
										  (g_vel.template get<velocity>(key.move(y,1))[y] - g_vel.template get<velocity>(key.move(y,-1))[y])+
										  fac6*g_vort.template get<vorticity>(key)[z]*
										  (g_vel.template get<velocity>(key.move(z,1))[y] - g_vel.template get<velocity>(key.move(z,-1))[y]);


528
		g_dwp.template get<rhs>(key)[z] = fac1*(g_vort.template get<vorticity>(key.move(x,1))[z]+g_vort.template get<vorticity>(key.move(x,-1))[z])+
incardon's avatar
incardon committed
529 530
				                          fac2*(g_vort.template get<vorticity>(key.move(y,1))[z]+g_vort.template get<vorticity>(key.move(y,-1))[z])+
										  fac3*(g_vort.template get<vorticity>(key.move(z,1))[z]+g_vort.template get<vorticity>(key.move(z,-1))[z])-
531
										  2.0f*(fac1+fac2+fac3)*g_vort.template get<vorticity>(key)[z]+
incardon's avatar
incardon committed
532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584
										  fac4*g_vort.template get<vorticity>(key)[x]*
										  (g_vel.template get<velocity>(key.move(x,1))[z] - g_vel.template get<velocity>(key.move(x,-1))[z])+
										  fac5*g_vort.template get<vorticity>(key)[y]*
										  (g_vel.template get<velocity>(key.move(y,1))[z] - g_vel.template get<velocity>(key.move(y,-1))[z])+
										  fac6*g_vort.template get<vorticity>(key)[z]*
										  (g_vel.template get<velocity>(key.move(z,1))[z] - g_vel.template get<velocity>(key.move(z,-1))[z]);

		++it;
	}
}


void rk_step1(particles_type & particles)
{
	auto it = particles.getDomainIterator();

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

		// Save the old vorticity
		particles.getProp<old_vort>(key)[x] = particles.getProp<vorticity>(key)[x];
		particles.getProp<old_vort>(key)[y] = particles.getProp<vorticity>(key)[y];
		particles.getProp<old_vort>(key)[z] = particles.getProp<vorticity>(key)[z];

		particles.getProp<vorticity>(key)[x] = particles.getProp<vorticity>(key)[x] + 0.5f * dt * particles.getProp<rhs_part>(key)[x];
		particles.getProp<vorticity>(key)[y] = particles.getProp<vorticity>(key)[y] + 0.5f * dt * particles.getProp<rhs_part>(key)[y];
		particles.getProp<vorticity>(key)[z] = particles.getProp<vorticity>(key)[z] + 0.5f * dt * particles.getProp<rhs_part>(key)[z];

		++it;
	}

	auto it2 = particles.getDomainIterator();

	while (it2.isNext())
	{
		auto key = it2.get();

		// Save the old position
		particles.getProp<old_pos>(key)[x] = particles.getPos(key)[x];
		particles.getProp<old_pos>(key)[y] = particles.getPos(key)[y];
		particles.getProp<old_pos>(key)[z] = particles.getPos(key)[z];

		particles.getPos(key)[x] = particles.getPos(key)[x] + 0.5f * dt * particles.getProp<p_vel>(key)[x];
		particles.getPos(key)[y] = particles.getPos(key)[y] + 0.5f * dt * particles.getProp<p_vel>(key)[y];
		particles.getPos(key)[z] = particles.getPos(key)[z] + 0.5f * dt * particles.getProp<p_vel>(key)[z];

		++it2;
	}

	particles.map();
}

585

incardon's avatar
incardon committed
586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616
void rk_step2(particles_type & particles)
{
	auto it = particles.getDomainIterator();

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

		particles.getProp<vorticity>(key)[x] = particles.getProp<old_vort>(key)[x] + dt * particles.getProp<rhs_part>(key)[x];
		particles.getProp<vorticity>(key)[y] = particles.getProp<old_vort>(key)[y] + dt * particles.getProp<rhs_part>(key)[y];
		particles.getProp<vorticity>(key)[z] = particles.getProp<old_vort>(key)[z] + dt * particles.getProp<rhs_part>(key)[z];

		++it;
	}

	auto it2 = particles.getDomainIterator();

	while (it2.isNext())
	{
		auto key = it2.get();

		particles.getPos(key)[x] = particles.getProp<old_pos>(key)[x] + dt * particles.getProp<p_vel>(key)[x];
		particles.getPos(key)[y] = particles.getProp<old_pos>(key)[y] + dt * particles.getProp<p_vel>(key)[y];
		particles.getPos(key)[z] = particles.getProp<old_pos>(key)[z] + dt * particles.getProp<p_vel>(key)[z];

		++it2;
	}

	particles.map();
}

617 618
template<typename grid, typename vector> void do_step(grid_type_s & psi,
													  grid_type & phi_v,
incardon's avatar
incardon committed
619 620
													  FDScheme<poisson_nn_helm> & fd,
													  vector & particles,
incardon's avatar
incardon committed
621 622 623
		                                              grid & g_vort,
													  grid & g_vel,
													  grid & g_dvort,
incardon's avatar
incardon committed
624 625
													  Box<3,float> & domain,
													  interpolate<particles_type,grid_type,mp4_kernel<float>> & inte,
incardon's avatar
incardon committed
626 627
													  petsc_solver<double>::return_type (& phi_s)[3],
													  petsc_solver<double> & solver)
incardon's avatar
incardon committed
628 629 630 631 632
{
	constexpr int rhs = 0;

	set_zero<vorticity>(g_vort);
	inte.template p2m<vorticity,vorticity>(particles,g_vort);
incardon's avatar
incardon committed
633

incardon's avatar
incardon committed
634 635 636
	g_vort.template ghost_put<add_,vorticity>();

	// Calculate velocity from vorticity
incardon's avatar
incardon committed
637
	comp_vel(psi,phi_v,fd,domain,g_vort,g_vel,phi_s,solver);
incardon's avatar
incardon committed
638 639 640 641 642 643 644 645
	calc_rhs(g_vort,g_vel,g_dvort);

	g_dvort.template ghost_get<rhs>();
	g_vel.template ghost_get<velocity>();
	set_zero<rhs_part>(particles);
	set_zero<p_vel>(particles);
	inte.template m2p<rhs,rhs_part>(g_dvort,particles);
	inte.template m2p<velocity,p_vel>(g_vel,particles);
646 647

	//! \cond [inte_m2p] \endcond
incardon's avatar
incardon committed
648 649
}

650 651 652 653 654 655 656 657 658 659 660 661 662
template<typename vector, typename grid> void check_point_and_save(vector & particles,
																   grid & g_vort,
																   grid & g_vel,
																   grid & g_dvort,
																   size_t i)
{
	Vcluster & v_cl = create_vcluster();

	if (v_cl.getProcessingUnits() < 24)
	{
		particles.write("part_out_" + std::to_string(i),VTK_WRITER | FORMAT_BINARY);
		g_vort.template ghost_get<vorticity>();
		g_vel.template ghost_get<velocity>();
incardon's avatar
incardon committed
663 664
		g_vel.write_frame("grid_velocity",i, VTK_WRITER | FORMAT_BINARY);
		g_vort.write_frame("grid_vorticity",i, VTK_WRITER | FORMAT_BINARY);
665 666 667 668
	}

	// In order to reduce the size of the saved data we apply a threshold.
	// We only save particles with vorticity higher than 0.1
669
	particles_type_s part_save(particles.getDecomposition(),0);
670 671 672 673 674 675 676

	auto it_s = particles.getDomainIterator();

	while (it_s.isNext())
	{
		auto p = it_s.get();

incardon's avatar
incardon committed
677
		float vort_magn = sqrt(particles.template getProp<vorticity>(p)[0] * particles.template getProp<vorticity>(p)[0] +
678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702
							   particles.template getProp<vorticity>(p)[1] * particles.template getProp<vorticity>(p)[1] +
							   particles.template getProp<vorticity>(p)[2] * particles.template getProp<vorticity>(p)[2]);

		if (vort_magn < 0.1)
		{
			++it_s;
			continue;
		}

		part_save.add();

		part_save.getLastPos()[0] = particles.getPos(p)[0];
		part_save.getLastPos()[1] = particles.getPos(p)[1];
		part_save.getLastPos()[2] = particles.getPos(p)[2];

		part_save.template getLastProp<0>() = vort_magn;

		++it_s;
	}

	part_save.map();

	particles.save("check_point");
}

incardon's avatar
incardon committed
703 704 705 706
int main(int argc, char* argv[])
{
	// Initialize
	openfpm_init(&argc,&argv);
707
	{
incardon's avatar
incardon committed
708
	// Domain, a rectangle
incardon's avatar
incardon committed
709 710
	// For the grid 1600x400x400 use
	// Box<3,float> domain({0.0,0.0,0.0},{22.0,5.57,5.57});
711
	Box<3,float> domain({0.0,0.0,0.0},{3.57,3.57,3.57});
incardon's avatar
incardon committed
712 713 714 715 716

	// Ghost (Not important in this case but required)
	Ghost<3,long int> g(2);

	// Grid points on x=128 y=64 z=64
incardon's avatar
incardon committed
717 718
	// if we use Re = 7500
	//  long int sz[] = {1600,400,400};
incardon's avatar
incardon committed
719
	long int sz[] = {128,128,128};
incardon's avatar
incardon committed
720 721 722 723
	size_t szu[] = {(size_t)sz[0],(size_t)sz[1],(size_t)sz[2]};

	periodicity<3> bc = {{PERIODIC,PERIODIC,PERIODIC}};

724 725
	//! \cond [construct grids] \endcond

incardon's avatar
incardon committed
726 727 728 729 730
	grid_type g_vort(szu,domain,g,bc);
	grid_type g_vel(g_vort.getDecomposition(),szu,g);
	grid_type g_dvort(g_vort.getDecomposition(),szu,g);
	particles_type particles(g_vort.getDecomposition(),0);

incardon's avatar
incardon committed
731 732 733 734 735 736
	// Construct an FDScheme is heavy so we construct it here
	// We also construct distributed temporal grids here because
	// they are expensive

	// Here we create temporal distributed grid, create a distributed grid is expensive we do it once outside
	// And the vectorial phi_v
737 738
	grid_type_s psi(g_vort.getDecomposition(),g_vort.getGridInfo().getSize(),g);
	grid_type phi_v(g_vort.getDecomposition(),g_vort.getGridInfo().getSize(),g);
incardon's avatar
incardon committed
739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757

	// In order to create a matrix that represent the poisson equation we have to indicate
	// we have to indicate the maximum extension of the stencil and we we need an extra layer
	// of points in case we have to impose particular boundary conditions.
	Ghost<3,long int> stencil_max(2);

	// We define a field phi_f
	typedef Field<phi,poisson_nn_helm> phi_f;

	// We assemble the poisson equation doing the
	// poisson of the Laplacian of phi using Laplacian
	// central scheme (where the both derivative use central scheme, not
	// forward composed backward like usually)
	typedef Lap<phi_f,poisson_nn_helm,CENTRAL_SYM> poisson;

	FDScheme<poisson_nn_helm> fd(stencil_max, domain, psi);

	fd.template impose_dit<0>(poisson(),psi,psi.getDomainIterator());

758
	//! \cond [construct grids] \endcond
incardon's avatar
incardon committed
759

760 761 762
	// It store the solution to compute velocity
	// It is used as initial guess every time we call the solver
	Vector<double,PETSC_BASE> phi_s[3];
incardon's avatar
incardon committed
763
	Vector<double,PETSC_BASE> x_;
764 765

	// Parallel object
incardon's avatar
incardon committed
766 767
	Vcluster & v_cl = create_vcluster();

768
	// print out the spacing of the grid
incardon's avatar
incardon committed
769 770
	if (v_cl.getProcessUnitID() == 0)
	{std::cout << "SPACING " << g_vort.spacing(0) << " " << g_vort.spacing(1) << " " << g_vort.spacing(2) << std::endl;}
incardon's avatar
incardon committed
771

772
	// initialize the ring step 1
incardon's avatar
incardon committed
773
	init_ring(g_vort,domain);
774

incardon's avatar
incardon committed
775 776 777 778 779 780
	x_.resize(g_vort.size(),g_vort.getLocalDomainSize());
	x_.setZero();

	// Create a PETSC solver to get the solution x
	petsc_solver<double> solver;

781
	// Do the helmotz projection step 2
incardon's avatar
incardon committed
782
	helmotz_hodge_projection(psi,fd,g_vort,domain,solver,x_,true);
incardon's avatar
incardon committed
783

784
	// initialize the particle on a mesh step 3
incardon's avatar
incardon committed
785 786
	remesh(particles,g_vort,domain);

787
	// Get the total number of particles
incardon's avatar
incardon committed
788 789 790 791
	size_t tot_part = particles.size_local();
	v_cl.sum(tot_part);
	v_cl.execute();

792
	// Now we set the size of phi_s
793 794 795
	phi_s[0].resize(g_vort.size(),g_vort.getLocalDomainSize());
	phi_s[1].resize(g_vort.size(),g_vort.getLocalDomainSize());
	phi_s[2].resize(g_vort.size(),g_vort.getLocalDomainSize());
incardon's avatar
incardon committed
796

797
	// And we initially set it to zero
incardon's avatar
incardon committed
798 799 800 801
	phi_s[0].setZero();
	phi_s[1].setZero();
	phi_s[2].setZero();

802 803 804 805 806
	// We create the interpolation object outside the loop cycles
	// to avoid to recreate it at every cycle. Internally this object
	// create some data-structure to optimize the conversion particle
	// position to sub-domain. If there are no re-balancing it is safe
	// to reuse-it
incardon's avatar
incardon committed
807
	interpolate<particles_type,grid_type,mp4_kernel<float>> inte(particles,g_vort);
incardon's avatar
incardon committed
808

incardon's avatar
incardon committed
809 810 811
	// With more than 24 core we rely on the HDF5 checkpoint restart file
	if (v_cl.getProcessingUnits() < 24)
	{g_vort.write("grid_vorticity_init", VTK_WRITER | FORMAT_BINARY);}
incardon's avatar
incardon committed
812

813

814 815
	// Before entering the loop we check if we want to restart from
	// a previous simulation
816 817
	size_t i = 0;

818
	// if we have an argument
819 820
	if (argc > 1)
	{
821
		// take that argument
822 823
		std::string restart(argv[1]);

824
		// convert it into number
825 826
		i = std::stoi(restart);

827
		// load the file
828 829
		particles.load("check_point_" + std::to_string(i));

830 831
		// Print to inform that we are restarting from a
		// particular position
incardon's avatar
incardon committed
832 833
		if (v_cl.getProcessUnitID() == 0)
		{std::cout << "Restarting from " << i << std::endl;}
834 835
	}

incardon's avatar
incardon committed
836
	// Time Integration
incardon's avatar
incardon committed
837
	for ( ; i < nsteps ; i++)
incardon's avatar
incardon committed
838
	{
839
		// do step 4-5-6-7
incardon's avatar
incardon committed
840
		do_step(psi,phi_v,fd,particles,g_vort,g_vel,g_dvort,domain,inte,phi_s,solver);
incardon's avatar
incardon committed
841

842
		// do step 8
incardon's avatar
incardon committed
843 844
		rk_step1(particles);

845
		// do step 9-10-11-12
incardon's avatar
incardon committed
846
		do_step(psi,phi_v,fd,particles,g_vort,g_vel,g_dvort,domain,inte,phi_s,solver);
incardon's avatar
incardon committed
847

848
		// do step 13
incardon's avatar
incardon committed
849 850
		rk_step2(particles);

851
		// so step 14
incardon's avatar
incardon committed
852 853 854 855
		set_zero<vorticity>(g_vort);
		inte.template p2m<vorticity,vorticity>(particles,g_vort);
		g_vort.template ghost_put<add_,vorticity>();

856
		// helmotz-hodge projection
incardon's avatar
incardon committed
857
		helmotz_hodge_projection(psi,fd,g_vort,domain,solver,x_,false);
858

incardon's avatar
incardon committed
859 860
		remesh(particles,g_vort,domain);

861 862 863
		// print the step number
		if (v_cl.getProcessUnitID() == 0)
		{std::cout << "Step " << i << std::endl;}
incardon's avatar
incardon committed
864

865 866
		// every 100 steps write the output
		if (i % 100 == 0)		{check_point_and_save(particles,g_vort,g_vel,g_dvort,i);}
867

incardon's avatar
incardon committed
868
	}
869
	}
incardon's avatar
incardon committed
870 871 872 873

	openfpm_finalize();
}