PETSc on 3D Cube with Periodic Boundary Conditions

General Discussion
1

Dear community,

I’m interested in using the HPDDM solver by PETSC on a 3D cube with periodic boundary conditions (PBC).

Following the 2D example diffusion-periodic-balanced-2d-PETSc.edp, I initially modified it on a basic 2D square so that the domain decomposition was balanced, and it works. Subsequently, I extended the same idea to the 3D cube, with the necessary modifications due to the increase in dimension. This is what I have implemented (I apologize in advance for the length of the code).

//  run with MPI:  ff-mpirun -np 4 script.edp

load "medit"
load "PETSc"                        // PETSc plugin
macro dimension()3// EOM            // 2D or 3D
include "macro_ddm.idp"             // additional DDM functions

int[int] labPeriodic = [1, 3, 2, 4, 5, 6];
int s = getARGV("-split", 1);                                   // refinement factor

string dataname = "part";
int[int] forder = [0];
macro Pk() P1, periodic=[[labPeriodic[0],x,z], [labPeriodic[1],x,z],[labPeriodic[2],y,z],[labPeriodic[3],y,z],[labPeriodic[4],x,y],[labPeriodic[5],x,y]]// EOM
macro PkNoPeriodic() P1// EOM

string Repo0 = "PBDDM_3D";
if(mpirank == 0){
    exec("mkdir "+Repo0);

    exec("mkdir "+Repo0+"/np_"+mpisize);
    exec("mkdir "+Repo0+"/np_"+mpisize+"/Mesh_0");
}

string Repo = Repo0+"/np_"+mpisize+"/Mesh_0";

real squaresizex = 10;
real squaresizey = 10;
real squaresizez = 10;

int n = getARGV("-global", 10);
mesh3 Th3 = cube(n, n, n, [squaresizex*x, squaresizey*y, squaresizez*z]);
mesh3 Th3Backup = Th3;

plot(Th3, cmm = "Starting Mesh");
Mat A;

