OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
Vortex in Cell 3D

Vortex in Cell 3D ring with ringlets

In this example we solve the Navier-Stokes equation in the vortex formulation in 3D for an incompressible fluid. (bold symbols are vectorial quantity)

Video 1

Video 2

Numerical method

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

\( \nabla \times \boldsymbol u = - \boldsymbol w \)

\( \frac{\displaystyle D \boldsymbol w}{\displaystyle dt} = ( \boldsymbol w \cdot \vec \nabla) \boldsymbol u + \nu \nabla^{2} \boldsymbol w \) (5)

Where \(w\) is the vorticity and \(u\) is the velocity of the fluid. With Reynold number defined as \(Re = \frac{uL}{\nu}\). The algorithm can be expressed with the following pseudo code.

 1) Initialize the vortex ring on grid
 2) Do an helmotz hodge projection to make the vorticity divergent free
 3) Initialize particles on the same position as the grid or remesh

 while (t < t_end) do
    4) Interpolate vorticity from the particles to mesh
    5) calculate velocity u from the vorticity w
    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 the vorticity into mesh and reinitialize the particles
       in a grid like position
 end while

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

 1) Initialize the vortex ring on grid
 2) Do an helmotz hodge projection to make the vorticity divergent free
 3) Initialize particles on the same position as the grid or remesh

 while (t < t_end) do
    4) Interpolate vorticity from the particles to mesh
    5) calculate velocity u from the vorticity w
    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 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) 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

In the following we explain how each step is implemented in the code

Inclusion

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. 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 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

For convenience we also define the particles type and the grid type and some convenient constants

#include "interpolation/interpolation.hpp"
#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"
#include "Solvers/petsc_solver_AMG_report.hpp"
constexpr int x = 0;
constexpr int y = 1;
constexpr int z = 2;
// The type of the grids
// The type of the particles
// radius of the torus
float ringr1 = 1.0;
// radius of the core of the torus
float sigma = 1.0/3.523;
// Reynold number (If you want to use 7500.0 you have to use a grid 1600x400x400)
//float tgtre = 7500.0;
float tgtre = 3000.0;
// Kinematic viscosity
float nu = 1.0/tgtre;
// Time step
// float dt = 0.0025 for Re 7500
float dt = 0.0125;
#ifdef TEST_RUN
const unsigned int nsteps = 10;
#else
const unsigned int nsteps = 10001;
#endif
// All the properties by index
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;
This is a distributed grid.
Distributed vector.

Step 1: Initialization of the vortex ring

In this function we initialize the vortex ring. The vortex ring is initialized accordingly to these formula.

\( w(t = 0) = \frac{\Gamma}{\pi \sigma^{2}} e^{-(s/ \sigma)^2} \)

\( s^2 = (z-z_c)^{2} + ((x-x_c)^2 + (y-y_c)^2 - R^2) \)

\( \Gamma = \nu Re \)

With this initialization the vortex ring look like the one in figure

Vortex ring initialization the arrow indicate the direction where the vortex point while the colors indicate the magnitude from blue (low) to red (high)
/*
* gr is the grid where we are initializing the vortex ring
* domain is the simulation domain
*
*/
void init_ring(grid_type & gr, const Box<3,float> & domain)
{
// To add some noise to the vortex ring we create two random
// vector
constexpr int nk = 32;
float ak[nk];
float bk[nk];
for (size_t i = 0 ; i < nk ; i++)
{
ak[i] = rand()/RAND_MAX;
bk[i] = rand()/RAND_MAX;
}
// We calculate the circuitation gamma
float gamma = nu * tgtre;
float rinv2 = 1.0f/(sigma*sigma);
float max_vorticity = gamma*rinv2/M_PI;
// We go through the grid to initialize the vortex
auto it = gr.getDomainIterator();
while (it.isNext())
{
auto key_d = it.get();
auto key = it.getGKey(key_d);
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);
float theta1 = atan2((ty-domain.getHigh(1)/2.0f),(tz-domain.getHigh(2)/2.0f));
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;
float rad1t = tx - 1.0f;
float rad1sq = rad1r*rad1r + rad1t*rad1t;
float 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);
// kill the axis term
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;
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);
++it;
}
}
This class represent an N-dimensional box.
Definition Box.hpp:61
St spacing(size_t i) const
Get the spacing of the grid in direction i.
grid_dist_iterator< dim, device_grid, decltype(device_grid::type_of_subiterator()), FREE > getDomainIterator() const
It return an iterator that span the full grid domain (each processor span its local domain)

