 |
v0.13.1 |
|
v0.13.1 |
In this tutorial, we show how to create sub-problems by using discrete manager (DM) and sub-problems DMs to create nested matrices and use them with field-split preconditioner. This example shows simplified use of field-split preconditioner, for more advanced example of use of sub discrete managers and field-split preconditioner see cell_forces.cpp
The field-split preconditioner is a block solver where on a block of matrices we apply Jacobi or Gauss–Seidel iterations. Alternatively, if a matrix has 2x2 blocks, a Schur complement can be approximated by preconditioner. We do not intend to explain mathematical details of block solver but merely show how to use PETSc preconditioner, i.e. PCFIELDSPLIT (see http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCFIELDSPLIT.html).
Block matrix is represented using nested matrix, see details MATNEST (see http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateNest.html).
We change only small part of the code described in tutorial COR-2: Solving the Poisson equation.
We solve the same problem to one shown in COR-2: Solving the Poisson equation. However, to remove some issues with field-split solver, we will make the upper diagonal block of matrix invertible, by adding stabilisation matrix, which in this case has the interpretation of the penalty;
\[ \left[ \begin{array}{cc} \mathbf{K}+\mathbf{S} & \mathbf{C}^\textrm{T}\\ \mathbf{C} & \mathbf{0} \end{array} \right] \left\{ \begin{array}{c} \mathbf{U}\\ \mathbf{L} \end{array} \right\} = \left[ \begin{array}{c} \mathbf{F}+\mathbf{S}\overline{U} \\ \mathbf{g} \end{array} \right],\\ \mathbf{K}= \int_\Omega (\nabla \boldsymbol\phi)^\textrm{T} \nabla \boldsymbol\phi \textrm{d}\Omega,\quad \mathbf{C} = \int_{\partial\Omega} \boldsymbol\psi^\textrm{T} \boldsymbol\phi \textrm{d}\partial\Omega,\\ \mathbf{F} = \int_\Omega \boldsymbol\phi^\textrm{T} f \textrm{d}\Omega,\quad \mathbf{g} = \int_{\partial\Omega} \boldsymbol\psi^\textrm{T} \overline{u} \textrm{d}\partial\Omega,\\ \mathbf{S} = \int_{\partial\Omega} \boldsymbol\phi^\textrm{T} \boldsymbol\phi \textrm{d}\partial\Omega. \]
This code is largely the same as in analytical_poisson.cpp, with only two places where we introduce changes. We will focus only on those parts. First we make pointer new penalty finite element instance,
and create finite element class instance itself
with that at hand we can add appropriate user data operators
Implementation of penalty operator OpS::iNtegrate is in file analytical_poisson_field_split.cpp and do not differ to what was shown in COR-3: Implementing operators for the Poisson equation.
We will start with declaration of data structure necessary for creation of nested matrix. Nested matrix is PETSc structure which is used to store block matrices and on which one can perform operations as on other matrices types.
Vector of IS (Index Set) is used to store global indices of block matrices.
Assembly of off-diagonal blocks is similar to diagonal term, note that of diagonal block is not square matrix and integration is only over finite elements entities on the boundary
Now we assemble global the right hand side vector, in the usual way using global DM
For details how to create nested matrix see PETSc manual (http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateNest.html)
In order to run the program, one should first go to the directory where the problem is located, compile the code and then run the executable file. Typically, this can be done as follows
where options in .petscrc file are
as result of this we get
0 KSP Residual norm 8.033746248668e-01 ** ParaSails Setup Pattern Statistics *********** symmetric : 1 thresh : 0.200000 num_levels : 1 Max cost (average) : 3.9e+01 (2.0e+01) Nnz (ratio) : 39 ( 0.07) Max setup pattern time: 0.0 ************************************************* ** ParaSails Setup Values Statistics ************ filter : 0.100000 loadbal : 0.000000 Final Nnz (ratio) : 39 ( 0.07) Max setup values time : 0.0 ************************************************* Setup (pattern and values) times: 0: 0.0 1: 0.0 ave: 0.0 ************************************************* ** ParaSails Setup Pattern Statistics *********** symmetric : 1 thresh : 0.200000 num_levels : 1 Max cost (average) : 6.0e+08 (4.9e+08) Nnz (ratio) : 103886 ( 1.48) Max setup pattern time: 0.0 ************************************************* ** ParaSails Setup Values Statistics ************ filter : 0.100000 loadbal : 0.000000 Final Nnz (ratio) : 91505 ( 1.31) Max setup values time : 0.4 ************************************************* Setup (pattern and values) times: 0: 0.3 1: 0.4 ave: 0.3 ************************************************* Residual norms for fieldsplit_1_ solve. 0 KSP Residual norm 8.569397744104e+03 1 KSP Residual norm 2.381604391507e+03 2 KSP Residual norm 1.122966950794e+03 3 KSP Residual norm 6.059895318115e+02 4 KSP Residual norm 2.785669509756e+02 5 KSP Residual norm 1.439683032222e+02 6 KSP Residual norm 6.002357292386e+01 1 KSP Residual norm 7.690377439031e-04 Residual norms for fieldsplit_1_ solve. 0 KSP Residual norm 2.956833059930e+05 1 KSP Residual norm 1.787320287092e+05 2 KSP Residual norm 9.550287792776e+04 3 KSP Residual norm 4.057604924174e+04 4 KSP Residual norm 2.237002256464e+04 5 KSP Residual norm 1.082742229824e+04 6 KSP Residual norm 3.839937493403e+03 7 KSP Residual norm 1.837737978611e+03 2 KSP Residual norm 4.494035527480e-06 Residual norms for fieldsplit_1_ solve. 0 KSP Residual norm 2.705597685965e+05 1 KSP Residual norm 1.587873121132e+05 2 KSP Residual norm 7.567285121798e+04 3 KSP Residual norm 4.027211131979e+04 4 KSP Residual norm 2.030643277176e+04 5 KSP Residual norm 9.996336070718e+03 6 KSP Residual norm 4.499209738782e+03 7 KSP Residual norm 2.298423470140e+03 3 KSP Residual norm 2.350518497486e-08 Residual norms for fieldsplit_1_ solve. 0 KSP Residual norm 3.304403201367e+05 1 KSP Residual norm 1.993083865911e+05 2 KSP Residual norm 9.397960969595e+04 3 KSP Residual norm 4.928650742565e+04 4 KSP Residual norm 2.451251367195e+04 5 KSP Residual norm 1.237756537851e+04 6 KSP Residual norm 5.480296653554e+03 7 KSP Residual norm 3.101516648503e+03 4 KSP Residual norm 9.568748302483e-11 Residual norms for fieldsplit_1_ solve. 0 KSP Residual norm 4.448880690454e+05 1 KSP Residual norm 2.628153049738e+05 2 KSP Residual norm 1.281287424988e+05 3 KSP Residual norm 7.292699051221e+04 4 KSP Residual norm 3.180368227288e+04 5 KSP Residual norm 1.546252803266e+04 6 KSP Residual norm 6.807281733103e+03 7 KSP Residual norm 3.144299396548e+03 5 KSP Residual norm 5.494536394249e-13 Approximation error 3.391e-11
Note that KPS solver make 5 iterations to converge, since in this case, Schur complaint is approximated by pre-conditioner. If from other hand we use full Schur complement with and use direct solver
we will get
0 KSP Residual norm 8.033746248668e-01 1 KSP Residual norm 6.498236020454e-01 2 KSP Residual norm 9.691076621610e-12
Note that we converged in three steps in that case.