fespace Wh(Th3, Pk);
fespace Ph(Th3, P0);
Ph part;
mesh3 Th3Extended = Th3;
{
    if(mpirank == 0)
    {
        partitionerSeq(part[], Th3, mpisize);
        fespace PhExtended(Th3Extended, P0);
        PhExtended partExtended;
        partExtended = part;

        for(int i = 0; i < mpisize; ++i) {
            Wh isBoundary;
            varf isBoundaryVarf([u], [v]) = on(labPeriodic, u = 1);
            isBoundary[] = isBoundaryVarf(0, Wh, tgv = 1);
            mesh3 Th3Local = trunc(Th3, abs(part - i) < 1e-6 && isBoundary > 0, label = 0);
            if(Th3Local.nt > 0){
                int j = 0;
                fespace Vh(Th3Local, P1);
                Vh u, uOn;
                int[int] touch(6);
                touch = 0;
                for(j = 0; j < labPeriodic.n; ++j) {
                    varf vOn(u, v) = on(labPeriodic[j], u = 1);
                    u[] = vOn(0, Vh, tgv = 1);

                    mesh3 Th3LocalMoved;
                    if(u[].linfty > 1e-12) {
                        Th3LocalMoved = trunc(Th3Local, u > 1e-12);
                        mesh3 Th3LocalMovedBackup = Th3LocalMoved;
                        // Face 1   --> 0
                        if(j == 0) {
                            Th3LocalMoved = movemesh(Th3LocalMoved, [x, y + squaresizey, z]);
                            touch[0] = 1;
                        }
                        // Face 3   --> 1
                        else if(j == 1) {
                            Th3LocalMoved = movemesh(Th3LocalMoved, [x, y - squaresizey, z]);
                            touch[1] = 1;
                        }
                        // Face 4   --> 2
                        else if(j == 2) {
                            Th3LocalMoved = movemesh(Th3LocalMoved, [x - squaresizex, y, z]);
                            touch[2] = 1;
                        }
                        // Face 2   --> 3
                        else if(j == 3) {
                            Th3LocalMoved = movemesh(Th3LocalMoved, [x + squaresizex, y, z]);
                            touch[3] = 1;
                        }
                        // Face 5   --> 4
                        else if(j == 4) {
                            Th3LocalMoved = movemesh(Th3LocalMoved, [x, y, z + squaresizez]);
                            touch[4] = 1;
                        }
                        // Face 6   --> 5
                        else if(j == 5) {
                            Th3LocalMoved = movemesh(Th3LocalMoved, [x, y, z - squaresizez]);
                            touch[5] = 1;
                        }
                    }
                    if(j == 5) // I think that this is for the corners // Touch[0] == 1 if there is the touch --> explain better
                    {
                        mesh3 Th3LocalMovedSecond;
                        mesh3 Th3LocalMod;
                        Vh u1, u2, u3;
                        // Edge 12  --> 02
                        if(touch[0] == 1 && touch[2] == 1) {
                            u1[] = 0;
                            u2[] = 0;
                            varf vOn0(u, v) = on(labPeriodic[0], u = 1);
                            varf vOn1(u, v) = on(labPeriodic[2], u = 1);
                            u1[] = vOn0(0, Vh, tgv = 1);
                            u2[] = vOn1(0, Vh, tgv = 1);
                            Th3LocalMod = trunc(Th3Local, u1*u2 > 1e-12);

                            Th3LocalMovedSecond = movemesh(Th3LocalMod, [x - squaresizex, y + squaresizey, z]);
                            Th3LocalMoved = Th3LocalMoved + Th3LocalMovedSecond;
                        }
                        // Edge 14  --> 03
                        if(touch[0] == 1 && touch[3] == 1) {
                            u1[] = 0;
                            u2[] = 0;
                            varf vOn0(u, v) = on(labPeriodic[0], u = 1);
                            varf vOn1(u, v) = on(labPeriodic[3], u = 1);
                            u1[] = vOn0(0, Vh, tgv = 1);
                            u2[] = vOn1(0, Vh, tgv = 1);
                            Th3LocalMod = trunc(Th3Local, u1*u2 > 1e-12);

                            Th3LocalMovedSecond = movemesh(Th3LocalMod, [x + squaresizex, y + squaresizey, z]);
                            Th3LocalMoved = Th3LocalMoved + Th3LocalMovedSecond;
                        }
                        // Edge 15  --> 04
                        if(touch[0] == 1 && touch[4] == 1) {
                            u1[] = 0;
                            u2[] = 0;
                            varf vOn0(u, v) = on(labPeriodic[0], u = 1);
                            varf vOn1(u, v) = on(labPeriodic[4], u = 1);
                            u1[] = vOn0(0, Vh, tgv = 1);
                            u2[] = vOn1(0, Vh, tgv = 1);
                            Th3LocalMod = trunc(Th3Local, u1*u2 > 1e-12);

                            Th3LocalMovedSecond = movemesh(Th3LocalMod, [x, y + squaresizey, z + squaresizez]);
                            Th3LocalMoved = Th3LocalMoved + Th3LocalMovedSecond;
                        }
                        // Edge 16  --> 05
                        if(touch[0] == 1 && touch[5] == 1) {
                            u1[] = 0;
                            u2[] = 0;
                            varf vOn0(u, v) = on(labPeriodic[0], u = 1);
                            varf vOn1(u, v) = on(labPeriodic[5], u = 1);
                            u1[] = vOn0(0, Vh, tgv = 1);
                            u2[] = vOn1(0, Vh, tgv = 1);
                            Th3LocalMod = trunc(Th3Local, u1*u2 > 1e-12);
                            
                            Th3LocalMovedSecond = movemesh(Th3LocalMod, [x, y + squaresizey, z - squaresizez]);
                            Th3LocalMoved = Th3LocalMoved + Th3LocalMovedSecond;
                        }
                        // Edge 32  --> 12
                        if(touch[1] == 1 && touch[2] == 1) {
                            u1[] = 0;
                            u2[] = 0;
                            varf vOn0(u, v) = on(labPeriodic[1], u = 1);
                            varf vOn1(u, v) = on(labPeriodic[2], u = 1);
                            u1[] = vOn0(0, Vh, tgv = 1);
                            u2[] = vOn1(0, Vh, tgv = 1);
                            Th3LocalMod = trunc(Th3Local, u1*u2 > 1e-12);
                            
                            Th3LocalMovedSecond = movemesh(Th3LocalMod, [x - squaresizex, y - squaresizey, z]);
                            Th3LocalMoved = Th3LocalMoved + Th3LocalMovedSecond;
                        }
                        // Edge 34  --> 13
                        if(touch[1] == 1 && touch[3] == 1) {
                            u1[] = 0;
                            u2[] = 0;
                            varf vOn0(u, v) = on(labPeriodic[1], u = 1);
                            varf vOn1(u, v) = on(labPeriodic[3], u = 1);
                            u1[] = vOn0(0, Vh, tgv = 1);
                            u2[] = vOn1(0, Vh, tgv = 1);
                            Th3LocalMod = trunc(Th3Local, u1*u2 > 1e-12);
                            
                            Th3LocalMovedSecond = movemesh(Th3LocalMod, [x + squaresizex, y - squaresizey, z]);
                            Th3LocalMoved = Th3LocalMoved + Th3LocalMovedSecond;
                        }
                        // Edge 35  --> 14
                        if(touch[1] == 1 && touch[4] == 1) {
                            u1[] = 0;
                            u2[] = 0;
                            varf vOn0(u, v) = on(labPeriodic[1], u = 1);
                            varf vOn1(u, v) = on(labPeriodic[4], u = 1);
                            u1[] = vOn0(0, Vh, tgv = 1);
                            u2[] = vOn1(0, Vh, tgv = 1);
                            Th3LocalMod = trunc(Th3Local, u1*u2 > 1e-12);
                            
                            Th3LocalMovedSecond = movemesh(Th3LocalMod, [x, y - squaresizey, z + squaresizez]);
                            Th3LocalMoved = Th3LocalMoved + Th3LocalMovedSecond;
                        }
                        // Edge 36  --> 15
                        if(touch[1] == 1 && touch[5] == 1) {
                            u1[] = 0;
                            u2[] = 0;
                            varf vOn0(u, v) = on(labPeriodic[1], u = 1);
                            varf vOn1(u, v) = on(labPeriodic[5], u = 1);
                            u1[] = vOn0(0, Vh, tgv = 1);
                            u2[] = vOn1(0, Vh, tgv = 1);
                            Th3LocalMod = trunc(Th3Local, u1*u2 > 1e-12);
                            
                            Th3LocalMovedSecond = movemesh(Th3LocalMod, [x, y - squaresizey, z - squaresizez]);
                            Th3LocalMoved = Th3LocalMoved + Th3LocalMovedSecond;
                        }
                        // Edge 25  --> 24
                        if(touch[2] == 1 && touch[4] == 1) {
                            u1[] = 0;
                            u2[] = 0;
                            varf vOn0(u, v) = on(labPeriodic[2], u = 1);
                            varf vOn1(u, v) = on(labPeriodic[4], u = 1);
                            u1[] = vOn0(0, Vh, tgv = 1);
                            u2[] = vOn1(0, Vh, tgv = 1);
                            Th3LocalMod = trunc(Th3Local, u1*u2 > 1e-12);
                            
                            Th3LocalMovedSecond = movemesh(Th3LocalMod, [x - squaresizex, y, z + squaresizez]);
                            Th3LocalMoved = Th3LocalMoved + Th3LocalMovedSecond;
                        }
                        // Edge 26  --> 25
                        if(touch[2] == 1 && touch[5] == 1) {
                            u1[] = 0;
                            u2[] = 0;
                            varf vOn0(u, v) = on(labPeriodic[2], u = 1);
                            varf vOn1(u, v) = on(labPeriodic[5], u = 1);
                            u1[] = vOn0(0, Vh, tgv = 1);
                            u2[] = vOn1(0, Vh, tgv = 1);
                            Th3LocalMod = trunc(Th3Local, u1*u2 > 1e-12);
                            
                            Th3LocalMovedSecond = movemesh(Th3LocalMod, [x - squaresizex, y, z - squaresizez]);
                            Th3LocalMoved = Th3LocalMoved + Th3LocalMovedSecond;
                        }
                        // Edge 45  --> 34
                        if(touch[3] == 1 && touch[4] == 1) {
                            u1[] = 0;
                            u2[] = 0;
                            varf vOn0(u, v) = on(labPeriodic[3], u = 1);
                            varf vOn1(u, v) = on(labPeriodic[4], u = 1);
                            u1[] = vOn0(0, Vh, tgv = 1);
                            u2[] = vOn1(0, Vh, tgv = 1);
                            Th3LocalMod = trunc(Th3Local, u1*u2 > 1e-12);

                            Th3LocalMovedSecond = movemesh(Th3LocalMod, [x + squaresizex, y, z + squaresizez]);
                            Th3LocalMoved = Th3LocalMoved + Th3LocalMovedSecond;
                        }
                        // Edge 46  --> 35
                        if(touch[3] == 1 && touch[5] == 1) {
                            u1[] = 0;
                            u2[] = 0;
                            varf vOn0(u, v) = on(labPeriodic[3], u = 1);
                            varf vOn1(u, v) = on(labPeriodic[5], u = 1);
                            u1[] = vOn0(0, Vh, tgv = 1);
                            u2[] = vOn1(0, Vh, tgv = 1);
                            Th3LocalMod = trunc(Th3Local, u1*u2 > 1e-12);
                            
                            Th3LocalMovedSecond = movemesh(Th3LocalMod, [x + squaresizex, y, z - squaresizez]);
                            Th3LocalMoved = Th3LocalMoved + Th3LocalMovedSecond;
                        }

                        mesh3 Th3LocalMovedThird;
                        // 125  --> 024
                        if(touch[0] == 1 && touch[2] == 1 && touch[4] == 1) {
                            u1[] = 0;
                            u2[] = 0;
                            u3[] = 0;
                            varf vOn0(u, v) = on(labPeriodic[0], u = 1);
                            varf vOn1(u, v) = on(labPeriodic[2], u = 1);
                            varf vOn2(u, v) = on(labPeriodic[4], u = 1);
                            u1[] = vOn0(0, Vh, tgv = 1);
                            u2[] = vOn1(0, Vh, tgv = 1);
                            u3[] = vOn2(0, Vh, tgv = 1);
                            Th3LocalMod = trunc(Th3Local, u1*u2*u3 > 1e-12);

                            Th3LocalMovedThird = movemesh(Th3LocalMod, [x - squaresizex, y + squaresizey, z + squaresizez]);
                            Th3LocalMoved = Th3LocalMoved + Th3LocalMovedThird;
                        }
                        // 145  --> 034
                        if(touch[0] == 1 && touch[3] == 1 && touch[4] == 1) {
                            u1[] = 0;
                            u2[] = 0;
                            u3[] = 0;
                            varf vOn0(u, v) = on(labPeriodic[0], u = 1);
                            varf vOn1(u, v) = on(labPeriodic[3], u = 1);
                            varf vOn2(u, v) = on(labPeriodic[4], u = 1);
                            u1[] = vOn0(0, Vh, tgv = 1);
                            u2[] = vOn1(0, Vh, tgv = 1);
                            u3[] = vOn2(0, Vh, tgv = 1);
                            Th3LocalMod = trunc(Th3Local, u1*u2*u3 > 1e-12);

                            Th3LocalMovedThird = movemesh(Th3LocalMod, [x + squaresizex, y + squaresizey, z + squaresizez]);
                            Th3LocalMoved = Th3LocalMoved + Th3LocalMovedThird;
                        }
                        // 345  --> 134
                        if(touch[1] == 1 && touch[3] == 1 && touch[4] == 1) {
                            u1[] = 0;
                            u2[] = 0;
                            u3[] = 0;
                            varf vOn0(u, v) = on(labPeriodic[1], u = 1);
                            varf vOn1(u, v) = on(labPeriodic[3], u = 1);
                            varf vOn2(u, v) = on(labPeriodic[4], u = 1);
                            u1[] = vOn0(0, Vh, tgv = 1);
                            u2[] = vOn1(0, Vh, tgv = 1);
                            u3[] = vOn2(0, Vh, tgv = 1);
                            Th3LocalMod = trunc(Th3Local, u1*u2*u3 > 1e-12);

                            Th3LocalMovedThird = movemesh(Th3LocalMod, [x + squaresizex, y - squaresizey, z + squaresizez]);
                            Th3LocalMoved = Th3LocalMoved + Th3LocalMovedThird;
                        }
                        //  235 --> 214
                        if(touch[2] == 1 && touch[1] == 1 && touch[4] == 1) {
                            u1[] = 0;
                            u2[] = 0;
                            u3[] = 0;
                            varf vOn0(u, v) = on(labPeriodic[2], u = 1);
                            varf vOn1(u, v) = on(labPeriodic[1], u = 1);
                            varf vOn2(u, v) = on(labPeriodic[4], u = 1);
                            u1[] = vOn0(0, Vh, tgv = 1);
                            u2[] = vOn1(0, Vh, tgv = 1);
                            u3[] = vOn2(0, Vh, tgv = 1);
                            Th3LocalMod = trunc(Th3Local, u1*u2*u3 > 1e-12);

                            Th3LocalMovedThird = movemesh(Th3LocalMod, [x - squaresizex, y - squaresizey, z + squaresizez]);
                            Th3LocalMoved = Th3LocalMoved + Th3LocalMovedThird;
                        }
                        // 126  --> 025
                        if(touch[0] == 1 && touch[2] == 1 && touch[5] == 1) {
                            u1[] = 0;
                            u2[] = 0;
                            u3[] = 0;
                            varf vOn0(u, v) = on(labPeriodic[0], u = 1);
                            varf vOn1(u, v) = on(labPeriodic[2], u = 1);
                            varf vOn2(u, v) = on(labPeriodic[5], u = 1);
                            u1[] = vOn0(0, Vh, tgv = 1);
                            u2[] = vOn1(0, Vh, tgv = 1);
                            u3[] = vOn2(0, Vh, tgv = 1);
                            Th3LocalMod = trunc(Th3Local, u1*u2*u3 > 1e-12);

                            Th3LocalMovedThird = movemesh(Th3LocalMod, [x - squaresizex, y + squaresizey, z - squaresizez]);
                            Th3LocalMoved = Th3LocalMoved + Th3LocalMovedThird;
                        }
                        // 146  --> 035
                        if(touch[0] == 1 && touch[3] == 1 && touch[5] == 1) {
                            u1[] = 0;
                            u2[] = 0;
                            u3[] = 0;
                            varf vOn0(u, v) = on(labPeriodic[0], u = 1);
                            varf vOn1(u, v) = on(labPeriodic[3], u = 1);
                            varf vOn2(u, v) = on(labPeriodic[5], u = 1);
                            u1[] = vOn0(0, Vh, tgv = 1);
                            u2[] = vOn1(0, Vh, tgv = 1);
                            u3[] = vOn2(0, Vh, tgv = 1);
                            Th3LocalMod = trunc(Th3Local, u1*u2*u3 > 1e-12);

                            Th3LocalMovedThird = movemesh(Th3LocalMod, [x + squaresizex, y + squaresizey, z - squaresizez]);
                            Th3LocalMoved = Th3LocalMoved + Th3LocalMovedThird;
                        }
                        // 346  --> 135
                        if(touch[1] == 1 && touch[3] == 1 && touch[5] == 1) {
                            u1[] = 0;
                            u2[] = 0;
                            u3[] = 0;
                            varf vOn0(u, v) = on(labPeriodic[1], u = 1);
                            varf vOn1(u, v) = on(labPeriodic[3], u = 1);
                            varf vOn2(u, v) = on(labPeriodic[5], u = 1);
                            u1[] = vOn0(0, Vh, tgv = 1);
                            u2[] = vOn1(0, Vh, tgv = 1);
                            u3[] = vOn2(0, Vh, tgv = 1);
                            Th3LocalMod = trunc(Th3Local, u1*u2*u3 > 1e-12);

                            Th3LocalMovedThird = movemesh(Th3LocalMod, [x + squaresizex, y - squaresizey, z - squaresizez]);
                            Th3LocalMoved = Th3LocalMoved + Th3LocalMovedThird;
                        }
                        // 263  --> 251
                        if(touch[2] == 1 && touch[5] == 1 && touch[1] == 1) {
                            u1[] = 0;
                            u2[] = 0;
                            u3[] = 0;
                            varf vOn0(u, v) = on(labPeriodic[2], u = 1);
                            varf vOn1(u, v) = on(labPeriodic[5], u = 1);
                            varf vOn2(u, v) = on(labPeriodic[1], u = 1);
                            u1[] = vOn0(0, Vh, tgv = 1);
                            u2[] = vOn1(0, Vh, tgv = 1);
                            u3[] = vOn2(0, Vh, tgv = 1);
                            Th3LocalMod = trunc(Th3Local, u1*u2*u3 > 1e-12);

                            Th3LocalMovedThird = movemesh(Th3LocalMod, [x - squaresizex, y - squaresizey, z - squaresizez]);
                            Th3LocalMoved = Th3LocalMoved + Th3LocalMovedThird;
                        }
                    }
                    if(Th3LocalMoved.nt > 0) {
                        real[int] backup = partExtended[];
                        fespace WhLocalMoved (Th3LocalMoved, P1);                    
                        WhLocalMoved isB;
                        isB = (z <= squaresizez + Th3LocalMoved.hmax && z >= -Th3LocalMoved.hmax && 
                        y <= squaresizey + Th3LocalMoved.hmax && y >= -Th3LocalMoved.hmax && 
                        x <= squaresizex+Th3LocalMoved.hmax && x >= -Th3LocalMoved.hmax);

                        try{
                            Th3LocalMoved = trunc(Th3LocalMoved, isB == 1);
                            cout << "On "+i+" we have trunced its mesh" << endl;
                        }
                        catch (...)
                        {
                            cout << "Prevention: on "+i+" its subdomain do not intersect boundary" << endl;
                        }
                        Th3LocalMoved = checkmesh(Th3LocalMoved, removeduplicate= 1, rebuildboundary = 1);
                        Th3Extended = Th3Extended + Th3LocalMoved;
                        cout << "Th3Extended done!" << endl;

                        partExtended = 0;
                        partExtended[](0:Th3Extended.nt - Th3LocalMoved.nt - 1) = backup;
                        partExtended[](Th3Extended.nt - Th3LocalMoved.nt:Th3Extended.nt - 1) = i;
                    }
                }
            }            
        }
        cout << "Changing labelling" << endl;
        Th3Extended = change(Th3Extended,
                                flabel = (x>=0 && x<=squaresizex && y>=0 && y<=squaresizey && z>=0 && z<=squaresizez)
                                            *((x == 0)*4 + (x == squaresizex)*2 +
                                                (y == 0)*1 + (y == squaresizey)*3 +
                                                (z == 0)*5 + (z == squaresizez)*6),
                                fregion = (x>=0 && x<=squaresizex && y>=0 && y<=squaresizey && z>=0 && z<=squaresizez)
                                            *((x == 0)*4 + (x == squaresizex)*2 +
                                                (y == 0)*1 + (y == squaresizey)*3 +
                                                (z == 0)*5 + (z == squaresizez)*6));
        
        savevtk(Repo+"/Th3Extended.vtu", Th3Extended, partExtended, dataname = dataname, order = forder);
        cout << "Here almost outside! #"+mpirank << endl;
        Th3 = Th3Extended;
        part = partExtended;
        plot(Th3, part, cmm="Th3 part After");
    }
}
cout << "*******************************" << endl << "Here! #"+mpirank << endl;
mpiBarrier(mpiCommWorld);

