1// run with MPI: ff-mpirun -np 4 script.edp 2// NBPROC 4 3 4load "PETSc" // PETSc plugin 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 = P1; // finite element space 12 13int s = getARGV("-split", 1); // refinement factor 14 15meshN Th = square(1, 1); 16fespace Wh(Th, Pk); // local finite element space 17int[int][int] intersection; // local-to-neighbors renumbering 18real[int] D; // partition of unity 19{ 20 int[int] l = [1, 2, 2, 2]; 21 Th = square(getARGV("-global", 40), getARGV("-global", 40), label = l); // global mesh 22 build(Th, s, intersection, D, Pk, mpiCommWorld) 23} 24 25real[int] rhs(Wh.ndof); // local right-hand side 26matrix<real> Loc; // local operator 27{ // local weak form 28 fespace Ph(Th, P0); 29 Ph kappa = x < 0.25 ? 10.0 : 1.0; 30 varf vPb(u, v) = intN(Th)(-1.0 * kappa * grad(u)' * grad(v)) + on(1, u = 0.0); 31 Loc = vPb(Wh, Wh, tgv = -2); 32 rhs = vPb(0, Wh, tgv = -2); 33} 34 35Mat T(Loc, intersection, D); 36func real[int] funcRes(real t, real[int]& in, real[int]& inT) { 37 real[int] out(in.n); 38 T = Loc; 39 MatMult(T, in, out); 40 out = inT - out; 41 return out; 42} 43real shift = 0; 44matrix Id; 45{ 46 real[int] D(Loc.n); 47 D = 1.0; 48 Id = D; 49} 50func int funcJ(real t, real[int]& in, real[int]& inT, real a) { 51 matrix B = (-1.0) * Loc + a * Id; 52 T = B; 53 shift = a; 54 return 0; 55} 56Wh w; 57func int funcM(int s, real t, real[int]& u) { 58 changeNumbering(T, w[], u, exchange = true, inverse = true); 59 plotMPI(Th, w, Pk, def, real, cmm = "Global solution at step " + s + " (time " + t + ")") 60} 61w = (0.5 - x)^2 + (0.5 - y)^2 < 0.2 ? 1.0 : 0.0; 62real[int] wPETSc; 63changeNumbering(T, w[], wPETSc); 64TSSolve(T, funcJ, funcRes, wPETSc, sparams = "-ts_type beuler -ts_dt 0.1 -ts_max_time 100 -ts_exact_final_time interpolate -ts_max_snes_failures -1 -ts_view -pc_type lu -ts_adapt_type basic -ts_rtol 1e-3", monitor = funcM); 65