What is the purpose of O = Loc in set

General Discussion
1

In the following code (Helmholtz-2d-PETSc-complex.edp) :
What is the purpose of the argument: " O = Loc" in the command set(…)?

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

load "PETSc-complex"                // PETSc plugin
include "macro_ddm.idp"             // additional DDM functions

macro grad(u)[dx(u), dy(u)]// EOM   // two-dimensional gradient
func Pk = P1;                       // finite element space

mesh Th = square(getARGV("-global", 40), getARGV("-global", 40)); // global mesh
Mat<complex> A;
int s = getARGV("-split", 1);
macro ThRefinementFactor()s//
MatCreate(Th, A, Pk);

fespace Wh(Th, Pk);                 // local finite element space

func real wedge(real a, real b) {
    if(y < 0.4 + 0.1 * 0.75 * x)
        return 2;
    else if(y < 0.8 - 0.2 * 0.75 * x)
        return 1.5;
    else
        return 3;
}
real omega = 2 * pi * 5;
real xi = 0, yi = 0;
func f = 80 * 100 * s * exp(-20 * 100 * s * ((x-0.5-xi)^2 + (y-0.25-yi)^2));
complex[int] rhs(Wh.ndof);                  // local right-hand side
matrix<complex> Loc;                        // local operator
{                                           // local weak form
    fespace Ph(Th, P0);
    Ph val = wedge(x, y);
    Ph k = omega / val;
    varf vPb(u, v) = int2d(Th)(-k^2 * u * v + (grad(u)' * grad(v))) + int1d(Th, -111111)(1i * k * u * v) + int1d(Th, 2)(1i * k * u * v) + int2d(Th)(f * v) + on(1, u = 0.0);
    Loc = vPb(Wh, Wh, tgv = -1);
    rhs = vPb(0, Wh, tgv = -1);
}

A = Loc;
if(usedARGV("-optimized_Schwarz") != -1)
    set(A, sparams = "-ksp_view -sub_pc_type lu", O = Loc);
else
    set(A, sparams = "-ksp_view -pc_type lu");
Wh<complex> u;                      // local solution

u[] = A^-1 * rhs;

macro def(u)u//
plotMPI(Th, u, Pk, def, complex, cmm = "Global solution");

int M = 4;
complex[int, int] RHS(rhs.n, M);
RHS(:, 0) = rhs;
for(int i = 1; i < M; ++i) {
    xi += 0.1;
    yi += 0.1;
    varf vPb(u, v) = int2d(Th)(f * v) + on(1, u = 0.0);
    RHS(:, i) = vPb(0, Wh, tgv = -1);
}
set(A, sparams = "-ksp_type hpddm -ksp_converged_reason -ksp_hpddm_type bgmres");
complex[int, int] B = A^-1 * RHS;
for(int i = 0; i < M; ++i) {
    u[] = B(:, i);
    macro params()cmm = "Global solution #" + i, wait = 1, fill = 1, value = 1, dim = 3//
    plotMPI(Th, u, Pk, def, complex, params);
}
complex[int] C(B.n * B.m);
C = A^-1 * RHS.asarray;
for(int i = 0; i < M; ++i) {
    u[] = B(:, i);
    u[] -= C(i * Wh.ndof:(i + 1) * Wh.ndof - 1);
    macro params()cmm = "Error #" + i, wait = 1, fill = 1, value = 1, dim = 3//
    plotMPI(Th, u, Pk, def, complex, params);
}
2

It is if you want to use an optimized Schwarz method instead of plain (restricted) additive Schwarz method.

3

I see. In the following code, at the set(...) line (right at the end), when I remove the “O = Loc” option and run the code on 2 processes, the two sub-matrices have sizes of 124x124 and 128x128, respectively. However, when I put the “O = Loc” option back in, both matrices end up with the same size: 128x128. So, I gather that the command increases the size of one of the matrices, but I don’t know the principle behind this addition.

// Version with ORAS
// K - w²M
// macro partitioner()scotch// EOM
load "PETSc"                // PETSc plugin
include "macro_ddm.idp"             // additional DDM functions
load "hpddm"                        // HPDDM plugin


mesh Th = square(60, 1, [0+300*x, 0+100*y], flags = 1); 

savemesh(Th,"Th.mesh");


macro dimension()2// EOM            // 2D or 3D
func Pk = [P1, P1];                 // finite element space

macro def(i)[i, i#B]// EOM          // vector field definition
macro init(i)[i, i]// EOM           // vector field initialization

real young = 1.0;
real nu = 0.0;
real omega = 1; 
real rho = 1;
real mu = young / (2.0 * (1.0 + nu));

real lambdax = (young * nu) / ((1.0 + nu) * (1.0 - 2.0 * nu));
macro epsilonx(u1,u2) [dx(u1),dy(u2),(dy(u1)+dx(u2))/sqrt(2.)]  // EOM 
macro divx(u,v) (dx(u) + dy(v)) // EOM 

real cp = sqrt((lambdax + 2.*mu) / rho);   // vitesse des ondes P
real cs = sqrt(mu / rho);                  // vitesse des ondes S
real fakeInterface = -111111;
real robinCoef = getARGV("-robinCoef", 1e0); //2.293119e1
// real label = (abs(fakeInterface) + 1) * 100;

fespace Qh(Th, Pk);
// Decompose the mesh



macro ThOverlap()1//

DmeshCreate(Th);


// plotDmesh(Th, cmm= "M");
Mat A;
MatCreate(Th, A, Pk);

real f = -1;
real[int] rhs(Qh.ndof);                  // local right-hand side

  

varf vaLoc0([u,v],[uu,vv]) = int2d(Th)(
            lambdax* ( divx(u,v) *divx(uu,vv))    
           + 2.*mu*(epsilonx(u,v)' * epsilonx(uu,vv)))
           - int2d(Th)( rho * (u * uu + v * vv) *omega*omega );

varf vaLoc1([u,v],[uu,vv]) = int2d(Th)(
            lambdax* ( divx(u,v) *divx(uu,vv))    
           + 2.*mu*(epsilonx(u,v)' * epsilonx(uu,vv)))
           - int2d(Th)( rho * (u * uu + v * vv) *omega*omega )
            + int1d(Th, fakeInterface)(
                robinCoef*rho*cp * (u*N.x + v*N.y) * (uu*N.x + vv*N.y)
            + robinCoef*rho*cs * (u*N.y - v*N.x) * (uu*N.y - vv*N.x)
            );

varf vaLoc2([u,v],[uu,vv]) = 
             int1d(Th, fakeInterface)(
                robinCoef*rho*cp * (u*N.x + v*N.y) * (uu*N.x + vv*N.y)
            + robinCoef*rho*cs * (u*N.y - v*N.x) * (uu*N.y - vv*N.x)
            );

matrix Loc1;   // local operator

Loc1 = vaLoc1(Qh, Qh, sym = 0, tgv = -2);

matrix Loc0;   // local operator

Loc0 = vaLoc0(Qh, Qh, sym = 0, tgv = -2);



varf vRHS([u,v],[uu,vv]) = int2d(Th)( f * (vv) );
rhs = vRHS(0, Qh, tgv = -2);



A = Loc0;   
// A = vaLoc0(Qh, Qh, sym = 0, tgv = -2);


string viewFile = "mat_view_rank_" + mpirank + ".txt";


exec("rm -f " + viewFile);


set(A, sparams = "-ksp_view -pc_type asm -pc_asm_type restrict -sub_pc_type cholesky -ksp_monitor_true_residual -ksp_rtol 1e-8 -mat_view ascii:" + viewFile, O = Loc1);

Qh<real> [u, v]; 

u[] = A^-1 * rhs;
4

In one case, the overlap is built geometrically, in the other, algebraically.