broadcast(processor(0), Th3);
broadcast(processor(0), Th3Extended);
broadcast(processor(0), part[]);
cout << "Broadcasting completed! #"+mpirank << endl;
plotD(Th3, part, cmm="Th3 before bwp");
mpiBarrier(mpiCommWorld);
cout << "Ready to buildWithPartitioning on #"+mpirank << endl;

// Partition of unity
real[int] DMat;
// Local-to-neighbors renumbering
int[int][int] intersectionMat;

buildWithPartitioning(Th3Extended, part[], s, intersectionMat, DMat, PkNoPeriodic, mpiCommWorld);
constructor(A, DMat.n, intersectionMat, DMat, communicator = mpiCommWorld);
savevtk(Repo+"/Th3ExtendedBuildWithPart.vtu", Th3Extended, part, dataname = dataname, order = forder);

cout << "Make buildWithPartitioning on #"+mpirank << endl;
mpiBarrier(mpiCommWorld);
fespace WhExtended(Th3Extended, PkNoPeriodic);
fespace WhNoPeriodic(Th3, PkNoPeriodic);
WhExtended isB2;
cout << "HEY! #"+mpirank << endl;
isB2 = (z <= squaresizez && z >= 0 && 
        y <= squaresizey && y >= 0 && 
        x <= squaresizex && x >= 0);
