MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
ex35p.cpp
Go to the documentation of this file.
1// MFEM Example 35 - Parallel Version
2//
3// Compile with: make ex35p
4//
5// Sample runs: mpirun -np 4 ex35p -p 0 -o 2
6// mpirun -np 4 ex35p -p 0 -o 2 -pbc '22 23 24' -em 0
7// mpirun -np 4 ex35p -p 1 -o 1 -rp 2
8// mpirun -np 4 ex35p -p 1 -o 2
9// mpirun -np 4 ex35p -p 2 -o 1 -rp 2 -c 15
10//
11// Device sample runs:
12//
13// Description: This example code demonstrates the use of MFEM to define and
14// solve simple complex-valued linear systems. It implements three
15// variants of a damped harmonic oscillator:
16//
17// 1) A scalar H1 field
18// -Div(a Grad u) - omega^2 b u + i omega c u = 0
19//
20// 2) A vector H(Curl) field
21// Curl(a Curl u) - omega^2 b u + i omega c u = 0
22//
23// 3) A vector H(Div) field
24// -Grad(a Div u) - omega^2 b u + i omega c u = 0
25//
26// In each case the field is driven by a forced oscillation, with
27// angular frequency omega, imposed at the boundary or a portion
28// of the boundary. The spatial variation of the boundary
29// condition is computed as an eigenmode of an appropriate
30// operator defined on a portion of the boundary i.e. a port
31// boundary condition.
32//
33// In electromagnetics the coefficients are typically named the
34// permeability, mu = 1/a, permittivity, epsilon = b, and
35// conductivity, sigma = c. The user can specify these constants
36// using either set of names.
37//
38// This example demonstrates how to transfer fields computed on a
39// boundary generated SubMesh to the full mesh and apply them as
40// boundary conditions. The default mesh and corresponding
41// boundary attributes were chosen to verify proper behavior on
42// both triangular and quadrilateral faces of tetrahedral,
43// wedge-shaped, and hexahedral elements.
44//
45// The example also demonstrates how to display a time-varying
46// solution as a sequence of fields sent to a single GLVis socket.
47//
48// We recommend viewing examples 11, 13, and 22 before viewing
49// this example.
50
51#include "mfem.hpp"
52#include <fstream>
53#include <iostream>
54
55using namespace std;
56using namespace mfem;
57
58static real_t mu_ = 1.0;
59static real_t epsilon_ = 1.0;
60static real_t sigma_ = 2.0;
61
62void SetPortBC(int prob, int dim, int mode, ParGridFunction &port_bc);
63
64int main(int argc, char *argv[])
65{
66 // 1. Initialize MPI and HYPRE.
67 Mpi::Init(argc, argv);
68 int num_procs = Mpi::WorldSize();
69 int myid = Mpi::WorldRank();
71
72 // 2. Parse command-line options.
73 const char *mesh_file = "../data/fichera-mixed.mesh";
74 int ser_ref_levels = 1;
75 int par_ref_levels = 1;
76 int order = 1;
77 Array<int> port_bc_attr;
78 int prob = 0;
79 int mode = 1;
80 real_t freq = -1.0;
81 real_t omega = 2.0 * M_PI;
82 real_t a_coef = 0.0;
83 bool herm_conv = true;
84 bool slu_solver = false;
85 bool visualization = 1;
86 bool mixed = true;
87 bool pa = false;
88 const char *device_config = "cpu";
89
90 OptionsParser args(argc, argv);
92 "Mesh file to use.");
94 "Number of times to refine the mesh uniformly in serial.");
96 "Number of times to refine the mesh uniformly in parallel.");
98 "Finite element order (polynomial degree).");
100 "Choose between 0: H_1, 1: H(Curl), or 2: H(Div) "
101 "damped harmonic oscillator.");
103 "Choose the index of the port eigenmode.");
105 "Stiffness coefficient (spring constant or 1/mu).");
107 "Mass coefficient (or epsilon).");
109 "Damping coefficient (or sigma).");
111 "Permeability of free space (or 1/(spring constant)).");
113 "Permittivity of free space (or mass constant).");
115 "Conductivity (or damping constant).");
117 "Frequency (in Hz).");
119 "Attributes of port boundary condition");
121 "--no-hermitian", "Use convention for Hermitian operators.");
122#ifdef MFEM_USE_SUPERLU
124 "--no-superlu", "Use the SuperLU Solver.");
125#endif
127 "--no-visualization",
128 "Enable or disable GLVis visualization.");
130 "--hex-mesh", "Mixed mesh of hexahedral mesh.");
132 "--no-partial-assembly", "Enable Partial Assembly.");
134 "Device configuration string, see Device::Configure().");
135 args.Parse();
136 if (!args.Good())
137 {
138 if (myid == 0)
139 {
140 args.PrintUsage(cout);
141 }
142 return 1;
143 }
144
145 if (!mixed || pa)
146 {
147 mesh_file = "../data/fichera.mesh";
148 }
149
150 if ( a_coef != 0.0 )
151 {
152 mu_ = 1.0 / a_coef;
153 }
154 if ( freq > 0.0 )
155 {
156 omega = 2.0 * M_PI * freq;
157 }
158 if (port_bc_attr.Size() == 0 &&
159 (strcmp(mesh_file, "../data/fichera-mixed.mesh") == 0 ||
160 strcmp(mesh_file, "../data/fichera.mesh") == 0))
161 {
162 port_bc_attr.SetSize(4);
163 port_bc_attr[0] = 7;
164 port_bc_attr[1] = 8;
165 port_bc_attr[2] = 11;
166 port_bc_attr[3] = 12;
167 }
168
169 if (myid == 0)
170 {
171 args.PrintOptions(cout);
172 }
173
174 MFEM_VERIFY(prob >= 0 && prob <=2,
175 "Unrecognized problem type: " << prob);
176
179
180 // 3. Enable hardware devices such as GPUs, and programming models such as
181 // CUDA, OCCA, RAJA and OpenMP based on command line options.
182 Device device(device_config);
183 if (myid == 0) { device.Print(); }
184
185 // 4. Read the (serial) mesh from the given mesh file on all processors. We
186 // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
187 // and volume meshes with the same code.
188 Mesh *mesh = new Mesh(mesh_file, 1, 1);
189 int dim = mesh->Dimension();
190
191 // 5. Refine the serial mesh on all processors to increase the resolution.
192 for (int l = 0; l < ser_ref_levels; l++)
193 {
194 mesh->UniformRefinement();
195 }
196
197 // 6a. Define a parallel mesh by a partitioning of the serial mesh. Refine
198 // this mesh further in parallel to increase the resolution. Once the
199 // parallel mesh is defined, the serial mesh can be deleted.
200 ParMesh pmesh(MPI_COMM_WORLD, *mesh);
201 delete mesh;
202 for (int l = 0; l < par_ref_levels; l++)
203 {
204 pmesh.UniformRefinement();
205 }
206
207 // 6b. Extract a submesh covering a portion of the boundary
208 ParSubMesh pmesh_port(ParSubMesh::CreateFromBoundary(pmesh, port_bc_attr));
209
210 // 7a. Define a parallel finite element space on the parallel mesh. Here we
211 // use continuous Lagrange, Nedelec, or Raviart-Thomas finite elements
212 // of the specified order.
213 if (dim == 1 && prob != 0 )
214 {
215 if (myid == 0)
216 {
217 cout << "Switching to problem type 0, H1 basis functions, "
218 << "for 1 dimensional mesh." << endl;
219 }
220 prob = 0;
221 }
222
223 FiniteElementCollection *fec = NULL;
224 switch (prob)
225 {
226 case 0: fec = new H1_FECollection(order, dim); break;
227 case 1: fec = new ND_FECollection(order, dim); break;
228 case 2: fec = new RT_FECollection(order - 1, dim); break;
229 default: break; // This should be unreachable
230 }
231 ParFiniteElementSpace fespace(&pmesh, fec);
232 HYPRE_BigInt size = fespace.GlobalTrueVSize();
233 if (myid == 0)
234 {
235 cout << "Number of finite element unknowns: " << size << endl;
236 }
237
238 // 7b. Define a parallel finite element space on the sub-mesh. Here we
239 // use continuous Lagrange, Nedelec, or L2 finite elements of
240 // the specified order.
241 FiniteElementCollection *fec_port = NULL;
242 switch (prob)
243 {
244 case 0: fec_port = new H1_FECollection(order, dim-1); break;
245 case 1:
246 if (dim == 3)
247 {
248 fec_port = new ND_FECollection(order, dim-1);
249 }
250 else
251 {
252 fec_port = new L2_FECollection(order - 1, dim-1,
255 }
256 break;
257 case 2: fec_port = new L2_FECollection(order - 1, dim-1,
260 default: break; // This should be unreachable
261 }
262 ParFiniteElementSpace fespace_port(&pmesh_port, fec_port);
263 HYPRE_BigInt size_port = fespace_port.GlobalTrueVSize();
264 if (myid == 0)
265 {
266 cout << "Number of finite element port BC unknowns: " << size_port
267 << endl;
268 }
269
270 // 8a. Define a parallel grid function on the SubMesh which will contain
271 // the field to be applied as a port boundary condition.
272 ParGridFunction port_bc(&fespace_port);
273 port_bc = 0.0;
274
275 SetPortBC(prob, dim, mode, port_bc);
276
277 // 8b. Save the SubMesh and associated port boundary condition in parallel.
278 // This output can be viewed later using GLVis:
279 // "glvis -np <np> -m port_mesh -g port_mode"
280 {
281 ostringstream mesh_name, port_name;
282 mesh_name << "port_mesh." << setfill('0') << setw(6) << myid;
283 port_name << "port_mode." << setfill('0') << setw(6) << myid;
284
285 ofstream mesh_ofs(mesh_name.str().c_str());
286 mesh_ofs.precision(8);
287 pmesh_port.Print(mesh_ofs);
288
289 ofstream port_ofs(port_name.str().c_str());
290 port_ofs.precision(8);
291 port_bc.Save(port_ofs);
292 }
293 // 8c. Send the port bc, computed on the SubMesh, to a GLVis server.
294 if (visualization && dim == 3)
295 {
296 char vishost[] = "localhost";
297 int visport = 19916;
298 socketstream port_sock(vishost, visport);
299 port_sock << "parallel " << num_procs << " " << myid << "\n";
300 port_sock.precision(8);
301 port_sock << "solution\n" << pmesh_port << port_bc
302 << "window_title 'Port BC'"
303 << "window_geometry 0 0 400 350" << flush;
304 }
305
306 // 9. Determine the list of true (i.e. parallel conforming) essential
307 // boundary dofs. In this example, the boundary conditions are defined
308 // using an eigenmode of the appropriate type computed on the SubMesh.
309 Array<int> ess_tdof_list;
310 Array<int> ess_bdr;
311 if (pmesh.bdr_attributes.Size())
312 {
313 ess_bdr.SetSize(pmesh.bdr_attributes.Max());
314 ess_bdr = 1;
315 fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
316 }
317
318 // 10. Set up the parallel linear form b(.) which corresponds to the
319 // right-hand side of the FEM linear system.
320 ParComplexLinearForm b(&fespace, conv);
321 b.Vector::operator=(0.0);
322
323 // 11a. Define the solution vector u as a parallel complex finite element
324 // grid function corresponding to fespace. Initialize u to equal zero.
325 ParComplexGridFunction u(&fespace);
326 u = 0.0;
327 pmesh_port.Transfer(port_bc, u.real());
328
329 // 11b. Send the transferred port bc field to a GLVis server.
330 {
331 ParGridFunction full_bc(&fespace);
332 ParTransferMap port_to_full(port_bc, full_bc);
333
334 full_bc = 0.0;
335 port_to_full.Transfer(port_bc, full_bc);
336
337 if (visualization)
338 {
339 char vishost[] = "localhost";
340 int visport = 19916;
341 socketstream full_sock(vishost, visport);
342 full_sock << "parallel " << num_procs << " " << myid << "\n";
343 full_sock.precision(8);
344 full_sock << "solution\n" << pmesh << full_bc
345 << "window_title 'Transferred BC'"
346 << "window_geometry 400 0 400 350"<< flush;
347 }
348 }
349
350 // 12. Set up the parallel sesquilinear form a(.,.) on the finite element
351 // space corresponding to the damped harmonic oscillator operator of the
352 // appropriate type:
353 //
354 // 0) A scalar H1 field
355 // -Div(a Grad) - omega^2 b + i omega c
356 //
357 // 1) A vector H(Curl) field
358 // Curl(a Curl) - omega^2 b + i omega c
359 //
360 // 2) A vector H(Div) field
361 // -Grad(a Div) - omega^2 b + i omega c
362 //
363 ConstantCoefficient stiffnessCoef(1.0/mu_);
364 ConstantCoefficient massCoef(-omega * omega * epsilon_);
365 ConstantCoefficient lossCoef(omega * sigma_);
366 ConstantCoefficient negMassCoef(omega * omega * epsilon_);
367
368 ParSesquilinearForm a(&fespace, conv);
369 if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
370 switch (prob)
371 {
372 case 0:
374 NULL);
376 new MassIntegrator(lossCoef));
377 break;
378 case 1:
380 NULL);
382 new VectorFEMassIntegrator(lossCoef));
383 break;
384 case 2:
386 NULL);
388 new VectorFEMassIntegrator(lossCoef));
389 break;
390 default: break; // This should be unreachable
391 }
392
393 // 13. Assemble the parallel bilinear form and the corresponding linear
394 // system, applying any necessary transformations such as: parallel
395 // assembly, eliminating boundary conditions, applying conforming
396 // constraints for non-conforming AMR, etc.
397 a.Assemble();
398
400 Vector B, U;
401
402 a.FormLinearSystem(ess_tdof_list, u, b, A, U, B);
403
404 if (myid == 0)
405 {
406 cout << "Size of linear system: "
407 << 2 * size << endl << endl;
408 }
409
410 if (!slu_solver)
411 {
412 // 14a. Set up the parallel bilinear form for the preconditioner
413 // corresponding to the appropriate operator
414 //
415 // 0) A scalar H1 field
416 // -Div(a Grad) - omega^2 b + i omega c
417 //
418 // 1) A vector H(Curl) field
419 // Curl(a Curl) + omega^2 b + i omega c
420 //
421 // 2) A vector H(Div) field
422 // -Grad(a Div) - omega^2 b + i omega c
423 ParBilinearForm pcOp(&fespace);
424 if (pa) { pcOp.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
425 switch (prob)
426 {
427 case 0:
431 break;
432 case 1:
436 break;
437 case 2:
441 break;
442 default: break; // This should be unreachable
443 }
444 pcOp.Assemble();
445
446 // 14b. Define and apply a parallel FGMRES solver for AU=B with a block
447 // diagonal preconditioner based on the appropriate multigrid
448 // preconditioner from hypre.
449 Array<int> blockTrueOffsets;
450 blockTrueOffsets.SetSize(3);
451 blockTrueOffsets[0] = 0;
452 blockTrueOffsets[1] = A->Height() / 2;
453 blockTrueOffsets[2] = A->Height() / 2;
454 blockTrueOffsets.PartialSum();
455
456 BlockDiagonalPreconditioner BDP(blockTrueOffsets);
457
458 Operator * pc_r = NULL;
459 Operator * pc_i = NULL;
460
461 if (pa)
462 {
463 pc_r = new OperatorJacobiSmoother(pcOp, ess_tdof_list);
464 }
465 else
466 {
467 OperatorHandle PCOp;
468 pcOp.FormSystemMatrix(ess_tdof_list, PCOp);
469
470 switch (prob)
471 {
472 case 0:
473 pc_r = new HypreBoomerAMG(*PCOp.As<HypreParMatrix>());
474 break;
475 case 1:
476 pc_r = new HypreAMS(*PCOp.As<HypreParMatrix>(), &fespace);
477 break;
478 case 2:
479 if (dim == 2 )
480 {
481 pc_r = new HypreAMS(*PCOp.As<HypreParMatrix>(), &fespace);
482 }
483 else
484 {
485 pc_r = new HypreADS(*PCOp.As<HypreParMatrix>(), &fespace);
486 }
487 break;
488 default: break; // This should be unreachable
489 }
490 }
491 pc_i = new ScaledOperator(pc_r,
493 -1.0:1.0);
494
495 BDP.SetDiagonalBlock(0, pc_r);
496 BDP.SetDiagonalBlock(1, pc_i);
497 BDP.owns_blocks = 1;
498
499 FGMRESSolver fgmres(MPI_COMM_WORLD);
500 fgmres.SetPreconditioner(BDP);
501 fgmres.SetOperator(*A.Ptr());
502 fgmres.SetRelTol(1e-6);
503 fgmres.SetMaxIter(1000);
504 fgmres.SetPrintLevel(1);
505 fgmres.Mult(B, U);
506 }
507#ifdef MFEM_USE_SUPERLU
508 else
509 {
510 // 14. Solve using a direct solver
511 // Transform to monolithic HypreParMatrix
512 HypreParMatrix *A_hyp = A.As<ComplexHypreParMatrix>()->GetSystemMatrix();
513 SuperLURowLocMatrix SA(*A_hyp);
514 SuperLUSolver superlu(MPI_COMM_WORLD);
515 superlu.SetPrintStatistics(true);
516 superlu.SetSymmetricPattern(false);
518 superlu.SetOperator(SA);
519 superlu.Mult(B, U);
520 delete A_hyp;
521 }
522#endif
523
524 // 15. Recover the parallel grid function corresponding to U. This is the
525 // local finite element solution on each processor.
526 a.RecoverFEMSolution(U, b, u);
527
528 // 16. Save the refined mesh and the solution in parallel. This output can be
529 // viewed later using GLVis: "glvis -np <np> -m mesh -g sol_r" or
530 // "glvis -np <np> -m mesh -g sol_i".
531 {
532 ostringstream mesh_name, sol_r_name, sol_i_name;
533 mesh_name << "mesh." << setfill('0') << setw(6) << myid;
534 sol_r_name << "sol_r." << setfill('0') << setw(6) << myid;
535 sol_i_name << "sol_i." << setfill('0') << setw(6) << myid;
536
537 ofstream mesh_ofs(mesh_name.str().c_str());
538 mesh_ofs.precision(8);
539 pmesh.Print(mesh_ofs);
540
541 ofstream sol_r_ofs(sol_r_name.str().c_str());
542 ofstream sol_i_ofs(sol_i_name.str().c_str());
543 sol_r_ofs.precision(8);
544 sol_i_ofs.precision(8);
545 u.real().Save(sol_r_ofs);
546 u.imag().Save(sol_i_ofs);
547 }
548
549 // 17. Send the solution by socket to a GLVis server.
550 if (visualization)
551 {
552 char vishost[] = "localhost";
553 int visport = 19916;
554 socketstream sol_sock_r(vishost, visport);
555 sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
556 sol_sock_r.precision(8);
557 sol_sock_r << "solution\n" << pmesh << u.real()
558 << "window_title 'Solution: Real Part'"
559 << "window_geometry 800 0 400 350" << flush;
560
561 MPI_Barrier(MPI_COMM_WORLD);
562
563 socketstream sol_sock_i(vishost, visport);
564 sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
565 sol_sock_i.precision(8);
566 sol_sock_i << "solution\n" << pmesh << u.imag()
567 << "window_title 'Solution: Imaginary Part'"
568 << "window_geometry 1200 0 400 350" << flush;
569 }
570 if (visualization)
571 {
572 ParGridFunction u_t(&fespace);
573 u_t = u.real();
574 char vishost[] = "localhost";
575 int visport = 19916;
576 socketstream sol_sock(vishost, visport);
577 sol_sock << "parallel " << num_procs << " " << myid << "\n";
578 sol_sock.precision(8);
579 sol_sock << "solution\n" << pmesh << u_t
580 << "window_title 'Harmonic Solution (t = 0.0 T)'"
581 << "window_geometry 0 432 600 450"
582 << "pause\n" << flush;
583 if (myid == 0)
584 cout << "GLVis visualization paused."
585 << " Press space (in the GLVis window) to resume it.\n";
586 int num_frames = 32;
587 int i = 0;
588 while (sol_sock)
589 {
590 real_t t = (real_t)(i % num_frames) / num_frames;
591 ostringstream oss;
592 oss << "Harmonic Solution (t = " << t << " T)";
593
594 add(cos( 2.0 * M_PI * t), u.real(),
595 sin(-2.0 * M_PI * t), u.imag(), u_t);
596 sol_sock << "parallel " << num_procs << " " << myid << "\n";
597 sol_sock << "solution\n" << pmesh << u_t
598 << "window_title '" << oss.str() << "'" << flush;
599 i++;
600 }
601 }
602
603 // 18. Free the used memory.
604 delete fec_port;
605 delete fec;
606
607 return 0;
608}
609
610/**
611 Solves the eigenvalue problem -Div(Grad x) = lambda x with homogeneous
612 Dirichlet boundary conditions on the boundary of the domain. Returns mode
613 number "mode" (counting from zero) in the ParGridFunction "x".
614*/
616{
617 int nev = std::max(mode + 2, 5);
618 int seed = 75;
619
620 ParFiniteElementSpace &fespace = *x.ParFESpace();
621 ParMesh &pmesh = *fespace.GetParMesh();
622
623 Array<int> ess_bdr;
624 if (pmesh.bdr_attributes.Size())
625 {
626 ess_bdr.SetSize(pmesh.bdr_attributes.Max());
627 ess_bdr = 1;
628 }
629
630 ParBilinearForm a(&fespace);
632 a.Assemble();
633 a.EliminateEssentialBCDiag(ess_bdr, 1.0);
634 a.Finalize();
635
636 ParBilinearForm m(&fespace);
638 m.Assemble();
639 // shift the eigenvalue corresponding to eliminated dofs to a large value
640 m.EliminateEssentialBCDiag(ess_bdr, numeric_limits<real_t>::min());
641 m.Finalize();
642
643 HypreParMatrix *A = a.ParallelAssemble();
645
646 HypreBoomerAMG amg(*A);
647 amg.SetPrintLevel(0);
648
649 HypreLOBPCG lobpcg(MPI_COMM_WORLD);
650 lobpcg.SetNumModes(nev);
651 lobpcg.SetRandomSeed(seed);
652 lobpcg.SetPreconditioner(amg);
653 lobpcg.SetMaxIter(200);
654 lobpcg.SetTol(1e-8);
655 lobpcg.SetPrecondUsageMode(1);
656 lobpcg.SetPrintLevel(1);
657 lobpcg.SetMassMatrix(*M);
658 lobpcg.SetOperator(*A);
659 lobpcg.Solve();
660
661 x = lobpcg.GetEigenvector(mode);
662
663 delete A;
664 delete M;
665}
666
667/**
668 Solves the eigenvalue problem -Curl(Curl x) = lambda x with homogeneous
669 Dirichlet boundary conditions, on the tangential component of x, on the
670 boundary of the domain. Returns mode number "mode" (counting from zero) in
671 the ParGridFunction "x".
672*/
674{
675 int nev = std::max(mode + 2, 5);
676
677 ParFiniteElementSpace &fespace = *x.ParFESpace();
678 ParMesh &pmesh = *fespace.GetParMesh();
679
680 Array<int> ess_bdr;
681 if (pmesh.bdr_attributes.Size())
682 {
683 ess_bdr.SetSize(pmesh.bdr_attributes.Max());
684 ess_bdr = 1;
685 }
686
687 ParBilinearForm a(&fespace);
689 a.Assemble();
690 a.EliminateEssentialBCDiag(ess_bdr, 1.0);
691 a.Finalize();
692
693 ParBilinearForm m(&fespace);
695 m.Assemble();
696 // shift the eigenvalue corresponding to eliminated dofs to a large value
697 m.EliminateEssentialBCDiag(ess_bdr, numeric_limits<real_t>::min());
698 m.Finalize();
699
700 HypreParMatrix *A = a.ParallelAssemble();
702
703 HypreAMS ams(*A,&fespace);
704 ams.SetPrintLevel(0);
705 ams.SetSingularProblem();
706
707 HypreAME ame(MPI_COMM_WORLD);
708 ame.SetNumModes(nev);
709 ame.SetPreconditioner(ams);
710 ame.SetMaxIter(100);
711 ame.SetTol(1e-8);
712 ame.SetPrintLevel(1);
713 ame.SetMassMatrix(*M);
714 ame.SetOperator(*A);
715 ame.Solve();
716
717 x = ame.GetEigenvector(mode);
718
719 delete A;
720 delete M;
721}
722
723/**
724 Solves the eigenvalue problem -Div(Grad x) = lambda x with homogeneous
725 Neumann boundary conditions on the boundary of the domain. Returns mode
726 number "mode" (counting from zero) in the ParGridFunction "x_l2". Note that
727 mode 0 is a constant field so higher mode numbers are often more
728 interesting. The eigenmode is solved using continuous H1 basis of the
729 appropriate order and then projected onto the L2 basis and returned.
730*/
732{
733 int nev = std::max(mode + 2, 5);
734 int seed = 75;
735
736 ParFiniteElementSpace &fespace_l2 = *x_l2.ParFESpace();
737 ParMesh &pmesh = *fespace_l2.GetParMesh();
738 int order_l2 = fespace_l2.FEColl()->GetOrder();
739
740 H1_FECollection fec(order_l2+1, pmesh.Dimension());
741 ParFiniteElementSpace fespace(&pmesh, &fec);
742 ParGridFunction x(&fespace);
743 x = 0.0;
744
745 GridFunctionCoefficient xCoef(&x);
746
747 if (mode == 0)
748 {
749 x = 1.0;
750 x_l2.ProjectCoefficient(xCoef);
751 return;
752 }
753
754 ParBilinearForm a(&fespace);
756 a.AddDomainIntegrator(new MassIntegrator); // Shift eigenvalues by 1
757 a.Assemble();
758 a.Finalize();
759
760 ParBilinearForm m(&fespace);
762 m.Assemble();
763 m.Finalize();
764
765 HypreParMatrix *A = a.ParallelAssemble();
767
768 HypreBoomerAMG amg(*A);
769 amg.SetPrintLevel(0);
770
771 HypreLOBPCG lobpcg(MPI_COMM_WORLD);
772 lobpcg.SetNumModes(nev);
773 lobpcg.SetRandomSeed(seed);
774 lobpcg.SetPreconditioner(amg);
775 lobpcg.SetMaxIter(200);
776 lobpcg.SetTol(1e-8);
777 lobpcg.SetPrecondUsageMode(1);
778 lobpcg.SetPrintLevel(1);
779 lobpcg.SetMassMatrix(*M);
780 lobpcg.SetOperator(*A);
781 lobpcg.Solve();
782
783 x = lobpcg.GetEigenvector(mode);
784
785 x_l2.ProjectCoefficient(xCoef);
786
787 delete A;
788 delete M;
789}
790
791// Compute eigenmode "mode" of either a Dirichlet or Neumann Laplacian or of a
792// Dirichlet curl curl operator based on the problem type and dimension of the
793// domain.
794void SetPortBC(int prob, int dim, int mode, ParGridFunction &port_bc)
795{
796 switch (prob)
797 {
798 case 0:
799 ScalarWaveGuide(mode, port_bc);
800 break;
801 case 1:
802 if (dim == 3)
803 {
804 VectorWaveGuide(mode, port_bc);
805 }
806 else
807 {
808 PseudoScalarWaveGuide(mode, port_bc);
809 }
810 break;
811 case 2:
812 PseudoScalarWaveGuide(mode, port_bc);
813 break;
814 }
815}
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
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
Definition array.cpp:103
@ GaussLegendre
Open type.
Definition fe_base.hpp:31
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY....
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 class to handle Block diagonal preconditioners in a matrix-free implementation.
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
Specialization of the ComplexOperator built from a pair of HypreParMatrices.
@ HERMITIAN
Native convention for Hermitian operators.
@ BLOCK_SYMMETRIC
Alternate convention for damping operators.
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
for Raviart-Thomas elements
FGMRES method.
Definition solvers.hpp:567
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the FGMRES method.
Definition solvers.cpp:1168
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
int GetOrder() const
Return the order (polynomial degree) of the FE collection, corresponding to the order/degree returned...
Definition fe_coll.hpp:225
const FiniteElementCollection * FEColl() const
Definition fespace.hpp:727
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
The Auxiliary-space Divergence Solver in hypre.
Definition hypre.hpp:1922
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 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
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
void SetPrintLevel(int print_level)
Definition hypre.hpp:1771
void SetMassMatrix(Operator &M)
Definition hypre.cpp:6343
void SetPrintLevel(int logging)
Definition hypre.cpp:6272
void SetPreconditioner(Solver &precond)
Definition hypre.cpp:6287
void SetTol(real_t tol)
Definition hypre.cpp:6250
void SetOperator(Operator &A)
Definition hypre.cpp:6296
void SetNumModes(int num_eigs)
Definition hypre.hpp:2103
void Solve()
Solve the eigenproblem.
Definition hypre.cpp:6410
void SetPrecondUsageMode(int pcg_mode)
Definition hypre.cpp:6281
void SetRandomSeed(int s)
Definition hypre.hpp:2105
void SetMaxIter(int max_iter)
Definition hypre.cpp:6266
const HypreParVector & GetEigenvector(unsigned int i) const
Extract a single eigenvector.
Definition hypre.cpp:6365
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
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.cpp:179
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
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
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 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
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
Operator * Ptr() const
Access the underlying Operator pointer.
Definition handle.hpp:87
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition solvers.hpp:313
Abstract operator.
Definition operator.hpp:25
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
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.
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.
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
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
ParMesh * GetParMesh() const
Definition pfespace.hpp:277
Class for parallel grid function.
Definition pgridfunc.hpp:33
void Save(std::ostream &out) const override
ParFiniteElementSpace * ParFESpace() const
void ProjectCoefficient(Coefficient &coeff) override
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
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
Subdomain representation of a topological parent in another ParMesh.
Definition psubmesh.hpp:52
static ParSubMesh CreateFromBoundary(const ParMesh &parent, Array< int > &boundary_attributes)
Create a surface ParSubMesh from it's parent.
Definition psubmesh.cpp:32
static void Transfer(const ParGridFunction &src, ParGridFunction &dst)
Transfer the dofs of a ParGridFunction.
ParTransferMap represents a mapping of degrees of freedom from a source ParGridFunction to a destinat...
void Transfer(const ParGridFunction &src, ParGridFunction &dst) const
Transfer the source ParGridFunction to the destination ParGridFunction.
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition fe_coll.hpp:386
Scaled Operator B: x -> a A(x).
Definition operator.hpp:727
void SetColumnPermutation(superlu::ColPerm col_perm)
Specify how to permute the columns of the matrix.
Definition superlu.cpp:399
void Mult(const Vector &x, Vector &y) const
Factor and solve the linear system .
Definition superlu.cpp:587
void SetSymmetricPattern(bool sym)
Specify whether the matrix has a symmetric pattern to avoid extra work (default false)
Definition superlu.cpp:454
void SetOperator(const Operator &op)
Set the operator/matrix.
Definition superlu.cpp:475
void SetPrintStatistics(bool print_stat)
Specify whether to print the solver statistics (default true)
Definition superlu.cpp:385
Vector data type.
Definition vector.hpp:80
int dim
Definition ex24.cpp:53
real_t freq
Definition ex24.cpp:54
real_t omega
Definition ex25.cpp:142
prob_type prob
Definition ex25.cpp:156
void SetPortBC(int prob, int dim, int mode, ParGridFunction &port_bc)
Definition ex35p.cpp:794
void ScalarWaveGuide(int mode, ParGridFunction &x)
Definition ex35p.cpp:615
void VectorWaveGuide(int mode, ParGridFunction &x)
Definition ex35p.cpp:673
void PseudoScalarWaveGuide(int mode, ParGridFunction &x_l2)
Definition ex35p.cpp:731
int main()
HYPRE_Int HYPRE_BigInt
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
@ PARMETIS
Sequential ordering on structure of using the PARMETIS package.
Definition superlu.hpp:71
const int visport
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition vector.cpp:316
float real_t
Definition config.hpp:43
const char vishost[]
RefCoord t[3]