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

Stokes incompressible 2D

In this example we try to solve the 2D 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 - \partial_{x} P = 0 \]

\[ \eta \partial^{2}_{x} v_y + \eta \partial^{2}_{y} v_y - \partial_{y} P = 0 \]

\[ \partial_{x} v_x + \partial_{y} v_y = 0 \]

\( v_x = 0 \quad v_y = 1 \) at x = L

\( v_x = 0 \quad v_y = 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 (2D problem 3D problem ...)
static const unsigned int dims = 2;
// number of fields in the system v_x, v_y, P so a total of 3
static const unsigned int nvar = 3;
// 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 2D 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};
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} \) 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

\( dx\_vx + dy\_vy \) Step 7. This is the third equation in the system

// Constant Field
struct eta
{
typedef void const_field;
static float val() {return 1.0;}
};
// Convenient constants
constexpr unsigned int v[] = {0,1};
constexpr unsigned int P = 2;
constexpr unsigned int ic = 2;
constexpr int x = 0;
constexpr int y = 1;
// Create field that we have v_x, v_y, P
typedef Field<v[x],lid_nn> v_x; // Definition 1 v_x
typedef Field<v[y],lid_nn> v_y; // Definition 2 v_y
typedef Field<P,lid_nn> Prs; // Definition 3 Prs
// Eq1 V_x
typedef D<x,Prs,lid_nn> p_x; // Step 2
typedef minus<p_x,lid_nn> m_p_x; // Step 3
// Eq2 V_y
// Eq3 Incompressibility
typedef D<x,v_x,lid_nn,FORWARD> dx_vx; // Step 5
typedef D<y,v_y,lid_nn,FORWARD> dy_vy; // Step 6
typedef sum<dx_vx,dy_vy,lid_nn> ic_eq; // 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

// Equation for boundary conditions
// Directional Avg
// Usefull constants (as MACRO)
#define EQ_1 0
#define EQ_2 1
#define EQ_3 2
Average.
Definition Average.hpp:27

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},{12.0,4.0});
// Ghost (Not important in this case but required)
Ghost<2,float> g(0.01);
// Grid points on x=96 and y=32
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<2,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.