Step 2: Helmotz-hodge projection

The Helnotz hodge projection is required in order to make the vorticity divergent free. The Helmotz-holde projection work in this way. A field can be divided into a curl-free part and a divergent-free part.

\( w = w_{rot} + w_{div} \)

with

\( \vec \nabla \times w_{rot} = 0 \)

\( \nabla \cdot w_{div} = 0 \)

To have a vorticity divergent free we have to get the component (3) \(w_{div} = w - w_{rot}\). In particular it hold

\( \nabla \cdot w = \nabla \cdot w_{rot} \)

Bacause \( \vec \nabla \times w_{rot} = 0 \) we can introduce a field \( \psi \) such that

(2) \( w_{rot} = \vec \nabla \psi \)

Doing the \( \nabla \cdot \vec \nabla \psi \) we obtain

\( \nabla \cdot \vec \nabla \psi = \nabla^{2} \psi = \nabla \cdot w_{rot} = \nabla \cdot w \)

so we lead to this equation

(1) \( \nabla^{2} \psi = \nabla \cdot w \)

Solving the equation for \( \psi \) we can obtain \( w_{rot} \) doing the gradient of \( \psi \) and finally correct \( w \) obtaining \( w_{div} \)

The helmotz_hodge_projection function do this correction to the vorticity

In particular it solve the equation (1) it calculate \( w_{rot} \) using (2) and correct the vorticity using using (3)

Poisson equation

To solve a poisson equation on a grid using finite-difference, we need to create an object that carry information about the system of equations

// Specification of the poisson equation for the helmotz-hodge projection
// 3D (dims = 3). The field is a scalar value (nvar = 1), bournary are periodic
// type of the the space is float. The grid type that store \psi
// The others indicate which kind of Matrix to use to construct the linear system and
// which kind of vector to construct for the right hand side. Here we use a PETSC Sparse Matrix
// and PETSC vector. NORMAL_GRID indicate that is a standard grid (non-staggered)
{
static const unsigned int dims = 3;
static const unsigned int nvar = 1;
static const bool boundary[];
typedef float stype;
static const int grid_type = NORMAL_GRID;
};
const bool poisson_nn_helm::boundary[] = {PERIODIC,PERIODIC,PERIODIC};
Sparse Matrix implementation.
Sparse Matrix implementation stub object when OpenFPM is compiled with no linear algebra support.
Definition Vector.hpp:40
Vector< double, PETSC_BASE > Vector_type
type of vector that store the solution
SparseMatrix< double, int, PETSC_BASE > SparseMatrix_type
type of sparse grid that store the Matrix A
static const bool boundary[]
specify the boundary conditions
static const unsigned int nvar
We have only one scalar unknown.
static const unsigned int dims
Number of dimansions of the problem.
grid_dist_id< 3, float, aggregate< float > > b_grid
Grid that store the solution.

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 \( \nabla^{2} \psi \)

// Convenient constant
constexpr unsigned int phi = 0;
// We define a field 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)
Definition eq.hpp:83
Laplacian second order on h (spacing)
Definition Laplacian.hpp:23

Before to construct the linear system we also calculate the divergence of the vorticity \( \nabla \cdot w \) that will be the right-hand-side of the equation

// ghost get
gr.template ghost_get<vorticity>();
// ghost size of the psi function
// Here we create a distributed grid to store the result of the helmotz projection
// Calculate the divergence of the vortex field
auto it = gr.getDomainIterator();
while (it.isNext())
{
auto key = it.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] );
++it;
}
Decomposition & getDecomposition()
Get the object that store the information about the decomposition.
const grid_sm< dim, T > & getGridInfo() const
Get an object containing the grid informations.
__device__ __host__ const size_t(& getSize() const)[N]
Return the size of the grid as an array.
Definition grid_sm.hpp:760

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 our SparseMatrix. In this example we are imposing the poisson equation with right hand side equal to the divergence of vorticity (note: to avoid to create another field we use \( \psi \) to preliminary store the divergence of the vorticity). Imposing the equations produce an invertible SparseMatrix A and a right-hand-side Vector b.

