MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
ex4p.cpp
Go to the documentation of this file.
1// MFEM Example 4 - Parallel Version
2// PETSc Modification
3//
4// Compile with: make ex4p
5//
6// Sample runs:
7// mpirun -np 4 ex4p -m ../../data/klein-bottle.mesh -o 2 --petscopts rc_ex4p
8// mpirun -np 4 ex4p -m ../../data/klein-bottle.mesh -o 2 --petscopts rc_ex4p_bddc --nonoverlapping
9//
10// Description: This example code solves a simple 2D/3D H(div) diffusion
11// problem corresponding to the second order definite equation
12// -grad(alpha div F) + beta F = f with boundary condition F dot n
13// = <given normal field>. Here, we use a given exact solution F
14// and compute the corresponding r.h.s. f. We discretize with
15// Raviart-Thomas finite elements.
16//
17// The example demonstrates the use of H(div) finite element
18// spaces with the grad-div and H(div) vector finite element mass
19// bilinear form, as well as the computation of discretization
20// error when the exact solution is known. Bilinear form
21// hybridization and static condensation are also illustrated.
22//
23// We recommend viewing examples 1-3 before viewing this example.
24
25#include "mfem.hpp"
26#include <fstream>
27#include <iostream>
28
29#ifndef MFEM_USE_PETSC
30#error This example requires that MFEM is built with MFEM_USE_PETSC=YES
31#endif
32
33using namespace std;
34using namespace mfem;
35
36// Exact solution, F, and r.h.s., f. See below for implementation.
37void F_exact(const Vector &, Vector &);
38void f_exact(const Vector &, Vector &);
40
41int main(int argc, char *argv[])
42{
43 // 1. Initialize MPI and HYPRE.
44 Mpi::Init(argc, argv);
45 int num_procs = Mpi::WorldSize();
46 int myid = Mpi::WorldRank();
48
49 // 2. Parse command-line options.
50 const char *mesh_file = "../../data/star.mesh";
51 int ser_ref_levels = -1;
52 int par_ref_levels = 2;
53 int order = 1;
54 bool set_bc = true;
55 bool static_cond = false;
56 bool hybridization = false;
57 bool visualization = 1;
58 bool use_petsc = true;
59 const char *petscrc_file = "";
60 bool use_nonoverlapping = false;
61
62 OptionsParser args(argc, argv);
64 "Mesh file to use.");
66 "Number of times to refine the mesh uniformly in serial.");
68 "Number of times to refine the mesh uniformly in parallel.");
70 "Finite element order (polynomial degree).");
71 args.AddOption(&set_bc, "-bc", "--impose-bc", "-no-bc", "--dont-impose-bc",
72 "Impose or not essential boundary conditions.");
73 args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
74 " solution.");
76 "--no-static-condensation", "Enable static condensation.");
78 "--no-hybridization", "Enable hybridization.");
80 "--no-visualization",
81 "Enable or disable GLVis visualization.");
83 "--no-petsc",
84 "Use or not PETSc to solve the linear system.");
86 "PetscOptions file to use.");
88 "-no-nonoverlapping", "--no-nonoverlapping",
89 "Use or not the block diagonal PETSc's matrix format "
90 "for non-overlapping domain decomposition.");
91 args.Parse();
92 if (!args.Good())
93 {
94 if (myid == 0)
95 {
96 args.PrintUsage(cout);
97 }
98 return 1;
99 }
100 if (myid == 0)
101 {
102 args.PrintOptions(cout);
103 }
104 // 2b. We initialize PETSc
105 if (use_petsc) { MFEMInitializePetsc(NULL,NULL,petscrc_file,NULL); }
106 kappa = freq * M_PI;
107
108 // 3. Read the (serial) mesh from the given mesh file on all processors. We
109 // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
110 // and volume, as well as periodic meshes with the same code.
111 Mesh *mesh = new Mesh(mesh_file, 1, 1);
112 int dim = mesh->Dimension();
113 int sdim = mesh->SpaceDimension();
114
115 // 4. Refine the serial mesh on all processors to increase the resolution. In
116 // this example we do 'ref_levels' of uniform refinement. We choose
117 // 'ref_levels' to be the largest number that gives a final mesh with no
118 // more than 1,000 elements.
119 {
120 if (ser_ref_levels < 0)
121 {
122 ser_ref_levels = (int)floor(log(1000./mesh->GetNE())/log(2.)/dim);
123 }
124 for (int l = 0; l < ser_ref_levels; l++)
125 {
126 mesh->UniformRefinement();
127 }
128 }
129
130 // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
131 // this mesh further in parallel to increase the resolution. Once the
132 // parallel mesh is defined, the serial mesh can be deleted.
133 ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
134 delete mesh;
135 {
136 for (int l = 0; l < par_ref_levels; l++)
137 {
138 pmesh->UniformRefinement();
139 }
140 }
141
142 // 6. Define a parallel finite element space on the parallel mesh. Here we
143 // use the Raviart-Thomas finite elements of the specified order.
144 FiniteElementCollection *fec = new RT_FECollection(order-1, dim);
145 ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
146 HYPRE_BigInt size = fespace->GlobalTrueVSize();
147 if (myid == 0)
148 {
149 cout << "Number of finite element unknowns: " << size << endl;
150 }
151
152 // 7. Determine the list of true (i.e. parallel conforming) essential
153 // boundary dofs. In this example, the boundary conditions are defined
154 // by marking all the boundary attributes from the mesh as essential
155 // (Dirichlet) and converting them to a list of true dofs.
156 Array<int> ess_tdof_list;
157 if (pmesh->bdr_attributes.Size())
158 {
159 Array<int> ess_bdr(pmesh->bdr_attributes.Max());
160 ess_bdr = set_bc ? 1 : 0;
161 fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
162 }
163
164 // 8. Set up the parallel linear form b(.) which corresponds to the
165 // right-hand side of the FEM linear system, which in this case is
166 // (f,phi_i) where f is given by the function f_exact and phi_i are the
167 // basis functions in the finite element fespace.
169 ParLinearForm *b = new ParLinearForm(fespace);
171 b->Assemble();
172
173 // 9. Define the solution vector x as a parallel finite element grid function
174 // corresponding to fespace. Initialize x by projecting the exact
175 // solution. Note that only values from the boundary faces will be used
176 // when eliminating the non-homogeneous boundary condition to modify the
177 // r.h.s. vector b.
178 ParGridFunction x(fespace);
181
182 // 10. Set up the parallel bilinear form corresponding to the H(div)
183 // diffusion operator grad alpha div + beta I, by adding the div-div and
184 // the mass domain integrators.
187 ParBilinearForm *a = new ParBilinearForm(fespace);
190
191 // 11. Assemble the parallel bilinear form and the corresponding linear
192 // system, applying any necessary transformations such as: parallel
193 // assembly, eliminating boundary conditions, applying conforming
194 // constraints for non-conforming AMR, static condensation,
195 // hybridization, etc.
196 FiniteElementCollection *hfec = NULL;
197 ParFiniteElementSpace *hfes = NULL;
198 if (static_cond)
199 {
200 a->EnableStaticCondensation();
201 }
202 else if (hybridization)
203 {
204 hfec = new DG_Interface_FECollection(order-1, dim);
205 hfes = new ParFiniteElementSpace(pmesh, hfec);
206 a->EnableHybridization(hfes, new NormalTraceJumpIntegrator(),
207 ess_tdof_list);
208 }
209 a->Assemble();
210
211 Vector B, X;
212 CGSolver *pcg = new CGSolver(MPI_COMM_WORLD);
213 pcg->SetRelTol(1e-12);
214 pcg->SetMaxIter(500);
215 pcg->SetPrintLevel(1);
216
217 if (!use_petsc)
218 {
220 a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
221
222 HYPRE_BigInt glob_size = A.GetGlobalNumRows();
223 if (myid == 0)
224 {
225 cout << "Size of linear system: " << glob_size << endl;
226 }
227
228 // 12. Define and apply a parallel PCG solver for A X = B with the 2D AMS or
229 // the 3D ADS preconditioners from hypre. If using hybridization, the
230 // system is preconditioned with hypre's BoomerAMG.
231 HypreSolver *prec = NULL;
232 pcg->SetOperator(A);
233 if (hybridization) { prec = new HypreBoomerAMG(A); }
234 else
235 {
236 ParFiniteElementSpace *prec_fespace =
237 (a->StaticCondensationIsEnabled() ? a->SCParFESpace() : fespace);
238 if (dim == 2) { prec = new HypreAMS(A, prec_fespace); }
239 else { prec = new HypreADS(A, prec_fespace); }
240 }
241 pcg->SetPreconditioner(*prec);
242 pcg->Mult(B, X);
243 delete prec;
244 }
245 else
246 {
248 PetscPreconditioner *prec = NULL;
249 a->SetOperatorType(use_nonoverlapping ?
251 a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
252
253 if (myid == 0)
254 {
255 cout << "Size of linear system: " << A.M() << endl;
256 }
257
258 pcg->SetOperator(A);
259 if (use_nonoverlapping)
260 {
261 ParFiniteElementSpace *prec_fespace =
262 (a->StaticCondensationIsEnabled() ? a->SCParFESpace() :
263 (hfes ? NULL : fespace));
264
265 // Auxiliary class for BDDC customization
267 // Inform the solver about the finite element space
268 opts.SetSpace(prec_fespace);
269 // Inform the solver about essential dofs
270 opts.SetEssBdrDofs(&ess_tdof_list);
271 // Create a BDDC solver with parameters
272 prec = new PetscBDDCSolver(A, opts);
273 }
274 else
275 {
276 // Create an empty preconditioner that can be customized at runtime.
277 prec = new PetscPreconditioner(A, "solver_");
278 }
279 pcg->SetPreconditioner(*prec);
280 pcg->Mult(B, X);
281 delete prec;
282 }
283 delete pcg;
284
285 // 13. Recover the parallel grid function corresponding to X. This is the
286 // local finite element solution on each processor.
287 a->RecoverFEMSolution(X, *b, x);
288
289 // 14. Compute and print the L^2 norm of the error.
290 {
292 if (myid == 0)
293 {
294 cout << "\n|| F_h - F ||_{L^2} = " << err << '\n' << endl;
295 }
296 }
297
298 // 15. Save the refined mesh and the solution in parallel. This output can
299 // be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
300 {
301 ostringstream mesh_name, sol_name;
302 mesh_name << "mesh." << setfill('0') << setw(6) << myid;
303 sol_name << "sol." << setfill('0') << setw(6) << myid;
304
305 ofstream mesh_ofs(mesh_name.str().c_str());
306 mesh_ofs.precision(8);
307 pmesh->Print(mesh_ofs);
308
309 ofstream sol_ofs(sol_name.str().c_str());
310 sol_ofs.precision(8);
311 x.Save(sol_ofs);
312 }
313
314 // 16. Send the solution by socket to a GLVis server.
315 if (visualization)
316 {
317 char vishost[] = "localhost";
318 int visport = 19916;
319 socketstream sol_sock(vishost, visport);
320 sol_sock << "parallel " << num_procs << " " << myid << "\n";
321 sol_sock.precision(8);
322 sol_sock << "solution\n" << *pmesh << x << flush;
323 }
324
325 // 17. Free the used memory.
326 delete hfes;
327 delete hfec;
328 delete a;
329 delete alpha;
330 delete beta;
331 delete b;
332 delete fespace;
333 delete fec;
334 delete pmesh;
335
336 // We finalize PETSc
337 if (use_petsc) { MFEMFinalizePetsc(); }
338
339 return 0;
340}
341
342
343// The exact solution (for non-surface meshes)
344void F_exact(const Vector &p, Vector &F)
345{
346 int dim = p.Size();
347
348 real_t x = p(0);
349 real_t y = p(1);
350 // real_t z = (dim == 3) ? p(2) : 0.0;
351
352 F(0) = cos(kappa*x)*sin(kappa*y);
353 F(1) = cos(kappa*y)*sin(kappa*x);
354 if (dim == 3)
355 {
356 F(2) = 0.0;
357 }
358}
359
360// The right hand side
361void f_exact(const Vector &p, Vector &f)
362{
363 int dim = p.Size();
364
365 real_t x = p(0);
366 real_t y = p(1);
367 // real_t z = (dim == 3) ? p(2) : 0.0;
368
369 real_t temp = 1 + 2*kappa*kappa;
370
371 f(0) = temp*cos(kappa*x)*sin(kappa*y);
372 f(1) = temp*cos(kappa*y)*sin(kappa*x);
373 if (dim == 3)
374 {
375 f(2) = 0;
376 }
377}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
int Size() const
Return the logical size of the array.
Definition array.hpp:144
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
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
A coefficient that is constant across space and time.
for Raviart-Thomas elements
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
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
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
HYPRE_BigInt GetGlobalNumRows() const
Return the global number of rows.
Definition hypre.hpp:679
Abstract class for hypre's solvers and preconditioners.
Definition hypre.hpp:1162
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
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
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1163
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static int WorldSize()
Return the size of MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, Operator *&A, Vector &X, Vector &B, int copy_interior=0)
Form a constrained linear system using a matrix-free approach.
Definition operator.cpp:114
@ PETSC_MATIS
ID for class PetscParMatrix, MATIS format.
Definition operator.hpp:289
@ PETSC_MATAIJ
ID for class PetscParMatrix, MATAIJ format.
Definition operator.hpp:288
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Reconstruct a solution vector x (e.g. a GridFunction) from the solution X of a constrained linear sys...
Definition operator.cpp:148
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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
bool Good() const
Return true if the command line options were parsed successfully.
Class for parallel bilinear form.
Abstract parallel finite element space.
Definition pfespace.hpp:29
void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const override
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
Class for parallel grid function.
Definition pgridfunc.hpp:33
void Save(std::ostream &out) const override
real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const override
void ProjectCoefficient(Coefficient &coeff) override
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Class for parallel linear form.
Class for parallel meshes.
Definition pmesh.hpp:34
void Print(std::ostream &out=mfem::out, const std::string &comments="") const override
Definition pmesh.cpp:4801
Auxiliary class for BDDC customization.
Definition petsc.hpp:831
void SetEssBdrDofs(const Array< int > *essdofs, bool loc=false)
Specify dofs on the essential boundary.
Definition petsc.hpp:850
void SetSpace(ParFiniteElementSpace *fe)
Definition petsc.hpp:845
Wrapper for PETSc's matrix class.
Definition petsc.hpp:319
PetscInt M() const
Returns the global number of rows.
Definition petsc.cpp:1005
Abstract class for PETSc's preconditioners.
Definition petsc.hpp:805
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition fe_coll.hpp:386
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Definition lininteg.hpp:347
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
Vector beta
const real_t alpha
Definition ex15.cpp:369
int dim
Definition ex24.cpp:53
int main()
HYPRE_Int HYPRE_BigInt
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
const int visport
void MFEMInitializePetsc()
Convenience functions to initialize/finalize PETSc.
Definition petsc.cpp:201
void MFEMFinalizePetsc()
Definition petsc.cpp:241
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition globals.hpp:71
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
const char vishost[]
real_t p(const Vector &x, real_t t)
real_t kappa
Definition ex4p.cpp:39
void f_exact(const Vector &, Vector &)
Definition ex4p.cpp:361
real_t freq
Definition ex4p.cpp:39
void F_exact(const Vector &, Vector &)
Definition ex4p.cpp:344