Macro ThN2O for mixed vectorial and scalar finite element space

Hello everyone!

I’ve recently watched the tutorials from professor Pierre Jolivet on Youtube about parallel computing in FreeFem. I have a code that works fine sequentially and I intend to swicht it to parallel. In the sequential script I have one scalar and one vectorial finite element space, both associated with an initial mesh Th. I need to go from local to global for both vectorial and scalar fields. Can it be obtained by using macro ThN2O for both cases all in the same code?

Thanks!

Yes, why would that matter?

As I’m using macro def(i) [i i#y] and macro init(i) [i i] for the vectorial field, the reduction from local to global works fine for the vectorial field. However, I couldn’t do the same for the scalar field.

Please share a minimal working example. N2O is a mesh information, it has nothing to do with a finite element space whatsoever, so as long as you use it appropriately afterwards, you can do whatever you want with it.

Ok, I see! Thanks for your reply and for the information. I’ll try to make a short code those ideas to share with you. Thanks again!

Here is an example. Of course, you are free to use any kind of vectorial finite elements, I just used Morley here.

load "Morley"
include "macro_ddm.idp"
load "PETSc"

mesh ThGlobal = square(100, 100);
mesh Th = square(100, 100);
fespace WhGlobal(ThGlobal, P2Morley);
fespace VhGlobal(ThGlobal, P1);
WhGlobal [uG, uGx, uGy] = [x^2, y^2, x^2+y^2];
VhGlobal pG = x + y;
int[int] n2o;
macro ThN2O()n2o//
buildDmesh(Th);
fespace Wh(Th, P2Morley);
fespace Vh(Th, P1);
Wh [u, ux, uy] = [x^2, y^2, x^2+y^2];
Vh p = x + y;
int[int] restWh = restrict(Wh, WhGlobal, n2o);
int[int] restVh = restrict(Vh, VhGlobal, n2o);
{
    Mat A;
    createMat(Th, A, P2); // trick => use P2, not P2Morley!
    real[int] tmp;
    ChangeNumbering(A, u[], tmp);
    ChangeNumbering(A, u[], tmp, inverse = true);
    WhGlobal [reducex, reducey, reducez];
    WhGlobal [beforex, beforey, beforez];
    for[i, v : restWh] beforex[][v] = u[][i];
    mpiReduce(beforex[], reducex[], processor(0, mpiCommWorld), mpiSUM);
    uG[] -= reducex[];
    if(mpirank == 0) assert(uG[].l2 < 1.0e-12);
}
{
    Mat A;
    createMat(Th, A, P1);
    real[int] tmp;
    ChangeNumbering(A, p[], tmp);
    ChangeNumbering(A, p[], tmp, inverse = true);
    VhGlobal reduce;
    VhGlobal before;
    before[] = 0;
    for[i, v : restVh] before[][v] = p[][i];
    mpiReduce(before[], reduce[], processor(0, mpiCommWorld), mpiSUM);
    pG[] -= reduce[];
    if(mpirank == 0) assert(pG[].l2 < 1.0e-12);
}

Thank you so much for sharing this example. It’ll certainly help me!

Best regards!

If I change Morley to Pk = [P1,P1,P1], I met Compile error : Invalide array size for vectorial fespace function. How do I solve this problem? Thanks!

load "PETSc"                        // PETSc plugin
macro dimension()2// EOM            // 2D or 3D
include "macro_ddm.idp"             // additional DDM functions
macro def(i)[i, i#B, i#C]//
macro init(i)[i, i, i]//

mesh ThGlobal = square(100, 100);
mesh Th = square(100, 100);
func Pk = [P1,P1,P1];
fespace WhGlobal(ThGlobal, Pk);
fespace VhGlobal(ThGlobal, P1);
WhGlobal [uG, uGx, uGy] = [x^2, y^2, x^2+y^2];
VhGlobal pG = x + y;
int[int] n2o;
macro ThN2O()n2o//
buildDmesh(Th);
fespace Wh(Th, Pk);
fespace Vh(Th, P1);
Wh [u, ux, uy] = [x^2, y^2, x^2+y^2];
Vh p = x + y;
int[int] restWh = restrict(Wh, WhGlobal, n2o);
int[int] restVh = restrict(Vh, VhGlobal, n2o);
{
    Mat A;
    createMat(Th, A, Pk); // trick => use P2, not P2Morley!
    real[int] tmp;
    ChangeNumbering(A, u[], tmp);
    ChangeNumbering(A, u[], tmp, inverse = true);
    WhGlobal [reducex, reducey, reducez];
    WhGlobal [beforex, beforey, beforez];
    for[i, v : restWh] beforex[][v] = u[][i];
    mpiReduce(beforex[], reducex[], processor(0, mpiCommWorld), mpiSUM);
    uG[] -= reducex[];
    if(mpirank == 0) assert(uG[].l2 < 1.0e-12);
}

{
    Mat A;
    createMat(Th, A, P1);
    real[int] tmp;
    ChangeNumbering(A, p[], tmp);
    ChangeNumbering(A, p[], tmp, inverse = true);
    VhGlobal reduce;
    VhGlobal before;
    before[] = 0;
    for[i, v : restVh] before[][v] = p[][i];
    mpiReduce(before[], reduce[], processor(0, mpiCommWorld), mpiSUM);
    pG[] -= reduce[];
    if(mpirank == 0) assert(pG[].l2 < 1.0e-12);
}```

402 @ fespace WhPartPrivate(ThTab[ 1 - 1], ThPkPart P1 );
403 @ WhPartPrivate def(func2vec) [func2vec, func2vecB, func2vecC]; the array size must be 1 not 3
Invalid array size for vectorial fespace function
current line = 403 mpirank 0 / 2
Compile error : Invalid array size for vectorial fespace function
line number :403, ;
error Compile error : Invalid array size for vectorial fespace function
line number :403, ;
code = 1 mpirank: 0

--- wrong.edp	2023-07-20 15:41:15
+++ right.edp	2023-07-20 15:41:26
@@ -1,8 +1,6 @@
 load "PETSc"                        // PETSc plugin
 macro dimension()2// EOM            // 2D or 3D
 include "macro_ddm.idp"             // additional DDM functions
-macro def(i)[i, i#B, i#C]//
-macro init(i)[i, i, i]//
 
 mesh ThGlobal = square(100, 100);
 mesh Th = square(100, 100);
@@ -22,6 +20,8 @@ int[int] restVh = restrict(Vh, VhGlobal, n2o);
 int[int] restVh = restrict(Vh, VhGlobal, n2o);
 {
     Mat A;
+    macro def(i)[i, i#B, i#C]//
+    macro init(i)[i, i, i]//
     createMat(Th, A, Pk); // trick => use P2, not P2Morley!
     real[int] tmp;
     ChangeNumbering(A, u[], tmp);

It works. Thank you very much!

I am writing to ask how to write a parallel solver for a transient problem with more than one finite space. For the case “A large fluid problem” in FreeFem documentation, we have Pk=[P2,P2,P1] for [u1,u2,p] and P1 for T, k. We can use macro def(i)[i, i#B, i#C]//and macro init(i)[i, i, i]// for [u1,u2,p]. It is similar to the above implementation. But now at each time step, we need to update T using u1and u2 in the previous time step, which are invisible in the code delimited by the second pair of curly braces.
{
Mat A;
createMat(Th, A, Pk); //for [u1,u2,p]

}
{
Mat B;
createMat(Th, B, P1); //for T

}

Thanks!

Why are you declaring A and B inside the brackets?

Here is the sample code of my simulation.
prll.edp (3.4 KB)
I do not know how to handle the following macros properly.

macro def(i)[i, i#B, i#C, i#D, i#E, i#F, i#G]//
macro init(i)[i, i, i, i, i, i, i]//
Mat A;
createMat(Th,A,Pk)
set(A, sparams = “-pc_type lu”);
A = vSt(Wh, Wh);
macro def(i)[i, i#B, i#C]//
macro init(i)[i, i, i]//
Mat B;
createMat(Th,B,Pk)
set(B, sparams = “-pc_type lu”);

Why did you remove all brackets?

I removed the brackets, I obtained the following errors.

Blockquote
current line = 96 mpirank 0 / 6
Compile error : The macro already exists, sorry
line number :96, )
error Compile error : The macro already exists, sorry
line number :96, )
code = 1 mpirank: 0

Then I used the brackets, such as

{
macro def(i)[i, i#B, i#C, i#D, i#E, i#F, i#G]//
macro init(i)[i, i, i, i, i, i, i]//
Mat A;
createMat(Th,A,Pk)
set(A, sparams = “-pc_type lu”);
A = vSt(Wh, Wh);
}
{
macro def(i)[i, i#B, i#C]//
macro init(i)[i, i, i]//
Mat B;
createMat(Th,B,Pk)
set(B, sparams = “-pc_type lu”);
}

I got the following errors

396 @ fespace WhPartPrivate(ThTab[ 1 - 1], ThPkPart Pk );
397 @ WhPartPrivate def(func2vec) [func2vec, func2vecB, func2vecC]; the array size must be 7 not 3
Error line number 397, in file macro: partitionPrivate in /opt/freefem/4.9/mpich3.3rc1-gnu9.4.0/lib/ff++/4.9/idp/macro_ddm.idp, before token ;
Invalide array size for vectorial fespace function
current line = 397 mpirank 0 / 6
Compile error : Invalide array size for vectorial fespace function
line number :397, ;
error Compile error : Invalide array size for vectorial fespace function
line number :397, ;
code = 1 mpirank: 0

Please post a minimal example reproducing the error, not something with 300+ lines of code.

Attached is a small example.
example.edp (3.0 KB)

You are using the same Pk in both calls to createMat(), but different def macros, why?
Also, you are still defining A and B inside brackets, which is wrong.

Yes. It should be

createMat(Th,B,Ps)

I commented

// macro def(i)[i, i#B, i#C]//
// macro init(i)[i, i, i]//

then I met problems

396 @ fespace WhPartPrivate(ThTab[ 1 - 1], ThPkPart Ps );
397 @ WhPartPrivate def(func2vec) [func2vec, func2vecB, func2vecC, func2vecD, func2vecE, func2vecF, func2vecG]; the array size must be 3 not 7
Error line number 397, in file macro: partitionPrivate in /opt/freefem/4.9/mpich3.3rc1-gnu9.4.0/lib/ff++/4.9/idp/macro_ddm.idp, before token ;
Invalide array size for vectorial fespace function
current line = 397 mpirank 0 / 6
Compile error : Invalide array size for vectorial fespace function
line number :397, ;
error Compile error : Invalide array size for vectorial fespace function
line number :397, ;
code = 1 mpirank: 0

I tried to use brackets to avoid the duplicate definition of macro def(i)[i, i#B, i#C]//and macro init(i)[i, i, i]// but it did not work.

Please send the updated code.