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: setting up PETSc preconditioners     *
 8  *******************************************************/
  
10 load "PETSc"
11 macro dimension()2// EOM
12 include "macro_ddm.idp"
 
14 macro grad(u)[dx(u), dy(u)]// EOM
15 load "Element_P3"
16 func Pk = P3;
 
18 border upper(t=0, pi)    { x=cos(t); y=sin(t); label=1; }
19 border lower(t=pi, 2*pi) { x=cos(t); y=sin(t); label=2; }
20 mesh Th = buildmesh(upper(100) + lower(100));
21 Mat A;
22 MatCreate(Th, A, Pk);
 
24 fespace Vh(Th, Pk);
25 varf vPb(u, v) = intN(Th)(grad(u)' * grad(v)) + intN(Th)(v) + on(2, u = 0);
26 matrix Loc = vPb(Vh, Vh, tgv = -2);
27 real[int] b = vPb(0, Vh, tgv = -1);
 
29 A = Loc;
 
31 Vh u;
32 set(A, sparams = "-pc_type cholesky -ksp_type bcgs -ksp_converged_reason");
33 ObjectView(A, object = "ksp");
34 u[] = A^-1 * b;
35 plotD(Th, u, cmm = "Solution with -pc_type cholesky");
 
37 set(A, sparams = "-pc_type gamg -ksp_type gmres");
 
39 u[] = 0;
40 u[] = A^-1 * b;
41 ObjectView(A, object = "ksp");
42 plotD(Th, u, cmm = "Solution with -pc_type gamg");