Petrov-Galerkin Least-Squares

General Discussion
1

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

2

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

(second script using matrix form with varf)

1 Like
3

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
image826×325 16.2 KB

4

Can you upload your full code?

5

// — 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);

}

6

Here is your code with corrections.
petrovgalerkin.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:

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;
1 Like
7

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.