MFEM  v4.6.0 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
77 int main(int argc, char *argv[])
78 {
79  const char *mesh_file = "../../data/star.mesh";
80  int ref_levels = 1;
81  int order = 3;
82  const char *fe = "h";
83  const char *device_config = "cpu";
84  bool visualization = true;
85
86  OptionsParser args(argc, argv);
87  args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
89  "Number of times to refine the mesh uniformly.");
90  args.AddOption(&order, "-o", "--order", "Polynomial degree.");
92  "FE type. h for H1, n for Hcurl, r for Hdiv, l for L2");
94  "--no-visualization",
95  "Enable or disable GLVis visualization.");
97  "Device configuration string, see Device::Configure().");
98  args.ParseCheck();
99
100  Device device(device_config);
101  device.Print();
102
103  bool H1 = false, ND = false, RT = false, L2 = false;
104  if (string(fe) == "h") { H1 = true; }
105  else if (string(fe) == "n") { ND = true; }
106  else if (string(fe) == "r") { RT = true; }
107  else if (string(fe) == "l") { L2 = true; }
108  else { MFEM_ABORT("Bad FE type. Must be 'h', 'n', 'r', or 'l'."); }
109
110  double kappa = (order+1)*(order+1); // Penalty used for DG discretizations
111
112  Mesh mesh(mesh_file, 1, 1);
113  int dim = mesh.Dimension();
114  MFEM_VERIFY(dim == 2 || dim == 3, "Spatial dimension must be 2 or 3.");
115  for (int l = 0; l < ref_levels; l++) { mesh.UniformRefinement(); }
116
117  FunctionCoefficient f_coeff(f(1.0)), u_coeff(u);
118  VectorFunctionCoefficient f_vec_coeff(dim, f_vec(RT)), u_vec_coeff(dim, u_vec);
119
120  int b1 = BasisType::GaussLobatto, b2 = BasisType::IntegratedGLL;
121  unique_ptr<FiniteElementCollection> fec;
122  if (H1) { fec.reset(new H1_FECollection(order, dim, b1)); }
123  else if (ND) { fec.reset(new ND_FECollection(order, dim, b1, b2)); }
124  else if (RT) { fec.reset(new RT_FECollection(order-1, dim, b1, b2)); }
125  else { fec.reset(new L2_FECollection(order, dim, b1)); }
126
127  FiniteElementSpace fes(&mesh, fec.get());
128  cout << "Number of DOFs: " << fes.GetTrueVSize() << endl;
129
130  Array<int> ess_dofs;
131  // In DG, boundary conditions are enforced weakly, so no essential DOFs.
132  if (!L2) { fes.GetBoundaryTrueDofs(ess_dofs); }
133
134  BilinearForm a(&fes);
135  if (H1 || L2)
136  {
139  }
140  else
141  {
143  }
144
145  if (ND) { a.AddDomainIntegrator(new CurlCurlIntegrator); }
146  else if (RT) { a.AddDomainIntegrator(new DivDivIntegrator); }
147  else if (L2)
148  {
151  }
152  // TODO: L2 diffusion not implemented with partial assembly
153  if (!L2) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
154  a.Assemble();
155
156  LinearForm b(&fes);
157  if (H1 || L2) { b.AddDomainIntegrator(new DomainLFIntegrator(f_coeff)); }
158  else { b.AddDomainIntegrator(new VectorFEDomainLFIntegrator(f_vec_coeff)); }
159  if (L2)
160  {
161  // DG boundary conditions are enforced weakly with this integrator.
163  }
164  b.Assemble();
165
166  GridFunction x(&fes);
167  if (H1 || L2) { x.ProjectCoefficient(u_coeff);}
168  else { x.ProjectCoefficient(u_vec_coeff); }
169
170  Vector X, B;
171  OperatorHandle A;
172  a.FormLinearSystem(ess_dofs, x, b, A, X, B);
173
174 #ifdef MFEM_USE_SUITESPARSE
175  LORSolver<UMFPackSolver> solv_lor(a, ess_dofs);
176 #else
177  LORSolver<GSSmoother> solv_lor(a, ess_dofs);
178 #endif
179
180  CGSolver cg;
181  cg.SetAbsTol(0.0);
182  cg.SetRelTol(1e-12);
183  cg.SetMaxIter(500);
184  cg.SetPrintLevel(1);
185  cg.SetOperator(*A);
186  cg.SetPreconditioner(solv_lor);
187  cg.Mult(B, X);
188
189  a.RecoverFEMSolution(X, b, x);
190
191  double er =
192  (H1 || L2) ? x.ComputeL2Error(u_coeff) : x.ComputeL2Error(u_vec_coeff);
193  cout << "L2 error: " << er << endl;
194
195  if (visualization)
196  {
197  // Save the solution and mesh to disk. The output can be viewed using
198  // GLVis as follows: "glvis -m mesh.mesh -g sol.gf"
199  x.Save("sol.gf");
200  mesh.Save("mesh.mesh");
201
202  // Also save the solution for visualization using ParaView
203  ParaViewDataCollection dc("LOR", &mesh);
204  dc.SetPrefixPath("ParaView");
205  dc.SetHighOrderOutput(true);
206  dc.SetLevelsOfDetail(order);
207  dc.RegisterField("u", &x);
208  dc.SetCycle(0);
209  dc.SetTime(0.0);
210  dc.Save();
211  }
212
213  return 0;
214 }
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
Definition: gridfunc.cpp:2786
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
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
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:77
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:718
virtual void Save(const std::string &fname, int precision=16) const
Definition: mesh.cpp:10705
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:10232
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:380
A general vector function coefficient.
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:712
void SetAbsTol(double atol)
Definition: solvers.hpp:200
std::function< void(const Vector &, Vector &)> f_vec(bool grad_div_problem)
Definition: lor_mms.hpp:68
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:219
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:2415
void SetLevelsOfDetail(int levels_of_detail_)
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:454
Vector data type.
Definition: vector.hpp:58
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:259
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:24
virtual void Save() override
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3696
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:255
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:327
double f(const Vector &p)