// 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);
// Here we define our Finite difference disctetization scheme object
FDScheme<poisson_nn_helm> fd(stencil_max, domain, psi);
// impose the poisson equation, using right hand side psi on the full domain (psi.getDomainIterator)
// the template paramereter instead define which property of phi carry the righ-hand-side
// in this case phi has only one property, so the property 0 carry the right hand side
fd.template impose_dit<0>(poisson(),psi,psi.getDomainIterator());
Finite Differences.
Definition FDScheme.hpp:127

Because we need \( x = A^{-1}b \). We have to invert and solve a linear system. In this case we use the Conjugate-gradient-Method an iterative solver. Such method is controlled by two parameters. One is the tollerance that determine when the method is converged, the second one is the maximum number of iterations to avoid that 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 \( \psi \).

timer tm_solve;
if (init == true)
{
// Set the conjugate-gradient as method to solve the linear solver
solver.setSolver(KSPBCGS);
// Set the absolute tolerance to determine that the method is converged
solver.setAbsTol(0.0001);
// Set the maximum number of iterations
solver.setMaxIter(500);
#ifdef USE_MULTIGRID
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");
#endif
tm_solve.start();
// Give to the solver A and b, return x, the solution
solver.solve(fd.getA(),x_,fd.getB());
tm_solve.stop();
}
else
{
tm_solve.start();
solver.solve(x_,fd.getB());
tm_solve.stop();
}
// copy the solution x to the grid psi
fd.template copy<phi>(x_,psi);
Class for cpu time benchmarking.
Definition timer.hpp:28
void stop()
Stop the timer.
Definition timer.hpp:119
void start()
Start the timer.
Definition timer.hpp:90

Note

Because we are solving the poisson equation in periodic boundary conditions the Matrix has determinat equal to zero. This mean that \( \psi \) 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 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 \( \nu = 0 \) and \( \nabla \cdot w = 0 \) we can rewrite (5) in a conservative way \( \frac{Dw}{dt} = div(w \otimes v) \) ). Is also good to notice that the solution that you get is the one with \( \int w = 0 \)

timer tm_solve;
if (init == true)
{
// Set the conjugate-gradient as method to solve the linear solver
solver.setSolver(KSPBCGS);
// Set the absolute tolerance to determine that the method is converged
solver.setAbsTol(0.0001);
// Set the maximum number of iterations
solver.setMaxIter(500);
#ifdef USE_MULTIGRID
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");
#endif
tm_solve.start();
// Give to the solver A and b, return x, the solution
solver.solve(fd.getA(),x_,fd.getB());
tm_solve.stop();
}
else
{
tm_solve.start();
solver.solve(x_,fd.getB());
tm_solve.stop();
}
// copy the solution x to the grid psi
fd.template copy<phi>(x_,psi);

Correction

After we got our solution for \( \psi \) we can calculate the correction of the vorticity doing the gradient of \( \psi \).

psi.template ghost_get<phi>();
// Correct the vorticity to make it divergence free
auto it2 = gr.getDomainIterator();
while (it2.isNext())
{
auto key = it2.get();
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)));
++it2;
}

We also do a sanity check and we control that the vorticity remain divergent-free. Getting the maximum value of the divergence and printing out its value

Step 3: Remeshing vorticity

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 in a grid like position and assign the property vorticity of the particles equal to the grid vorticity.

void remesh(particles_type & vd, grid_type & gr,Box<3,float> & domain)
{
// 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();
}
grid_dist_id_iterator_dec< Decomposition > getGridIterator(const size_t(&sz)[dim])
void map(size_t opt=NONE)
It move all the particles that does not belong to the local processor to the respective processor.
void clear()
remove all the elements
auto getLastPos() -> decltype(v_pos.template get< 0 >(0))
Get the position of the last element.
void add()
Add local particle.

Step 5: Compute velocity from vorticity

Computing the velocity from vorticity is done in the following way. Given

\( \vec \nabla \times u = -w \)

We intrododuce the stream line function defined as

\( \nabla \times \phi = u \) (7)

\( \nabla \cdot \phi = 0 \)

We obtain

\( \nabla \times \nabla \times \phi = -w = \vec \nabla (\nabla \cdot \phi) - \nabla^{2} \phi \)

Because the divergence of \( \phi \) is defined to be zero we have

\( \nabla^{2} \phi = w \)

The velocity can be recovered by the equation (7)

Putting into code what explained before, we again generate a poisson object

