1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 //
12 // ---------------------------------
13 // Low-Order Refined Solvers Miniapp
14 // ---------------------------------
15 //
16 // This miniapp illustrates the use of low-order refined preconditioners for
17 // finite element problems defined using H1, H(curl), H(div), or L2 finite
18 // element spaces. The following problems are solved, depending on the chosen
19 // finite element space:
20 //
21 // H1 and L2: definite Helmholtz problem, u - Delta u = f
22 // (in L2 discretized using the symmetric interior penalty DG method)
23 //
24 // H(curl): definite Maxwell problem, u + curl curl u = f
25 //
26 // H(div): grad-div problem, u - grad(div u) = f
27 //
28 // In each case, the high-order finite element problem is preconditioned using a
29 // low-order finite element discretization defined on a Gauss-Lobatto refined
30 // mesh. The low-order problem is solved using a direct solver if MFEM is
31 // compiled with SuiteSparse enabled, or with one iteration of symmetric
32 // Gauss-Seidel smoothing otherwise.
33 //
34 // For vector finite element spaces, the special "Integrated" basis type is used
35 // to obtain spectral equivalence between the high-order and low-order refined
36 // discretizations. This basis is defined in reference [1] and spectral
37 // equivalence is shown in [2]:
38 //
39 // [1]. M. Gerritsma. Edge functions for spectral element methods. Spectral and
40 // High Order Methods for Partial Differential Equations. (2010)
41 // [2]. C. Dohrmann. Spectral equivalence properties of higher-order tensor
42 // product finite elements and applications to preconditioning. (2021)
43 //
44 // The action of the high-order operator is computed using MFEM's partial
45 // assembly/matrix-free algorithms (except in the case of L2, which remains
46 // future work).
47 //
48 // Compile with: make lor_solvers
49 //
50 // Sample runs:
51 //
52 // lor_solvers -fe h
53 // lor_solvers -fe n
54 // lor_solvers -fe r
55 // lor_solvers -fe l
56 // lor_solvers -m ../../data/amr-quad.mesh -fe h
57 // lor_solvers -m ../../data/amr-quad.mesh -fe n
58 // lor_solvers -m ../../data/amr-quad.mesh -fe r
59 // lor_solvers -m ../../data/amr-quad.mesh -fe l
60
61 #include "mfem.hpp"
62 #include <fstream>
63 #include <iostream>
64 #include <memory>
65
66 #include "lor_mms.hpp"
67
68 using namespace std;
69 using namespace mfem;
70
71 bool grad_div_problem = false;
72
main(int argc,char * argv[])73 int main(int argc, char *argv[])
74 {
75 const char *mesh_file = "../../data/star.mesh";
76 int ref_levels = 1;
77 int order = 3;
78 const char *fe = "h";
79 bool visualization = true;
80
81 OptionsParser args(argc, argv);
82 args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
83 args.AddOption(&ref_levels, "-r", "--refine",
84 "Number of times to refine the mesh uniformly.");
85 args.AddOption(&order, "-o", "--order", "Polynomial degree.");
86 args.AddOption(&fe, "-fe", "--fe-type",
87 "FE type. h for H1, n for Hcurl, r for Hdiv, l for L2");
88 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
89 "--no-visualization",
90 "Enable or disable GLVis visualization.");
91 args.ParseCheck();
92
93 bool H1 = false, ND = false, RT = false, L2 = false;
94 if (string(fe) == "h") { H1 = true; }
95 else if (string(fe) == "n") { ND = true; }
96 else if (string(fe) == "r") { RT = true; }
97 else if (string(fe) == "l") { L2 = true; }
98 else { MFEM_ABORT("Bad FE type. Must be 'h', 'n', 'r', or 'l'."); }
99
100 if (RT) { grad_div_problem = true; }
101 double kappa = (order+1)*(order+1); // Penalty used for DG discretizations
102
103 Mesh mesh(mesh_file, 1, 1);
104 int dim = mesh.Dimension();
105 MFEM_VERIFY(dim == 2 || dim == 3, "Spatial dimension must be 2 or 3.");
106 for (int l = 0; l < ref_levels; l++) { mesh.UniformRefinement(); }
107
108 FunctionCoefficient f_coeff(f), u_coeff(u);
109 VectorFunctionCoefficient f_vec_coeff(dim, f_vec), u_vec_coeff(dim, u_vec);
110
111 int b1 = BasisType::GaussLobatto, b2 = BasisType::IntegratedGLL;
112 unique_ptr<FiniteElementCollection> fec;
113 if (H1) { fec.reset(new H1_FECollection(order, dim, b1)); }
114 else if (ND) { fec.reset(new ND_FECollection(order, dim, b1, b2)); }
115 else if (RT) { fec.reset(new RT_FECollection(order-1, dim, b1, b2)); }
116 else { fec.reset(new L2_FECollection(order, dim, b1)); }
117
118 FiniteElementSpace fes(&mesh, fec.get());
119 cout << "Number of DOFs: " << fes.GetTrueVSize() << endl;
120
121 Array<int> ess_dofs;
122 // In DG, boundary conditions are enforced weakly, so no essential DOFs.
123 if (!L2) { fes.GetBoundaryTrueDofs(ess_dofs); }
124
125 BilinearForm a(&fes);
126 if (H1 || L2)
127 {
128 a.AddDomainIntegrator(new MassIntegrator);
129 a.AddDomainIntegrator(new DiffusionIntegrator);
130 }
131 else
132 {
133 a.AddDomainIntegrator(new VectorFEMassIntegrator);
134 }
135
136 if (ND) { a.AddDomainIntegrator(new CurlCurlIntegrator); }
137 else if (RT) { a.AddDomainIntegrator(new DivDivIntegrator); }
138 else if (L2)
139 {
140 a.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(-1.0, kappa));
141 a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(-1.0, kappa));
142 }
143 // TODO: L2 diffusion not implemented with partial assembly
144 if (!L2) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
145 a.Assemble();
146
147 LinearForm b(&fes);
148 if (H1 || L2) { b.AddDomainIntegrator(new DomainLFIntegrator(f_coeff)); }
149 else { b.AddDomainIntegrator(new VectorFEDomainLFIntegrator(f_vec_coeff)); }
150 if (L2)
151 {
152 // DG boundary conditions are enforced weakly with this integrator.
153 b.AddBdrFaceIntegrator(new DGDirichletLFIntegrator(u_coeff, -1.0, kappa));
154 }
155 b.Assemble();
156
157 GridFunction x(&fes);
158 if (H1 || L2) { x.ProjectCoefficient(u_coeff);}
159 else { x.ProjectCoefficient(u_vec_coeff); }
160
161 Vector X, B;
162 OperatorHandle A;
163 a.FormLinearSystem(ess_dofs, x, b, A, X, B);
164
165 #ifdef MFEM_USE_SUITESPARSE
166 LORSolver<UMFPackSolver> solv_lor(a, ess_dofs);
167 #else
168 LORSolver<GSSmoother> solv_lor(a, ess_dofs);
169 #endif
170
171 CGSolver cg;
172 cg.SetAbsTol(0.0);
173 cg.SetRelTol(1e-12);
174 cg.SetMaxIter(500);
175 cg.SetPrintLevel(1);
176 cg.SetOperator(*A);
177 cg.SetPreconditioner(solv_lor);
178 cg.Mult(B, X);
179
180 a.RecoverFEMSolution(X, b, x);
181
182 double er =
183 (H1 || L2) ? x.ComputeL2Error(u_coeff) : x.ComputeL2Error(u_vec_coeff);
184 cout << "L2 error: " << er << endl;
185
186 if (visualization)
187 {
188 // Save the solution and mesh to disk. The output can be viewed using
189 // GLVis as follows: "glvis -m mesh.mesh -g sol.gf"
190 x.Save("sol.gf");
191 mesh.Save("mesh.mesh");
192
193 // Also save the solution for visualization using ParaView
194 ParaViewDataCollection dc("LOR", &mesh);
195 dc.SetPrefixPath("ParaView");
196 dc.SetHighOrderOutput(true);
197 dc.SetLevelsOfDetail(order);
198 dc.RegisterField("u", &x);
199 dc.SetCycle(0);
200 dc.SetTime(0.0);
201 dc.Save();
202 }
203
204 return 0;
205 }
206