# Petrov-Galerkin Least-Squares

**URL:** https://community.freefem.org/t/petrov-galerkin-least-squares/3979
**Category:** General Discussion
**Created:** [June 21, 2025, 11:30am UTC](https://community.freefem.org/t/petrov-galerkin-least-squares/3979 "2025-06-21T11:30:05Z")
**Posts on this page:** 7
**Page:** 1

<div class="post-metadata">

### Author: ![lipesantosa](https://avatars.discourse-cdn.com/v4/letter/l/71c47a/32.png) [@lipesantosa](https://community.freefem.org/u/lipesantosa)
#### Post date: [June 21, 2025, 11:30am UTC](https://community.freefem.org/t/petrov-galerkin-least-squares/3979/1 "2025-06-21T11:30:05Z")

</div>

I’m having trouble with an Incompressible Navier Stokes script using the Petrov Galerkin method with Least Squares.

```
    matrix J = Jac(Vh2, Ph); // Matrix Jacobian
    matrix JGLS = GLS(Vh2, Ph);  

    J = J + JGLS; // Soma as matrizes

    // --- Vetor de resíduos ---
    real[int] R = Res(0, Vh2, Ph);

    // --- Resolver sistema linear ---
    set(J, solver=UMFPACK);
    real[int] Delta = J^-1 * (-R);

```

I’m having trouble with an Incompressible Navier Stokes script using the Petrov Galerkin method with Least Squares. In line real[int] R = Res(0, Vh2, Ph) the errors appears:

Error line number 82, in file teste.edp, before token )

current line = 82  
Compile error :  
line number :82, )  
error Compile error :  
line number :82, )  
code = 1 mpirank: 0

---

<div class="post-metadata">

### Author: ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [June 21, 2025, 1:28pm UTC](https://community.freefem.org/t/petrov-galerkin-least-squares/3979/2 "2025-06-21T13:28:44Z")

</div>

When you write the right-hand side you need to write `R=res(0,Xh)` where Xh is the space of velocity and pressure together.  
For that you can look at the example

> **[Tutorial to write a transient Stokes solver in matrix form](https://doc.freefem.org/tutorials/timeDependentStokes.html)**

(second script using matrix form with `varf`)

---

<div class="post-metadata">

### Author: ![lipesantosa](https://avatars.discourse-cdn.com/v4/letter/l/71c47a/32.png) [@lipesantosa](https://community.freefem.org/u/lipesantosa)
#### Post date: [June 21, 2025, 5:18pm UTC](https://community.freefem.org/t/petrov-galerkin-least-squares/3979/3 "2025-06-21T17:18:14Z")

</div>

Então, eu consegui resolver esse problema e arrumei a formatação das outras matrizes mas agora aparece o seguinte erro:

current line = 76  
Exec error : Check BilinearOperator N M  
– number :1  
Exec error : Check BilinearOperator N M  
– number :1  
err code 8 , mpirank 0

l

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/c/c51b60e3a8678e5fcc28f0d071f1b190e063b366.png)

---

<div class="post-metadata">

### Author: ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [June 21, 2025, 5:46pm UTC](https://community.freefem.org/t/petrov-galerkin-least-squares/3979/4 "2025-06-21T17:46:43Z")

</div>

Can you upload your full code?

---

<div class="post-metadata">

### Author: ![lipesantosa](https://avatars.discourse-cdn.com/v4/letter/l/71c47a/32.png) [@lipesantosa](https://community.freefem.org/u/lipesantosa)
#### Post date: [June 21, 2025, 5:51pm UTC](https://community.freefem.org/t/petrov-galerkin-least-squares/3979/5 "2025-06-21T17:51:09Z")

</div>

// — Parâmetros físicos e numéricos —  
real nu = 0.01; // Viscosidade cinemática  
real dt = 0.001; // Passo de tempo  
int nmax = 50; // Máx passos de tempo  
int imax = 20; // Máx iterações Newton por passo  
real tol = 1e-6; // Tolerância de Newton  
real deltaGLS = 1e-2; // Parâmetro de estabilização GLS (τ)

