1// run with MPI: ff-mpirun -np 4 script.edp 2// NBPROC 4 3 4load "PETSc" 5macro dimension()2// EOM 6include "macro_ddm.idp" 7 8real c = 6.25; 9int N = 80; 10 11// domain: unit square 12border aa(t=0,1) { x=t; y=0; }; 13border bb(t=0,1) { x=1; y=t; }; 14border cc(t=0,1) { x=1-t; y=1; }; 15border dd(t=0,1) { x=0; y=1-t; }; 16 17mesh M = buildmesh(aa(N)+bb(N)+cc(N)+dd(N)); 18 19load "Element_P3" 20func Pk = P1; 21Mat A; 22buildMat(M, getARGV("-split", 1), A, Pk, mpiCommWorld) 23fespace Vh(M, Pk); 24Vh u; 25func BC = cos(pi*x)*cos(pi*y); 26varf vInit(w, v) = on(aa, bb, cc, dd, w = BC); 27varf vJ(w, v) = int2d(M)(dx(w)*dx(v) + dy(w)*dy(v) - c*exp(u)*w*v) + on(aa, bb, cc, dd, w = 0); 28varf vRes(w, v) = int2d(M)(dx(u)*dx(v) + dy(u)*dy(v) - c*exp(u)*v) + on(aa, bb, cc, dd, w = u); 29func real[int] funcRes(real[int]& inPETSc) { 30 changeNumbering(A, u[], inPETSc, inverse = true, exchange = true); 31 real[int] out(Vh.ndof); 32 out = vRes(0, Vh, tgv = -2); 33 real[int] outPETSc; 34 changeNumbering(A, out, outPETSc); 35 return outPETSc; 36} 37func int funcJ(real[int]& inPETSc) { 38 changeNumbering(A, u[], inPETSc, inverse = true, exchange = true); 39 A = vJ(Vh, Vh, tgv = -2); 40 return 0; 41} 42if(hasType("PC", "hpddm")) { 43 bool stage = usedARGV("-stage") != -1; 44 for(int i = 0; i < 3; ++i) { 45 string params; 46 if(i == 0) { 47 params = ""; 48 if(stage) 49 PetscLogStagePush("Everything"); 50 } 51 else if(i == 1) { 52 params = "-pc_hpddm_coarse_ksp_reuse_preconditioner"; 53 if(stage) 54 PetscLogStagePush("Fine level"); 55 } 56 else if(i == 2) { 57 params = "-pc_hpddm_levels_1_ksp_reuse_preconditioner"; 58 if(stage) 59 PetscLogStagePush("Nothing"); 60 } 61 set(A, sparams = "-pc_type hpddm -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_levels_1_eps_nev 10 -pc_hpddm_define_subdomains -pc_hpddm_has_neumann -ksp_converged_reason -pc_hpddm_levels_1_pc_type asm " + params); 62 real[int] bPETSc; 63 u[] = vInit(0, Vh, tgv = -2); 64 changeNumbering(A, u[], bPETSc); 65 real[int] xPETSc = bPETSc; 66 SNESSolve(A, funcJ, funcRes, bPETSc, xPETSc, sparams = "-snes_monitor -snes_type newtonls -snes_converged_reason -ksp_converged_reason -snes_view"); 67 changeNumbering(A, u[], xPETSc, inverse = true, exchange = true); 68 macro def(u)u// 69 plotMPI(M, u, Pk, def, real, cmm = "Global solution") 70 set(A, sparams = "-pc_type none"); 71 if(stage) 72 PetscLogStagePop(); 73 } 74} 75