MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
ex13p.cpp
Go to the documentation of this file.
1// MFEM Example 13 - Parallel Version
2//
3// Compile with: make ex13p
4//
5// Sample runs: mpirun -np 4 ex13p -m ../data/star.mesh
6// mpirun -np 4 ex13p -m ../data/square-disc.mesh -o 2 -n 4
7// mpirun -np 4 ex13p -m ../data/beam-tet.mesh
8// mpirun -np 4 ex13p -m ../data/beam-tet.mesh -nc -o 2 -rs 1
9// mpirun -np 4 ex13p -m ../data/beam-hex.mesh
10// mpirun -np 4 ex13p -m ../data/escher.mesh
11// mpirun -np 4 ex13p -m ../data/fichera.mesh
12// mpirun -np 4 ex13p -m ../data/fichera-q2.vtk
13// mpirun -np 4 ex13p -m ../data/fichera-q3.mesh
14// mpirun -np 4 ex13p -m ../data/square-disc-nurbs.mesh
15// mpirun -np 4 ex13p -m ../data/beam-hex-nurbs.mesh
16// mpirun -np 4 ex13p -m ../data/amr-quad.mesh -o 2
17// mpirun -np 4 ex13p -m ../data/amr-hex.mesh
18// mpirun -np 4 ex13p -m ../data/mobius-strip.mesh -n 8 -o 2
19// mpirun -np 4 ex13p -m ../data/klein-bottle.mesh -n 10 -o 2
20//
21// Description: This example code solves the Maxwell (electromagnetic)
22// eigenvalue problem curl curl E = lambda E with homogeneous
23// Dirichlet boundary conditions E x n = 0.
24//
25// We compute a number of the lowest nonzero eigenmodes by
26// discretizing the curl curl operator using a Nedelec FE space of
27// the specified order in 2D or 3D.
28//
29// The example highlights the use of the AME subspace eigenvalue
30// solver from HYPRE, which uses LOBPCG and AMS internally.
31// Reusing a single GLVis visualization window for multiple
32// eigenfunctions is also illustrated.
33//
34// We recommend viewing examples 3 and 11 before viewing this
35// example.
36
37#include "mfem.hpp"
38#include <fstream>
39#include <iostream>
40
41using namespace std;
42using namespace mfem;
43
44int main(int argc, char *argv[])
45{
46 // 1. Initialize MPI and HYPRE.
47 Mpi::Init(argc, argv);
48 int num_procs = Mpi::WorldSize();
49 int myid = Mpi::WorldRank();
51
52 // 2. Parse command-line options.
53 const char *mesh_file = "../data/beam-tet.mesh";
54 int ser_ref_levels = 2;
55 int par_ref_levels = 1;
56 int order = 1;
57 int nev = 5;
58 bool nc = false;
59 bool visualization = 1;
60 const char *device_config = "cpu";
61
62 OptionsParser args(argc, argv);
63 args.AddOption(&mesh_file, "-m", "--mesh",
64 "Mesh file to use.");
65 args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
66 "Number of times to refine the mesh uniformly in serial.");
67 args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
68 "Number of times to refine the mesh uniformly in parallel.");
69 args.AddOption(&order, "-o", "--order",
70 "Finite element order (polynomial degree) or -1 for"
71 " isoparametric space.");
72 args.AddOption(&nev, "-n", "--num-eigs",
73 "Number of desired eigenmodes.");
74 args.AddOption(&nc, "-nc", "--non-conforming", "-c",
75 "--conforming",
76 "Mark the mesh as nonconforming before partitioning.");
77 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
78 "--no-visualization",
79 "Enable or disable GLVis visualization.");
80 args.AddOption(&device_config, "-d", "--device",
81 "Device configuration string, see Device::Configure().");
82 args.Parse();
83 if (!args.Good())
84 {
85 if (myid == 0)
86 {
87 args.PrintUsage(cout);
88 }
89 return 1;
90 }
91 if (myid == 0)
92 {
93 args.PrintOptions(cout);
94 }
95
96 // 3. Enable hardware devices such as GPUs, and programming models such as
97 // CUDA, OCCA, RAJA and OpenMP based on command line options.
98 Device device(device_config);
99 if (myid == 0) { device.Print(); }
100
101 // 4. Read the (serial) mesh from the given mesh file on all processors. We
102 // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
103 // and volume meshes with the same code.
104 Mesh *mesh = new Mesh(mesh_file, 1, 1);
105 int dim = mesh->Dimension();
106 if (nc)
107 {
108 mesh->EnsureNCMesh(true);
109 }
110
111 // 5. Refine the serial mesh on all processors to increase the resolution. In
112 // this example we do 'ref_levels' of uniform refinement (2 by default, or
113 // specified on the command line with -rs).
114 for (int lev = 0; lev < ser_ref_levels; lev++)
115 {
116 mesh->UniformRefinement();
117 }
118
119 // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
120 // this mesh further in parallel to increase the resolution (1 time by
121 // default, or specified on the command line with -rp). Once the parallel
122 // mesh is defined, the serial mesh can be deleted.
123 ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
124 delete mesh;
125 for (int lev = 0; lev < par_ref_levels; lev++)
126 {
127 pmesh->UniformRefinement();
128 }
129
130 // 7. Define a parallel finite element space on the parallel mesh. Here we
131 // use the Nedelec finite elements of the specified order.
133 ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
134 HYPRE_BigInt size = fespace->GlobalTrueVSize();
135 if (myid == 0)
136 {
137 cout << "Number of unknowns: " << size << endl;
138 }
139
140 // 8. Set up the parallel bilinear forms a(.,.) and m(.,.) on the finite
141 // element space. The first corresponds to the curl curl, while the second
142 // is a simple mass matrix needed on the right hand side of the
143 // generalized eigenvalue problem below. The boundary conditions are
144 // implemented by marking all the boundary attributes from the mesh as
145 // essential. The corresponding degrees of freedom are eliminated with
146 // special values on the diagonal to shift the Dirichlet eigenvalues out
147 // of the computational range. After serial and parallel assembly we
148 // extract the corresponding parallel matrices A and M.
149 ConstantCoefficient one(1.0);
150 Array<int> ess_bdr;
151 if (pmesh->bdr_attributes.Size())
152 {
153 ess_bdr.SetSize(pmesh->bdr_attributes.Max());
154 ess_bdr = 1;
155 }
156
157 ParBilinearForm *a = new ParBilinearForm(fespace);
158 a->AddDomainIntegrator(new CurlCurlIntegrator(one));
159 if (pmesh->bdr_attributes.Size() == 0)
160 {
161 // Add a mass term if the mesh has no boundary, e.g. periodic mesh or
162 // closed surface.
163 a->AddDomainIntegrator(new VectorFEMassIntegrator(one));
164 }
165 a->Assemble();
166 a->EliminateEssentialBCDiag(ess_bdr, 1.0);
167 a->Finalize();
168
169 ParBilinearForm *m = new ParBilinearForm(fespace);
171 m->Assemble();
172 // shift the eigenvalue corresponding to eliminated dofs to a large value
173 m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<real_t>::min());
174 m->Finalize();
175
176 HypreParMatrix *A = a->ParallelAssemble();
178
179 delete a;
180 delete m;
181
182 // 9. Define and configure the AME eigensolver and the AMS preconditioner for
183 // A to be used within the solver. Set the matrices which define the
184 // generalized eigenproblem A x = lambda M x.
185 HypreAMS *ams = new HypreAMS(*A,fespace);
186 ams->SetPrintLevel(0);
187 ams->SetSingularProblem();
188
189 HypreAME *ame = new HypreAME(MPI_COMM_WORLD);
190 ame->SetNumModes(nev);
191 ame->SetPreconditioner(*ams);
192 ame->SetMaxIter(100);
193 ame->SetTol(1e-8);
194 ame->SetPrintLevel(1);
195 ame->SetMassMatrix(*M);
196 ame->SetOperator(*A);
197
198 // 10. Compute the eigenmodes and extract the array of eigenvalues. Define a
199 // parallel grid function to represent each of the eigenmodes returned by
200 // the solver.
201 Array<real_t> eigenvalues;
202 ame->Solve();
203 ame->GetEigenvalues(eigenvalues);
204 ParGridFunction x(fespace);
205
206 // 11. Save the refined mesh and the modes in parallel. This output can be
207 // viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
208 {
209 ostringstream mesh_name, mode_name;
210 mesh_name << "mesh." << setfill('0') << setw(6) << myid;
211
212 ofstream mesh_ofs(mesh_name.str().c_str());
213 mesh_ofs.precision(8);
214 pmesh->Print(mesh_ofs);
215
216 for (int i=0; i<nev; i++)
217 {
218 // convert eigenvector from HypreParVector to ParGridFunction
219 x = ame->GetEigenvector(i);
220
221 mode_name << "mode_" << setfill('0') << setw(2) << i << "."
222 << setfill('0') << setw(6) << myid;
223
224 ofstream mode_ofs(mode_name.str().c_str());
225 mode_ofs.precision(8);
226 x.Save(mode_ofs);
227 mode_name.str("");
228 }
229 }
230
231 // 12. Send the solution by socket to a GLVis server.
232 if (visualization)
233 {
234 char vishost[] = "localhost";
235 int visport = 19916;
236 socketstream mode_sock(vishost, visport);
237 mode_sock.precision(8);
238
239 for (int i=0; i<nev; i++)
240 {
241 if ( myid == 0 )
242 {
243 cout << "Eigenmode " << i+1 << '/' << nev
244 << ", Lambda = " << eigenvalues[i] << endl;
245 }
246
247 // convert eigenvector from HypreParVector to ParGridFunction
248 x = ame->GetEigenvector(i);
249
250 mode_sock << "parallel " << num_procs << " " << myid << "\n"
251 << "solution\n" << *pmesh << x << flush
252 << "window_title 'Eigenmode " << i+1 << '/' << nev
253 << ", Lambda = " << eigenvalues[i] << "'" << endl;
254
255 char c;
256 if (myid == 0)
257 {
258 cout << "press (q)uit or (c)ontinue --> " << flush;
259 cin >> c;
260 }
261 MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
262
263 if (c != 'c')
264 {
265 break;
266 }
267 }
268 mode_sock.close();
269 }
270
271 // 13. Free the used memory.
272 delete ame;
273 delete ams;
274 delete M;
275 delete A;
276
277 delete fespace;
278 delete fec;
279 delete pmesh;
280
281 return 0;
282}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Size() const
Return the logical size of the array.
Definition array.hpp:144
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY....
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void EliminateEssentialBCDiag(const Array< int > &bdr_attr_is_ess, real_t value)
Perform elimination and set the diagonal entry to the given value.
A coefficient that is constant across space and time.
Integrator for for Nedelec elements.
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
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
void SetTol(real_t tol)
Definition hypre.cpp:6570
void SetNumModes(int num_eigs)
Definition hypre.cpp:6562
void SetPreconditioner(HypreSolver &precond)
Definition hypre.cpp:6601
void SetPrintLevel(int logging)
Definition hypre.cpp:6592
void SetMassMatrix(const HypreParMatrix &M)
Definition hypre.cpp:6622
const HypreParVector & GetEigenvector(unsigned int i) const
Extract a single eigenvector.
Definition hypre.cpp:6665
void GetEigenvalues(Array< real_t > &eigenvalues) const
Collect the converged eigenvalues.
Definition hypre.cpp:6641
void SetOperator(const HypreParMatrix &A)
Definition hypre.cpp:6607
void Solve()
Solve the eigenproblem.
Definition hypre.cpp:6629
void SetMaxIter(int max_iter)
Definition hypre.cpp:6586
The Auxiliary-space Maxwell Solver in hypre.
Definition hypre.hpp:1845
void SetPrintLevel(int print_lvl)
Definition hypre.cpp:5805
void SetSingularProblem()
Set this option when solving a curl-curl problem with zero mass term.
Definition hypre.hpp:1903
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
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 Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition mesh.cpp:10626
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).
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition fe_coll.hpp:465
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.
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
void Assemble(int skip_zeros=1)
Assemble the local matrix.
Abstract parallel finite element space.
Definition pfespace.hpp:29
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
Class for parallel grid function.
Definition pgridfunc.hpp:33
void Save(std::ostream &out) const override
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
int close()
Close the socketstream.
int dim
Definition ex24.cpp:53
int main()
HYPRE_Int HYPRE_BigInt
real_t a
Definition lissajous.cpp:41
const int visport
const char vishost[]