Th3 = trunc(Th3Extended, isB2 == 1);

savevtk(Repo+"/Th3AlmostFinal.vtu", Th3);
savemesh(Th3, "Th3AlmostFinal_"+mpirank+".mesh");
cout << "Th3 trunced #"+mpirank << endl;
mpiBarrier(mpiCommWorld);

matrix R = interpolate(WhNoPeriodic, WhExtended);
cout << "R made #"+mpirank << endl;
mpiBarrier(mpiCommWorld);
periodicity(R, intersectionMat, DMat);
mpiBarrier(mpiCommWorld);
cout << "Command 'periodicity' done on #"+mpirank << endl;
cout << "*******************************" << endl;

//Mat B;
//constructor(B, DMat.n, intersectionMat, DMat, communicator = mpiCommWorld);
mpiBarrier(mpiCommWorld);
cout << "Command 'constructor' done on #"+mpirank << endl;
mpiBarrier(mpiCommWorld);

plotDmesh(Th3, cmm="Th3 after bwp");

cout << "END! #"+mpirank << endl;
savevtk(Repo+"/Th3Final.vtu", Th3);

cout << "DOFS on #"+mpirank+"in P1-periodic "+Wh.ndof << endl;

Attempting to run the code, it seems to work intermittently, with instances of success and failure (e.g., with -np 8 -global 10). The output is as follows:

