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: geometric multigrid                  *
 8  *******************************************************/
  
10 load "PETSc-complex"
11 load "Element_Mixte3d"
12 macro dimension()3// EOM
13 include "macro_ddm.idp"
14 include "cube.idp"
 
16 macro Curl(ux, uy, uz)[dy(uz)-dz(uy), dz(ux)-dx(uz), dx(uy)-dy(ux)]// EOM
17 macro CrossN(ux, uy, uz)[uy*N.z-uz*N.y, uz*N.x-ux*N.z, ux*N.y-uy*N.x]// EOM
18 macro def(i)[i, i#y, i#z]// EOM
19 macro init(i)[i, i, i]// EOM
20 macro defPart(i)i// EOM
21 macro initPart(i)i// EOM
 
23 int level = 2;
24 int s = 3;
 
26 Mat<complex>[int] MG(level);
27 meshN[int] Th(level);
28 matrix[int] P(level - 1);
29 int Robin = 2;
30 int[int] LL = [Robin, Robin, Robin, Robin, Robin, Robin];
31 real n = getARGV("-n", 40) / s;
32 Th[level - 1] = cube(n, n, n, [x, y, z], label = LL);
33 func Pk = Edge03d;
34 func PkPart = Edge03ds0;
35 buildMatEdgeRecursive(Th, s, level, P, MG, Pk, mpiCommWorld, PkPart, defPart, initPart);
 
37 func k = 6 * pi;
38 func f = exp(-60 * ((x-0.5)^2 + (y-0.5)^2 + (z-0.5)^2));
 
40 complex[int] rhs;
41 matrix<complex>[int] Opt(level);
42 for(int i = 0; i < level; ++i) {
43     fespace Vh(Th[i], Pk);
44     varf vPb([Ex,Ey,Ez], [vx,vy,vz]) =
45         intN(Th[i])(Curl(vx,vy,vz)'*Curl(Ex,Ey,Ez))
46       + intN(Th[i])(-k^2*[vx,vy,vz]'*[Ex,Ey,Ez])
47       + intN1(Th[i],Robin)(1i*k*CrossN(vx,vy,vz)'*CrossN(Ex,Ey,Ez))
48       + intN1(Th[i],-111111)(1i*k*CrossN(vx,vy,vz)'*CrossN(Ex,Ey,Ez));
49     Opt[i] = vPb(Vh, Vh, sym = 1);
50     MG[i] = Opt[i];
51     if(i == 0) {
52         varf vRHS([Ex,Ey,Ez],[vx,vy,vz]) = -intN(Th[i])([vx,vy,vz]'*[0,0,f]);
53         rhs.resize(Vh.ndof);
54         rhs = vRHS(0, Vh);
55     }
56 }
57 set(MG, P, sparams = "-pc_type mg -ksp_monitor -ksp_view_final_residual -ksp_type fgmres -ksp_gmres_restart 200 -ksp_max_it 200");
58 set(MG, 0, sparams = "-mg_coarse_ksp_type gmres -mg_coarse_ksp_rtol 1e-1 -mg_coarse_ksp_pc_side right " + " -mg_coarse_ksp_max_it 50 -mg_coarse_ksp_converged_reason -mg_coarse_pc_type asm -mg_coarse_sub_pc_factor_mat_solver_type mumps -mg_coarse_sub_pc_type cholesky -mg_coarse_pc_asm_type restrict", O = Opt[level - 1]);
59 set(MG, level - 1, sparams = "-mg_levels_ksp_type richardson -mg_levels_ksp_pc_side left -mg_levels_pc_type asm -mg_levels_sub_pc_type cholesky -mg_levels_sub_pc_factor_mat_solver_type mumps -mg_levels_pc_asm_type restrict", O = Opt[0]);
60 fespace Vh(Th[0], Pk);
61 Vh<complex> def(u);
62 u[] = MG[0]^-1 * rhs;
63 mesh3 ThFine = Th[0];
64 fespace VhFine(ThFine, P1);
65 VhFine plt = sqrt(real(u)^2 + real(uy)^2 + real(uz)^2);
66 plotD(ThFine, plt, cmm = "Solution");
67 int[int] fforder = [1, 1];
68 savevtk("maxwell.vtu", ThFine, [real(u), real(uy), real(uz)], [imag(u), imag(uy), imag(uz)], order = fforder, bin = 1);