OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
Stokes incompressible 3D eigen

Stokes incompressible 3D

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

General Model properties

In order to do this we have to define a structure that contain the main information about the system

#include "Grid/grid_dist_id.hpp"
#include "Matrix/SparseMatrix.hpp"
#include "Vector/Vector.hpp"
#include "FiniteDifference/FDScheme.hpp"
#include "FiniteDifference/util/common.hpp"
#include "FiniteDifference/eq.hpp"
#include "Solvers/umfpack_solver.hpp"
#include "Solvers/petsc_solver.hpp"
struct lid_nn
{
// dimensionaly of the equation ( 3D problem ...)
static const unsigned int dims = 3;
// number of fields in the system v_x, v_y, v_z, P so a total of 4
static const unsigned int nvar = 4;
// boundary conditions PERIODIC OR NON_PERIODIC
static const bool boundary[];
// type of space float, double, ...
typedef float stype;
// type of base grid, it is the distributed grid that will store the result
// Note the first property is a 3D vector (velocity), the second is a scalar (Pressure)
// type of SparseMatrix, for the linear system, this parameter is bounded by the solver
// that you are using, in case of umfpack using <double,int> it is the only possible choice
// By default SparseMatrix is EIGEN based
typedef SparseMatrix<double,int> SparseMatrix_type;
// type of Vector for the linear system, this parameter is bounded by the solver
// that you are using, in case of umfpack using <double> it is the only possible choice
typedef Vector<double> Vector_type;
// Define that the underline grid where we discretize the system of equation is staggered
static const int grid_type = STAGGERED_GRID;
};
const bool lid_nn::boundary[] = {NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
This class decompose a space into sub-sub-domains and distribute them across processors.
Sparse Matrix implementation.
Sparse Matrix implementation stub object when OpenFPM is compiled with no linear algebra support.
Definition Vector.hpp:40
This is a distributed grid.
[Definition of the system]

System of equation modeling

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

// Constant Field
struct eta
{
typedef void const_field;
static float val() {return 1.0;}
};
// Model the equations
constexpr unsigned int v[] = {0,1,2};
constexpr unsigned int P = 3;
constexpr unsigned int ic = 3;
constexpr int x = 0;
constexpr int y = 1;
constexpr int z = 2;
// Eq1 V_x
typedef D<x,Prs,lid_nn> p_x; // Step 2
typedef minus<p_x,lid_nn> m_p_x; // Step3
// Eq2 V_y
// Eq3 V_z
// Eq4 Incompressibility
typedef D<x,v_x,lid_nn,FORWARD> dx_vx; // Step5
typedef D<y,v_y,lid_nn,FORWARD> dy_vy; // Step6
typedef D<z,v_z,lid_nn,FORWARD> dz_vz; // Step 7
Derivative second order on h (spacing)
Definition eq.hpp:83
Test structure used for several test.
[Definition of the system]
It ancapsulate the minus operation.
Definition sum.hpp:142
It model an expression expr1 * expr2.
Definition mul.hpp:120
It model an expression expr1 + ... exprn.
Definition sum.hpp:93

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(x)

$ = velocity(y)

  • = pressure

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

Initialization

After model our equation we:

  • Initialize the library
  • Define some useful constants
  • define Ghost size
  • Non-periodic boundary conditions
  • Padding domain expansion

Padding and Ghost differ in the fact the padding extend the domain. Ghost is an extension for each sub-domain

// Initialize
openfpm_init(&argc,&argv);
// 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<2,float> domain({0.0,0.0},{3.0,1.0});
// Ghost (Not important in this case but required)
Ghost<2,float> g(0.01);
// Grid points on x=256 and y=64
long int sz[] = {96,32};
size_t szu[2];
szu[0] = (size_t)sz[0];
szu[1] = (size_t)sz[1];
// We need one more point on the left and down part of the domain
// This is given by the boundary conditions that we impose.
//
Padding<2> pd({1,1},{0,0});
This class represent an N-dimensional box.
Definition Box.hpp:61
Class that contain Padding information on each direction positive and Negative direction.
Definition Ghost.hpp:127

Distributed grid that store the solution

See also
Grid instantiation

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

// It is the maximum extension of the stencil (order 2 laplacian stencil has extension 1)
Ghost<3,long int> stencil_max(1);
// Finite difference scheme
FDScheme<lid_nn> fd(pd, stencil_max, domain, g_dist);
Finite Differences.
Definition FDScheme.hpp:127

Impose the equation on the domain

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.

// start and end of the bulk
fd.impose(ic_eq(),0.0, EQ_4, {0,0,0},{sz[0]-2,sz[1]-2,sz[2]-2},true);
fd.impose(Prs(), 0.0, EQ_4, {0,0,0},{0,0,0});
fd.impose(vx_eq(),0.0, EQ_1, {1,0},{sz[0]-2,sz[1]-2,sz[2]-2});
fd.impose(vy_eq(),0.0, EQ_2, {0,1},{sz[0]-2,sz[1]-2,sz[2]-2});
fd.impose(vz_eq(),0.0, EQ_3, {0,0,1},{sz[0]-2,sz[1]-2,sz[2]-2});
// v_x
// R L (Right,Left)
fd.impose(v_x(),0.0, EQ_1, {0,0,0}, {0,sz[1]-2,sz[2]-2});
fd.impose(v_x(),0.0, EQ_1, {sz[0]-1,0,0},{sz[0]-1,sz[1]-2,sz[2]-2});
// T B (Top,Bottom)
fd.impose(avg_y_vx_f(),0.0, EQ_1, {0,-1,0}, {sz[0]-1,-1,sz[2]-2});
fd.impose(avg_y_vx(),0.0, EQ_1, {0,sz[1]-1,0},{sz[0]-1,sz[1]-1,sz[2]-2});
// A F (Forward,Backward)
fd.impose(avg_z_vx_f(),0.0, EQ_1, {0,-1,-1}, {sz[0]-1,sz[1]-1,-1});
fd.impose(avg_z_vx(),0.0, EQ_1, {0,-1,sz[2]-1},{sz[0]-1,sz[1]-1,sz[2]-1});
// v_y
// R L
fd.impose(avg_x_vy_f(),0.0, EQ_2, {-1,0,0}, {-1,sz[1]-1,sz[2]-2});
fd.impose(avg_x_vy(),1.0, EQ_2, {sz[0]-1,0,0},{sz[0]-1,sz[1]-1,sz[2]-2});
// T B
fd.impose(v_y(), 0.0, EQ_2, {0,0,0}, {sz[0]-2,0,sz[2]-2});
fd.impose(v_y(), 0.0, EQ_2, {0,sz[1]-1,0},{sz[0]-2,sz[1]-1,sz[2]-2});
// F A
fd.impose(avg_z_vy(),0.0, EQ_2, {-1,0,sz[2]-1}, {sz[0]-1,sz[1]-1,sz[2]-1});
fd.impose(avg_z_vy_f(),0.0, EQ_2, {-1,0,-1}, {sz[0]-1,sz[1]-1,-1});
// v_z
// R L
fd.impose(avg_x_vz_f(),0.0, EQ_3, {-1,0,0}, {-1,sz[1]-2,sz[2]-1});
fd.impose(avg_x_vz(),1.0, EQ_3, {sz[0]-1,0,0},{sz[0]-1,sz[1]-2,sz[2]-1});
// T B
fd.impose(avg_y_vz(),0.0, EQ_3, {-1,sz[1]-1,0},{sz[0]-1,sz[1]-1,sz[2]-1});
fd.impose(avg_y_vz_f(),0.0, EQ_3, {-1,-1,0}, {sz[0]-1,-1,sz[2]-1});
// F A
fd.impose(v_z(),0.0, EQ_3, {0,0,0}, {sz[0]-2,sz[1]-2,0});
fd.impose(v_z(),0.0, EQ_3, {0,0,sz[2]-1},{sz[0]-2,sz[1]-2,sz[2]-1});
// When we pad the grid, there are points of the grid that are not
// touched by the previous condition. Mathematically this lead
// to have too many variables for the conditions that we are imposing.
// Here we are imposing variables that we do not touch to zero
//
// L R
fd.impose(Prs(), 0.0, EQ_4, {-1,-1,-1},{-1,sz[1]-1,sz[2]-1});
fd.impose(Prs(), 0.0, EQ_4, {sz[0]-1,-1,-1},{sz[0]-1,sz[1]-1,sz[2]-1});
// T B
fd.impose(Prs(), 0.0, EQ_4, {0,sz[1]-1,-1}, {sz[0]-2,sz[1]-1,sz[2]-1});
fd.impose(Prs(), 0.0, EQ_4, {0,-1 ,-1}, {sz[0]-2,-1, sz[2]-1});
// F A
fd.impose(Prs(), 0.0, EQ_4, {0,0,sz[2]-1}, {sz[0]-2,sz[1]-2,sz[2]-1});
fd.impose(Prs(), 0.0, EQ_4, {0,0,-1}, {sz[0]-2,sz[1]-2,-1});
// Impose v_x v_y v_z padding
fd.impose(v_x(), 0.0, EQ_1, {-1,-1,-1},{-1,sz[1]-1,sz[2]-1});
fd.impose(v_y(), 0.0, EQ_2, {-1,-1,-1},{sz[0]-1,-1,sz[2]-1});
fd.impose(v_z(), 0.0, EQ_3, {-1,-1,-1},{sz[0]-1,sz[1]-1,-1});

Solve the system of equation

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 a serial direct solver Umfpack.

// Create an UMFPACK solver
// Give to the solver A and b, return x, the solution
auto x = solver.solve(fd.getA(),fd.getB());
stub when library compiled without eigen
static Vector< T > solve(const SparseMatrix< T, id_type, impl > &A, const Vector< T, impl > &b)
stub solve

Copy the solution on the grid and write on VTK

Once we have the solution we copy it on the grid

// Bring the solution to grid
fd.template copy<velocity,pressure>(x,{0,0},{sz[0]-1,sz[1]-1,sz[2]-1},g_dist);
g_dist.write("lid_driven_cavity_p_umfpack");

Finalize

At the very end of the program we have always to de-initialize the library

openfpm_finalize();

Full code

#include "Grid/grid_dist_id.hpp"
#include "Matrix/SparseMatrix.hpp"
#include "Vector/Vector.hpp"
#include "FiniteDifference/FDScheme.hpp"
#include "FiniteDifference/util/common.hpp"
#include "FiniteDifference/eq.hpp"
#include "Solvers/umfpack_solver.hpp"
#include "Solvers/petsc_solver.hpp"
struct lid_nn
{
// dimensionaly of the equation ( 3D problem ...)
static const unsigned int dims = 3;
// number of fields in the system v_x, v_y, v_z, P so a total of 4
static const unsigned int nvar = 4;
// boundary conditions PERIODIC OR NON_PERIODIC
static const bool boundary[];
// type of space float, double, ...
typedef float stype;
// type of base grid, it is the distributed grid that will store the result
// Note the first property is a 3D vector (velocity), the second is a scalar (Pressure)
// type of SparseMatrix, for the linear system, this parameter is bounded by the solver
// that you are using, in case of umfpack using <double,int> it is the only possible choice
// By default SparseMatrix is EIGEN based
typedef SparseMatrix<double,int> SparseMatrix_type;
// type of Vector for the linear system, this parameter is bounded by the solver
// that you are using, in case of umfpack using <double> it is the only possible choice
typedef Vector<double> Vector_type;
// Define that the underline grid where we discretize the system of equation is staggered
static const int grid_type = STAGGERED_GRID;
};
const bool lid_nn::boundary[] = {NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
// Constant Field
struct eta
{
typedef void const_field;
static float val() {return 1.0;}
};
// Model the equations
constexpr unsigned int v[] = {0,1,2};
constexpr unsigned int P = 3;
constexpr unsigned int ic = 3;
constexpr int x = 0;
constexpr int y = 1;
constexpr int z = 2;
// Eq1 V_x
typedef D<x,Prs,lid_nn> p_x; // Step 2
typedef minus<p_x,lid_nn> m_p_x; // Step3
// Eq2 V_y
// Eq3 V_z
// Eq4 Incompressibility
typedef D<x,v_x,lid_nn,FORWARD> dx_vx; // Step5
typedef D<y,v_y,lid_nn,FORWARD> dy_vy; // Step6
typedef D<z,v_z,lid_nn,FORWARD> dz_vz; // Step 7
// Directional Avg
#define EQ_1 0
#define EQ_2 1
#define EQ_3 2
#define EQ_4 3
#include "Vector/vector_dist.hpp"
#include "data_type/aggregate.hpp"
int main(int argc, char* argv[])
{
// Initialize
openfpm_init(&argc,&argv);
// velocity in the grid is the property 0, pressure is the property 1
constexpr int velocity = 0;
constexpr int pressure = 1;
// Domain
Box<3,float> domain({0.0,0.0,0.0},{3.0,1.0,1.0});
// Ghost (Not important in this case but required)
Ghost<3,float> g(0.01);
// Grid points on x=36 and y=12 z=12
long int sz[] = {36,12,12};
size_t szu[3];
szu[0] = (size_t)sz[0];
szu[1] = (size_t)sz[1];
szu[2] = (size_t)sz[2];
// We need one more point on the left and down part of the domain
// This is given by the boundary conditions that we impose.
//
Padding<3> pd({1,1,1},{0,0,0});
// It is the maximum extension of the stencil (order 2 laplacian stencil has extension 1)
Ghost<3,long int> stencil_max(1);
// Finite difference scheme
FDScheme<lid_nn> fd(pd, stencil_max, domain, g_dist);
// start and end of the bulk
fd.impose(ic_eq(),0.0, EQ_4, {0,0,0},{sz[0]-2,sz[1]-2,sz[2]-2},true);
fd.impose(Prs(), 0.0, EQ_4, {0,0,0},{0,0,0});
fd.impose(vx_eq(),0.0, EQ_1, {1,0},{sz[0]-2,sz[1]-2,sz[2]-2});
fd.impose(vy_eq(),0.0, EQ_2, {0,1},{sz[0]-2,sz[1]-2,sz[2]-2});
fd.impose(vz_eq(),0.0, EQ_3, {0,0,1},{sz[0]-2,sz[1]-2,sz[2]-2});
// v_x
// R L (Right,Left)
fd.impose(v_x(),0.0, EQ_1, {0,0,0}, {0,sz[1]-2,sz[2]-2});
fd.impose(v_x(),0.0, EQ_1, {sz[0]-1,0,0},{sz[0]-1,sz[1]-2,sz[2]-2});
// T B (Top,Bottom)
fd.impose(avg_y_vx_f(),0.0, EQ_1, {0,-1,0}, {sz[0]-1,-1,sz[2]-2});
fd.impose(avg_y_vx(),0.0, EQ_1, {0,sz[1]-1,0},{sz[0]-1,sz[1]-1,sz[2]-2});
// A F (Forward,Backward)
fd.impose(avg_z_vx_f(),0.0, EQ_1, {0,-1,-1}, {sz[0]-1,sz[1]-1,-1});
fd.impose(avg_z_vx(),0.0, EQ_1, {0,-1,sz[2]-1},{sz[0]-1,sz[1]-1,sz[2]-1});
// v_y
// R L
fd.impose(avg_x_vy_f(),0.0, EQ_2, {-1,0,0}, {-1,sz[1]-1,sz[2]-2});
fd.impose(avg_x_vy(),1.0, EQ_2, {sz[0]-1,0,0},{sz[0]-1,sz[1]-1,sz[2]-2});
// T B
fd.impose(v_y(), 0.0, EQ_2, {0,0,0}, {sz[0]-2,0,sz[2]-2});
fd.impose(v_y(), 0.0, EQ_2, {0,sz[1]-1,0},{sz[0]-2,sz[1]-1,sz[2]-2});
// F A
fd.impose(avg_z_vy(),0.0, EQ_2, {-1,0,sz[2]-1}, {sz[0]-1,sz[1]-1,sz[2]-1});
fd.impose(avg_z_vy_f(),0.0, EQ_2, {-1,0,-1}, {sz[0]-1,sz[1]-1,-1});
// v_z
// R L
fd.impose(avg_x_vz_f(),0.0, EQ_3, {-1,0,0}, {-1,sz[1]-2,sz[2]-1});
fd.impose(avg_x_vz(),1.0, EQ_3, {sz[0]-1,0,0},{sz[0]-1,sz[1]-2,sz[2]-1});
// T B
fd.impose(avg_y_vz(),0.0, EQ_3, {-1,sz[1]-1,0},{sz[0]-1,sz[1]-1,sz[2]-1});
fd.impose(avg_y_vz_f(),0.0, EQ_3, {-1,-1,0}, {sz[0]-1,-1,sz[2]-1});
// F A
fd.impose(v_z(),0.0, EQ_3, {0,0,0}, {sz[0]-2,sz[1]-2,0});
fd.impose(v_z(),0.0, EQ_3, {0,0,sz[2]-1},{sz[0]-2,sz[1]-2,sz[2]-1});
// When we pad the grid, there are points of the grid that are not
// touched by the previous condition. Mathematically this lead
// to have too many variables for the conditions that we are imposing.
// Here we are imposing variables that we do not touch to zero
//
// L R
fd.impose(Prs(), 0.0, EQ_4, {-1,-1,-1},{-1,sz[1]-1,sz[2]-1});
fd.impose(Prs(), 0.0, EQ_4, {sz[0]-1,-1,-1},{sz[0]-1,sz[1]-1,sz[2]-1});
// T B
fd.impose(Prs(), 0.0, EQ_4, {0,sz[1]-1,-1}, {sz[0]-2,sz[1]-1,sz[2]-1});
fd.impose(Prs(), 0.0, EQ_4, {0,-1 ,-1}, {sz[0]-2,-1, sz[2]-1});
// F A
fd.impose(Prs(), 0.0, EQ_4, {0,0,sz[2]-1}, {sz[0]-2,sz[1]-2,sz[2]-1});
fd.impose(Prs(), 0.0, EQ_4, {0,0,-1}, {sz[0]-2,sz[1]-2,-1});
// Impose v_x v_y v_z padding
fd.impose(v_x(), 0.0, EQ_1, {-1,-1,-1},{-1,sz[1]-1,sz[2]-1});
fd.impose(v_y(), 0.0, EQ_2, {-1,-1,-1},{sz[0]-1,-1,sz[2]-1});
fd.impose(v_z(), 0.0, EQ_3, {-1,-1,-1},{sz[0]-1,sz[1]-1,-1});
// Create an UMFPACK solver
// Give to the solver A and b, return x, the solution
auto x = solver.solve(fd.getA(),fd.getB());
// Bring the solution to grid
fd.template copy<velocity,pressure>(x,{0,0},{sz[0]-1,sz[1]-1,sz[2]-1},g_dist);
g_dist.write("lid_driven_cavity_p_umfpack");
openfpm_finalize();
}