|
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 try to solve the 3D stokes equation for incompressible fluid with Reynold number = 0
Such equation require the inversion of a system. We will show how to produce such system and solve it using Finite differences with staggered grid. The system of equation to solve is the following
\[ \eta \partial^{2}_{x} v_x + \eta \partial^{2}_{y} v_x + \eta \partial^{2}_{z} v_x - \partial_{x} P = 0 \]
\[ \eta \partial^{2}_{x} v_y + \eta \partial^{2}_{y} v_y + \eta \partial^{2}_{z} v_y - \partial_{y} P = 0 \]
\[ \eta \partial^{2}_{x} v_z + \eta \partial^{2}_{y} v_z + \eta \partial^{2}_{z} v_z - \partial_{y} P = 0 \]
\[ \partial_{x} v_x + \partial_{y} v_y + \partial_{z} v_z = 0 \]
\( v_x = 0 \quad v_y = 1 \quad v_z = 1 \) at x = L
\( v_x = 0 \quad v_y = 0 \quad v_z = 0 \) otherwise
In order to do this we have to define a structure that contain the main information about the system
We model the equations. This part will change in the near future to use template expression parsing and simplify the process of defining equations.
\( \eta v_x \nabla v_x = eta\_lap\_vx \quad \nabla = \partial^{2}_{x} + \partial^{2}_{y} + \partial^{2}_{z}\) Step1
\( \partial_{x} P = p\_x \) Step 2
\( -p\_x = m\_ p\_ x \) Step 3
\( eta\_lap\_vx - m\_p\_x \) Step4. This is the first equation in the system
The second equation definition is similar to the first one
\( \partial^{forward}_{x} v_x = dx\_vx \) Step5
\( \partial^{forward}_{y} v_y = dy\_vy \) Step6
\( \partial^{forward}_{z} v_z = dz\_vz \) Step7
\( dx\_vx + dy\_vy + dz_vz \) Step 8. This is the third equation in the system
In case of boundary conditions and staggered grid we need a particular set of equations at boundaries. Explain in detail is out of the scope of this example, but to qualitatively have an idea consider the following staggered cell
\verbatim +--$--+ | | # * # | | 0--$--+
$ = velocity(y)
As we can see several properties has different position in the cell. Consider we want to impose \(v_y = 0\) on \(x=0\). Because there are not points at \(x = 0\) we have to interpolate between $ of this cell and $ of the previous cell on y Averaging using the Avg operator. The same apply for \(v_x\) on \(y=0\). Similar arguments can be done for other boundaries in order to finally reach a well defined system. Such observation can be extended In 3D leading to the following equations
After model our equation we:
Padding and Ghost differ in the fact the padding extend the domain. Ghost is an extension for each sub-domain
Distributed grid that store the solution
Solving the system above require the solution of a system like that
\( Ax = b \quad x = A^{-1}b\)
where A is the system the discretize the left hand side of the equations + boundary conditions and b discretize the right hand size + boundary conditions
FDScheme is the object that we use to produce the Matrix A and the vector b. Such object require the maximum extension of the stencil
Here we impose the system of equation, we start from the incompressibility Eq imposed in the bulk with the exception of the first point {0,0} and than we set P = 0 in {0,0}, why we are doing this is again mathematical to have a well defined system, an intuitive explanation is that P and P + c are both solution for the incompressibility equation, this produce an ill-posed problem to make it well posed we set one point in this case {0,0} the pressure to a fixed constant for convenience P = 0
The best way to understand what we are doing is to draw a smaller example like 8x8. Considering that we have one additional point on the left for padding we have a grid 9x9. If on each point we have v_x v_y and P unknown we have 9x9x3 = 243 unknown. In order to fully determine and unique solution we have to impose 243 condition. The code under impose (in the case of 9x9) between domain and bulk 243 conditions.
Once we imposed all the equations we can retrieve the Matrix A and the vector b and pass these two element to the solver. In this example we are using PETSC solvers direct/Iterative solvers. While Umfpack has only one solver, PETSC wrap several solvers. The function best_solve set the solver in the modality to try multiple solvers to solve your system. The subsequent call to solve produce a report of all the solvers tried comparing them in error-convergence and speed. If you do not use best_solve try to solve your system with the default solver GMRES (That is the most robust iterative solver method)
Once we have the solution we copy it on the grid
At the very end of the program we have always to de-initialize the library
1.9.8