MFEM  v4.6.0 Finite element discretization library
plor_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 // Parallel 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 hypre's AMG preconditioners:
31 // BoomerAMG is used for H1 and L2 problems, AMS is used for H(curl) and 2D
32 // H(div) problems, and ADS is used for 3D H(div) problems.
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 plor_solvers
49 //
50 // Sample runs:
51 //
52 // mpirun -np 4 plor_solvers -m ../../data/fichera.mesh -fe h
53 // mpirun -np 4 plor_solvers -m ../../data/fichera.mesh -fe n
54 // mpirun -np 4 plor_solvers -m ../../data/fichera.mesh -fe r
55 // mpirun -np 4 plor_solvers -m ../../data/fichera.mesh -fe l
56 // mpirun -np 4 plor_solvers -m ../../data/amr-hex.mesh -fe h -rs 0 -o 2
57 // mpirun -np 4 plor_solvers -m ../../data/amr-hex.mesh -fe l -rs 0 -o 2
58 //
59 // Device sample runs:
60 // mpirun -np 4 plor_solvers -m ../../data/fichera.mesh -fe h -d cuda
61 // mpirun -np 4 plor_solvers -m ../../data/fichera.mesh -fe n -d cuda
62 // mpirun -np 4 plor_solvers -m ../../data/fichera.mesh -fe r -d cuda
63 // mpirun -np 4 plor_solvers -m ../../data/fichera.mesh -fe l -d cuda
64
65 #include "mfem.hpp"
66 #include <fstream>
67 #include <iostream>
68 #include <memory>
69
70 #include "lor_mms.hpp"
71
72 using namespace std;
73 using namespace mfem;
74
75 int main(int argc, char *argv[])
76 {
77  Mpi::Init();
78  Hypre::Init();
79
80  const char *mesh_file = "../../data/star.mesh";
81  int ser_ref_levels = 1, par_ref_levels = 1;
82  int order = 3;
83  const char *fe = "h";
84  const char *device_config = "cpu";
85  bool visualization = true;
86
87  OptionsParser args(argc, argv);
88  args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
90  "Number of times to refine the mesh uniformly in serial.");
92  "Number of times to refine the mesh uniformly in parallel.");
93  args.AddOption(&order, "-o", "--order", "Polynomial degree.");
95  "FE type. h for H1, n for Hcurl, r for Hdiv, l for L2");
97  "--no-visualization",
98  "Enable or disable GLVis visualization.");
100  "Device configuration string, see Device::Configure().");
101  args.ParseCheck();
102
103  Device device(device_config);
104  if (Mpi::Root()) { device.Print(); }
105
106  bool H1 = false, ND = false, RT = false, L2 = false;
107  if (string(fe) == "h") { H1 = true; }
108  else if (string(fe) == "n") { ND = true; }
109  else if (string(fe) == "r") { RT = true; }
110  else if (string(fe) == "l") { L2 = true; }
111  else { MFEM_ABORT("Bad FE type. Must be 'h', 'n', 'r', or 'l'."); }
112
113  double kappa = (order+1)*(order+1); // Penalty used for DG discretizations
114
115  Mesh serial_mesh(mesh_file, 1, 1);
116  int dim = serial_mesh.Dimension();
117  MFEM_VERIFY(dim == 2 || dim == 3, "Spatial dimension must be 2 or 3.");
118  for (int l = 0; l < ser_ref_levels; l++) { serial_mesh.UniformRefinement(); }
119  ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
120  for (int l = 0; l < par_ref_levels; l++) { mesh.UniformRefinement(); }
121  serial_mesh.Clear();
122
123  if (mesh.ncmesh && (RT || ND))
124  { MFEM_ABORT("LOR AMS and ADS solvers are not supported with AMR meshes."); }
125
126  FunctionCoefficient f_coeff(f(1.0)), u_coeff(u);
127  VectorFunctionCoefficient f_vec_coeff(dim, f_vec(RT)), u_vec_coeff(dim, u_vec);
128
129  int b1 = BasisType::GaussLobatto, b2 = BasisType::IntegratedGLL;
130  unique_ptr<FiniteElementCollection> fec;
131  if (H1) { fec.reset(new H1_FECollection(order, dim, b1)); }
132  else if (ND) { fec.reset(new ND_FECollection(order, dim, b1, b2)); }
133  else if (RT) { fec.reset(new RT_FECollection(order-1, dim, b1, b2)); }
134  else { fec.reset(new L2_FECollection(order, dim, b1)); }
135
136  ParFiniteElementSpace fes(&mesh, fec.get());
137  HYPRE_Int ndofs = fes.GlobalTrueVSize();
138  if (Mpi::Root()) { cout << "Number of DOFs: " << ndofs << endl; }
139
140  Array<int> ess_dofs;
141  // In DG, boundary conditions are enforced weakly, so no essential DOFs.
142  if (!L2) { fes.GetBoundaryTrueDofs(ess_dofs); }
143
144  ParBilinearForm a(&fes);
145  if (H1 || L2)
146  {
149  }
150  else
151  {
153  }
154
155  if (ND) { a.AddDomainIntegrator(new CurlCurlIntegrator); }
156  else if (RT) { a.AddDomainIntegrator(new DivDivIntegrator); }
157  else if (L2)
158  {
161  }
162  // TODO: L2 diffusion not implemented with partial assembly
163  if (!L2) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
164  a.Assemble();
165
166  ParLinearForm b(&fes);
167  if (H1 || L2) { b.AddDomainIntegrator(new DomainLFIntegrator(f_coeff)); }
168  else { b.AddDomainIntegrator(new VectorFEDomainLFIntegrator(f_vec_coeff)); }
169  if (L2)
170  {
171  // DG boundary conditions are enforced weakly with this integrator.
173  }
174  b.Assemble();
175
176  ParGridFunction x(&fes);
177  if (H1 || L2) { x.ProjectCoefficient(u_coeff);}
178  else { x.ProjectCoefficient(u_vec_coeff); }
179
180  Vector X, B;
181  OperatorHandle A;
182  a.FormLinearSystem(ess_dofs, x, b, A, X, B);
183
184
185  unique_ptr<Solver> solv_lor;
186  if (H1 || L2)
187  {
188  solv_lor.reset(new LORSolver<HypreBoomerAMG>(a, ess_dofs));
189  }
190  else if (RT && dim == 3)
191  {
193  }
194  else
195  {
196  solv_lor.reset(new LORSolver<HypreAMS>(a, ess_dofs));
197  }
198
199  CGSolver cg(MPI_COMM_WORLD);
200  cg.SetAbsTol(0.0);
201  cg.SetRelTol(1e-12);
202  cg.SetMaxIter(500);
203  cg.SetPrintLevel(1);
204  cg.SetOperator(*A);
205  cg.SetPreconditioner(*solv_lor);
206  cg.Mult(B, X);
207
208  a.RecoverFEMSolution(X, b, x);
209
210  double er =
211  (H1 || L2) ? x.ComputeL2Error(u_coeff) : x.ComputeL2Error(u_vec_coeff);
212  if (Mpi::Root()) { cout << "L2 error: " << er << endl; }
213
214  if (visualization)
215  {
216  // Save the solution and mesh to disk. The output can be viewed using
217  // GLVis as follows: "glvis -np <np> -m mesh -g sol"
218  x.Save("sol");
219  mesh.Save("mesh");
220
221  // Also save the solution for visualization using ParaView
222  ParaViewDataCollection dc("PLOR", &mesh);
223  dc.SetPrefixPath("ParaView");
224  dc.SetHighOrderOutput(true);
225  dc.SetLevelsOfDetail(order);
226  dc.RegisterField("u", &x);
227  dc.SetCycle(0);
228  dc.SetTime(0.0);
229  dc.Save();
230  }
231
232  return 0;
233 }
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
Definition: solvers.hpp:493
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
int main(int argc, char *argv[])
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
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:718
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
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:561
STL namespace.
(Q div u, div v) for RT elements
Class for parallel linear form.
Definition: plinearform.hpp:26
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 Save(const std::string &fname, int precision=16) const override
Definition: pmesh.cpp:4941
void SetHighOrderOutput(bool high_order_output_)
void SetTime(double t)
Set physical time (for time-dependent simulations)
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
Definition: pgridfunc.hpp:285
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:281
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:380
A general vector function coefficient.
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
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
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:909
double a
Definition: lissajous.cpp:41
int dim
Definition: ex24.cpp:53
Class for parallel bilinear form.
NCMesh * ncmesh
Optional nonconforming mesh extension.
Definition: mesh.hpp:278
void SetLevelsOfDetail(int levels_of_detail_)
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:678
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
virtual void Save() override
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
Class for parallel grid function.
Definition: pgridfunc.hpp:32
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
Class for parallel meshes.
Definition: pmesh.hpp:32
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)