# Adapting the mAL preconditioner to stokes-fieldsplit-2d

**URL:** <https://community.freefem.org/t/adapting-the-mal-preconditioner-to-stokes-fieldsplit-2d/4041>\
**Category:** General Discussion\
**Created:** [August 15, 2025, 1:56pm UTC](https://community.freefem.org/t/adapting-the-mal-preconditioner-to-stokes-fieldsplit-2d/4041 "2025-08-15T13:56:36Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![AzizTakhirov](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/aziztakhirov/32/26_2.png) [@AzizTakhirov](https://community.freefem.org/u/AzizTakhirov)\
**Post date:** [August 15, 2025, 1:56pm UTC](https://community.freefem.org/t/adapting-the-mal-preconditioner-to-stokes-fieldsplit-2d/4041/1 "2025-08-15T13:56:36Z")

</div>

Hello,

I want to adapt the implementation of the mAL-preconditioner explained at [https://github.com/prj-/moulin2019al](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,

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [August 15, 2025, 6:17pm UTC](https://community.freefem.org/t/adapting-the-mal-preconditioner-to-stokes-fieldsplit-2d/4041/2 "2025-08-15T18:17:19Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![AzizTakhirov](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/aziztakhirov/32/26_2.png) [@AzizTakhirov](https://community.freefem.org/u/AzizTakhirov)\
**Post date:** [January 26, 2026, 8:27pm UTC](https://community.freefem.org/t/adapting-the-mal-preconditioner-to-stokes-fieldsplit-2d/4041/3 "2026-01-26T20:27:11Z")

</div>

I would like to understand the correct syntax for implementing this preconditioner. Of course, the end goal isn’t 2D Stokes, but 3D problems. Here is the code attached.

[stokes-fieldsplit-2d-PETSc\_bug.edp](https://community.freefem.org/uploads/short-url/1d57GMkFmZWPBaXWwlhRz9dWJH8.edp) (4.2 KB)

---

<div class="post-metadata">

**Author:** ![AzizTakhirov](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/aziztakhirov/32/26_2.png) [@AzizTakhirov](https://community.freefem.org/u/AzizTakhirov)\
**Post date:** [January 29, 2026, 4:28pm UTC](https://community.freefem.org/t/adapting-the-mal-preconditioner-to-stokes-fieldsplit-2d/4041/4 "2026-01-29T16:28:21Z")

</div>

I have found the issue, it was simply that rtol was 1e-1, changing it to 1e-5 gives the right solution.