|
OpenFPM_pdata
4.1.0
Project that contain the implementation of distributed structures
|
|
|
OpenFPM_pdata
4.1.0
Project that contain the implementation of distributed structures
|
|
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 2In 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 whileThis 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 whileIn the following we explain how each step is implemented in the code
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
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
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)
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
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 \)
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
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.
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 \).
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 \)
After we got our solution for \( \psi \) we can calculate the correction of the vorticity doing the gradient of \( \psi \).
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
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.
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
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)
We save the component \( i \) of \( \phi \) into phi_v
Once we filled phi_v we can implement (7) and calculate the curl of phi_v to recover the velocity v
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
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 \)
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 \)
The do step function assemble multiple steps some of them already explained. First we interpolate the vorticity from particles to mesh
than we calculate velocity out of vorticity and the right-hand-side recalling step 5 and 6
Finally we interpolate velocity and right-hand-side back to particles
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.
In this example we solve the Navier-Stokes equation in the vortex formulation in 3D for an incompressible fluid. (bold symbols are vectorial quantity)
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 whileThis 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 whileIn the following we explain how each step is implemented in the code
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
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)
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
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 \)
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
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.
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 \).
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 \)
After we got our solution for \( \psi \) we can calculate the correction of the vorticity doing the gradient of \( \psi \).
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
1.9.8