// — Malha: Cavidade quadrada —  
border a1(t=0,1){x=t; y=0; label=3;}  
border a2(t=0,1){x=1; y=t; label=4;}  
border a3(t=0,1){x=1-t; y=1; label=1;}  
border a4(t=0,1){x=0; y=1-t; label=2;}  
int n = 40;  
mesh Th = buildmesh(a1(n)+a2(n)+a3(n)+a4(n));

// — Espaços de funções —  
fespace Vh(Th, P1); // Velocidade escalar  
fespace Ph(Th, P1); // Pressão  
fespace Vh2(Th, [P1, P1]); // Velocidade vetorial (2 componentes)

// — Variáveis —  
Vh u1n=0, u2n=0, u1np=0, u2np=0, du1=0, du2=0, v1, v2;  
Ph pn=0, dp=0, q;  
func f1 = 0;  
func f2 = 0;

// — Loop no tempo —  
for (int nt=0; nt\<nmax; nt++) {  
cout \<\< "Tempo " \<\< nt \<\< endl;  
u1np = u1n; u2np = u2n; pn = pn;

```
int iter = 0;
bool convergiu = false;

while (iter < imax && !convergiu) {

    // --- Resíduo Total (R) ---
    varf Res([du1, du2, dp], [v1, v2, q], optimize=0) =
        int2d(Th)(
            (u1np/dt)*v1 + (u2np/dt)*v2
            + nu*(dx(u1np)*dx(v1) + dy(u1np)*dy(v1)
                + dx(u2np)*dx(v2) + dy(u2np)*dy(v2))
            + (u1np*dx(u1np) + u2np*dy(u1np))*v1
            + (u1np*dx(u2np) + u2np*dy(u2np))*v2
            - pn*(dx(v1) + dy(v2))
            - f1*v1 - f2*v2
        )
        + int2d(Th)( -(dx(u1np) + dy(u2np)) * q )
        + on(1, du1=1, du2=0) + on(2,3,4, du1=0, du2=0);

    // --- Jacobiano sem estabilização ---
    varf Jac([du1, du2, dp], [v1, v2, q]) = 
        int2d(Th)(
            (du1/dt)*v1 + (du2/dt)*v2
            + nu*(dx(du1)*dx(v1) + dy(du1)*dy(v1) + dx(du2)*dx(v2) + dy(du2)*dy(v2))
            + (u1np*dx(du1) + u2np*dy(du1))*v1
            + (u1np*dx(du2) + u2np*dy(du2))*v2
            - dp*(dx(v1) + dy(v2))
            + q*(dx(du1) + dy(du2))
        )+ on(1, du1=1, du2=0) + on(2,3,4, du1=0, du2=0);

    // --- Estabilização GLS (Petrov-Galerkin Least Squares) ---
    varf GLS([du1, du2, dp], [v1, v2, q]) = 
        int2d(Th)(
            deltaGLS * (
                ((du1/dt) + u1np*dx(du1) + u2np*dy(du1) + dx(dp)) * 
                ((v1/dt) + u1np*dx(v1) + u2np*dy(v1) + dx(q)) +
                ((du2/dt) + u1np*dx(du2) + u2np*dy(du2) + dy(dp)) * 
                ((v2/dt) + u1np*dx(v2) + u2np*dy(v2) + dy(q)) +
                (dx(du1)+dy(du2)) * (dx(v1)+dy(v2))
            )
        )+ on(1, du1=1, du2=0) + on(2,3,4, du1=0, du2=0);

    // --- Montagem da matriz global ---
    matrix J = Jac(Vh2, Vh2, solver=UMFPACK); // Corrigido para usar as variáveis corretas
    matrix JGLS = GLS(Vh2, Vh2);  

    J = J + JGLS; // Soma as matrizes

    // --- Vetor de resíduos ---
    real[int] R = Res(0, Vh2); // Vetor de resíduos 

    // --- Resolver sistema linear ---
    real[int] minusR = -1 * R;
    set(J, solver=UMFPACK);
    real[int] Delta = J^-1 * minusR;

    // --- Atualizar solução ---
    [du1[], du2[], dp[]] = Delta;
    u1np[] += du1[];
    u2np[] += du2[];
    pn[] += dp[];

    // --- Convergência ---
    real normDelta = sqrt(du1[].l2^2 + du2[].l2^2 + dp[].l2^2);
    cout << " Iter Newton " << iter << ", ||Delta|| = " << normDelta << endl;
    if (normDelta < tol) convergiu = true;

    iter++;
}

if (!convergiu)
    cout << "Newton NAO convergiu no tempo " << nt << endl;

// Atualiza para próximo passo
u1n[] = u1np[];
u2n[] = u2np[];

// Visualização
plot([u1n, u2n], cmm="Velocidade - t=" + nt, wait=0, value=true, fill=true);
//plot(pn, cmm="Pressao - t=" + nt, wait=0, value=true, fill=true);

```

}

