MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
ex3p.cpp
Go to the documentation of this file.
1// MFEM Example 3 - Parallel Version
2// PETSc Modification
3//
4// Compile with: make ex3p
5//
6// Sample runs:
7// mpirun -np 4 ex3p -m ../../data/klein-bottle.mesh -o 2 -f 0.1 --petscopts rc_ex3p
8// mpirun -np 4 ex3p -m ../../data/klein-bottle.mesh -o 2 -f 0.1 --petscopts rc_ex3p_bddc --nonoverlapping
9//
10// Description: This example code solves a simple electromagnetic diffusion
11// problem corresponding to the second order definite Maxwell
12// equation curl curl E + E = f with boundary condition
13// E x n = <given tangential field>. Here, we use a given exact
14// solution E and compute the corresponding r.h.s. f.
15// We discretize with Nedelec finite elements in 2D or 3D.
16//
17// The example demonstrates the use of H(curl) finite element
18// spaces with the curl-curl and the (vector finite element) mass
19// bilinear form, as well as the computation of discretization
20// error when the exact solution is known. Static condensation is
21// also illustrated.
22//
23// The example also show how to use the non-overlapping feature of
24// the ParBilinearForm class to obtain the linear operator in
25// a format suitable for the BDDC preconditioner in PETSc.
26//
27// We recommend viewing examples 1-2 before viewing this example.
28
29#include "mfem.hpp"
30#include <fstream>
31#include <iostream>
32
33#ifndef MFEM_USE_PETSC
34#error This example requires that MFEM is built with MFEM_USE_PETSC=YES
35#endif
36
37using namespace std;
38using namespace mfem;
39
40// Exact solution, E, and r.h.s., f. See below for implementation.
41void E_exact(const Vector &, Vector &);
42void f_exact(const Vector &, Vector &);
44int dim;
45
46int main(int argc, char *argv[])
47{
48 // 1. Initialize MPI and HYPRE.
49 Mpi::Init(argc, argv);
50 int num_procs = Mpi::WorldSize();
51 int myid = Mpi::WorldRank();
53
54 // 2. Parse command-line options.
55 const char *mesh_file = "../../data/beam-tet.mesh";
56 int ser_ref_levels = -1;
57 int par_ref_levels = 2;
58 int order = 1;
59 bool static_cond = false;
60 bool visualization = 1;
61 bool use_petsc = true;
62 const char *petscrc_file = "";
63 bool use_nonoverlapping = false;
64
65 OptionsParser args(argc, argv);
67 "Mesh file to use.");
69 "Number of times to refine the mesh uniformly in serial.");
71 "Number of times to refine the mesh uniformly in parallel.");
73 "Finite element order (polynomial degree).");
74 args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
75 " solution.");
77 "--no-static-condensation", "Enable static condensation.");
79 "--no-visualization",
80 "Enable or disable GLVis visualization.");
82 "--no-petsc",
83 "Use or not PETSc to solve the linear system.");
85 "PetscOptions file to use.");
87 "-no-nonoverlapping", "--no-nonoverlapping",
88 "Use or not the block diagonal PETSc's matrix format "
89 "for non-overlapping domain decomposition.");
90 args.Parse();
91 if (!args.Good())
92 {
93 if (myid == 0)
94 {
95 args.PrintUsage(cout);
96 }
97 return 1;
98 }
99 if (myid == 0)
100 {
101 args.PrintOptions(cout);
102 }
103 // 2b. We initialize PETSc
104 if (use_petsc) { MFEMInitializePetsc(NULL,NULL,petscrc_file,NULL); }
105 kappa = freq * M_PI;
106
107 // 3. Read the (serial) mesh from the given mesh file on all processors. We
108 // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
109 // and volume meshes with the same code.
110 Mesh *mesh = new Mesh(mesh_file, 1, 1);
111 dim = mesh->Dimension();
112#if PETSC_VERSION_LT(3,21,0)
113 if (dim == 3 && use_petsc && use_nonoverlapping)
114 {
115 cout << "\nFor three-dimensional runs you need a version of PETSc greater or equal 3.21.\n\n";
116 delete mesh;
119 return MFEM_SKIP_RETURN_VALUE;
120 }
121#endif
122 int sdim = mesh->SpaceDimension();
123
124 // 4. Refine the serial mesh on all processors to increase the resolution. In
125 // this example we do 'ref_levels' of uniform refinement. We choose
126 // 'ref_levels' to be the largest number that gives a final mesh with no
127 // more than 1,000 elements.
128 {
129 if (ser_ref_levels < 0)
130 {
131 ser_ref_levels = (int)floor(log(1000./mesh->GetNE())/log(2.)/dim);
132 }
133 for (int l = 0; l < ser_ref_levels; l++)
134 {
135 mesh->UniformRefinement();
136 }
137 }
138
139 // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
140 // this mesh further in parallel to increase the resolution. Once the
141 // parallel mesh is defined, the serial mesh can be deleted.
142 ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
143 delete mesh;
144 {
145 for (int l = 0; l < par_ref_levels; l++)
146 {
147 pmesh->UniformRefinement();
148 }
149 }
150
151 // 6. Define a parallel finite element space on the parallel mesh. Here we
152 // use the Nedelec finite elements of the specified order.
154 ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
155 HYPRE_BigInt size = fespace->GlobalTrueVSize();
156 if (myid == 0)
157 {
158 cout << "Number of finite element unknowns: " << size << endl;
159 }
160
161 // 7. Determine the list of true (i.e. parallel conforming) essential
162 // boundary dofs. In this example, the boundary conditions are defined
163 // by marking all the boundary attributes from the mesh as essential
164 // (Dirichlet) and converting them to a list of true dofs.
165 Array<int> ess_tdof_list;
166 if (pmesh->bdr_attributes.Size())
167 {
168 Array<int> ess_bdr(pmesh->bdr_attributes.Max());
169 ess_bdr = 1;
170 fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
171 }
172
173 // 8. Set up the parallel linear form b(.) which corresponds to the
174 // right-hand side of the FEM linear system, which in this case is
175 // (f,phi_i) where f is given by the function f_exact and phi_i are the
176 // basis functions in the finite element fespace.
178 ParLinearForm *b = new ParLinearForm(fespace);
180 b->Assemble();
181
182 // 9. Define the solution vector x as a parallel finite element grid function
183 // corresponding to fespace. Initialize x by projecting the exact
184 // solution. Note that only values from the boundary edges will be used
185 // when eliminating the non-homogeneous boundary condition to modify the
186 // r.h.s. vector b.
187 ParGridFunction x(fespace);
190
191 // 10. Set up the parallel bilinear form corresponding to the EM diffusion
192 // operator curl muinv curl + sigma I, by adding the curl-curl and the
193 // mass domain integrators.
194 Coefficient *muinv = new ConstantCoefficient(1.0);
196 ParBilinearForm *a = new ParBilinearForm(fespace);
199
200 // 11. Assemble the parallel bilinear form and the corresponding linear
201 // system, applying any necessary transformations such as: parallel
202 // assembly, eliminating boundary conditions, applying conforming
203 // constraints for non-conforming AMR, static condensation, etc.
204 if (static_cond) { a->EnableStaticCondensation(); }
205 a->Assemble();
206
207 Vector B, X;
208 if (!use_petsc)
209 {
211 a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
212
213 if (myid == 0)
214 {
215 cout << "Size of linear system: " << A.GetGlobalNumRows() << endl;
216 }
217
218 // 12. Define and apply a parallel PCG solver for AX=B with the AMS
219 // preconditioner from hypre.
220 ParFiniteElementSpace *prec_fespace =
221 (a->StaticCondensationIsEnabled() ? a->SCParFESpace() : fespace);
222 HypreSolver *ams = new HypreAMS(A, prec_fespace);
223 HyprePCG *pcg = new HyprePCG(A);
224 pcg->SetTol(1e-10);
225 pcg->SetMaxIter(500);
226 pcg->SetPrintLevel(2);
227 pcg->SetPreconditioner(*ams);
228 pcg->Mult(B, X);
229 delete pcg;
230 delete ams;
231 }
232 else
233 {
235 a->SetOperatorType(use_nonoverlapping ?
237 a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
238
239 if (myid == 0)
240 {
241 cout << "Size of linear system: " << A.M() << endl;
242 }
243
244 // 12. Define and apply a parallel PCG solver.
245 ParFiniteElementSpace *prec_fespace =
246 (a->StaticCondensationIsEnabled() ? a->SCParFESpace() : fespace);
247 PetscPCGSolver *pcg = new PetscPCGSolver(A);
248 PetscPreconditioner *prec = NULL;
249 pcg->SetTol(1e-10);
250 pcg->SetMaxIter(500);
251 pcg->SetPrintLevel(2);
252 if (use_nonoverlapping)
253 {
254 // Auxiliary class for BDDC customization
256 // Inform the solver about the finite element space
257 opts.SetSpace(prec_fespace);
258 // Inform the solver about essential dofs
259 opts.SetEssBdrDofs(&ess_tdof_list);
260 // Create a BDDC solver with parameters
261 prec = new PetscBDDCSolver(A,opts);
262 }
263 else
264 {
265 // Create an empty preconditioner object that can
266 // be customized at runtime
267 prec = new PetscPreconditioner(A,"solver_");
268 }
269 pcg->SetPreconditioner(*prec);
270 pcg->Mult(B, X);
271 delete pcg;
272 delete prec;
273 }
274
275 // 13. Recover the parallel grid function corresponding to X. This is the
276 // local finite element solution on each processor.
277 a->RecoverFEMSolution(X, *b, x);
278
279 // 14. Compute and print the L^2 norm of the error.
280 {
282 if (myid == 0)
283 {
284 cout << "\n|| E_h - E ||_{L^2} = " << err << '\n' << endl;
285 }
286 }
287
288 // 15. Save the refined mesh and the solution in parallel. This output can
289 // be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
290 {
291 ostringstream mesh_name, sol_name;
292 mesh_name << "mesh." << setfill('0') << setw(6) << myid;
293 sol_name << "sol." << setfill('0') << setw(6) << myid;
294
295 ofstream mesh_ofs(mesh_name.str().c_str());
296 mesh_ofs.precision(8);
297 pmesh->Print(mesh_ofs);
298
299 ofstream sol_ofs(sol_name.str().c_str());
300 sol_ofs.precision(8);
301 x.Save(sol_ofs);
302 }
303
304 // 16. Send the solution by socket to a GLVis server.
305 if (visualization)
306 {
307 char vishost[] = "localhost";
308 int visport = 19916;
309 socketstream sol_sock(vishost, visport);
310 sol_sock << "parallel " << num_procs << " " << myid << "\n";
311 sol_sock.precision(8);
312 sol_sock << "solution\n" << *pmesh << x << flush;
313 }
314
315 // 17. Free the used memory.
316 delete a;
317 delete sigma;
318 delete muinv;
319 delete b;
320 delete fespace;
321 delete fec;
322 delete pmesh;
323
324 // We finalize PETSc
325 if (use_petsc) { MFEMFinalizePetsc(); }
326
327 return 0;
328}
329
330
331void E_exact(const Vector &x, Vector &E)
332{
333 if (dim == 3)
334 {
335 E(0) = sin(kappa * x(1));
336 E(1) = sin(kappa * x(2));
337 E(2) = sin(kappa * x(0));
338 }
339 else
340 {
341 E(0) = sin(kappa * x(1));
342 E(1) = sin(kappa * x(0));
343 if (x.Size() == 3) { E(2) = 0.0; }
344 }
345}
346
347void f_exact(const Vector &x, Vector &f)
348{
349 if (dim == 3)
350 {
351 f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
352 f(1) = (1. + kappa * kappa) * sin(kappa * x(2));
353 f(2) = (1. + kappa * kappa) * sin(kappa * x(0));
354 }
355 else
356 {
357 f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
358 f(1) = (1. + kappa * kappa) * sin(kappa * x(0));
359 if (x.Size() == 3) { f(2) = 0.0; }
360 }
361}
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
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.
Integrator for for Nedelec 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 Maxwell Solver in hypre.
Definition hypre.hpp:1845
PCG solver in hypre.
Definition hypre.hpp:1275
void SetPrintLevel(int print_lvl)
Definition hypre.cpp:4156
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition hypre.cpp:4161
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre's PCG.
Definition hypre.cpp:4184
void SetMaxIter(int max_iter)
Definition hypre.cpp:4146
void SetTol(real_t tol)
Definition hypre.cpp:4136
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
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 void Finalize()
Finalize MPI (if it has been initialized and not yet already finalized).
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).
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition fe_coll.hpp:465
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
void SetPreconditioner(Solver &precond)
Definition petsc.cpp:3110
virtual void Mult(const Vector &b, Vector &x) const
Application of the solver.
Definition petsc.cpp:3190
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
void SetTol(real_t tol)
Definition petsc.cpp:2371
void SetPrintLevel(int plev)
Definition petsc.cpp:2453
void SetMaxIter(int max_iter)
Definition petsc.cpp:2426
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
real_t sigma(const Vector &x)
Definition maxwell.cpp:102
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 kappa
Definition ex3p.cpp:43
int dim
Definition ex3p.cpp:44
void E_exact(const Vector &, Vector &)
Definition ex3p.cpp:331
void f_exact(const Vector &, Vector &)
Definition ex3p.cpp:347
real_t freq
Definition ex3p.cpp:43