MFEM  v4.5.2 Finite element discretization library
lor_solvers.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
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 //
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 // Device sample runs:
62 // lor_solvers -fe h -d cuda
63 // lor_solvers -fe n -d cuda
64 // lor_solvers -fe r -d cuda
65 // lor_solvers -fe l -d cuda
66
67 #include "mfem.hpp"
68 #include <fstream>
69 #include <iostream>
70 #include <memory>
71
72 #include "lor_mms.hpp"
73
74 using namespace std;
75 using namespace mfem;
76
78
79 int main(int argc, char *argv[])
80 {
81  const char *mesh_file = "../../data/star.mesh";
82  int ref_levels = 1;
83  int order = 3;
84  const char *fe = "h";
85  const char *device_config = "cpu";
86  bool visualization = true;
87
88  OptionsParser args(argc, argv);
89  args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
91  "Number of times to refine the mesh uniformly.");
92  args.AddOption(&order, "-o", "--order", "Polynomial degree.");
94  "FE type. h for H1, n for Hcurl, r for Hdiv, l for L2");
96  "--no-visualization",
97  "Enable or disable GLVis visualization.");
99  "Device configuration string, see Device::Configure().");
100  args.ParseCheck();
101
102  Device device(device_config);
103  device.Print();
104
105  bool H1 = false, ND = false, RT = false, L2 = false;
106  if (string(fe) == "h") { H1 = true; }
107  else if (string(fe) == "n") { ND = true; }
108  else if (string(fe) == "r") { RT = true; }
109  else if (string(fe) == "l") { L2 = true; }
110  else { MFEM_ABORT("Bad FE type. Must be 'h', 'n', 'r', or 'l'."); }
111
112  if (RT) { grad_div_problem = true; }
113  double kappa = (order+1)*(order+1); // Penalty used for DG discretizations
114
115  Mesh mesh(mesh_file, 1, 1);
116  int dim = mesh.Dimension();
117  MFEM_VERIFY(dim == 2 || dim == 3, "Spatial dimension must be 2 or 3.");
118  for (int l = 0; l < ref_levels; l++) { mesh.UniformRefinement(); }
119
120  FunctionCoefficient f_coeff(f), u_coeff(u);
121  VectorFunctionCoefficient f_vec_coeff(dim, f_vec), u_vec_coeff(dim, u_vec);
122
123  int b1 = BasisType::GaussLobatto, b2 = BasisType::IntegratedGLL;
124  unique_ptr<FiniteElementCollection> fec;
125  if (H1) { fec.reset(new H1_FECollection(order, dim, b1)); }
126  else if (ND) { fec.reset(new ND_FECollection(order, dim, b1, b2)); }
127  else if (RT) { fec.reset(new RT_FECollection(order-1, dim, b1, b2)); }
128  else { fec.reset(new L2_FECollection(order, dim, b1)); }
129
130  FiniteElementSpace fes(&mesh, fec.get());
131  cout << "Number of DOFs: " << fes.GetTrueVSize() << endl;
132
133  Array<int> ess_dofs;
134  // In DG, boundary conditions are enforced weakly, so no essential DOFs.
135  if (!L2) { fes.GetBoundaryTrueDofs(ess_dofs); }
136
137  BilinearForm a(&fes);
138  if (H1 || L2)
139  {
142  }
143  else
144  {
146  }
147
148  if (ND) { a.AddDomainIntegrator(new CurlCurlIntegrator); }
149  else if (RT) { a.AddDomainIntegrator(new DivDivIntegrator); }
150  else if (L2)
151  {
154  }
155  // TODO: L2 diffusion not implemented with partial assembly
156  if (!L2) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
157  a.Assemble();
158
159  LinearForm b(&fes);
160  if (H1 || L2) { b.AddDomainIntegrator(new DomainLFIntegrator(f_coeff)); }
161  else { b.AddDomainIntegrator(new VectorFEDomainLFIntegrator(f_vec_coeff)); }
162  if (L2)
163  {
164  // DG boundary conditions are enforced weakly with this integrator.
166  }
167  b.Assemble();
168
169  GridFunction x(&fes);
170  if (H1 || L2) { x.ProjectCoefficient(u_coeff);}
171  else { x.ProjectCoefficient(u_vec_coeff); }
172
173  Vector X, B;
174  OperatorHandle A;
175  a.FormLinearSystem(ess_dofs, x, b, A, X, B);
176
177 #ifdef MFEM_USE_SUITESPARSE
178  LORSolver<UMFPackSolver> solv_lor(a, ess_dofs);
179 #else
180  LORSolver<GSSmoother> solv_lor(a, ess_dofs);
181 #endif
182
183  CGSolver cg;
184  cg.SetAbsTol(0.0);
185  cg.SetRelTol(1e-12);
186  cg.SetMaxIter(500);
187  cg.SetPrintLevel(1);
188  cg.SetOperator(*A);
189  cg.SetPreconditioner(solv_lor);
190  cg.Mult(B, X);
191
192  a.RecoverFEMSolution(X, b, x);
193
194  double er =
195  (H1 || L2) ? x.ComputeL2Error(u_coeff) : x.ComputeL2Error(u_vec_coeff);
196  cout << "L2 error: " << er << endl;
197
198  if (visualization)
199  {
200  // Save the solution and mesh to disk. The output can be viewed using
201  // GLVis as follows: "glvis -m mesh.mesh -g sol.gf"
202  x.Save("sol.gf");
203  mesh.Save("mesh.mesh");
204
205  // Also save the solution for visualization using ParaView
206  ParaViewDataCollection dc("LOR", &mesh);
207  dc.SetPrefixPath("ParaView");
208  dc.SetHighOrderOutput(true);
209  dc.SetLevelsOfDetail(order);
210  dc.RegisterField("u", &x);
211  dc.SetCycle(0);
212  dc.SetTime(0.0);
213  dc.Save();
214  }
215
216  return 0;
217 }
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
Definition: gridfunc.cpp:2763
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
Definition: solvers.hpp:493
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
int Dimension() const
Definition: mesh.hpp:1047
Helper class for ParaView visualization data.
Integrator for (curl u, curl v) for Nedelec elements.
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
int main(int argc, char *argv[])
Definition: lor_solvers.cpp:79
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:712
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:279
double kappa
Definition: ex24.cpp:54
STL namespace.
(Q div u, div v) for RT elements
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9878
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
void u_vec(const Vector &xvec, Vector &u)
Definition: lor_mms.hpp:50
void SetHighOrderOutput(bool high_order_output_)
void SetTime(double t)
Set physical time (for time-dependent simulations)
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:373
A general vector function coefficient.
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:590
void SetAbsTol(double atol)
Definition: solvers.hpp:200
void SetRelTol(double rtol)
Definition: solvers.hpp:199
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
double a
Definition: lissajous.cpp:41
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
int dim
Definition: ex24.cpp:53
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2396
void SetLevelsOfDetail(int levels_of_detail_)
Definition: lor_solvers.cpp:77
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Definition: lininteg.hpp:346
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:507
A general function coefficient.
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:447
Vector data type.
Definition: vector.hpp:60
void f_vec(const Vector &xvec, Vector &f)
Definition: lor_mms.hpp:68
virtual void Save(const char *fname, int precision=16) const
Definition: mesh.cpp:10351
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:252
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:24
virtual void Save() override
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3673
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
void ParseCheck(std::ostream &out=mfem::out)
Definition: optparser.cpp:252
void SetPrefixPath(const std::string &prefix)
Set the path where the DataCollection will be saved.
Represents a solver of type SolverType created using the low-order refined version of the given Bilin...
Definition: lor.hpp:203
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:320
double f(const Vector &p)