Exec error : Periodic 3d:  the both number of face is not the same 
   -- number :1
 err code 8 ,  mpirank 3

I tried visualizing the results on Paraview, and indeed, for some processors, the domain partition is not periodic (usually a missing triangle), even though the overall process seems correct.

How should I proceed?

Thank you in advance for your help.

Best regards,

2

Is there a particular reason why you don’t use the simpler approach of defining macro ThPeriodicity(), as in https://github.com/FreeFem/FreeFem-sources/blob/develop/examples/hpddm/diffusion-periodic-2d-PETSc.edp#L22-L24?

3

Dear @prj,
Thanks for the reply. I’ve already seen this option, but in 3D it does a non-balanced Domain Decomposition, since it considers an “external shell” and then the inside is decomposed with the typical overlapping DDM. Certainly, this works, but in terms of balancing the DOFs is not convenient.

To see this, I do a very quick modification of the previous code, applying this DDM instead of mine with -np 4 -global 25 and printing out the final DOFs of the periodic space it gives me:

DOFS on #1in P1-periodic 4224
DOFS on #3in P1-periodic 4377
DOFS on #2in P1-periodic 4296
DOFS on #0in P1-periodic 11529

Lastly, increasing the -np, the DOFS on #0in P1-periodic remains the same.

4

