1 /*******************************************************
 2  *  This file is part of the FreeFEM tutorial made by  *
 3  *      Pierre Jolivet <pierre@joliv.et>               *
 4  *                                                     *
 5  *  See https://joliv.et/FreeFem-tutorial for more     *
 6  *                                                     *
 7  *   Description: systems with multiple right-hand     *
 8  *                sides                                *
 9  *******************************************************/
 
11 load "PETSc-complex"
12 macro dimension()2// EOM
13 include "macro_ddm.idp"
 
15 border a(t=0, 2*pi) { x=cos(t); y=sin(t); label=1; }
16 mesh Th;
17 if(!mpirank)
18     Th = buildmesh(a(300));
19 broadcast(processor(0), Th);
 
21 Mat<complex> A;
22 func Pk = P2;
23 MatCreate(Th, A, Pk);
 
25 real k = 20.0;
26 macro grad(u)[dx(u), dy(u)]// EOM
27 varf vPb(u, v) = int2d(Th)(- (grad(u)' * grad(v)) + k^2 * u * v);
28 fespace Vh(Th, Pk);
29 A = vPb(Vh, Vh);
 
31 complex[int, int] rhs(Vh.ndof, 5), sol(Vh.ndof, 5);
32 for(int i = -2; i < 3; ++i) {
33     func f = 100 * exp(-20 * ((x-i/3.0)^2 + (y-i/3.0)^2));
34     varf vRhs(u, v) = int2d(Th)(f * v);
35     rhs(:, i + 2) = vRhs(0, Vh);
36 }
 
38 complex[int, int] B(A.n, rhs.m), X(A.n, rhs.m);
39 ChangeNumbering(A, rhs, B);
40 set(A, sparams = "-ksp_type preonly -pc_type lu");
41 KSPSolve(A, B, X);
42 ChangeNumbering(A, sol, X, inverse = true, exchange = true);
43 for(int i = 0; i < rhs.m; ++i) {
44     Vh u;
45     u[] = sol(:, i).re;
46     plotD(Th, u, cmm = "-ksp_type preonly: u_" + i);
47 }
 
49 X = 0;
50 set(A, sparams = "-ksp_type hpddm -ksp_hpddm_type bgmres -pc_type asm -sub_pc_type lu");
51 KSPSolve(A, B, X);
52 ChangeNumbering(A, sol, X, inverse = true, exchange = true);
53 for(int i = 0; i < rhs.m; ++i) {
54     Vh u;
55     u[] = sol(:, i).re;
56     plotD(Th, u, cmm = "-ksp_type hpddm: u_" + i);
57 }