diff --git a/example/Numerics/Vortex_in_cell/main_vic_petsc.cpp b/example/Numerics/Vortex_in_cell/main_vic_petsc.cpp
index c26dce08fee1fd51a8b2a6d86bb02bbf9de3d07f..495ab6cf4f6039bac31c9c5bb0f7df82cbd3eb23 100644
--- a/example/Numerics/Vortex_in_cell/main_vic_petsc.cpp
+++ b/example/Numerics/Vortex_in_cell/main_vic_petsc.cpp
@@ -3,16 +3,16 @@
* # Vortex in Cell 3D ring with ringlets # {#vic_ringlets}
*
* In this example we solve the Navier-Stokes equation in the vortex formulation in 3D
- * for an incompressible fluid.
+ * for an incompressible fluid. (bold symbols are vectorial quantity)
*
* ## Numerical method ## {#num_vic_mt}
*
* In this code we solve the Navier-stokes equation for incompressible fluid in the
* vorticity formulation. We first recall the Navier-stokes equation in vorticity formulation
*
- * \f$ \nabla \times u = -w \f$
+ * \f$ \nabla \times \boldsymbol u = - \boldsymbol w \f$
*
- * \f$ \frac{Dw}{dt} = (w \cdot \vec \nabla) u + \nu \nabla^{2} w \f$ (5)
+ * \f$ \frac{\displaystyle D \boldsymbol w}{\displaystyle dt} = ( \boldsymbol w \cdot \vec \nabla) \boldsymbol u + \nu \nabla^{2} \boldsymbol w \f$ (5)
*
* Where \f$w\f$ is the vorticity and \f$u\f$ is the velocity of the fluid.
* With high Reynold number \f$Re = \frac{uL}{\nu}\f$ and the term \f$uL\f$ significantly
@@ -26,17 +26,18 @@
3) Initialize particles on the same position as the grid or remesh
while (t < t_end) do
- 4) Interpolate mesh vorticity on particles
+ 4) Interpolate vorticity from the particles to mesh
5) calculate velocity u from the vorticity w
- 6) interpolate velocity u to particles
- 7) calculate the right-hand-side on grid and interpolate on particles
+ 6) calculate the right-hand-side on grid and interpolate on particles
+ 7) interpolate velocity u to particles
8) move particles accordingly to the velocity
- 9) interpolate vorticity w back to grid
+ 9) interpolate the vorticity into mesh and reinitialize the particles
+ in a grid like position
end while
\endverbatim
*
- * This pseudo code show how to solve the equation above using euler integration
+ * This pseudo code show how to solve the equation above using euler integration.
* In case of Runge-kutta of order two the pseudo code change into
*
*
@@ -47,18 +48,19 @@
3) Initialize particles on the same position as the grid or remesh
while (t < t_end) do
- 4) Interpolate mesh vorticity on particles
+ 4) 4) Interpolate vorticity from the particles to mesh
5) calculate velocity u from the vorticity w
- 6) interpolate velocity u to particles
- 7) calculate the right-hand-side on grid and interpolate on particles
+ 6) calculate the right-hand-side on grid and interpolate on particles
+ 7) interpolate velocity u to particles
8) move particles accordingly to the velocity and save the old position in x_old
- 9) interpolate vorticity w back to grid
- 10) Interpolate mesh vorticity on particles
- 11) calculate velocity u from the vorticity w
+ 9) Interpolate vorticity on mesh on the particles
+ 10) calculate velocity u from the vorticity w
+ 11) calculate the right-hand-side on grid and interpolate on particles
12) interpolate velocity u to particles
- 13) calculate the right-hand-side on grid and interpolate on particles
- 14) move particles accordingly to the velocity starting from x_old
+ 13) move particles accordingly to the velocity starting from x_old
+ 14) interpolate the vorticity into mesh and reinitialize the particles
+ in a grid like position
end while
\endverbatim
@@ -68,16 +70,17 @@
* ## Inclusion ## {#num_vic_inc}
*
* This example code need several components. First because is a particle
- * mesh example we have to activate "grid_dist_id.hpp" and "vector_dist_id.hpp".
+ * mesh example we have to activate **grid_dist_id.hpp** and **vector_dist_id.hpp**.
* Because we use a finite-difference scheme and linear-algebra to calculate the
- * velocity out of the vorticity, we have to include "FDScheme.hpp" to produce
- * out of the finite difference scheme a matrix that represent the linear-system
- * to solve. "SparseMatrix.hpp" is the Sparse-Matrix that contain the linear
- * system. "Vector.hpp" is the data-structure that contain the solution of the
- * linear system. "petsc_solver.hpp" is the library to invert the linear system.
+ * velocity out of the vorticity, we have to include **FDScheme.hpp** to produce
+ * from the finite difference scheme a matrix that represent the linear-system
+ * to solve. **SparseMatrix.hpp** is the Sparse-Matrix that will contain the linear
+ * system to solve in order to get the velocity out of the vorticity.
+ * **Vector.hpp** is the data-structure that contain the solution of the
+ * linear system. **petsc_solver.hpp** is the library to use in order invert the linear system.
* Because we have to interpolate between particles and grid we the to include
- * "interpolate.hpp" as interpolation kernel we use the mp4, so we include the
- * "mp4_kernel.hpp"
+ * **interpolate.hpp** as interpolation kernel we use the mp4, so we include the
+ * **mp4_kernel.hpp**
*
* For convenience we also define the particles type and the grid type and some
* convenient constants
@@ -106,9 +109,10 @@ constexpr int x = 0;
constexpr int y = 1;
constexpr int z = 2;
-////////////////// Simulation Parameters //////////////////////
-
+// The type of the grids
typedef grid_dist_id<3,float,aggregate<float[3]>> grid_type;
+
+// The type of the particles
typedef vector_dist<3,float,aggregate<float[3],float[3],float[3],float[3],float[3]>> particles_type;
// radius of the torus
@@ -127,6 +131,7 @@ float nu = 0.0001535;
float dt = 0.006;
+// All the properties by index
constexpr unsigned int vorticity = 0;
constexpr unsigned int velocity = 0;
constexpr unsigned int p_vel = 1;
@@ -136,6 +141,50 @@ constexpr unsigned int old_pos = 4;
//! \cond [inclusion] \endcond
+template<typename grid> void calc_and_print_max_div_and_int(grid & g_vort)
+{
+ //! \cond [sanity_int_div] \endcond
+
+ 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;}
+
+ //! \cond [sanity_int_div] \endcond
+}
+
/*! \page Vortex_in_cell_petsc Vortex in Cell 3D
*
* # Step 1: Initialization of the vortex ring # {#vic_ring_init}
@@ -302,21 +351,17 @@ const bool poisson_nn_helm::boundary[] = {PERIODIC,PERIODIC,PERIODIC};
*
* \snippet Numerics/Vortex_in_cell/main_vic_petsc.cpp poisson_syseq
*
- * Once created this object we can define the equation we are trying to solve
- * the left-hand-side of the equation \f$ \nabla^{2} \psi \f$
+ * Once created this object we can define the equation we are trying to solve.
+ * In particular the code below define the left-hand-side of the equation \f$ \nabla^{2} \psi \f$
*
* \snippet Numerics/Vortex_in_cell/main_vic_petsc.cpp poisson_obj_eq
*
* Before to construct the linear system we also calculate the divergence of the
* vorticity \f$ \nabla \cdot w \f$ that will be the right-hand-side
+ * of the equation
*
* \snippet Numerics/Vortex_in_cell/main_vic_petsc.cpp calc_div_vort
*
- * For statistic purpose we also calculate the integral of the vorticity and
- * the maximum divergence of the vorticity
- *
- * \snippet Numerics/Vortex_in_cell/main_vic_petsc.cpp sanity_div_integral
- *
* Finally we can create the object FDScheme using the object **poisson_nn_helm**
* as template variable. In addition to the constructor we have to specify the maximum extension of the stencil, the domain and the
* grid that will store the result. At this point we can impose an equation to construct
@@ -334,17 +379,17 @@ const bool poisson_nn_helm::boundary[] = {PERIODIC,PERIODIC,PERIODIC};
* the method go into infinite loop. After we set the parameters of the solver we can the
* the solution **x**. Finaly we copy back the solution **x** into the grid \f$ \psi \f$.
*
- * \snippet Numerics/Vortex_in_cell/main_vic_petsc.cpp div_free_check
+ * \snippet Numerics/Vortex_in_cell/main_vic_petsc.cpp solve_petsc
*
* ### Note ###
*
* Because we are solving the poisson equation in periodic boundary conditions the Matrix has
* determinat equal to zero. This mean that \f$ \psi \f$ has no unique solution (if it has one).
- * In order to recover one we have to ensure that the integral of the righ hand side or vorticity
+ * In order to recover one, we have to ensure that the integral of the righ hand side or vorticity
* is zero. (In our case is the case). We have to ensure that across time the integral of the
* vorticity is conserved. (In our case is the case if we consider the \f$ \nu = 0 \f$ and \f$
* \nabla \cdot w = 0 \f$ we can rewrite (5) in a conservative way \f$ \frac{Dw}{dt} = div(w \otimes v) \f$ ).
- * Is also good to notice that the solution that you get is the one with \f$ \int v = 0 \f$
+ * Is also good to notice that the solution that you get is the one with \f$ \int w = 0 \f$
*
*
*
@@ -416,52 +461,9 @@ void helmotz_hodge_projection(grid_type & gr, const Box<3,float> & domain)
//! \cond [calc_div_vort] \endcond
- //! \cond [sanity_div_integral] \endcond
-
- // Get iterator across the domain points
- auto it5 = gr.getDomainIterator();
-
- // maximum vorticity
- double max_vort = 0.0;
-
- // Total integral
- double int_vort[3] = {0.0,0.0,0.0};
-
- // For each grid point
- while (it5.isNext())
- {
- auto key = it5.get();
-
- double tmp;
-
- // Calculate the vorticity
- tmp = 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] ) +
- 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] );
-
- // Calculate the integral of the vorticity on each direction
- int_vort[x] += gr.template get<vorticity>(key)[x];
- int_vort[y] += gr.template get<vorticity>(key)[y];
- int_vort[z] += gr.template get<vorticity>(key)[z];
-
- // Store the maximum vorticity
- 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;}
+ calc_and_print_max_div_and_int(gr);
- //! \cond [sanity_div_integral] \endcond
//! \cond [create_fdscheme] \endcond
@@ -523,46 +525,15 @@ void helmotz_hodge_projection(grid_type & gr, const Box<3,float> & domain)
//! \cond [vort_correction] \endcond
- //! \cond [div_free_check] \endcond
-
- // Here we check the maximum divergence of the ring
- // We also copy the solution to the original grid
- double max = 0.0;
-
- gr.template ghost_get<vorticity>();
- auto it3 = gr.getDomainIterator();
-
- while (it3.isNext())
- {
- auto key = it3.get();
-
- 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] ) +
- 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] );
-
-
-
- if (psi.template get<phi>(key) > max)
- max = psi.template get<phi>(key);
-
- ++it3;
- }
-
- v_cl.max(max);
- v_cl.execute();
-
- if (v_cl.getProcessUnitID() == 0)
- {std::cout << "Maximum divergence of the ring MAX " << max << std::endl;}
-
- //! \cond [div_free_check] \endcond
-
+ calc_and_print_max_div_and_int(gr);
}
+
/*! \page Vortex_in_cell_petsc Vortex in Cell 3D
*
* # Step 3: Remeshing vorticity # {#vic_remesh_vort}
*
- * After that we initialized the vorticity on the grid, initialize the particles
+ * After that we initialized the vorticity on the grid, we initialize the particles
* in a grid like position and we interpolate the vorticity on particles. Because
* of the particles position being in a grid-like position and the symmetry of the
* interpolation kernels, the re-mesh step simply reduce to initialize the particle
@@ -642,34 +613,29 @@ void remesh(particles_type & vd, grid_type & gr,Box<3,float> & domain)
*
* \snippet Numerics/Vortex_in_cell/main_vic_petsc.cpp poisson_obj_eq
*
- * For statistic reason we do a check of the total integral of the vorticity
- * and maximum divergence of the vorticity
- *
- * \snippet Numerics/Vortex_in_cell/main_vic_petsc.cpp sanity_int_div
- *
* In order to calculate the velocity out of the vorticity, we solve a poisson
* equation like we did in helmotz-projection equation, but we do it for each
- * component \f$ i \f$ of the vorticity.
+ * component \f$ i \f$ of the vorticity. Qnce we have the solution in **psi_s**
+ * we copy the result back into the grid **gr_ps**. We than calculate the
+ * quality of the solution printing the norm infinity of the residual and
+ * finally we save in the grid vector vield **phi_v** the compinent \f$ i \f$
+ * (Copy from phi_s to phi_v is necessary because in phi_s is not a grid
+ * and cannot be used as a grid like object)
*
* \snippet Numerics/Vortex_in_cell/main_vic_petsc.cpp solve_poisson_comp
*
- * We save the \f$ \phi \f$ component into **phi_s**
+ * We save the component \f$ i \f$ of \f$ \phi \f$ into **phi_v**
*
- * \snippet Numerics/Vortex_in_cell/main_vic_petsc.cpp solve_poisson_comp
+ * \snippet Numerics/Vortex_in_cell/main_vic_petsc.cpp copy_to_phi_v
*
- * Once we filled phi_s we can implement (7) doing the curl of **phi_s**
- * and recovering the velocity v
+ * Once we filled phi_v we can implement (7) and calculate the curl of **phi_v**
+ * to recover the velocity v
*
* \snippet Numerics/Vortex_in_cell/main_vic_petsc.cpp curl_phi_v
*
- * \snippet Numerics/Vortex_in_cell/main_vic_petsc.cpp solve_poisson_comp
- *
- * \snippet Numerics/Vortex_in_cell/main_vic_petsc.cpp cross_check_curl_u_vort
*
*/
-//! \cond [comp_vel] \endcond
-
void comp_vel(Box<3,float> & domain, grid_type & g_vort,grid_type & g_vel, petsc_solver<double>::return_type (& phi_s)[3])
{
//! \cond [poisson_obj_eq] \endcond
@@ -688,61 +654,23 @@ void comp_vel(Box<3,float> & domain, grid_type & g_vort,grid_type & g_vel, petsc
//! \cond [poisson_obj_eq] \endcond
+ // Maximum size of the stencil
Ghost<3,long int> g(2);
- // Here we create a distributed grid to store the result of the helmotz projection
+ // Here we create a distributed grid to store the result
grid_dist_id<3,float,aggregate<float>> gr_ps(g_vort.getDecomposition(),g_vort.getGridInfo().getSize(),g);
grid_dist_id<3,float,aggregate<float[3]>> phi_v(g_vort.getDecomposition(),g_vort.getGridInfo().getSize(),g);
- // Check one we control the the divergence of vorticity is zero
-
- //! \cond [sanity_int_div] \endcond
-
- 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;}
-
- //! \cond [sanity_int_div] \endcond
+ // We calculate and print the maximum divergence of the vorticity
+ // and the integral of it
+ calc_and_print_max_div_and_int(g_vort);
// For each component solve a poisson
for (size_t i = 0 ; i < 3 ; i++)
{
//! \cond [solve_poisson_comp] \endcond
- // Copy the vorticity in gr_ps
+ // Copy the vorticity component i in gr_ps
auto it2 = gr_ps.getDomainIterator();
// calculate the velocity from the curl of phi
@@ -750,19 +678,20 @@ void comp_vel(Box<3,float> & domain, grid_type & g_vort,grid_type & g_vel, petsc
{
auto key = it2.get();
+ // copy
gr_ps.get<vorticity>(key) = g_vort.template get<vorticity>(key)[i];
++it2;
}
+ // Maximum size of the stencil
Ghost<3,long int> stencil_max(2);
// Finite difference scheme
FDScheme<poisson_nn_helm> fd(stencil_max, domain, gr_ps);
- poisson ps;
-
- fd.template impose_dit<phi>(ps,gr_ps,gr_ps.getDomainIterator());
+ // impose the poisson equation using gr_ps = vorticity for the right-hand-side (on the full domain)
+ fd.template impose_dit<phi>(poisson(),gr_ps,gr_ps.getDomainIterator());
// Create an PETSC solver
petsc_solver<double> solver;
@@ -771,29 +700,34 @@ void comp_vel(Box<3,float> & domain, grid_type & g_vort,grid_type & g_vel, petsc
solver.setAbsTol(0.001);
solver.setMaxIter(500);
- PetscBool flg;
+ // Get the sparse matrix that represent the left-hand-side
+ // of the equation
SparseMatrix<double,int,PETSC_BASE> & A = fd.getA();
+ // Get the vector that represent the right-hand-side
Vector<double,PETSC_BASE> & b = fd.getB();
- // Give to the solver A and b, return x, the solution
+ // Give to the solver A and b in this case we are giving
+ // an intitial guess phi_s calculated in the previous
+ // time step
solver.solve(A,phi_s[i],b);
- //! \cond [solve_poisson_comp] \endcond
-
- //! \cond [print_residual_copy] \endcond
-
// Calculate the residual
solError serr;
serr = solver.get_residual_error(A,phi_s[i],b);
+ Vcluster & v_cl = create_vcluster();
if (v_cl.getProcessUnitID() == 0)
{std::cout << "Solved component " << i << " Error: " << serr.err_inf << std::endl;}
// copy the solution to grid
fd.template copy<phi>(phi_s[i],gr_ps);
+ //! \cond [solve_poisson_comp] \endcond
+
+ //! \cond [copy_to_phi_v] \endcond
+
auto it3 = gr_ps.getDomainIterator();
// calculate the velocity from the curl of phi
@@ -806,7 +740,7 @@ void comp_vel(Box<3,float> & domain, grid_type & g_vort,grid_type & g_vel, petsc
++it3;
}
- //! \cond [print_residual_copy] \endcond
+ //! \cond [copy_to_phi_v] \endcond
}
//! \cond [curl_phi_v] \endcond
@@ -839,89 +773,8 @@ void comp_vel(Box<3,float> & domain, grid_type & g_vort,grid_type & g_vel, petsc
}
//! \cond [curl_phi_v] \endcond
-
- //! \cond [cross_check_curl_u_vort] \endcond
-
- g_vel.template ghost_get<velocity>();
-
- // We check that curl u = vorticity
-
- auto it4 = phi_v.getDomainIterator();
-
- double norm_max = 0.0;
-
- // calculate the velocity from the curl of phi
- while (it4.isNext())
- {
- auto key = it4.get();
-
- float phi_xy = (g_vel.template get<phi>(key.move(y,1))[x] - g_vel.template get<phi>(key.move(y,-1))[x])*inv_dy;
- float phi_xz = (g_vel.template get<phi>(key.move(z,1))[x] - g_vel.template get<phi>(key.move(z,-1))[x])*inv_dz;
- float phi_yx = (g_vel.template get<phi>(key.move(x,1))[y] - g_vel.template get<phi>(key.move(x,-1))[y])*inv_dx;
- float phi_yz = (g_vel.template get<phi>(key.move(z,1))[y] - g_vel.template get<phi>(key.move(z,-1))[y])*inv_dz;
- float phi_zx = (g_vel.template get<phi>(key.move(x,1))[z] - g_vel.template get<phi>(key.move(x,-1))[z])*inv_dx;
- float phi_zy = (g_vel.template get<phi>(key.move(y,1))[z] - g_vel.template get<phi>(key.move(y,-1))[z])*inv_dy;
-
- phi_v.template get<velocity>(key)[x] = phi_zy - phi_yz + g_vort.template get<vorticity>(key)[x];
- phi_v.template get<velocity>(key)[y] = phi_xz - phi_zx + g_vort.template get<vorticity>(key)[y];
- phi_v.template get<velocity>(key)[z] = phi_yx - phi_xy + g_vort.template get<vorticity>(key)[z];
-
- if (phi_v.template get<velocity>(key)[x] > norm_max)
- norm_max = phi_v.template get<velocity>(key)[x];
-
- if (phi_v.template get<velocity>(key)[y] > norm_max)
- norm_max = phi_v.template get<velocity>(key)[y];
-
- if (phi_v.template get<velocity>(key)[z] > norm_max)
- norm_max = phi_v.template get<velocity>(key)[z];
-
- ++it4;
- }
-
- v_cl.max(norm_max);
- v_cl.execute();
-
- if (v_cl.getProcessUnitID() == 0)
- {std::cout << "Norm max: " << norm_max << std::endl;}
-
- //! \cond [cross_check_curl_u_vort] \endcond
}
-//! \cond [comp_vel] \endcond
-
-void remesh(particles_type & vd, grid_type & g_rhs, grid_type & g_vel,Box<3,float> & domain)
-{
- constexpr int rhs_g = 0;
-
- // Remove all particles
- vd.clear();
-
- auto git = vd.getGridIterator(g_rhs.getGridInfo().getSize());
-
- while (git.isNext())
- {
- auto p = git.get();
- auto key = git.get_dist();
-
- vd.add();
-
- vd.getLastPos()[x] = g_rhs.spacing(x)*p.get(x) + domain.getLow(x);
- vd.getLastPos()[y] = g_rhs.spacing(y)*p.get(y) + domain.getLow(y);
- vd.getLastPos()[z] = g_rhs.spacing(z)*p.get(z) + domain.getLow(z);
-
- vd.template getLastProp<rhs_part>()[x] = g_rhs.template get<rhs_g>(key)[x];
- vd.template getLastProp<rhs_part>()[y] = g_rhs.template get<rhs_g>(key)[y];
- vd.template getLastProp<rhs_part>()[z] = g_rhs.template get<rhs_g>(key)[z];
-
- vd.template getLastProp<p_vel>()[x] = g_vel.template get<velocity>(key)[x];
- vd.template getLastProp<p_vel>()[y] = g_vel.template get<velocity>(key)[y];
- vd.template getLastProp<p_vel>()[z] = g_vel.template get<velocity>(key)[z];
-
- ++git;
- }
-
- vd.map();
-}
template<unsigned int prp> void set_zero(grid_type & gr)
{
@@ -958,12 +811,29 @@ template<unsigned int prp> void set_zero(particles_type & vd)
}
}
+/*! \page Vortex_in_cell_petsc Vortex in Cell 3D
+ *
+ * # Step 6: Compute right hand side # {#vic_rhs_calc}
+ *
+ * Computing the right hand side is performed calculating the term
+ * \f$ (w \cdot \nabla) u \f$. For the nabla operator we use second
+ * order finite difference central scheme. The commented part is the
+ * term \f$ \nu \nabla^{2} w \f$ that we said to neglect
+ *
+ * \snippet Numerics/Vortex_in_cell/main_vic_petsc.cpp calc_rhs
+ *
+ */
+
+//! \cond [calc_rhs] \endcond
// Calculate the right hand side of the vorticity formulation
template<typename grid> void calc_rhs(grid & g_vort, grid & g_vel, grid & g_dwp)
{
+ // usefull constant
constexpr int rhs = 0;
+ // calculate several pre-factors for the stencil finite
+ // difference
float fac1 = 2.0f*nu/(g_vort.spacing(0));
float fac2 = 2.0f*nu/(g_vort.spacing(1));
float fac3 = 2.0f*nu/(g_vort.spacing(2));
@@ -1018,6 +888,27 @@ template<typename grid> void calc_rhs(grid & g_vort, grid & g_vel, grid & g_dwp)
}
}
+//! \cond [calc_rhs] \endcond
+
+/*! \page Vortex_in_cell_petsc Vortex in Cell 3D
+ *
+ * # Step 8: Runge-Kutta # {#vic_runge_kutta1}
+ *
+ * Here we do the first step of the runge kutta update. In particular we
+ * update the vorticity and position of the particles. The right-hand-side
+ * of the vorticity update is calculated on the grid and interpolated
+ * on the particles. The Runge-Kutta of order two
+ * require the following update for the vorticity and position as first step
+ *
+ * \f$ \boldsymbol w = \boldsymbol w + \frac{1}{2} \boldsymbol {rhs} \delta t \f$
+ *
+ * \f$ \boldsymbol x = \boldsymbol x + \frac{1}{2} \boldsymbol u \delta t \f$
+ *
+ * \snippet Numerics/Vortex_in_cell/main_vic_petsc.cpp runge_kutta_1
+ *
+ */
+
+//! \cond [runge_kutta_1] \endcond
void rk_step1(particles_type & particles)
{
@@ -1060,6 +951,28 @@ void rk_step1(particles_type & particles)
particles.map();
}
+//! \cond [runge_kutta_1] \endcond
+
+/*! \page Vortex_in_cell_petsc Vortex in Cell 3D
+ *
+ * # Step 13: Runge-Kutta # {#vic_runge_kutta2}
+ *
+ * Here we do the second step of the Runge-Kutta update. In particular we
+ * update the vorticity and position of the particles. The right-hand-side
+ * of the vorticity update is calculated on the grid and interpolated
+ * on the particles. The Runge-Kutta of order two
+ * require the following update for the vorticity and position as first step
+ *
+ * \f$ \boldsymbol w = \boldsymbol w + \frac{1}{2} \boldsymbol {rhs} \delta t \f$
+ *
+ * \f$ \boldsymbol x = \boldsymbol x + \frac{1}{2} \boldsymbol u \delta t \f$
+ *
+ * \snippet Numerics/Vortex_in_cell/main_vic_petsc.cpp runge_kutta_2
+ *
+ */
+
+//! \cond [runge_kutta_2] \endcond
+
void rk_step2(particles_type & particles)
{
auto it = particles.getDomainIterator();
@@ -1091,6 +1004,27 @@ void rk_step2(particles_type & particles)
particles.map();
}
+//! \cond [runge_kutta_2] \endcond
+
+/*! \page Vortex_in_cell_petsc Vortex in Cell 3D
+ *
+ * # Step 4-5-6-7: Do step # {#vic_do_step}
+ *
+ * The do step function assemble multiple steps some of them already explained.
+ * First we interpolate the vorticity from particles to mesh
+ *
+ * \snippet Numerics/Vortex_in_cell/main_vic_petsc.cpp do_step_p2m
+ *
+ * than we calculate velocity out of vorticity and the right-hand-side
+ * recalling step 5 and 6
+ *
+ * \snippet Numerics/Vortex_in_cell/main_vic_petsc.cpp step_56
+ *
+ * Finally we interpolate velocity and right-hand-side back to particles
+ *
+ * \snippet Numerics/Vortex_in_cell/main_vic_petsc.cpp inte_m2p
+ *
+ */
template<typename grid, typename vector> void do_step(vector & particles,
grid & g_vort,
@@ -1102,36 +1036,111 @@ template<typename grid, typename vector> void do_step(vector & particles,
{
constexpr int rhs = 0;
+ //! \cond [do_step_p2m] \endcond
+
set_zero<vorticity>(g_vort);
inte.template p2m<vorticity,vorticity>(particles,g_vort);
g_vort.template ghost_put<add_,vorticity>();
+ //! \cond [do_step_p2m] \endcond
+
+ //! \cond [step_56] \endcond
+
// Calculate velocity from vorticity
comp_vel(domain,g_vort,g_vel,phi_s);
-
calc_rhs(g_vort,g_vel,g_dvort);
+ //! \cond [step_56] \endcond
+
+ //! \cond [inte_m2p] \endcond
+
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);
+
+ //! \cond [inte_m2p] \endcond
}
+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>();
+ g_vel.write("grid_velocity_" + std::to_string(i), VTK_WRITER | FORMAT_BINARY);
+ g_vort.write("grid_vorticity_" + std::to_string(i), VTK_WRITER | FORMAT_BINARY);
+ }
+
+ // In order to reduce the size of the saved data we apply a threshold.
+ // We only save particles with vorticity higher than 0.1
+ vector_dist<3,float,aggregate<float>> part_save(particles.getDecomposition(),0);
+
+ auto it_s = particles.getDomainIterator();
+
+ while (it_s.isNext())
+ {
+ auto p = it_s.get();
+
+ float vort_magn = sqrt(particles.template getProp<vorticity>(p)[0] * particles.template getProp<vorticity>(p)[0] +
+ 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();
+
+ // Save and HDF5 file for checkpoint restart
+ particles.save("check_point");
+}
+
+/*! \page Vortex_in_cell_petsc Vortex in Cell 3D
+ *
+ * # Main # {#vic_main}
+ *
+ * The main function as usual call the function **openfpm_init** it define
+ * the domain where the simulation take place. The ghost size of the grid
+ * the size of the grid in grid units on each direction, the periodicity
+ * of the domain, in this case PERIODIC in each dimension and we create
+ * our basic data structure for the simulation. A grid for the vorticity
+ * **g_vort** a grid for the velocity **g_vel** a grid for the right-hand-side
+ * of the vorticity update, and the **particles** vector. Additionally
+ * we define the data structure **phi_s[3]** that store the velocity solution
+ * in the previous time-step. keeping track of the previous solution for
+ * the velocity help the interative-solver to find the solution more quickly.
+ * Using the old velocity configuration as initial guess the solver will
+ * converge in few iterations refining the old one.
+ *
+ */
int main(int argc, char* argv[])
{
// Initialize
openfpm_init(&argc,&argv);
- // It store the solution to compute velocity
- // It is used as initial guess every time we call the solver
- petsc_solver<double>::return_type phi_s[3];
-
- // velocity in the grid is the property 0, pressure is the property 1
- constexpr int velocity = 0;
- constexpr int pressure = 1;
-
// Domain, a rectangle
Box<3,float> domain({0.0,0.0,0.0},{22.0,5.0,5.0});
@@ -1149,29 +1158,46 @@ int main(int argc, char* argv[])
grid_type g_dvort(g_vort.getDecomposition(),szu,g);
particles_type particles(g_vort.getDecomposition(),0);
+ // 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];
+
+ // Parallel object
Vcluster & v_cl = create_vcluster();
+ // print out the spacing of the grid
if (v_cl.getProcessUnitID() == 0)
{std::cout << "SPACING " << g_vort.spacing(0) << " " << g_vort.spacing(1) << " " << g_vort.spacing(2) << std::endl;}
+ // initialize the ring step 1
init_ring(g_vort,domain);
+
+ // Do the helmotz projection step 2
helmotz_hodge_projection(g_vort,domain);
+ // initialize the particle on a mesh step 3
remesh(particles,g_vort,domain);
+ // Get the total number of particles
size_t tot_part = particles.size_local();
v_cl.sum(tot_part);
v_cl.execute();
- // Now we set phi_s
+ // Now we set the size of phi_s
phi_s[0].resize(tot_part,particles.size_local());
phi_s[1].resize(tot_part,particles.size_local());
phi_s[2].resize(tot_part,particles.size_local());
+ // And we initially set it to zero
phi_s[0].setZero();
phi_s[1].setZero();
phi_s[2].setZero();
+ // 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
interpolate<particles_type,grid_type,mp4_kernel<float>> inte(particles,g_vort);
// With more than 24 core we rely on the HDF5 checkpoint restart file
@@ -1179,59 +1205,56 @@ int main(int argc, char* argv[])
{g_vort.write("grid_vorticity_init", VTK_WRITER | FORMAT_BINARY);}
- // Before entering the loop we check if we have to restart a previous simulation
-
+ // Before entering the loop we check if we want to restart from
+ // a previous simulation
size_t i = 0;
+ // if we have an argument
if (argc > 1)
{
+ // take that argument
std::string restart(argv[1]);
+ // convert it into number
i = std::stoi(restart);
+ // load the file
particles.load("check_point_" + std::to_string(i));
- particles.map();
-
+ // Print to inform that we are restarting from a
+ // particular position
if (v_cl.getProcessUnitID() == 0)
{std::cout << "Restarting from " << i << std::endl;}
}
// Time Integration
-
for ( ; i < 10001 ; i++)
{
+ // do step 4-5-6-7
do_step(particles,g_vort,g_vel,g_dvort,domain,inte,phi_s);
+ // do step 8
rk_step1(particles);
+ // do step 9-10-11-12
do_step(particles,g_vort,g_vel,g_dvort,domain,inte,phi_s);
+ // do step 13
rk_step2(particles);
+ // so step 14
set_zero<vorticity>(g_vort);
inte.template p2m<vorticity,vorticity>(particles,g_vort);
g_vort.template ghost_put<add_,vorticity>();
remesh(particles,g_vort,domain);
- std::cout << "Step " << i << std::endl;
+ // print the step number
+ if (v_cl.getProcessUnitID() == 0)
+ {std::cout << "Step " << i << std::endl;}
- if (i % 100 == 0)
- {
- if (v_cl.getProcessingUnits() < 24)
- {
- particles.write("part_out_" + std::to_string(i),VTK_WRITER | FORMAT_BINARY);
- g_vort.ghost_get<vorticity>();
- g_vel.ghost_get<velocity>();
- g_vel.write("grid_velocity_" + std::to_string(i), VTK_WRITER | FORMAT_BINARY);
- g_vort.write("grid_vorticity_" + std::to_string(i), VTK_WRITER | FORMAT_BINARY);
- }
-
- // Save also for checkpoint restart
-
- particles.save("check_point_" + std::to_string(i+1));
- }
+ // every 100 steps write the output
+ if (i % 100 == 0) {check_point_and_save(particles,g_vort,g_vel,g_dvort,i);}
}
diff --git a/example/Vector/0_simple/main.cpp b/example/Vector/0_simple/main.cpp
index bcb7903fee72d4ee890abc8642501512efa2140f..994a04262897faadcf1467be8ecebdcd4002676f 100644
--- a/example/Vector/0_simple/main.cpp
+++ b/example/Vector/0_simple/main.cpp
@@ -3,6 +3,7 @@
* \subpage Vector_0_simple
* \subpage Vector_1_celllist
* \subpage Vector_1_ghost_get
+ * \subpage Vector_1_HDF5
* \subpage Vector_2_expression
* \subpage Vector_3_md
* \subpage Vector_4_reo_root
diff --git a/example/Vector/1_HDF5_save_load/Makefile b/example/Vector/1_HDF5_save_load/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..703912dfbb69414b088bbaa7fa532977f2ba03c5
--- /dev/null
+++ b/example/Vector/1_HDF5_save_load/Makefile
@@ -0,0 +1,27 @@
+include ../../example.mk
+
+CC=mpic++
+
+LDIR =
+OPT=
+
+OBJ = main.o
+
+all: hdf5
+
+
+%.o: %.cpp
+ $(CC) -O3 $(OPT) -g -c --std=c++11 -o $@ $< $(INCLUDE_PATH)
+
+hdf5: $(OBJ)
+ $(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
+
+
+run: hdf5
+ mpirun -np 2 ./hdf5
+
+.PHONY: clean all run
+
+clean:
+ rm -f *.o *~ core hdf5
+
diff --git a/example/Vector/1_HDF5_save_load/main.cpp b/example/Vector/1_HDF5_save_load/main.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..89535613254b0bd6dcee21fd5157821c4be72ee3
--- /dev/null
+++ b/example/Vector/1_HDF5_save_load/main.cpp
@@ -0,0 +1,368 @@
+ /*! \page Vector Vector
+ *
+ * \subpage Vector_0_simple
+ * \subpage Vector_1_celllist
+ * \subpage Vector_1_ghost_get
+ * \subpage Vector_2_expression
+ * \subpage Vector_3_md
+ * \subpage Vector_4_reo_root
+ * \subpage Vector_4_cp
+ * \subpage Vector_4_mp_cl
+ * \subpage Vector_5_md_vl_sym
+ * \subpage Vector_5_md_vl_sym_crs
+ * \subpage Vector_6_complex_usage
+ * \subpage Vector_7_sph_dlb
+ * \subpage Vector_7_sph_dlb_opt
+ *
+ */
+
+
+/*!
+ * \page Vector_1_HDF5 HDF5 save and load
+ *
+ *
+ * [TOC]
+ *
+ *
+ * # HDF5 Save and load # {#HDF5_save_and_load_parallel}
+ *
+ *
+ * This example show how to save and load a vector in the parallel
+ * format HDF5.
+ *
+ * ## inclusion ## {#e0_v_inclusion}
+ *
+ * In order to use distributed vectors in our code we have to include the file Vector/vector_dist.hpp
+ *
+ * \snippet Vector/1_HDF5_save_load/main.cpp inclusion
+ *
+ */
+
+//! \cond [inclusion] \endcond
+
+#include "Vector/vector_dist.hpp"
+
+//! \cond [inclusion] \endcond
+
+int main(int argc, char* argv[])
+{
+
+ /*!
+ * \page Vector_1_HDF5 HDF5 save and load
+ *
+ * ## Initialization ## {#HDF5_initialization}
+ *
+ * Here we
+ * * Initialize the library
+ * * we create a Box that define our domain
+ * * An array that define our boundary conditions
+ * * A Ghost object that will define the extension of the ghost part in physical units
+ *
+ * \snippet Vector/1_HDF5_save_load/main.cpp Initialization and parameters
+ *
+ *
+ */
+
+ //! \cond [Initialization and parameters] \endcond
+
+ // initialize the library
+ openfpm_init(&argc,&argv);
+
+ // Here we define our domain a 2D box with internals from 0 to 1.0 for x and y
+ Box<3,float> domain({0.0,0.0,0.0},{22.0,5.0,5.0});
+
+ // Here we define the boundary conditions of our problem
+ size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+
+ // extended boundary around the domain, and the processor domain
+ Ghost<3,float> g(0.01);
+
+ //! \cond [Initialization and parameters] \endcond
+
+ /*!
+ * \page Vector_1_HDF5 HDF5 save and load
+ *
+ * ## Vector instantiation ## {#HDF5_vector_inst}
+ *
+ * Here we are creating a distributed vector defined by the following parameters
+ *
+ * * **3** is the Dimensionality of the space where the objects live
+ * * **float** is the type used for the spatial coordinate of the particles
+ * * the information stored by each particle **float[3],float[3],float[3],float[3],float[3]**.
+ * The list of properties must be put into an aggregate data structure aggregate<prop1,prop2,prop3, ... >
+ *
+ * vd is the instantiation of the object
+ *
+ * The Constructor instead require:
+ *
+ * * Number of particles, 4096 in this case
+ * * Domain where is defined this structure
+ * * bc boundary conditions
+ * * g Ghost
+ *
+ *
+ * \snippet Vector/1_HDF5_save_load/main.cpp vector instantiation
+ *
+ */
+
+ //! \cond [vector instantiation] \endcond
+
+ vector_dist<3,float, aggregate<float[3],float[3],float[3],float[3],float[3]> > vd(4096,domain,bc,g);
+
+ //! \cond [vector instantiation] \endcond
+
+ /*!
+ * \page Vector_1_HDF5 HDF5 save and load
+ *
+ * ## Assign position ## {#HDF5_vector_assign}
+ *
+ * Get an iterator that go through the 4096 particles. Initially all the particles
+ * has an undefined position state. In this cycle we define its position. In this
+ * example we use iterators. Iterators are convenient way to explore/iterate data-structures in an
+ * convenient and easy way
+ *
+ * \snippet Vector/1_HDF5_save_load/main.cpp assign position
+ *
+ *
+ */
+
+ //! \cond [assign position] \endcond
+
+ auto it = vd.getDomainIterator();
+
+ while (it.isNext())
+ {
+ auto key = it.get();
+
+ // we define x, assign a random position between 0.0 and 1.0
+ vd.getPos(key)[0] = 22.0*((float)rand() / RAND_MAX);
+
+ // we define y, assign a random position between 0.0 and 1.0
+ vd.getPos(key)[1] = 5.0*((float)rand() / RAND_MAX);
+
+ // we define y, assign a random position between 0.0 and 1.0
+ vd.getPos(key)[1] = 5.0*((float)rand() / RAND_MAX);
+
+ // next particle
+ ++it;
+ }
+
+ //! \cond [assign position] \endcond
+
+ /*!
+ * \page Vector_1_HDF5 HDF5 save and load
+ *
+ * ## Mapping particles ## {#e0_s_map}
+ *
+ * On a parallel program, once we define the position, we distribute the particles according to the underlying space decomposition
+ * The default decomposition is created even before assigning the position to the object, and is calculated
+ * giving to each processor an equal portion of space minimizing the surface to reduce communication.
+ *
+ * \snippet Vector/1_HDF5_save_load/main.cpp map
+ *
+ *
+ */
+
+ //! \cond [map] \endcond
+
+ vd.map();
+
+ //! \cond [map] \endcond
+
+ /*!
+ * \page Vector_1_HDF5 HDF5 save and load
+ *
+ * ## Assign values to particles property ## {#assign_prop}
+ *
+ * We Iterate across all the particles, we assign some value
+ * to all the particles properties. Each particle has a scalar,
+ * vector and tensor property.
+ *
+ * \snippet Vector/1_HDF5_save_load/main.cpp assign property
+ *
+ *
+ */
+
+ //! \cond [assign property] \endcond
+
+ // Get a particle iterator
+ it = vd.getDomainIterator();
+
+ // For each particle ...
+ while (it.isNext())
+ {
+ // ... p
+ auto p = it.get();
+
+ // we set the properties of the particle p
+
+ vd.template getProp<0>(p)[0] = 1.0;
+ vd.template getProp<0>(p)[1] = 1.0;
+ vd.template getProp<0>(p)[2] = 1.0;
+
+ vd.template getProp<1>(p)[0] = 2.0;
+ vd.template getProp<1>(p)[1] = 2.0;
+ vd.template getProp<1>(p)[2] = 2.0;
+
+ vd.template getProp<2>(p)[0] = 3.0;
+ vd.template getProp<2>(p)[1] = 3.0;
+ vd.template getProp<2>(p)[2] = 3.0;
+
+ vd.template getProp<3>(p)[0] = 4.0;
+ vd.template getProp<3>(p)[1] = 4.0;
+ vd.template getProp<3>(p)[2] = 4.0;
+
+ vd.template getProp<4>(p)[0] = 5.0;
+ vd.template getProp<4>(p)[1] = 5.0;
+ vd.template getProp<4>(p)[2] = 5.0;
+
+
+ // next particle
+ ++it;
+ }
+
+ //! \cond [assign property] \endcond
+
+ /*!
+ * \page Vector_1_HDF5 HDF5 save and load
+ *
+ * ## Parallel IO save ## {#HDF5_par_save}
+ *
+ * To save the file we use the function **save**. The system save
+ * all the information about the particles in the file (positions and
+ * all the properties, even complex properties)
+ *
+ * \see \ref vector_example_cp
+ *
+ * It is good to notice that independently from the number of
+ * processor this function produce only one file.
+ *
+ *
+ * \note The saved file can be used for checkpoint-restart. the status
+ * of the particles or the application in general can be saved periodically
+ * and can be used later-on to restart from the latest save
+ *
+ * \snippet Vector/1_HDF5_save_load/main.cpp save_part
+ *
+ */
+
+ //! \cond [save_part] \endcond
+
+ vd.save("particles_save.hdf5");
+
+ //! \cond [save_part] \endcond
+
+ /*!
+ * \page Vector_1_HDF5 HDF5 save and load
+ *
+ * ## Parallel IO load ## {#HDF5_par_load}
+ *
+ * To load the file we use the function **load**. The system load
+ * all the information about the particles (position and
+ * all the properties, even complex)
+ *
+ * \see \ref vector_example_cp
+ *
+ * It is good to notice that this function work independently
+ * from the number of processors used to save the file
+ *
+ * \warning Despite the fact that the function works for any number of processors
+ * the **domain** parameter, the **dimensionality**, the **space type**
+ * and the **properties** of the particles must match the saved one.
+ * At the moment there is not check or control, so in case the
+ * parameters does not match the load will produce error or
+ * corrupted data.
+ *
+ *
+ * \snippet Vector/1_HDF5_save_load/main.cpp hdf5_load
+ *
+ */
+
+ //! \cond [hdf5_load] \endcond
+
+ vector_dist<3,float, aggregate<float[3],float[3],float[3],float[3],float[3]> > vd2(vd.getDecomposition(),0);
+
+ // load the particles on another vector
+ vd2.load("particles_load.hdf5");
+
+ //! \cond [hdf5_load] \endcond
+
+ /*!
+ * \page Vector_1_HDF5 HDF5 save and load
+ *
+ * ## Process loaded data ## {#HDF5_par_}
+ *
+ * Because the function **load** works independently from the number
+ * of processors the file created with save can be used also to post-process
+ * the saved data.
+ * Once loaded the particles on the distributed vector **vd2** we can
+ * use **vd2** to post-process the data. In this case we calculate
+ * the magnitude of the property 0 and we write it to a vtk file.
+ *
+ * \snippet Vector/1_HDF5_save_load/main.cpp hdf5_post_process
+ *
+ */
+
+ //! \cond [hdf5_post_process] \endcond
+
+ vector_dist<3,float, aggregate<float> > vd3(vd.getDecomposition(),0);
+
+ auto it2 = vd2.getDomainIterator();
+
+ while (it2.isNext())
+ {
+ auto p = it2.get();
+
+ float magn_vort = sqrt(vd2.getProp<0>(p)[0]*vd2.getProp<0>(p)[0] +
+ vd2.getProp<0>(p)[1]*vd2.getProp<0>(p)[1] +
+ vd2.getProp<0>(p)[2]*vd2.getProp<0>(p)[2]);
+
+ if (magn_vort < 0.1)
+ {
+ ++it2;
+ continue;
+ }
+
+ vd3.add();
+
+ vd3.getLastPos()[0] = vd2.getPos(p)[0];
+ vd3.getLastPos()[1] = vd2.getPos(p)[1];
+ vd3.getLastPos()[2] = vd2.getPos(p)[2];
+
+ vd3.template getLastProp<0>() = magn_vort;
+
+ ++it2;
+ }
+
+ // We write the vtk file out from vd3
+ vd3.write("output", VTK_WRITER | FORMAT_BINARY );
+ vd3.save("particles_post_process_2");
+
+ //! \cond [hdf5_post_process] \endcond
+
+ /*!
+ * \page Vector_1_HDF5 HDF5 save and load
+ *
+ * ## Finalize ## {#finalize_e0_sim}
+ *
+ * At the very end of the program we have always de-initialize the library
+ *
+ * \snippet Vector/0_simple/main.cpp finalize
+ *
+ */
+
+ //! \cond [finalize] \endcond
+
+ openfpm_finalize();
+
+ //! \cond [finalize] \endcond
+
+ /*!
+ * \page Vector_0_simple Vector 0 simple
+ *
+ * ## Full code ## {#code_e0_sim}
+ *
+ * \include Vector/0_simple/main.cpp
+ *
+ */
+}