1// run with MPI: ff-mpirun -np 4 script.edp 2// NBPROC 4 3 4load "hpddm" 5macro dimension()2// EOM // 2D or 3D 6include "macro_ddm.idp" // additional DDM functions 7 8macro def(i)i// EOM // scalar field definition 9macro init(i)i// EOM // scalar field initialization 10macro grad(u)[dx(u), dy(u)]// EOM // two-dimensional gradient 11func Pk = P2; // finite element space 12 13real k = getARGV("-waven", 40.0); 14func f = 100 * exp(-10 * k * ((x-0.5)^2 + (y-0.5)^2)); 15 16real lambda = 2 * pi / k; 17 18real epsilonA = 0; 19real epsilonE = 0; 20 21int Dirichlet = 1; 22int Robin = 2; 23 24int s = getARGV("-split", 3); 25real nw = 15.0 / s; 26meshN ThCoarse = square(nw * 1/lambda, nw * 1/lambda); 27 28{ 29 int[int] chlab = [1, Dirichlet, 2, Dirichlet, 3, Dirichlet, 4, Dirichlet]; 30 ThCoarse = change(ThCoarse, refe = chlab); 31} 32 33meshN Th; 34fespace Wh(Th, Pk); 35real[int] part; 36fespace PhCoarse(ThCoarse, P0); 37PhCoarse partCoarse; 38if(mpirank == 0) 39 partitionerSeq(partCoarse[], ThCoarse, mpisize); 40partitionerPar(partCoarse[], ThCoarse, mpiCommWorld, mpisize); 41int[int][int] intersection; // local-to-neighbors renumbering 42real[int] D; // partition of unity 43{ 44 Th = trunc(ThCoarse, 1, split = s); 45 fespace Ph(Th, P0); 46 Ph part; 47 part = partCoarse; 48 buildWithPartitioning(Th, part[], 1, intersection, D, Pk, mpiCommWorld); 49} 50 51real epsilon = epsilonA; 52varf vPb(u, v) = intN(Th)(-(k^2 - 1i*epsilon)*u*v + grad(u)'*grad(v)) 53 + intN1(Th, Robin)(1i*k*u*v) 54 - intN(Th)(f*v) 55 + on(Dirichlet, u = 0); 56matrix<complex> Mat = vPb(Wh, Wh, sym = 1); 57fespace WhCoarse(ThCoarse, Pk); 58matrix<complex> MatCoarse; 59{ 60 meshN ThBackup = Th; 61 Th = ThCoarse; 62 epsilon = epsilonE; 63 MatCoarse = vPb(Wh, Wh, sym = 1); 64 epsilon = epsilonA; 65 Th = ThBackup; 66} 67set(MatCoarse, solver = sparsesolverSym, sym = 1); 68schwarz<complex> ACoarseSeq(MatCoarse); 69set(ACoarseSeq, prefix = "level_2_"); 70matrix R = interpolate(Wh, WhCoarse); 71complex[int] rhs = vPb(0, Wh); 72 73schwarz<complex> A(Mat, intersection, D); 74set(A, sparams = "-hpddm_verbosity " + (!mpirank ? "1" : "0") + " -hpddm_schwarz_method ras -hpddm_variant right -hpddm_gmres_restart 200 -hpddm_max_it 200"); 75 76Wh<complex> u; 77u[] = A^-1 * rhs; 78plotMPI(Th, real(u), Pk, def, real, cmm = "Global solution with a one-level method") 79 80set(A, sparams = "-hpddm_verbosity " + (!mpirank ? "3" : "0") + " -hpddm_schwarz_coarse_correction deflated -hpddm_variant flexible -hpddm_level_2_verbosity 0"); 81func complex[int] correctionExact(complex[int]& in) { 82 complex[int] out(in.n), tmp(WhCoarse.ndof), tmpReduced(WhCoarse.ndof); 83 for[i, di: D] in[i] *= di; 84 tmp = R' * in; 85 mpiAllReduce(tmp, tmpReduced, mpiCommWorld, mpiSUM); 86 tmp = ACoarseSeq^-1 * tmpReduced; 87 out = R * tmp; 88 return out; 89} 90if(mpisize > 1 && isSetOption("schwarz_coarse_correction")) // two-level Schwarz methods 91 attachCoarseOperator(mpiCommWorld, A, correctionExact); 92 93u[] = 0; 94u[] = A^-1 * rhs; 95plotMPI(Th, real(u), Pk, def, real, cmm = "Global solution with a redundant coarse correction") 96 97int[int][int] intersectionCoarse; 98real[int] DCoarse; 99{ 100 buildWithPartitioning(ThCoarse, partCoarse[], 1, intersectionCoarse, DCoarse, Pk, mpiCommWorld); 101 meshN ThBackup = Th; 102 Th = ThCoarse; 103 epsilon = epsilonE; 104 MatCoarse = vPb(Wh, Wh, sym = 1); 105 epsilon = epsilonA; 106 Th = ThBackup; 107} 108R = interpolate(Wh, WhCoarse); 109schwarz<complex> ACoarse(MatCoarse, intersectionCoarse, DCoarse); 110set(ACoarse, sparams = "-hpddm_level_2_verbosity " + (!mpirank ? "1" : "0") + " -hpddm_level_2_schwarz_method ras -hpddm_level_2_tol 1e-1 -hpddm_level_2_gmres_restart 200 -hpddm_level_2_max_it 200 -hpddm_level_2_krylov_method gcrodr -hpddm_level_2_recycle 10", prefix = "level_2_"); 111macro mplot()wait = 1, cmm = "Transfers on the correction at iteration " + ijk// 112int ijk = 1; 113func complex[int] correctionInexact(complex[int]& in) { 114 complex[int] out(in.n), tmp(WhCoarse.ndof), tmpReduced(WhCoarse.ndof); 115 if(!NoGraphicWindow && ijk > 0 && ijk < 10) { 116 Wh outW; 117 for[i, xi: in] outW[][i] = real(xi); 118 plotMPI(Th, outW, Pk, def, real, mplot) 119 } 120 tmpReduced = R' * in; 121 exchange(ACoarse, tmpReduced, scaled = true); 122 if(!NoGraphicWindow && ijk > 0 && ijk < 10) { 123 WhCoarse outWCoarse; 124 for[i, xi: tmpReduced] outWCoarse[][i] = real(tmpReduced[i]); 125 plotMPI(ThCoarse, outWCoarse, Pk, def, real, mplot) 126 ++ijk; 127 } 128 tmp = ACoarse^-1 * tmpReduced; 129 out = R * tmp; 130 exchange(A, out); 131 return out; 132} 133if(mpisize > 1 && isSetOption("schwarz_coarse_correction")) // two-level Schwarz methods 134 attachCoarseOperator(mpiCommWorld, A, correctionInexact); 135else 136 exit(0); 137 138u[] = 0; 139u[] = A^-1 * rhs; 140plotMPI(Th, real(u), Pk, def, real, cmm = "Global solution with a one-level inner coarse correction (two levels in total)") 141