Right, so, as you increase the size of the problem and the number of processes, this should even out. Also, if you are only interested in standard geometries, you can still do a somewhat balanced decomposition in a much simpler manner. Here is one example.

assert(mpisize > 2);
int n = 5;
mesh Th = square(4 * n, 4 * n);
fespace Vh(Th, P1);
varf vPb(u, v) = on(1, 2, 3, 4, u = 1.0);
Vh u;
u[] = vPb(0, Vh, tgv = -1);
fespace Ph(Th, P0);
Ph part = u;
AddLayers(Th, part[], 1, u[]);
plot(part);
plot(u);
int[int] n2oInner, n2oOuter;
mesh ThInner = trunc(Th, u < 1.0E-2, new2old = n2oInner);
mesh ThOuter = trunc(Th, u > 1.0E-2, new2old = n2oOuter);
include "macro_ddm.idp"
fespace PhInner(ThInner, P0);
PhInner partInner;
partitionerSeq(partInner[], ThInner, mpisize - 2);
partitionerPar(partInner[], ThInner, mpiCommWorld, mpisize - 2);
fespace PhOuter(ThOuter, P0);
PhOuter partOuter = (x < 0.25 && y < 0.25) || (x > 0.75 && y > 0.75) || (x < 0.25 && y > 0.75) || (x > 0.75 && y < 0.25) ? 0.0 : 1.0;
part[] = 0;
partInner[] += 2.0;
part[](n2oInner) = partInner[];
part[](n2oOuter) = partOuter[];
plot(part);

You could do the same for the cube. Then, instead of having only a single process in charge of the periodic faces, you could have two (or more, for example four: one for all corners, one for the bottom and top faces, one for left and right faces, and one for front and back faces).

5

Thanks a lot @prj,
Now I will work on your code in order to understand all the passages.

This is theoretically true, but the balancing is strictly dependent on the number of processor, the dimension (obviously), but also on the regularity of the mesh.
Indeed, my goal was to find a consistent method for performing a balanced DDM on a cube without relying on the regularity of the underlying mesh, as I intend to apply this method to an adapted periodic mesh.

6

Early optimisation is the root of all evil. The current load balancing may not be the best, but as you change your mesh, it may be perfectly fine. So I would not bother with this until it becomes a clear bottleneck.

7

Thank you so much, I understand.
I will work with this defined DDM and hope to obtain good results in a low computational time.