v_y{bb}; // Definition 2 v_y
Field Prs{cc}; // Definition 3 Prs
coeff eta{1.0,bl}; // Coefficient
// Create the derivatives and Laplacians that are needed
Laplacian lap;
Der dx;
Der dy;
Der dx_f;
Der dy_f;
// Equations
auto vx_eq = eta*lap(v_x) - dx(Prs); // Eq1 V_x
auto vy_eq = eta*lap(v_y) - dy(Prs); // Eq2 V_y
auto ic_eq = dx_f(v_x) + dy_f(v_y); // Eq3 Incom
Avg avg_vy{v_y};
Avg avg_vx{v_x};
Avg avg_vy_f{v_y};
Avg avg_vx_f{v_x};
//! \cond [def equation] \endcond
// 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});
//! \cond [init] \endcond
/*!
* \page Stokes_0_2D_petsc Stokes incompressible 2D petsc
*
* Distributed grid that store the solution
*
* \see \ref e0_s_grid_inst
*
* \snippet Numerics/Stoke_flow/0_2D_incompressible/main_petsc.cpp grid inst
*
*/
//! \cond [grid inst] \endcond
grid_dist_id<2,float,aggregate> g_dist(szu,domain,g);
//! \cond [grid inst] \endcond
/*!
* \page Stokes_0_2D_petsc Stokes incompressible 2D petsc
*
* Solving the system above require the solution of a system like that
*
* \f$ Ax = b \quad x = A^{-1}b\f$
*
* 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
*
* \snippet Numerics/Stoke_flow/0_2D_incompressible/main_petsc.cpp fd scheme
*
*/
//! \cond [fd scheme] \endcond
// 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 fd(pd, stencil_max, domain, g_dist);
//! \cond [fd scheme] \endcond
/*!
* \page Stokes_0_2D_petsc Stokes incompressible 2D petsc
*
* ## Impose the equation on the domain ## {#num_sk_inc_2D_ps_ied}
*
* 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.
*
* \snippet Numerics/Stoke_flow/0_2D_incompressible/main_petsc.cpp impose eq dom
*
*
*/
//! \cond [impose eq dom] \endcond
fd.impose(ic_eq,0.0,{0,0},{sz[0]-2,sz[1]-2},cc,true);
fd.impose(Prs,0.0,{0,0},{0,0},cc);
// Here we impose the Eq1 and Eq2
fd.impose(vx_eq,0.0,{1,0},{sz[0]-2,sz[1]-2},ll);
fd.impose(vy_eq,0.0,{0,1},{sz[0]-2,sz[1]-2},bb);
// v_x and v_y
// Imposing B1
fd.impose(v_x,0.0,{0,0},{0,sz[1]-2},ll);
fd.impose(avg_vy_f,0.0,{-1,0},{-1,sz[1]-1},bb);
// Imposing B2
fd.impose(v_x,0.0,{sz[0]-1,0},{sz[0]-1,sz[1]-2},ll);
fd.impose(avg_vy,1.0,{sz[0]-1,0},{sz[0]-1,sz[1]-1},bb);
// Imposing B3
fd.impose(avg_vx_f,0.0,{0,-1},{sz[0]-1,-1},ll);
fd.impose(v_y,0.0,{0,0},{sz[0]-2,0},bb);
// Imposing B4
fd.impose(avg_vx,0.0,{0,sz[1]-1},{sz[0]-1,sz[1]-1},ll);
fd.impose(v_y,0.0,{0,sz[1]-1},{sz[0]-2,sz[1]-1},bb);
// 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,{-1,-1},{sz[0]-1,-1},cc);
fd.impose(Prs,0.0,{-1,sz[1]-1},{sz[0]-1,sz[1]-1},cc);
fd.impose(Prs,0.0,{-1,0},{-1,sz[1]-2},cc);
fd.impose(Prs,0.0,{sz[0]-1,0},{sz[0]-1,sz[1]-2},cc);
// Impose v_x Padding Impose v_y padding
fd.impose(v_x,0.0,{-1,-1},{-1,sz[1]-1},cc);
fd.impose(v_y,0.0,{-1,-1},{sz[0]-1,-1},cc);
//! \cond [impose eq dom] \endcond
/*!
* \page Stokes_0_2D_petsc Stokes incompressible 2D petsc
*
* ## Solve the system of equation ## {#num_sk_inc_2D_sse}
*
* 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)
*
* \snippet Numerics/Stoke_flow/0_2D_incompressible/main_petsc.cpp solver
*
*/
//! \cond [solver] \endcond
// Create a PETSC solver
petsc_solver 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());
//! \cond [solver] \endcond
/*!
* \page Stokes_0_2D_petsc Stokes incompressible 2D petsc
*
* ## Copy the solution on the grid and write on VTK ## {#num_sk_inc_2D_ps_csg}
*
* Once we have the solution we copy it on the grid
*
* \snippet Numerics/Stoke_flow/0_2D_incompressible/main_petsc.cpp copy write
*
*/
//! \cond [copy write] \endcond
fd.template copy(x,{0,0},{sz[0]-1,sz[1]-1},g_dist);
g_dist.write("lid_driven_cavity_p_petsc");
//! \cond [cony write] \endcond
/*!
* \page Stokes_0_2D_petsc Stokes incompressible 2D petsc
*
* ## Finalize ## {#num_sk_inc_2D_ps_fin}
*
* At the very end of the program we have always to de-initialize the library
*
* \snippet Numerics/Stoke_flow/0_2D_incompressible/main_petsc.cpp fin lib
*
*/
//! \cond [fin lib] \endcond
openfpm_finalize();
//! \cond [fin lib] \endcond
/*!
* \page Stokes_0_2D_petsc Stokes incompressible 2D petsc
*
* # Full code # {#num_sk_inc_2D_ps_code}
*
* \include Numerics/Stoke_flow/0_2D_incompressible/main_petsc.cpp
*
*/
}
#else
int main(int argc, char* argv[])
{
return 0;
}
#endif