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.
Block matrix
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. \]
Code dissection
static char help[] =
"...\n\n";
static const bool debug =
false;
double operator()(
const double x,
const double y,
const double z)
const {
return 1 + x * x + y * y + z * z * z;
}
};
const double z) const {
grad(0) = 2 * x;
grad(1) = 2 * y;
grad(2) = 3 * z * z;
return grad;
}
};
double operator()(
const double x,
const double y,
const double z)
const {
return 4 + 6 * z;
}
};
true),
EntityType col_type,
nbRows = row_data.getIndices().size();
nbCols = col_data.getIndices().size();
if (row_side == col_side && row_type == col_type) {
} else {
}
}
private:
auto t_row_base = row_data.getFTensor0N();
const double alpha = t_w * area;
for (
int rr = 0; rr !=
nbRows; rr++) {
auto t_col_base = col_data.getFTensor0N(gg, 0);
for (
int cc = 0; cc !=
nbCols; cc++) {
a += alpha * t_row_base * t_col_base;
++t_col_base;
}
++t_row_base;
}
++t_w;
}
}
const int *row_indices = &*row_data.getIndices().data().begin();
const int *col_indices = &*col_data.getIndices().data().begin();
&*
locMat.data().begin(), ADD_VALUES);
&*
locMat.data().begin(), ADD_VALUES);
}
}
};
int main(
int argc,
char *argv[]) {
try {
PetscBool flg_test = PETSC_FALSE;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Poisson's problem options",
"none");
PETSC_NULL);
CHKERR PetscOptionsBool(
"-test",
"if true is ctest",
"", flg_test,
&flg_test, PETSC_NULL);
ierr = PetscOptionsEnd();
boost::shared_ptr<ForcesAndSourcesCore>
domain_lhs_fe;
boost::shared_ptr<ForcesAndSourcesCore>
boundary_lhs_fe;
boost::shared_ptr<ForcesAndSourcesCore>
domain_rhs_fe;
boost::shared_ptr<ForcesAndSourcesCore>
boundary_rhs_fe;
boost::shared_ptr<ForcesAndSourcesCore>
domain_error;
boost::shared_ptr<PoissonExample::PostProcFE>
post_proc_volume;
boost::shared_ptr<ForcesAndSourcesCore> null;
boost::shared_ptr<ForcesAndSourcesCore> boundary_penalty_lhs_fe;
{
boundary_lhs_fe, domain_rhs_fe, boundary_rhs_fe, false);
global_error, domain_error);
const double beta = 1;
boundary_penalty_lhs_fe = boost::shared_ptr<ForcesAndSourcesCore>(
boundary_penalty_lhs_fe->getOpPtrVector().push_back(
new OpS(beta));
boundary_rhs_fe->getOpPtrVector().push_back(
}
Simple *simple_interface;
{
CHKERR simple_interface->getOptions();
CHKERR simple_interface->loadFile();
1);
CHKERR simple_interface->addBoundaryField(
"L",
H1,
CHKERR simple_interface->addDataField(
"ERROR",
L2,
CHKERR simple_interface->setFieldOrder(
"U",
CHKERR simple_interface->setFieldOrder(
"L",
CHKERR simple_interface->setFieldOrder(
"ERROR", 0);
CHKERR simple_interface->setUp();
}
DM dm;
CHKERR simple_interface->getDM(&dm);
{
CHKERR DMCreateGlobalVector(dm, &
F);
DM dm_sub_KK, dm_sub_LU;
ublas::matrix<Mat> nested_matrices(2, 2);
ublas::vector<IS> nested_is(2);
CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_KK);
CHKERR DMSetType(dm_sub_KK,
"DMMOFEM");
CHKERR DMSetFromOptions(dm_sub_KK);
simple_interface->getDomainFEName().c_str());
simple_interface->getBoundaryFEName().c_str());
CHKERR DMCreateMatrix(dm_sub_KK, &nested_matrices(0, 0));
domain_lhs_fe->ksp_B = nested_matrices(0, 0);
dm_sub_KK, simple_interface->getDomainFEName(), domain_lhs_fe);
boundary_penalty_lhs_fe->ksp_B = nested_matrices(0, 0);
simple_interface->getBoundaryFEName(),
boundary_penalty_lhs_fe);
CHKERR MatAssemblyBegin(nested_matrices(0, 0), MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(nested_matrices(0, 0), MAT_FINAL_ASSEMBLY);
CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_LU);
CHKERR DMSetType(dm_sub_LU,
"DMMOFEM");
CHKERR DMSetFromOptions(dm_sub_LU);
simple_interface->getBoundaryFEName().c_str());
CHKERR DMCreateMatrix(dm_sub_LU, &nested_matrices(1, 0));
boundary_lhs_fe->ksp_B = nested_matrices(1, 0);
dm_sub_LU, simple_interface->getBoundaryFEName(), boundary_lhs_fe);
CHKERR MatAssemblyBegin(nested_matrices(1, 0), MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(nested_matrices(1, 0), MAT_FINAL_ASSEMBLY);
CHKERR MatTranspose(nested_matrices(1, 0), MAT_INITIAL_MATRIX,
&nested_matrices(0, 1));
domain_rhs_fe->ksp_f =
F;
domain_rhs_fe);
boundary_rhs_fe->ksp_f =
F;
boundary_rhs_fe);
nested_matrices(1, 1) = PETSC_NULL;
MatGetType(nested_matrices(0, 0), &
type);
cerr <<
"K " <<
type << endl;
MatGetType(nested_matrices(1, 0), &
type);
cerr <<
"C " <<
type << endl;
MatGetType(nested_matrices(0, 1), &
type);
cerr <<
"CT " <<
type << endl;
std::string wait;
cerr << "UU" << endl;
MatView(nested_matrices(0, 0), PETSC_VIEWER_DRAW_WORLD);
std::cin >> wait;
cerr << "LU" << endl;
MatView(nested_matrices(1, 0), PETSC_VIEWER_DRAW_WORLD);
std::cin >> wait;
cerr << "UL" << endl;
MatView(nested_matrices(0, 1), PETSC_VIEWER_DRAW_WORLD);
std::cin >> wait;
}
CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &nested_is[0], 2, &nested_is[0],
&nested_matrices(0, 0), &
A);
KSP solver;
CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
CHKERR KSPSetFromOptions(solver);
PC pc;
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs) {
CHKERR PCFieldSplitSetIS(pc, NULL, nested_is[0]);
CHKERR PCFieldSplitSetIS(pc, NULL, nested_is[1]);
} else {
"Only works with pre-conditioner PCFIELDSPLIT");
}
for (
int i = 0;
i != 2;
i++) {
for (
int j = 0;
j != 2;
j++) {
CHKERR MatDestroy(&nested_matrices(
i,
j));
}
}
}
{
domain_error);
global_error);
if (flg_test == PETSC_TRUE) {
}
}
{
post_proc_volume);
post_proc_volume->writeFile("out_vol.h5m");
}
CHKERR VecDestroy(&global_error);
}
return 0;
}
Creating penalty finite element instance
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,
boost::shared_ptr<ForcesAndSourcesCore> boundary_penalty_lhs_fe;
and create finite element class instance itself
with that at hand we can add appropriate user data operators
boundary_penalty_lhs_fe->getOpPtrVector().push_back(new PoissonExample::OpS(beta));
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.
Code dissection
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.
DM dm_sub_KK,dm_sub_LU;
ublas::matrix<Mat> nested_matrices(2,2);
ublas::vector<IS> nested_is(2);
Vector of IS (Index Set) is used to store global indices of block matrices.
Assembly of K and S matrices
- We create sub DM from DM
CHKERR DMCreate(PETSC_COMM_WORLD,&dm_sub_KK);
CHKERR DMSetType(dm_sub_KK,
"DMMOFEM");
- Configure and set-up DM
- Get global indices of DM for sub DM for the upper diagonal block
- Create matrix and assemble matrix \(\mathbf{K}\) and penalty matrix \(\mathbf{S}\) by iteration over domain finite elements entities and boundary finite elements entities
CHKERR DMCreateMatrix(dm_sub_KK,&nested_matrices(0,0));
domain_lhs_fe->ksp_B = nested_matrices(0,0);
boundary_penalty_lhs_fe->ksp_B = nested_matrices(0,0);
CHKERR MatAssemblyBegin(nested_matrices(0,0),MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(nested_matrices(0,0),MAT_FINAL_ASSEMBLY);
Assembly of C matrix
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
CHKERR DMCreate(PETSC_COMM_WORLD,&dm_sub_LU);
CHKERR DMSetType(dm_sub_LU,
"DMMOFEM");
CHKERR DMCreateMatrix(dm_sub_LU,&nested_matrices(1,0));
boundary_lhs_fe->ksp_B = nested_matrices(1,0);
CHKERR MatAssemblyBegin(nested_matrices(1,0),MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(nested_matrices(1,0),MAT_FINAL_ASSEMBLY);
CHKERR MatTranspose(nested_matrices(1,0),MAT_INITIAL_MATRIX,&nested_matrices(0,1));
Assembly of the right hand side vector
Now we assemble global the right hand side vector, in the usual way using global DM
domain_rhs_fe->ksp_f =
F;
boundary_rhs_fe->ksp_f =
F;
Creation of global nested matrix
For details how to create nested matrix see PETSc manual (http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateNest.html)
PETSC_COMM_WORLD,
2,&nested_is[0],2,&nested_is[0],&nested_matrices(0,0),&
A
);
Solving problem
- Create solver instance, and set-up from command line
KSP solver;
CHKERR KSPCreate(PETSC_COMM_WORLD,&solver);
CHKERR KSPSetFromOptions(solver);
- Get access to solver pre-conditioner and set-up PCFIELDSPLIT (see http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCFIELDSPLIT.html)
PC pc;
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&is_pcfs);
if(is_pcfs) {
CHKERR PCFieldSplitSetIS(pc,NULL,nested_is[0]);
CHKERR PCFieldSplitSetIS(pc,NULL,nested_is[1]);
} else {
SETERRQ(
PETSC_COMM_WORLD,
"Only works with pre-conditioner PCFIELDSPLIT"
);
}
Note that we set nested matrix as operator for pre-conditioner, and using PCFieldSplitSetIS we setting elements of the field to PCFIELDSPLIT pre-conditioner.
- Finally set-up solver and solve system of equations
Running the program
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
cd mofem_install/um/build/basic_finite_elements/poisson
make -j2
mpirun -np 2 ./analytical_poisson_field_split -file_name cube_2part.h5m -
order 3 -pc_type fieldsplit
where options in .petscrc file are
-ksp_type fgmres
-ksp_atol 1e-12
-ksp_rtol 0
-ksp_monitor ascii
-pc_fieldsplit_type schur
-pc_fieldsplit_schur_precondition self
-fieldsplit_0_ksp_type gmres
#-fieldsplit_0_ksp_monitor
-fieldsplit_0_ksp_max_it 100
-fieldsplit_0_ksp_atol 1e-6
-fieldsplit_0_ksp_rtol 1e-6
-fieldsplit_0_pc_type lu
-fieldsplit_0_pc_factor_mat_solver_type mumps
-fieldsplit_0_pc_type gamg
-fieldsplit_0_pc_mg_smoothup 5
-fieldsplit_0_pc_mg_smoothdown 5
-fieldsplit_0_mg_coarse_ksp_type preonly
-fieldsplit_0_mg_coarse_pc_type hypre
-fieldsplit_0_mg_coarse_pc_hypre_type parasails
-fieldsplit_0_mg_coarse_pc_hypre_parasails_sym SPD
-fieldsplit_0_mg_coarse_pc_hypre_parasails_thresh 0.2
-fieldsplit_0_mg_coarse_pc_hypre_parasails_logging 1
-fieldsplit_1_ksp_type minres
-fieldsplit_1_ksp_monitor
-fieldsplit_1_ksp_max_it 100
-fieldsplit_1_ksp_atol 1e-2
-fieldsplit_1_ksp_rtol 1e-2
-fieldsplit_1_pc_type lsc
-fieldsplit_1_lsc_ksp_type cg
-fieldsplit_1_lsc_ksp_atol 1e-5
-fieldsplit_1_lsc_ksp_rtol 1e-5
-fieldsplit_1_lsc_ksp_max_it 100
-fieldsplit_1_lsc_pc_type hypre
-fieldsplit_1_lsc_pc_hypre_type parasails
-fieldsplit_1_lsc_pc_hypre_parasails_sym SPD
-fieldsplit_1_lsc_pc_hypre_parasails_nlevels 1
-fieldsplit_1_lsc_pc_hypre_parasails_thresh 0.2
-fieldsplit_1_lsc_pc_hypre_parasails_logging 1
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
mpirun -np 1 ./analytical_poisson_field_split \
-file_name cube_1part.h5m -
order 3 \
-pc_type fieldsplit \
-pc_fieldsplit_schur_precondition full \
-fieldsplit_1_pc_type lu -fieldsplit_1_pc_factor_mat_solver_packag mumps \
-dm_mat_type aij
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.
Exercises
- Use direct solver for pre-conditioner for two blocks and check how it changes convergence