PCD preconditioned

zxss · February 19, 2026, 11:58pm

Dear all,

I am new to preconditioners, I learnt from the oseen-2d-PETSc.edp the pressure convection diffusion preconditioner. And I have modified it to solve an unsteady Navier-Stokes equation with semi-implicit scheme for the nonlinear cenvective term. (u^n \cdot \nabla) u^{n+1}. It works well for lid driven cavity problem with Re = 100, both 2D and 3D. Here is the 2d code
2D-lid-driven.edp (4.8 KB)

However, when I change for another configuration, flow past a cylinder, it doesn’t work. It’s not the problem PCD preconditioner because it can be run but the results are wrong from paraview, not correct even when I choose a directly method with -pc_type lu, it may be caused by the block structure because if use a whole monolithic matrix(

varf vF([u, v], [uu, vv]) = int2d(Th) ( [u,v]'*[uu,vv]/dt + nu*Grad(u, v)'*Grad(uu, vv) + UgradV(oldu, oldv, u, v)' * [uu, vv] )
   + int2d(Th) ( [oldu,oldv]'*[uu,vv]/dt )

), the results are normal, the corresponding codes are
2D-flow-cylinder.edp (4.1 KB)
and
2D-full-cylinder.edp (2.2 KB).
However, for the whole matrix code, I don’t know how to set PCD preconditioner.
There is an unfinished code which I tried to implement PCD preconditioner for the whole matrix not block matrix.
2D-full-lid-driven.edp (3.4 KB)

All in all, I have two problems, one is that the block structure is good for implement PCD preconditioner(modifying from oseen code), but it doesn’t work for flow past a cylinder, maybe there is a bug I don’t find in my code(but I have checked for a long time). The second problem is can PCD preconditioner be implemented for the whole matrix instead of a block matrix, I didn’t find corresponding examples.

zxss · February 20, 2026, 8:10pm

I change with a 3*3 block matrix (2D problem with 2 blocks for velocity and 1 for pressure), still find that this code can work for lid driven cavity but fail for flow past a cylinder, but I only change the mesh and Dirichlet boundary condition. The both codes can be run, but for the flow past a cylinder, the results are weird in paraview.

// This is a FreeFEM script to solve lid driven cavity
// Du/Dt -nu * Delta u + u \cdot grad(u) + grad(p) = 0
// div(u) = 0

load "PETSc"
include "macro_ddm.idp"

//   Parameters
int n = 20;
real Re = 100;
real nu = 1./Re;
real T = 4, t = 0, dt = 0.05; 
int Nt = T/dt; // number of time steps
int[int] fforder = [1, 1, 1];

//   mesh
mesh Th = square(n, n);
broadcast(processor(0), Th);
DmeshCreate(Th);

Mat dF1, dF2, dC;

MatCreate(Th, dF1, P2);
MatCreate(Th, dF2, P2);
MatCreate(Th, dC, P1);

//   fespace
fespace Vh(Th, P2 );
fespace Qh(Th, P1 );
Vh uh, vh;
Qh ph;

Vh oldu = 0, oldv = 0;

// Macros for variational forms
macro grad(u) [dx(u), dy(u)] //
macro Ugradv(u1, u2, v) ( u1*dx(v) + u2*dy(v) ) //

varf vF1(u, uu) = int2d(Th) ( u*uu/dt + nu*grad(u)'*grad(uu) + Ugradv(oldu, oldv, u) * uu)
   + int2d(Th) ( oldu*uu/dt )+ on(3, u=1) + on(1,4, u=0);
varf vF2(v, vv) = int2d(Th) ( v*vv/dt + nu*grad(v)'*grad(vv) + Ugradv(oldu, oldv, v) * vv)
   + int2d(Th) ( oldv*vv/dt )+ on(3, v=0) + on(1,4, v=0);

varf vB1(p, uu) = int2d(Th) ( - p * dx(uu) );
varf vB2(p, vv) = int2d(Th) ( - p * dy(vv) );

varf vC(p, pp) = int2d(Th) (0.0 * p * pp);  // pressure stabilization

varf vonBorder1(u, uu) = on(1,2,4, u=1);
varf vonBorder2(v, vv) = on(1,2,4, v=0);

matrix B1, B2;
for(int i=0; i<=Nt; i++){
    t = i*dt;
    dF1 = vF1(Vh, Vh, tgv=-1);
    dF2 = vF2(Vh, Vh, tgv=-1);
    dC = vC(Qh, Qh);

    B1 = vB1(Qh, Vh); // test function determines row, trial function determines column
    B2 = vB2(Qh, Vh);
    Mat dB1(dF1, dC, B1); // dB is a Mat with the same distribution as the row of dF and the column of dC, but with values>
    Mat dB2(dF2, dC, B2);
    { // for Dirichlet BC on dF
        real[int] onBorder = vonBorder1(0, Vh, tgv=-1);
        real[int] onBorderPETSc;
        ChangeNumbering(dF1, onBorder, onBorderPETSc);
        MatZeroRows(dB1, onBorderPETSc);
    }
    { // for Dirichlet BC on dF2
        real[int] onBorder = vonBorder2(0, Vh, tgv=-1);
        real[int] onBorderPETSc;
        ChangeNumbering(dF2, onBorder, onBorderPETSc);
        MatZeroRows(dB2, onBorderPETSc);
    }
    Mat dA = [[dF1, 0, dB1],
              [0, dF2, dB2],
              [dB1', dB2', dC]];

    real[int] rhsV1 = vF1(0, Vh, tgv=-1); // why VF, set dirichlet boundary conditions
    real[int] rhsV2 = vF2(0, Vh, tgv=-1);
    real[int] rhsP(Qh.ndof);
    rhsP = 0;
    real[int] rhsPETSc;
    ChangeNumbering([dF1, dF2, dC], [rhsV1, rhsV2, rhsP], rhsPETSc);

    set(dA, sparams="-pc_type lu");

    real[int] uhPETSc;
    KSPSolve(dA, rhsPETSc, uhPETSc);
    ChangeNumbering([dF1, dF2, dC], [uh[], vh[], ph[]], uhPETSc, inverse=true, exchange=true);

    oldu[] = uh[];
    oldv[] = vh[];

    savevtk("seperate-lid-2d/Out", bin=1, Th, uh, vh, ph, dataname="u v p", order=fforder, append = i == 0 ? false : true);

    if(mpirank == 0){
        cout<<"Computing Time: "<<t<<endl;
    }
}

and

// This is a FreeFEM script to solve flow past a cylinder
// Du/Dt -nu * Delta u + u \cdot grad(u) + grad(p) = 0
// div(u) = 0

load "PETSc"
include "macro_ddm.idp"

//   Parameters
int n = 20;
real nu = 1e-3;
real T = 4, t = 0, dt = 0.05; 
int Nt = T/dt; // number of time steps
int[int] fforder = [1, 1, 1];

//   mesh
real D=0.1; // diameter of circle
real H=0.41; // height
real W=2.5; // weight
real cx0 = 0.5, cy0 = 0.2; // center of cyl.
real Um = 1;
mesh Th;
if (mpirank == 0){
    border fr1(t=0,W){x=t; y=0; label=1;}
    border fr2(t=0,H){x=W; y=t; label=2;}
    border fr3(t=W,0){x=t; y=H; label=1;}
    border fr4(t=H,0){x=0; y=t; label=3;}
    border fr5(t=0,2*pi){x=cx0+D*cos(t)/2; y=cy0+D*sin(t)/2; label=4;}
    Th=buildmesh(fr1(5*n)+fr2(n)+fr3(5*n)+fr4(n)+fr5(-n*10));
}
broadcast(processor(0), Th);
DmeshCreate(Th);

Mat dF1, dF2, dC;

MatCreate(Th, dF1, P2);
MatCreate(Th, dF2, P2);
MatCreate(Th, dC, P1);

//   fespace
fespace Vh(Th, P2 );
fespace Qh(Th, P1 );
Vh uh, vh;
Qh ph;

Vh oldu = 0, oldv = 0;

// Macros for variational forms
macro grad(u) [dx(u), dy(u)] //
macro Ugradv(u1, u2, v) ( u1*dx(v) + u2*dy(v) ) //

varf vF1(u, uu) = int2d(Th) ( u*uu/dt + nu*grad(u)'*grad(uu) + Ugradv(oldu, oldv, u) * uu)
   + int2d(Th) ( oldu*uu/dt )+ on(3, u=6*Um*y*(H-y)/H/H) + on(1,4, u=0);
varf vF2(v, vv) = int2d(Th) ( v*vv/dt + nu*grad(v)'*grad(vv) + Ugradv(oldu, oldv, v) * vv)
   + int2d(Th) ( oldv*vv/dt )+ on(3, v=0) + on(1,4, v=0);

varf vB1(p, uu) = int2d(Th) ( - p * dx(uu) );
varf vB2(p, vv) = int2d(Th) ( - p * dy(vv) );

varf vC(p, pp) = int2d(Th) (0.0 * p * pp);  // pressure stabilization

varf vonBorder1(u, uu) = on(1,3,4, u=6*Um*y*(H-y)/H/H);
varf vonBorder2(v, vv) = on(1,3,4, v=0);

matrix B1, B2;
for(int i=0; i<=Nt; i++){
    t = i*dt;
    dF1 = vF1(Vh, Vh, tgv=-1);
    dF2 = vF2(Vh, Vh, tgv=-1);
    dC = vC(Qh, Qh);

    B1 = vB1(Qh, Vh); // test function determines row, trial function determines column
    B2 = vB2(Qh, Vh);
    Mat dB1(dF1, dC, B1); // dB is a Mat with the same distribution as the row of dF and the column of dC, but with values>
    Mat dB2(dF2, dC, B2);
    { // for Dirichlet BC on dF
        real[int] onBorder = vonBorder1(0, Vh, tgv=-1);
        real[int] onBorderPETSc;
        ChangeNumbering(dF1, onBorder, onBorderPETSc);
        MatZeroRows(dB1, onBorderPETSc);
    }
    { // for Dirichlet BC on dF2
        real[int] onBorder = vonBorder2(0, Vh, tgv=-1);
        real[int] onBorderPETSc;
        ChangeNumbering(dF2, onBorder, onBorderPETSc);
        MatZeroRows(dB2, onBorderPETSc);
    }
    Mat dA = [[dF1, 0, dB1],
              [0, dF2, dB2],
              [dB1', dB2', dC]];

    real[int] rhsV1 = vF1(0, Vh, tgv=-1); // why VF, set dirichlet boundary conditions
    real[int] rhsV2 = vF2(0, Vh, tgv=-1);
    real[int] rhsP(Qh.ndof);
    rhsP = 0;
    real[int] rhsPETSc;
    ChangeNumbering([dF1, dF2, dC], [rhsV1, rhsV2, rhsP], rhsPETSc);

    set(dA, sparams="-pc_type lu");

    real[int] uhPETSc;
    KSPSolve(dA, rhsPETSc, uhPETSc);
    ChangeNumbering([dF1, dF2, dC], [uh[], vh[], ph[]], uhPETSc, inverse=true, exchange=true);

    oldu[] = uh[];
    oldv[] = vh[];

    savevtk("seperate-cylinder-2d/Out", bin=1, Th, uh, vh, ph, dataname="u v p", order=fforder, append = i == 0 ? false : true);

    if(mpirank == 0){
        cout<<"Computing Time: "<<t<<endl;
    }
}

Topic		Replies	Views
Solver for unsteady N-S equation General Discussion	33	572	January 29, 2026
Supply PETSc numbering for shared-memory block matrix General Discussion	30	1651	September 19, 2022
Iterative solution of the Navier-Stokes equation General Discussion	4	749	July 26, 2021
Poisson equation with lagrange multiplier using PETSc General Discussion	5	1004	February 27, 2023
Parallel fluid code General Discussion	10	751	June 9, 2021

PCD preconditioned

Related topics