fd.impose(ic_eq(),0.0, EQ_3, {0,0},{sz[0]-2,sz[1]-2},true);
fd.impose(Prs(), 0.0, EQ_3, {0,0},{0,0});
// Here we impose the Eq1 and Eq2
fd.impose(vx_eq(),0.0, EQ_1, {1,0},{sz[0]-2,sz[1]-2});
fd.impose(vy_eq(),0.0, EQ_2, {0,1},{sz[0]-2,sz[1]-2});
// v_x and v_y
// Imposing B1
fd.impose(v_x(),0.0, EQ_1, {0,0},{0,sz[1]-2});
fd.impose(avg_vy_f(),0.0, EQ_2 , {-1,0},{-1,sz[1]-1});
// Imposing B2
fd.impose(v_x(),0.0, EQ_1, {sz[0]-1,0},{sz[0]-1,sz[1]-2});
fd.impose(avg_vy(),1.0, EQ_2, {sz[0]-1,0},{sz[0]-1,sz[1]-1});
// Imposing B3
fd.impose(avg_vx_f(),0.0, EQ_1, {0,-1},{sz[0]-1,-1});
fd.impose(v_y(), 0.0, EQ_2, {0,0},{sz[0]-2,0});
// Imposing B4
fd.impose(avg_vx(),0.0, EQ_1, {0,sz[1]-1},{sz[0]-1,sz[1]-1});
fd.impose(v_y(), 0.0, EQ_2, {0,sz[1]-1},{sz[0]-2,sz[1]-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
//
// Padding pressure
fd.impose(Prs(), 0.0, EQ_3, {-1,-1},{sz[0]-1,-1});
fd.impose(Prs(), 0.0, EQ_3, {-1,sz[1]-1},{sz[0]-1,sz[1]-1});
fd.impose(Prs(), 0.0, EQ_3, {-1,0},{-1,sz[1]-2});
fd.impose(Prs(), 0.0, EQ_3, {sz[0]-1,0},{sz[0]-1,sz[1]-2});
// Impose v_x Padding Impose v_y padding
fd.impose(v_x(), 0.0, EQ_1, {-1,-1},{-1,sz[1]-1});
fd.impose(v_y(), 0.0, EQ_2, {-1,-1},{sz[0]-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 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)

// Create a PETSC solver
// Set the maxumum nunber if iterations
solver.setMaxIter(1000);
solver.setRestart(200);
// Give to the solver A and b, return x, the solution
auto x = solver.try_solve(fd.getA(),fd.getB());
In case T does not match the PETSC precision compilation create a stub structure.

Copy the solution on the grid and write on VTK

Once we have the solution we copy it on the grid

Finalize

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

openfpm_finalize();

Full code

#include "config.h"
#ifdef HAVE_PETSC
#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/petsc_solver.hpp"
#include "Solvers/petsc_solver.hpp"
struct lid_nn
{
// dimensionaly of the equation (2D problem 3D problem ...)
static const unsigned int dims = 2;
// number of fields in the system v_x, v_y, P so a total of 3
static const unsigned int nvar = 3;
// 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 2D 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
// Here we choose PETSC implementation
typedef SparseMatrix<double,int,PETSC_BASE> SparseMatrix_type;
// type of Vector for the linear system, this parameter is bounded by the solver
// that you are using
typedef Vector<double,PETSC_BASE> 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};
// Constant Field
struct eta
{
typedef void const_field;
static float val() {return 1.0;}
};
// Convenient constants
constexpr unsigned int v[] = {0,1};
constexpr unsigned int P = 2;
constexpr unsigned int ic = 2;
constexpr int x = 0;
constexpr int y = 1;
// Create field that we have v_x, v_y, P
typedef Field<v[x],lid_nn> v_x; // Definition 1 v_x
typedef Field<v[y],lid_nn> v_y; // Definition 2 v_y
typedef Field<P,lid_nn> Prs; // Definition 3 Prs
// Eq1 V_x
typedef D<x,Prs,lid_nn> p_x; // Step 2
typedef minus<p_x,lid_nn> m_p_x; // Step 3
// Eq2 V_y
// Eq3 Incompressibility
typedef D<x,v_x,lid_nn,FORWARD> dx_vx; // Step 5
typedef D<y,v_y,lid_nn,FORWARD> dy_vy; // Step 6
typedef sum<dx_vx,dy_vy,lid_nn> ic_eq; // Step 7
// Equation for boundary conditions
// Directional Avg
// Usefull constants (as MACRO)
#define EQ_1 0
#define EQ_2 1
#define EQ_3 2
#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, a rectangle
Box<2,float> domain({0.0,0.0},{12.0,4.0});
// Ghost (Not important in this case but required)
Ghost<2,float> g(0.01);
// Grid points on x=96 and y=32
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});
// It is the maximum extension of the stencil (order 2 laplacian stencil has extension 1)
Ghost<2,long int> stencil_max(1);
// Finite difference scheme
FDScheme<lid_nn> fd(pd, stencil_max, domain, g_dist);
fd.impose(ic_eq(),0.0, EQ_3, {0,0},{sz[0]-2,sz[1]-2},true);
fd.impose(Prs(), 0.0, EQ_3, {0,0},{0,0});
// Here we impose the Eq1 and Eq2
fd.impose(vx_eq(),0.0, EQ_1, {1,0},{sz[0]-2,sz[1]-2});
fd.impose(vy_eq(),0.0, EQ_2, {0,1},{sz[0]-2,sz[1]-2});
// v_x and v_y
// Imposing B1
fd.impose(v_x(),0.0, EQ_1, {0,0},{0,sz[1]-2});
fd.impose(avg_vy_f(),0.0, EQ_2 , {-1,0},{-1,sz[1]-1});
// Imposing B2
fd.impose(v_x(),0.0, EQ_1, {sz[0]-1,0},{sz[0]-1,sz[1]-2});
fd.impose(avg_vy(),1.0, EQ_2, {sz[0]-1,0},{sz[0]-1,sz[1]-1});
// Imposing B3
fd.impose(avg_vx_f(),0.0, EQ_1, {0,-1},{sz[0]-1,-1});
fd.impose(v_y(), 0.0, EQ_2, {0,0},{sz[0]-2,0});
// Imposing B4
fd.impose(avg_vx(),0.0, EQ_1, {0,sz[1]-1},{sz[0]-1,sz[1]-1});
fd.impose(v_y(), 0.0, EQ_2, {0,sz[1]-1},{sz[0]-2,sz[1]-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
//
// Padding pressure
fd.impose(Prs(), 0.0, EQ_3, {-1,-1},{sz[0]-1,-1});
fd.impose(Prs(), 0.0, EQ_3, {-1,sz[1]-1},{sz[0]-1,sz[1]-1});
fd.impose(Prs(), 0.0, EQ_3, {-1,0},{-1,sz[1]-2});
fd.impose(Prs(), 0.0, EQ_3, {sz[0]-1,0},{sz[0]-1,sz[1]-2});
// Impose v_x Padding Impose v_y padding
fd.impose(v_x(), 0.0, EQ_1, {-1,-1},{-1,sz[1]-1});
fd.impose(v_y(), 0.0, EQ_2, {-1,-1},{sz[0]-1,-1});
// Create a PETSC solver
// Set the maxumum nunber if iterations
solver.setMaxIter(1000);
solver.setRestart(200);
// Give to the solver A and b, return x, the solution
auto x = solver.try_solve(fd.getA(),fd.getB());
fd.template copy<velocity,pressure>(x,{0,0},{sz[0]-1,sz[1]-1},g_dist);
g_dist.write("lid_driven_cavity_p_petsc");
openfpm_finalize();
}
#else
int main(int argc, char* argv[])
{
return 0;
}
#endif