MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
grad_div.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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// H(div) saddle-point system solver
14// ---------------------------------
15//
16// Solves the grad-div problem u - grad(div(u)) = f using a variety of solver
17// techniques. This miniapp supports solving this problem using a variety of
18// matrix-free and matrix-based preconditioning methods, including:
19//
20// * Matrix-free block-diagonal preconditioning for the saddle-point system.
21// * ADS-AMG preconditioning.
22// * Low-order-refined ADS-AMG preconditioning (matrix-free).
23// * Hybridization with AMG preconditioning.
24//
25// The problem setup is the same as in the LOR solvers miniapps (in the
26// miniapps/solvers directory). Dirichlet conditions are enforced on the normal
27// component of u.
28//
29// Sample runs:
30//
31// grad_div -sp -ams -lor -hb
32// mpirun -np 4 grad_div -sp -ams -lor -hb -m ../../data/fichera-q2.mesh -rp 0
33
34#include "mfem.hpp"
35#include <iostream>
36#include <memory>
39
40using namespace std;
41using namespace mfem;
42
43ParMesh LoadParMesh(const char *mesh_file, int ser_ref = 0, int par_ref = 0);
44void SolveCG(Operator &A, Solver &P, const Vector &B, Vector &X);
45
46int main(int argc, char *argv[])
47{
48 Mpi::Init(argc, argv);
50
51 const char *mesh_file = "../../data/star.mesh";
52 const char *device_config = "cpu";
53 int ser_ref = 1;
54 int par_ref = 1;
55 int order = 3;
56 bool use_saddle_point = false;
57 bool use_ams = false;
58 bool use_lor_ams = false;
59 bool use_hybridization = false;
60
61 OptionsParser args(argc, argv);
62 args.AddOption(&device_config, "-d", "--device",
63 "Device configuration string, see Device::Configure().");
64 args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
65 args.AddOption(&ser_ref, "-rs", "--serial-refine",
66 "Number of times to refine the mesh in serial.");
67 args.AddOption(&par_ref, "-rp", "--parallel-refine",
68 "Number of times to refine the mesh in parallel.");
69 args.AddOption(&order, "-o", "--order", "Polynomial degree.");
70 args.AddOption(&use_saddle_point,
71 "-sp", "--saddle-point", "-no-sp", "--no-saddle-point",
72 "Enable or disable saddle-point solver.");
73 args.AddOption(&use_ams, "-ams", "--ams", "-no-ams", "--no-ams",
74 "Enable or disable AMS solver.");
75 args.AddOption(&use_lor_ams, "-lor", "--lor-ams", "-no-lor", "--no-lor-ams",
76 "Enable or disable LOR-AMS solver.");
77 args.AddOption(&use_hybridization,
78 "-hb", "--hybridization", "-no-hb", "--no-hybridization",
79 "Enable or disable hybridization solver.");
80 args.ParseCheck();
81
82 if (!use_saddle_point && !use_ams && !use_lor_ams && !use_hybridization)
83 {
84 if (Mpi::Root()) { cout << "No solver enabled. Exiting.\n"; }
85 return 0;
86 }
87
88 Device device(device_config);
89 if (Mpi::Root()) { device.Print(); }
90
91 ParMesh mesh = LoadParMesh(mesh_file, ser_ref, par_ref);
92 const int dim = mesh.Dimension();
93 MFEM_VERIFY(dim == 2 || dim == 3, "Spatial dimension must be 2 or 3.");
94
96 RT_FECollection fec_rt(order-1, dim, b1, b2);
97 ParFiniteElementSpace fes_rt(&mesh, &fec_rt);
98
99 Array<int> ess_rt_dofs;
100 fes_rt.GetBoundaryTrueDofs(ess_rt_dofs);
101
102 VectorFunctionCoefficient f_vec_coeff(dim, f_vec(true)), u_vec_coeff(dim,
103 u_vec);
104
105 ParLinearForm b(&fes_rt);
106 b.AddDomainIntegrator(new VectorFEDomainLFIntegrator(f_vec_coeff));
107 b.UseFastAssembly(true);
108 b.Assemble();
109
110 ConstantCoefficient alpha_coeff(1.0);
111 ConstantCoefficient beta_coeff(1.0);
112
113 ParGridFunction x(&fes_rt);
114 x.ProjectCoefficient(u_vec_coeff);
115
116 cout.precision(4);
117 cout << scientific;
118
119 if (use_saddle_point)
120 {
121 if (Mpi::Root()) { cout << "\nSaddle point solver... " << flush; }
123
124 const int mt = FiniteElement::INTEGRAL;
125 L2_FECollection fec_l2(order-1, dim, b2, mt);
126 ParFiniteElementSpace fes_l2(&mesh, &fec_l2);
127
128 HdivSaddlePointSolver saddle_point_solver(
129 mesh, fes_rt, fes_l2, alpha_coeff, beta_coeff, ess_rt_dofs,
130 HdivSaddlePointSolver::Mode::GRAD_DIV);
131
132 const Array<int> &offsets = saddle_point_solver.GetOffsets();
133
134 BlockVector X_block(offsets), B_block(offsets);
135 B_block.GetBlock(0) = 0.0;
136 b.ParallelAssemble(B_block.GetBlock(1));
137 B_block.GetBlock(1) *= -1.0;
138 B_block.SyncFromBlocks();
139
140 x.ParallelProject(X_block.GetBlock(1));
141 saddle_point_solver.SetBC(X_block.GetBlock(1));
142
143 X_block = 0.0;
144 saddle_point_solver.Mult(B_block, X_block);
145
146 if (Mpi::Root())
147 {
148 cout << "Done.\nIterations: "
149 << saddle_point_solver.GetNumIterations()
150 << "\nElapsed: " << tic_toc.RealTime() << endl;
151 }
152
153 X_block.SyncToBlocks();
154 x.SetFromTrueDofs(X_block.GetBlock(1));
155 const real_t error = x.ComputeL2Error(u_vec_coeff);
156 if (Mpi::Root()) { cout << "L2 error: " << error << endl; }
157 }
158
159 if (use_ams)
160 {
161 if (Mpi::Root()) { cout << "\nAMS solver... " << flush; }
163
164 ParBilinearForm a(&fes_rt);
165 a.AddDomainIntegrator(new DivDivIntegrator(alpha_coeff));
166 a.AddDomainIntegrator(new VectorFEMassIntegrator(beta_coeff));
167 a.Assemble();
168
170 Vector B, X;
171 b.Assemble();
172 x.ProjectCoefficient(u_vec_coeff);
173 a.FormLinearSystem(ess_rt_dofs, x, b, A, X, B);
174 HypreParMatrix &Ah = *A.As<HypreParMatrix>();
175
176 std::unique_ptr<Solver> prec;
177 if (dim == 2) { prec.reset(new HypreAMS(Ah, &fes_rt)); }
178 else { prec.reset(new HypreADS(Ah, &fes_rt)); }
179
180 SolveCG(Ah, *prec, B, X);
181 x.SetFromTrueDofs(X);
182 const real_t error = x.ComputeL2Error(u_vec_coeff);
183 if (Mpi::Root()) { cout << "L2 error: " << error << endl; }
184 }
185
186 if (use_lor_ams)
187 {
188 const int b2_lor = BasisType::IntegratedGLL;
189 RT_FECollection fec_rt_lor(order-1, dim, b1, b2_lor);
190 ParFiniteElementSpace fes_rt_lor(&mesh, &fec_rt_lor);
191
192 ParLinearForm b_lor(&fes_rt_lor);
193 b_lor.AddDomainIntegrator(new VectorFEDomainLFIntegrator(f_vec_coeff));
194 b_lor.UseFastAssembly(true);
195 b_lor.Assemble();
196
197 if (Mpi::Root()) { cout << "\nLOR-AMS solver... " << flush; }
199
200 ParBilinearForm a(&fes_rt_lor);
201 a.SetAssemblyLevel(AssemblyLevel::PARTIAL);
202 a.AddDomainIntegrator(new DivDivIntegrator(alpha_coeff));
203 a.AddDomainIntegrator(new VectorFEMassIntegrator(beta_coeff));
204 a.Assemble();
205
206 ParGridFunction x_lor(&fes_rt_lor);
207 x_lor.ProjectCoefficient(u_vec_coeff);
208
210 Vector B, X;
211 a.FormLinearSystem(ess_rt_dofs, x_lor, b_lor, A, X, B);
212
213 std::unique_ptr<Solver> prec;
214 if (dim == 2) { prec.reset(new LORSolver<HypreAMS>(a, ess_rt_dofs)); }
215 else { prec.reset(new LORSolver<HypreADS>(a, ess_rt_dofs)); }
216
217 SolveCG(*A, *prec, B, X);
218 a.RecoverFEMSolution(X, b_lor, x_lor);
219 const real_t error = x_lor.ComputeL2Error(u_vec_coeff);
220 if (Mpi::Root()) { cout << "L2 error: " << error << endl; }
221 }
222
223 if (use_hybridization)
224 {
225 if (Mpi::Root()) { cout << "\nHybridization solver... " << flush; }
227
228 DG_Interface_FECollection fec_hb(order-1, dim);
229 ParFiniteElementSpace fes_hb(&mesh, &fec_hb);
230
231 ParBilinearForm a(&fes_rt);
232 a.AddDomainIntegrator(new DivDivIntegrator(alpha_coeff));
233 a.AddDomainIntegrator(new VectorFEMassIntegrator(beta_coeff));
234 a.EnableHybridization(&fes_hb, new NormalTraceJumpIntegrator, ess_rt_dofs);
235 a.Assemble();
236
238 Vector B, X;
239 b.Assemble();
240 x.ProjectCoefficient(u_vec_coeff);
241 a.FormLinearSystem(ess_rt_dofs, x, b, A, X, B);
242
243 HypreBoomerAMG amg_hb(*A.As<HypreParMatrix>());
244 amg_hb.SetPrintLevel(0);
245
246 SolveCG(*A, amg_hb, B, X);
247 a.RecoverFEMSolution(X, b, x);
248 const real_t error = x.ComputeL2Error(u_vec_coeff);
249 if (Mpi::Root()) { cout << "L2 error: " << error << endl; }
250 }
251
252 return 0;
253}
254
255ParMesh LoadParMesh(const char *mesh_file, int ser_ref, int par_ref)
256{
257 Mesh serial_mesh = Mesh::LoadFromFile(mesh_file);
258 for (int i = 0; i < ser_ref; ++i) { serial_mesh.UniformRefinement(); }
259 ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
260 serial_mesh.Clear();
261 for (int i = 0; i < par_ref; ++i) { mesh.UniformRefinement(); }
262 return mesh;
263}
264
265void SolveCG(Operator &A, Solver &P, const Vector &B, Vector &X)
266{
267 CGSolver cg(MPI_COMM_WORLD);
268 cg.SetAbsTol(0.0);
269 cg.SetRelTol(1e-12);
270 cg.SetMaxIter(500);
271 cg.SetPrintLevel(0);
272 cg.SetOperator(A);
273 cg.SetPreconditioner(P);
274 X = 0.0;
275 cg.Mult(B, X);
276 if (Mpi::Root())
277 {
278 cout << "Done.\nIterations: " << cg.GetNumIterations()
279 << "\nElapsed: " << tic_toc.RealTime() << endl;
280 }
281}
@ GaussLobatto
Closed type.
Definition fe_base.hpp:32
@ GaussLegendre
Open type.
Definition fe_base.hpp:31
@ IntegratedGLL
Integrated GLL indicator functions.
Definition fe_base.hpp:39
A class to handle Vectors in a block fashion.
void SyncFromBlocks() const
Synchronize the memory location flags (i.e. the memory validity flags) of the big/monolithic block-ve...
Vector & GetBlock(int i)
Get the i-th vector in the block.
Conjugate gradient method.
Definition solvers.hpp:513
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.hpp:526
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the Conjugate Gradient method.
Definition solvers.cpp:718
A coefficient that is constant across space and time.
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition device.hpp:123
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition device.cpp:286
for Raviart-Thomas elements
void GetBoundaryTrueDofs(Array< int > &boundary_dofs, int component=-1)
Get a list of all boundary true dofs, boundary_dofs. For spaces with 'vdim' > 1, the 'component' para...
Definition fespace.cpp:632
Solve the H(div) saddle-point system using MINRES with matrix-free block-diagonal preconditioning.
const Array< int > & GetOffsets() const
Return the offsets of the block system.
void Mult(const Vector &b, Vector &x) const override
Solve the linear system for L2 (scalar) and RT (flux) unknowns.
void SetBC(const Vector &x_rt)
Sets the Dirichlet boundary conditions at the RT essential DOFs.
int GetNumIterations() const
Get the number of MINRES iterations.
The Auxiliary-space Divergence Solver in hypre.
Definition hypre.hpp:1922
The Auxiliary-space Maxwell Solver in hypre.
Definition hypre.hpp:1845
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
void SetPrintLevel(int print_level)
Definition hypre.hpp:1771
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:173
int GetNumIterations() const
Returns the number of iterations taken during the last call to Mult()
Definition solvers.hpp:260
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
void SetMaxIter(int max_it)
Definition solvers.hpp:211
void SetAbsTol(real_t atol)
Definition solvers.hpp:210
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
Represents a solver of type SolverType created using the low-order refined version of the given Bilin...
Definition lor.hpp:204
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
void UseFastAssembly(bool use_fa)
Which assembly algorithm to use: the new device-compatible fast assembly (true), or the legacy CPU-on...
Mesh data type.
Definition mesh.hpp:56
void Clear()
Clear the contents of the Mesh.
Definition mesh.hpp:730
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
static Mesh LoadFromFile(const std::string &filename, int generate_edges=0, int refine=1, bool fix_orientation=true)
Definition mesh.cpp:4225
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
static bool Root()
Return true if the rank in MPI_COMM_WORLD is zero.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
Pointer to an Operator of a specified type.
Definition handle.hpp:34
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get().
Definition handle.hpp:104
Abstract operator.
Definition operator.hpp:25
void ParseCheck(std::ostream &out=mfem::out)
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
Class for parallel bilinear form.
Abstract parallel finite element space.
Definition pfespace.hpp:29
Class for parallel grid function.
Definition pgridfunc.hpp:33
real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const override
void SetFromTrueDofs(const Vector &tv) override
Set the GridFunction from the given true-dof vector.
void ProjectCoefficient(Coefficient &coeff) override
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
void ParallelProject(Vector &tv) const
Returns the vector restricted to the true dofs.
Class for parallel linear form.
void Assemble()
Assembles the ParLinearForm i.e. sums over all domain/bdr integrators.
Class for parallel meshes.
Definition pmesh.hpp:34
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition fe_coll.hpp:386
Base class for solvers.
Definition operator.hpp:683
double RealTime()
Return the number of real seconds elapsed since the stopwatch was started.
Definition tic_toc.cpp:432
void Start()
Start the stopwatch. The elapsed time is not cleared.
Definition tic_toc.cpp:411
void Clear()
Clear the elapsed time on the stopwatch and restart it if it's running.
Definition tic_toc.cpp:406
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Definition lininteg.hpp:347
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
int dim
Definition ex24.cpp:53
int main()
void SolveCG(Operator &A, Solver &P, const Vector &B, Vector &X)
Definition grad_div.cpp:265
ParMesh LoadParMesh(const char *mesh_file, int ser_ref=0, int par_ref=0)
Definition grad_div.cpp:255
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
void u_vec(const Vector &xvec, Vector &u)
Definition lor_mms.hpp:50
std::function< void(const Vector &, Vector &)> f_vec(bool grad_div_problem)
Definition lor_mms.hpp:68
StopWatch tic_toc
Definition tic_toc.cpp:450
float real_t
Definition config.hpp:43