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