MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
plor_solvers.cpp
Go to the documentation of this file.
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 // 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 //
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 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 #include "mfem.hpp"
60 #include <fstream>
61 #include <iostream>
62 #include <memory>
63 
64 #include "lor_mms.hpp"
65 
66 using namespace std;
67 using namespace mfem;
68 
69 bool grad_div_problem = false;
70 
71 int main(int argc, char *argv[])
72 {
73  MPI_Session mpi;
74 
75  const char *mesh_file = "../../data/star.mesh";
76  int ser_ref_levels = 1, par_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(&ser_ref_levels, "-rs", "--refine-serial",
84  "Number of times to refine the mesh uniformly in serial.");
85  args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
86  "Number of times to refine the mesh uniformly in parallel.");
87  args.AddOption(&order, "-o", "--order", "Polynomial degree.");
88  args.AddOption(&fe, "-fe", "--fe-type",
89  "FE type. h for H1, n for Hcurl, r for Hdiv, l for L2");
90  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
91  "--no-visualization",
92  "Enable or disable GLVis visualization.");
93  args.ParseCheck();
94 
95  bool H1 = false, ND = false, RT = false, L2 = false;
96  if (string(fe) == "h") { H1 = true; }
97  else if (string(fe) == "n") { ND = true; }
98  else if (string(fe) == "r") { RT = true; }
99  else if (string(fe) == "l") { L2 = true; }
100  else { MFEM_ABORT("Bad FE type. Must be 'h', 'n', 'r', or 'l'."); }
101 
102  if (RT) { grad_div_problem = true; }
103  double kappa = (order+1)*(order+1); // Penalty used for DG discretizations
104 
105  Mesh serial_mesh(mesh_file, 1, 1);
106  int dim = serial_mesh.Dimension();
107  MFEM_VERIFY(dim == 2 || dim == 3, "Spatial dimension must be 2 or 3.");
108  for (int l = 0; l < ser_ref_levels; l++) { serial_mesh.UniformRefinement(); }
109  ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
110  for (int l = 0; l < par_ref_levels; l++) { mesh.UniformRefinement(); }
111  serial_mesh.Clear();
112 
113  if (mesh.ncmesh && (RT || ND))
114  { MFEM_ABORT("LOR AMS and ADS solvers are not supported with AMR meshes."); }
115 
116  FunctionCoefficient f_coeff(f), u_coeff(u);
117  VectorFunctionCoefficient f_vec_coeff(dim, f_vec), u_vec_coeff(dim, u_vec);
118 
119  int b1 = BasisType::GaussLobatto, b2 = BasisType::IntegratedGLL;
120  unique_ptr<FiniteElementCollection> fec;
121  if (H1) { fec.reset(new H1_FECollection(order, dim, b1)); }
122  else if (ND) { fec.reset(new ND_FECollection(order, dim, b1, b2)); }
123  else if (RT) { fec.reset(new RT_FECollection(order-1, dim, b1, b2)); }
124  else { fec.reset(new L2_FECollection(order, dim, b1)); }
125 
126  ParFiniteElementSpace fes(&mesh, fec.get());
127  HYPRE_Int ndofs = fes.GlobalTrueVSize();
128  if (mpi.Root()) { cout << "Number of DOFs: " << ndofs << 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  ParBilinearForm 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  {
150  a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(-1.0, kappa));
151  }
152  // TODO: L2 diffusion not implemented with partial assembly
153  if (!L2) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
154  a.Assemble();
155 
156  ParLinearForm 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.
162  b.AddBdrFaceIntegrator(new DGDirichletLFIntegrator(u_coeff, -1.0, kappa));
163  }
164  b.Assemble();
165 
166  ParGridFunction 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  ParLORDiscretization lor(a, ess_dofs);
175  ParFiniteElementSpace &fes_lor = lor.GetParFESpace();
176 
177  unique_ptr<Solver> solv_lor;
178  if (H1 || L2)
179  {
180  solv_lor.reset(new LORSolver<HypreBoomerAMG>(lor));
181  }
182  else if (RT && dim == 3)
183  {
184  solv_lor.reset(new LORSolver<HypreADS>(lor, &fes_lor));
185  }
186  else
187  {
188  solv_lor.reset(new LORSolver<HypreAMS>(lor, &fes_lor));
189  }
190 
191  CGSolver cg(MPI_COMM_WORLD);
192  cg.SetAbsTol(0.0);
193  cg.SetRelTol(1e-12);
194  cg.SetMaxIter(500);
195  cg.SetPrintLevel(1);
196  cg.SetOperator(*A);
197  cg.SetPreconditioner(*solv_lor);
198  cg.Mult(B, X);
199 
200  a.RecoverFEMSolution(X, b, x);
201 
202  double er =
203  (H1 || L2) ? x.ComputeL2Error(u_coeff) : x.ComputeL2Error(u_vec_coeff);
204  if (mpi.Root()) { cout << "L2 error: " << er << endl; }
205 
206  if (visualization)
207  {
208  // Save the solution and mesh to disk. The output can be viewed using
209  // GLVis as follows: "glvis -np <np> -m mesh -g sol"
210  x.Save("sol");
211  mesh.Save("mesh");
212 
213  // Also save the solution for visualization using ParaView
214  ParaViewDataCollection dc("PLOR", &mesh);
215  dc.SetPrefixPath("ParaView");
216  dc.SetHighOrderOutput(true);
217  dc.SetLevelsOfDetail(order);
218  dc.RegisterField("u", &x);
219  dc.SetCycle(0);
220  dc.SetTime(0.0);
221  dc.Save();
222  }
223 
224  return 0;
225 }
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:97
Conjugate gradient method.
Definition: solvers.hpp:316
Create and assemble a low-order refined version of a ParBilinearForm.
Definition: lor.hpp:154
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(...
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:616
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
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:275
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:841
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
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:493
virtual void Save(const char *fname, int precision=16) const
Definition: pmesh.cpp:4499
(Q div u, div v) for RT elements
bool grad_div_problem
Definition: lor_solvers.cpp:71
Class for parallel linear form.
Definition: plinearform.hpp:26
A simple convenience class that calls MPI_Init() at construction and MPI_Finalize() at destruction...
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:71
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9143
void SetMaxIter(int max_it)
Definition: solvers.hpp:100
void Assemble(int skip_zeros=1)
Assemble the local matrix.
void u_vec(const Vector &xvec, Vector &u)
Definition: lor_mms.hpp:50
bool Root() const
Return true if WorldRank() == 0.
void SetHighOrderOutput(bool high_order_output_)
void AddBdrFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:83
int Dimension() const
Definition: mesh.hpp:911
void SetTime(double t)
Set physical time (for time-dependent simulations)
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:341
A general vector function coefficient.
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:39
void SetAbsTol(double atol)
Definition: solvers.hpp:99
void SetRelTol(double rtol)
Definition: solvers.hpp:98
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 RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
double a
Definition: lissajous.cpp:41
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const
Definition: pgridfunc.hpp:273
ParFiniteElementSpace & GetParFESpace() const
Return the LOR ParFiniteElementSpace.
Definition: lor.cpp:448
int dim
Definition: ex24.cpp:53
void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new interior Face Integrator. Assumes ownership of bfi.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
Class for parallel bilinear form.
NCMesh * ncmesh
Optional non-conforming mesh extension.
Definition: mesh.hpp:207
void SetLevelsOfDetail(int levels_of_detail_)
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:828
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Definition: lininteg.hpp:266
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:330
A general function coefficient.
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:416
Vector data type.
Definition: vector.hpp:60
void f_vec(const Vector &xvec, Vector &f)
Definition: lor_mms.hpp:68
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:92
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
virtual void Save() override
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Class for parallel meshes.
Definition: pmesh.hpp:32
int main()
void ParseCheck(std::ostream &out=mfem::out)
Definition: optparser.cpp:251
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:190
void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new boundary Face Integrator. Assumes ownership of bfi.
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:285
double f(const Vector &p)