// Convenient constant
constexpr unsigned int phi = 0;
// We define a field 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)

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 \( i \) 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 \( i \) (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)

// Copy the vorticity component i in gr_ps
auto it2 = gr_ps.getDomainIterator();
// calculate the velocity from the curl of phi
while (it2.isNext())
{
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);
// 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());
solver.setAbsTol(0.01);
// Get the vector that represent the right-hand-side
Vector<double,PETSC_BASE> & b = fd.getB();
timer tm_solve;
tm_solve.start();
// 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(phi_s[i],b);
tm_solve.stop();
// Calculate the residual
solError serr;
serr = solver.get_residual_error(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);
size_t getProcessUnitID()
Get the process unit id.
Implementation of VCluster class.
Definition VCluster.hpp:59
It contain statistic of the error of the calculated solution.
PetscReal err_inf
infinity norm of the error

We save the component \( i \) of \( \phi \) into phi_v

auto it3 = gr_ps.getDomainIterator();
// calculate the velocity from the curl of phi
while (it3.isNext())
{
auto key = it3.get();
phi_v.get<velocity>(key)[i] = gr_ps.get<phi>(key);
++it3;
}

Once we filled phi_v we can implement (7) and calculate the curl of phi_v to recover the velocity v

phi_v.ghost_get<phi>();
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);
auto it3 = phi_v.getDomainIterator();
// calculate the velocity from the curl of phi
while (it3.isNext())
{
auto key = it3.get();
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;
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;
}

Step 6: Compute right hand side

Computing the right hand side is performed calculating the term \( (w \cdot \nabla) u \). For the nabla operator we use second order finite difference central scheme. The commented part is the term \( \nu \nabla^{2} w \) that we said to neglect

// 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)*g_vort.spacing(0));
float fac2 = 2.0f*nu/(g_vort.spacing(1)*g_vort.spacing(1));
float fac3 = 2.0f*nu/(g_vort.spacing(2)*g_vort.spacing(2));
float fac4 = 0.5f/(g_vort.spacing(0));
float fac5 = 0.5f/(g_vort.spacing(1));
float fac6 = 0.5f/(g_vort.spacing(2));
auto it = g_dwp.getDomainIterator();
g_vort.template ghost_get<vorticity>();
while (it.isNext())
{
auto key = it.get();
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])+
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])-
2.0f*(fac1+fac2+fac3)*g_vort.template get<vorticity>(key)[x]+
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]);
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])+
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])-
2.0f*(fac1+fac2+fac3)*g_vort.template get<vorticity>(key)[y]+
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]);
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])+
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])-
2.0f*(fac1+fac2+fac3)*g_vort.template get<vorticity>(key)[z]+
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;
}
}

Step 8: Runge-Kutta

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

\( \boldsymbol w = \boldsymbol w + \frac{1}{2} \boldsymbol {rhs} \delta t \)

\( \boldsymbol x = \boldsymbol x + \frac{1}{2} \boldsymbol u \delta t \)

void rk_step1(particles_type & particles)
{
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;
}
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;
}
}
vect_dist_key_dx get()
Get the actual key.
auto getProp(vect_dist_key_dx vec_key) -> decltype(v_prp.template get< id >(vec_key.getKey()))
Get the property of an element.
auto getPos(vect_dist_key_dx vec_key) -> decltype(v_pos.template get< 0 >(vec_key.getKey()))
Get the position of an element.
vector_dist_iterator getDomainIterator() const
Get an iterator that traverse the particles in the domain.

Step 13: Runge-Kutta

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

\( \boldsymbol w = \boldsymbol w + \frac{1}{2} \boldsymbol {rhs} \delta t \)

\( \boldsymbol x = \boldsymbol x + \frac{1}{2} \boldsymbol u \delta t \)

void rk_step2(particles_type & particles)
{
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;
}
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;
}
}

Step 4-5-6-7: Do step

The do step function assemble multiple steps some of them already explained. First we interpolate the vorticity from particles to mesh

inte.template p2m<vorticity,vorticity>(particles,g_vort);
g_vort.template ghost_put<add_,vorticity>();
Meta-function to return a compatible zero-element.

than we calculate velocity out of vorticity and the right-hand-side recalling step 5 and 6

// Calculate velocity from vorticity
comp_vel(domain,g_vort,g_vel,phi_s,solver);
calc_rhs(g_vort,g_vel,g_dvort);

