Adapting the mAL preconditioner to stokes-fieldsplit-2d

General Discussion
1

Hello,

I want to adapt the implementation of the mAL-preconditioner explained at https://github.com/prj-/moulin2019al to solving the 2D Stokes problem (stokes-fieldsplit-2d-PETSc.edp).

load “PETSc”
macro trueRestrict()true//
macro removeZeros()true//
macro dimension()2//
include “macro_ddm.idp”

macro def(i)[i, i#B, i#C]//
macro init(i)[i, i, i]//

real Re = 1.;
real gamma = 0.3;
real nu = 1.0/Re;
func Pk = [P2, P2, P1];

real time = mpiWtime();
mesh ThGlobal = square(getARGV(“-global”, 40), getARGV(“-global”, 40), [x, y]); // global mesh
ThGlobal = trunc(ThGlobal, (x < 0.5) || (y < 0.5), label = 5);
mesh th = movemesh(ThGlobal, [-x, y]);
th = ThGlobal + th;

fespace Wh(th, Pk); // complete space [u, v, p]
fespace Qh(th, P1); // pressure space for Schur complement
int[int][int] intersection; // local-to-neighbors renumbering
real[int] D; // partition of unity
Wh [uc1, uc2, pc];

int split = getARGV(“-split”, 1); // refinement parameter
build(th, split, intersection, D, Pk, mpiCommWorld)
[uc1, uc2, pc] = [1, 0, 0];

time = mpiWtime() - time;
if(mpirank == 0) cout << " ### Building mesh done in " << time << “s” << endl;

macro grad(u)[dx(u), dy(u)]//
macro div(u)(dx(u#1) + dy(u#2))//

varf vRes([u1, u2, p], [v1, v2, q]) = on(1, 3, 5, u1 = 0, u2 = 0) + on(2, u1 = y*(0.5-y), u2 = 0);
varf vJ([u1, u2, p], [v1, v2, q]) = int2d(th)( nu*( grad(u1)’ * grad(v1) + grad(u2)’ * grad(v2) ) + gamma*div(u)*div(v)

  • div(u) * q - div(v) * p) + on(1, 3, 5, u1 = 0, u2 = 0) + on(2, u1 = y*(0.5-y), u2 = 0);

verbosity = 1;
Mat A(Wh.ndof, intersection, D);
verbosity = 0;

/# Fields #/
Wh [vX, vY, p] = [1, 2, 3]; // numbering of each field
string[int] names(3); // prefix of each field
names[0] = “vX”; // %\color{DarkGreen}{x})-velocity
names[1] = “vY”; // %\color{DarkGreen}{y})-velocity
names[2] = “p”; // pressure
/# EndFields #/
/# Correspondance #/
Qh pIdx; // function from the pressure space
pIdx = 1:pIdx.n; // numbering of the unknowns of Qh
// renumbering into the complete space by doing an interpolation on Wh
Wh [listX, listY, listP] = [0, 0, pIdx];
/# EndCorrespondance #/

/# Schur #/
matrix[int] S(1); // array with a single matrix
varf vSchur(p, q) = int2d(th)
(-1.0/(gamma + 1.0/Re) * p * q); // %\color{DarkGreen}{\cref{eq:approximatedshurcomplement}}) with %\color{DarkGreen}{s=0})
S[0] = vSchur(Qh, Qh); // matrix assembly
/# EndSchur #/

if(mpirank == 0) cout << “PETSc…” << endl;
time = mpiWtime();

/# V #/
real tolV = getARGV(“-velocity_tol”, 1.0e-1); // default to %\color{DarkGreen}{10^{-1}})
// monodimensional velocity solver
string paramsV = “-ksp_type gmres -ksp_pc_side right " +
“-ksp_rtol " + tolV + " -ksp_gmres_restart 50 -pc_type asm " +
“-pc_asm_overlap 1 -sub_pc_type lu -sub_pc_factor_mat_solver_type mumps”;
if(usedARGV(”-st_ksp_converged_reason”) == -1)
paramsV = paramsV + " -ksp_converged_reason";
/# EndV #/
/# XY #/
// each velocity component gets the same monodimensional solver
// defined by paramsV
string paramsXY = “-prefix_push fieldsplit_vX_ " + paramsV + " -prefix_pop”

  • " -prefix_push fieldsplit_vY_ " + paramsV + " -prefix_pop";
    /# EndXY #/

/# P #/
string paramsP = "-prefix_push fieldsplit_p_ " +
“-ksp_type cg -ksp_max_it 5 -pc_type jacobi -prefix_pop”;
/# EndP #/

/# Krylov #/
string paramsKrylov = “-ksp_type fgmres -ksp_monitor”

  • " -ksp_rtol 1.0e-1 -ksp_gmres_restart 200";
    /# EndKrylov #/

/# AllParams #/
string params = paramsXY + " " + paramsP + " " + paramsKrylov +
" -pc_type fieldsplit -pc_fieldsplit_type multiplicative"; // “multiplicative” gives the desired lower block-triangular structure of the preconditioner
/# EndAllParams #/

set(A, sparams = params, fields = vX, names = names, schurPreconditioner = S, schurList = listX);

real[int] out(Wh.ndof);
out = vRes(0, Wh, tgv = -1);
matrix J = vJ(Wh, Wh, tgv = -1);
A = J;
real[int] xPETSc;
changeNumbering(A, uc1, xPETSc);
uc1 = 0.0;
uc1 = A^-1 * out;
changeNumbering(A, uc1, xPETSc, inverse = true, exchange = true);

time = mpiWtime() - time;
if(mpirank == 0) cout << " ### PETSc done in " << time << “s” << endl;

// macro def2(u)[u#1, u#2]// EOM
macro def1(u)u// EOM
// plotMPI(th, def2(uc), [P2, P2], def2, real, cmm = “Global velocity with fieldsplit preconditioner”);
plotMPI(th, pc, P1, def1, real, cmm = “Global pressure with fieldsplit preconditioner”);

However, I am getting a wrong result and am unable to find the issue. If anyone could point me in the right direction, that would be great.

Thanks,

2

There are other perfectly functioning preconditioners for Stokes, do you have a specific need for using augmented Lagrangian?

Also, your code is unusable, so it’s not possible to help you.