1// run with MPI: ff-mpirun -np 4 script.edp 2// NBPROC 4 3 4load "PETSc" 5macro dimension()2// 6include "macro_ddm.idp" 7 8macro def(u)u// 9macro init(u)u// 10border C(t = 0, 2*pi){ x=cos(t); y=sin(t); } 11mesh Th = buildmesh(C(getARGV("-global", 100))); 12 13func Pk = P1dc; 14fespace Vh(Th, Pk); 15real[int] D; 16int[int][int] intersection; 17build(Th, 1, intersection, D, Pk, mpiCommWorld); 18Vh v1 = y, v2 = -x; 19 20macro diff(n)(abs(n) - (n))// 21 22varf vK(c, w) = int2d(Th)((v1 * dx(c) + v2 * dy(c)) * w) + intalledges(Th)((1 - nTonEdge) * w * diff(v1 * N.x + v2 * N.y) * 0.5 * jump(c)); 23varf vM(c, w) = int2d(Th)(c * w); 24matrix M = vM(Vh, Vh); 25matrix K = vK(Vh, Vh); 26 27Mat dK(K, intersection, D); 28Mat dM(dK, M); 29set(M, sparams = "-pc_type jacobi"); 30Mat dT(Vh.ndof, intersection, D); 31 32func real[int] funcExplicit(real t, real[int]& in) { 33 real[int] temp(in.n), out(in.n); 34 MatMult(dK, in, temp); 35 temp *= -1; 36 KSPSolve(dM, temp, out); 37 return out; 38} 39func real[int] funcRes(real t, real[int]& in, real[int]& inT) { 40 real[int] temp(in.n), out(in.n); 41 MatMult(dM, inT, out); 42 MatMult(dK, in, temp); 43 out += temp; 44 return out; 45} 46real shift = 0; 47func int funcJ(real t, real[int]& in, real[int]& inT, real a) { 48 if(abs(a - shift) > 1.0e-10) { 49 matrix B = a * M + K; 50 dT = B; 51 shift = a; 52 } 53 return 0; 54} 55Vh w; 56string type; 57func int funcM(int s, real t, real[int]& u) { 58 changeNumbering(dT, w[], u, inverse = true, exchange = true); 59 plotMPI(Th, w, Pk, def, real, cmm = "Global solution at step " + s + " (time " + t + ") with " + type + " method") 60} 61for(int j = 0; j < 3; ++j) { 62 Vh c = exp(-10 * ((x - 0.3)^2 + (y - 0.3)^2)); 63 real[int] cPETSc; 64 changeNumbering(dT, c[], cPETSc); 65 real T = 3 * pi; 66 string deflt = "-ts_view -ts_max_snes_failures -1 -ts_exact_final_time interpolate -ts_adapt_type basic -ts_max_time " + string(T); 67 if(j == 0) { 68 type = "explicit"; 69 TSSolve(dT, funcJ, funcExplicit, cPETSc, sparams = "-ts_type rk -ts_rk_type 5f " + deflt, monitor = funcM); 70 } 71 else if(j == 1) { 72 type = "implicit"; 73 TSSolve(dT, funcJ, funcRes, cPETSc, sparams = "-ts_type rosw " + deflt, monitor = funcM); 74 } 75 else { 76 type = "semi-implicit"; 77 TSSolve(dT, funcJ, funcRes, cPETSc, sparams = "-ts_type arkimex -ts_arkimex_type 5 " + deflt, monitor = funcM); 78 } 79 changeNumbering(dT, c[], cPETSc, inverse = true, exchange = true); 80 macro params()cmm = "Global solution with " + type + " method", wait = 1, dim = 3, fill = 1// EOM 81 plotMPI(Th, c, Pk, def, real, params) 82} 83