Finally we interpolate velocity and right-hand-side back to particles

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

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.

Vortex in Cell 3D ring with ringlets

In this example we solve the Navier-Stokes equation in the vortex formulation in 3D for an incompressible fluid. (bold symbols are vectorial quantity)

Numerical method

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

\( \nabla \times \boldsymbol u = - \boldsymbol w \)

\( \frac{\displaystyle D \boldsymbol w}{\displaystyle dt} = ( \boldsymbol w \cdot \vec \nabla) \boldsymbol u + \nu \nabla^{2} \boldsymbol w \) (5)

Where \(w\) is the vorticity and \(u\) is the velocity of the fluid. With high Reynold number \(Re = \frac{uL}{\nu}\) and the term \(uL\) significantly smaller than the reynold number we have that \(\nu\) is small and the term \(\nu \nabla^{2} w\) is negligible. The algorithm can be expressed with the following pseudo code.

 1) Initialize the vortex ring on grid
 2) Do an helmotz hodge projection to make the vorticity divergent free
 3) Initialize particles on the same position as the grid or remesh

 while (t < t_end) do
    4) Interpolate vorticity from the particles to mesh
    5) calculate velocity u from the vorticity w
    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 the vorticity into mesh and reinitialize the particles
       in a grid like position
 end while

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

 1) Initialize the vortex ring on grid
 2) Do an helmotz hodge projection to make the vorticity divergent free
 3) Initialize particles on the same position as the grid or remesh

 while (t < t_end) do
    4) 4) Interpolate vorticity from the particles to mesh
    5) calculate velocity u from the vorticity w
    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 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) 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

In the following we explain how each step is implemented in the code

Inclusion

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. 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 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

For convenience we also define the particles type and the grid type and some convenient constants

#include "interpolation/interpolation.hpp"
#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"
#include "Solvers/petsc_solver_AMG_report.hpp"
constexpr int x = 0;
constexpr int y = 1;
constexpr int z = 2;
// The type of the grids
// The type of the particles
// radius of the torus
float ringr1 = 1.0;
// radius of the core of the torus
float sigma = 1.0/3.523;
// Reynold number (If you want to use 7500.0 you have to use a grid 1600x400x400)
//float tgtre = 7500.0;
float tgtre = 3000.0;
// Kinematic viscosity
float nu = 1.0/tgtre;
// Time step
// float dt = 0.0025 for Re 7500
float dt = 0.0125;
#ifdef TEST_RUN
const unsigned int nsteps = 10;
#else
const unsigned int nsteps = 10001;
#endif
// All the properties by index
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;

Step 2: Helmotz-hodge projection

The Helnotz hodge projection is required in order to make the vorticity divergent free. The Helmotz-holde projection work in this way. A field can be divided into a curl-free part and a divergent-free part.

\( w = w_{rot} + w_{div} \)

with

\( \vec \nabla \times w_{rot} = 0 \)

\( \nabla \cdot w_{div} = 0 \)

To have a vorticity divergent free we have to get the component (3) \(w_{div} = w - w_{rot}\). In particular it hold

\( \nabla \cdot w = \nabla \cdot w_{rot} \)

Bacause \( \vec \nabla \times w_{rot} = 0 \) we can introduce a field \( \psi \) such that

(2) \( w_{rot} = \vec \nabla \psi \)

Doing the \( \nabla \cdot \vec \nabla \psi \) we obtain

\( \nabla \cdot \vec \nabla \psi = \nabla^{2} \psi = \nabla \cdot w_{rot} = \nabla \cdot w \)

so we lead to this equation

(1) \( \nabla^{2} \psi = \nabla \cdot w \)

Solving the equation for \( \psi \) we can obtain \( w_{rot} \) doing the gradient of \( \psi \) and finally correct \( w \) obtaining \( w_{div} \)

The helmotz_hodge_projection function do this correction to the vorticity

In particular it solve the equation (1) it calculate \( w_{rot} \) using (2) and correct the vorticity using using (3)

Poisson equation

To solve a poisson equation on a grid using finite-difference, we need to create an object that carry information about the system of equations