---

<div class="post-metadata">

### Author: ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [June 21, 2025, 7:32pm UTC](https://community.freefem.org/t/petrov-galerkin-least-squares/3979/6 "2025-06-21T19:32:40Z")

</div>

Here is your code with corrections.  
[petrovgalerkin.edp](https://community.freefem.org/uploads/short-url/6CBZrMMuKf6vEtLY5K7hlmqTwfD.edp) (3.7 KB)  
I don’t apply the JGLS because I don’t know what you want to do with it.  
The corrections are as follows:

```auto
18c18
< fespace Vh(Th, P1); // Velocidade escalar
---
> fespace Vh(Th, P2); // Velocidade escalar
20c20,21
< fespace Vh2(Th, [P1, P1]); // Velocidade vetorial (2 componentes)
---
> fespace Xh(Th,[P2,P2,P1]);// space of velocity (2 components) + pressure
> Xh [wx,wy,wp];
46a48
> - 1.e-8*pn*q
50c52
< + on(1, du1=1, du2=0) + on(2,3,4, du1=0, du2=0);
---
> + on(1, du1=u1np-1., du2=u2np) + on(2,3,4, du1=u1np, du2=u2np);
60c62,63
< + q*(dx(du1) + dy(du2))
---
> - 1.e-8*dp*q
> - q*(dx(du1) + dy(du2))
76,77c79,80
< matrix J = Jac(Vh2, Vh2, solver=UMFPACK); // Corrigido para usar as variáveis corretas
< matrix JGLS = GLS(Vh2, Vh2);  
---
> matrix J = Jac(Xh, Xh, solver=UMFPACK); // Corrigido para usar as variáveis corretas
> matrix JGLS = GLS(Xh, Xh);  
79c82
< J = J + JGLS; // Soma as matrizes
---
> // J = J + JGLS; // Soma as matrizes
82c85
< real[int] R = Res(0, Vh2); // Vetor de resíduos 
---
> real[int] R = Res(0, Xh); // Vetor de resíduos 
85d87
< real[int] minusR = -1 * R;
87c89
< real[int] Delta = J^-1 * minusR;
---
> wx[] = J^-1 * R;
90c92,94
< [du1[], du2[], dp[]] = Delta;
---
> du1=-wx;
> du2=-wy;
> dp=-wp;

```

---

<div class="post-metadata">

### Author: ![lipesantosa](https://avatars.discourse-cdn.com/v4/letter/l/71c47a/32.png) [@lipesantosa](https://community.freefem.org/u/lipesantosa)
#### Post date: [June 21, 2025, 9:33pm UTC](https://community.freefem.org/t/petrov-galerkin-least-squares/3979/7 "2025-06-21T21:33:12Z")

</div>

Thank you very much for the corrections!! The term JGLS would be the Jacobian of the stabilization that comes from GLS. I’m trying to add it to the code to complete the method.