Error operator in Petrov-Galerkin Least-Squares

General Discussion
1

I am getting an error in my script, “error operator,” whenever I try to subtract with a variable or function. In this case, the variables would be Vh fx, fy, fT. This happens whenever I have a problem to solve.

        Vh fx = (u1exact*dx(u1exact) + u2exact*dy(u1exact)) + dx(pexact)
                - nu*(dxx(u1exact) + dyy(u1exact));
        Vh fy = (u1exact*dx(u2exact) + u2exact*dy(u2exact)) + dy(pexact)
                - nu*(dxx(u2exact) + dyy(u2exact));
        ThT fT = (u1exact*dx(Texact) + u2exact*dy(Texact))
                - kappa*(dxx(Texact) + dyy(Texact));          

// ===== GLS Parte 2: convecção (u·∇u , u·∇v) =====

          + int2d(Th)(

              (tau) \* (

                  (u1old\*dx(u1) + u2old\*dy(u1))\*(u1old\*dx(v1) + u2old\*dy(v1))

              + (u1old\*dx(u2) + u2old\*dy(u2))\*(u1old\*dx(v2) + u2old\*dy(v2)) 

              - fx\*(u1old\*dx(v1) + u2old\*dy(v1)) 

              - fy\*(u1old\*dx(v2) + u2old\*dy(v2))

            )

        )

Erro:Error line number 126, in file PGLS.edp, before token -

current line = 125
Compile error :
line number :126, -
error Compile error :
line number :126, -
code = 1 mpirank: 0

2

To help you, I need the full script

3 
// ---- Malha ----

*border* a1(*t*=0,1){x=t; y=0; label=2;}

*border* a2(*t*=0,1){x=1; y=t; label=2;}

*border* a3(*t*=0,1){x=1-t; y=1; label=1;}

*border* a4(*t*=0,1){x=0; y=1-t; label=2;}

*mesh* Th = buildmesh(a1(n)+a2(n)+a3(n)+a4(n));



// ---- Espaços ----

*fespace* Vh(Th, P1);

*fespace* Qh(Th, P1);

*fespace* ThT(Th, P1);



Vh u1, u2, u1old, u2old, v1, v2, ht;

Qh p, q;

Vh du1, du2;

ThT T, Told, theta, dTn, Pe;



// ---- Parâmetros ----

*real* Re = 100;      // Reynolds

*real* nu = 1.0/Re;   // Viscosidade

*real* dt = 0.001;     // Passo de tempo

*int*  nmax = 10;     // Passos de tempo

*int*  itmax = 20;    // Iterações de Picard

*real* kappa = 0.01;  // Condutividade

*real* pi = acos(-1);



// ---- Solução exata (MMS) ----

*func* u1exactFunc =  pi \* cos(pi\*y) \* sin(pi\*x);

*func* u2exactFunc = -pi \* cos(pi\*x) \* sin(pi\*y);

*func* pexactFunc  =      cos(pi\*x) \* cos(pi\*y);

*func* TexactFunc  =      sin(pi\*x) \* sin(pi\*y);



// ---- Inicialização no estado exato ----

u1old = u1exactFunc;

u2old = u2exactFunc;

Told  = TexactFunc;



// ---- Estabilização ----

*macro* umod() sqrt(u1old^2 + u2old^2)  //

ht = hTriangle;

*macro* tau() (1.0e-3 / sqrt((2/dt)^2 + (2.0\*umod/ht)^2 + (4.0\*nu/(ht^2))^2)) //

*macro* Peclet() max(1e-8, 1./3.\*umod\*ht/(2\*kappa)) //

*macro* tauT() ht^2/(2\*Peclet\*kappa) \* sqrt(Peclet^2/(36 + Peclet^2)) //



// ---- Loop no tempo ----

for (*int* nTimeStep = 0; nTimeStep < nmax; ++nTimeStep) {

    cout << "\\n>>> Passo = " << nTimeStep << ", t = " << nTimeStep\*dt << endl;



    // Projeta solução exata

    Vh u1exact = u1exactFunc;

    Vh u2exact = u2exactFunc;

    Qh pexact  = pexactFunc;

    ThT Texact = TexactFunc;



    // Termos-fonte consistentes

    Vh fx = (u1exact\*dx(u1exact) + u2exact\*dy(u1exact)) + dx(pexact)

            - nu\*(dxx(u1exact) + dyy(u1exact));

    Vh fy = (u1exact\*dx(u2exact) + u2exact\*dy(u2exact)) + dy(pexact)

            - nu\*(dxx(u2exact) + dyy(u2exact));

    ThT fT = (u1exact\*dx(Texact) + u2exact\*dy(Texact))

            - kappa\*(dxx(Texact) + dyy(Texact));



    *func* fxfunc = fx;

    *func* fyfunc = fy;

    *func* fTfunc = fT;



    // estado inicial Picard

    u1\[\] = u1old\[\];

    u2\[\] = u2old\[\];



    // auxiliares (u/dt)

    Vh u1DivDt = u1/dt, u2DivDt = u2/dt;

    Vh u1oldDivDt = u1old/dt, u2oldDivDt = u2old/dt;



    // ===== Navier–Stokes com GLS EXPANDIDO =====

    for (*int* it=0; it<itmax; ++it) {

        cout << "   Picard it = " << it << endl;



        *problem* NSE(\[u1, u2, p\], \[v1, v2, q\]) =

        // --- Galerkin padrão ---

            int2d(Th)(

                    (1.0/dt)\*(u1\*v1 + u2\*v2)

                + nu\*(dx(u1)\*dx(v1) + dy(u1)\*dy(v1) + dx(u2)\*dx(v2) + dy(u2)\*dy(v2))

                + u1old\*dx(u1)\*v1 + u2old\*dy(u1)\*v1

                + u1old\*dx(u2)\*v2 + u2old\*dy(u2)\*v2

                - p\*(dx(v1) + dy(v2))

                - q\*(dx(u1) + dy(u2))

                - 1e-6\*p\*q

            )



            // --- GLS (resíduo convectivo + grad p - f) x (u·∇v - ∇q) ---

            + int2d(Th)(

                (tau)

                \* (

                    // Resíduo componente 1 (R1)

                    ( (u1)/dt + u1old\*dx(u1) + u2old\*dy(u1) + dx(p))

                    -nu \* ( u1old\*dxx(v1) + u2old\*dyy(v1) - dx(q) )



                // Resíduo componente 2 (R2)

                    + ( (u2)/dt + u1old\*dx(u2) + u2old\*dy(u2) + dy(p))

                    - nu \* ( u1old\*dxx(v2) + u2old\*dyy(v2) - dy(q) )

                )

            )

            - int2d(Th)(

                (tau) \* (u1old/dt + u2/dt)

            )



            // --- BCs: ajuste para MMS (u_exact em toda a borda recomendado) ---

                        // RHS do passo anterior + forçantes (Galerkin)

        - int2d(Th)((1.0/dt)\*(u1old\*v1 + u2old\*v2))

        - int2d(Th)( fx\*v1 + fy\*v2 )

        + on(1,2, u1 = u1exact, u2 = u2exact)

        ;

    }