// Specification of the poisson equation for the helmotz-hodge projection
// 3D (dims = 3). The field is a scalar value (nvar = 1), bournary are periodic
// type of the the space is float. The grid type that store \psi
// The others indicate which kind of Matrix to use to construct the linear system and
// which kind of vector to construct for the right hand side. Here we use a PETSC Sparse Matrix
// and PETSC vector. NORMAL_GRID indicate that is a standard grid (non-staggered)
{
static const unsigned int dims = 3;
static const unsigned int nvar = 1;
static const bool boundary[];
typedef float stype;
static const int grid_type = NORMAL_GRID;
};
const bool poisson_nn_helm::boundary[] = {PERIODIC,PERIODIC,PERIODIC};

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 \( \nabla^{2} \psi \)

// Convenient constant
constexpr unsigned int phi = 0;
// We define a field 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)

Before to construct the linear system we also calculate the divergence of the vorticity \( \nabla \cdot w \) that will be the right-hand-side of the equation

// ghost get
gr.template ghost_get<vorticity>();
// ghost size of the psi function
// Here we create a distributed grid to store the result of the helmotz projection
// Calculate the divergence of the vortex field
auto it = gr.getDomainIterator();
while (it.isNext())
{
auto key = it.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] );
++it;
}

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 our SparseMatrix. In this example we are imposing the poisson equation with right hand side equal to the divergence of vorticity (note: to avoid to create another field we use \( \psi \) to preliminary store the divergence of the vorticity). Imposing the equations produce an invertible SparseMatrix A and a right-hand-side Vector b.

// 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);
// Here we define our Finite difference disctetization scheme object
FDScheme<poisson_nn_helm> fd(stencil_max, domain, psi);
// impose the poisson equation, using right hand side psi on the full domain (psi.getDomainIterator)
// the template paramereter instead define which property of phi carry the righ-hand-side
// in this case phi has only one property, so the property 0 carry the right hand side
fd.template impose_dit<0>(poisson(),psi,psi.getDomainIterator());

Because we need \( x = A^{-1}b \). We have to invert and solve a linear system. In this case we use the Conjugate-gradient-Method an iterative solver. Such method is controlled by two parameters. One is the tollerance that determine when the method is converged, the second one is the maximum number of iterations to avoid that 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 \( \psi \).

timer tm_solve;
if (init == true)
{
// Set the conjugate-gradient as method to solve the linear solver
solver.setSolver(KSPBCGS);
// Set the absolute tolerance to determine that the method is converged
solver.setAbsTol(0.0001);
// Set the maximum number of iterations
solver.setMaxIter(500);
#ifdef USE_MULTIGRID
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");
#endif
tm_solve.start();
// Give to the solver A and b, return x, the solution
solver.solve(fd.getA(),x_,fd.getB());
tm_solve.stop();
}
else
{
tm_solve.start();
solver.solve(x_,fd.getB());
tm_solve.stop();
}
// copy the solution x to the grid psi
fd.template copy<phi>(x_,psi);

Note

Because we are solving the poisson equation in periodic boundary conditions the Matrix has determinat equal to zero. This mean that \( \psi \) 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 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 \( \nu = 0 \) and \( \nabla \cdot w = 0 \) we can rewrite (5) in a conservative way \( \frac{Dw}{dt} = div(w \otimes v) \) ). Is also good to notice that the solution that you get is the one with \( \int w = 0 \)

timer tm_solve;
if (init == true)
{
// Set the conjugate-gradient as method to solve the linear solver
solver.setSolver(KSPBCGS);
// Set the absolute tolerance to determine that the method is converged
solver.setAbsTol(0.0001);
// Set the maximum number of iterations
solver.setMaxIter(500);
#ifdef USE_MULTIGRID
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");
#endif
tm_solve.start();
// Give to the solver A and b, return x, the solution
solver.solve(fd.getA(),x_,fd.getB());
tm_solve.stop();
}
else
{
tm_solve.start();
solver.solve(x_,fd.getB());
tm_solve.stop();
}
// copy the solution x to the grid psi
fd.template copy<phi>(x_,psi);

Correction

After we got our solution for \( \psi \) we can calculate the correction of the vorticity doing the gradient of \( \psi \).

psi.template ghost_get<phi>();
// Correct the vorticity to make it divergence free
auto it2 = gr.getDomainIterator();
while (it2.isNext())
{
auto key = it2.get();
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)));
++it2;
}

We also do a sanity check and we control that the vorticity remain divergent-free. Getting the maximum value of the divergence and printing out its value