MFEM  v4.6.0 Finite element discretization library
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
55 using namespace std;
56 using namespace mfem;
57
58 static double mu_ = 1.0;
59 static double epsilon_ = 1.0;
60 static double sigma_ = 2.0;
61
62 void SetPortBC(int prob, int dim, int mode, ParGridFunction &port_bc);
63
64 int 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();
70  Hypre::Init();
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  double freq = -1.0;
81  double omega = 2.0 * M_PI;
82  double 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
178  herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
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,
253  BasisType::GaussLegendre,
254  FiniteElement::INTEGRAL);
255  }
256  break;
257  case 2: fec_port = new L2_FECollection(order - 1, dim-1,
258  BasisType::GaussLegendre,
259  FiniteElement::INTEGRAL); break;
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
399  OperatorHandle A;
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,
492  (conv == ComplexOperator::HERMITIAN) ?
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  double t = (double)(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<double>::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<double>::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.
794 void 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 }
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:1031
int visport
void ScalarWaveGuide(int mode, ParGridFunction &x)
Definition: ex35p.cpp:615
The Auxiliary-space Maxwell Solver in hypre.
Definition: hypre.hpp:1743
The Auxiliary-space Divergence Solver in hypre.
Definition: hypre.hpp:1820
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
void SetMassMatrix(const HypreParMatrix &M)
Definition: hypre.cpp:6471
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
void SetPortBC(int prob, int dim, int mode, ParGridFunction &port_bc)
Definition: ex35p.cpp:794
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
const HypreParVector & GetEigenvector(unsigned int i) const
Extract a single eigenvector.
Definition: hypre.cpp:6214
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
Integrator for (curl u, curl v) for Nedelec elements.
int GetOrder() const
Return the order (polynomial degree) of the FE collection, corresponding to the order/degree returned...
Definition: fe_coll.hpp:225
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1168
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
FGMRES method.
Definition: solvers.hpp:544
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
void SetPreconditioner(HypreSolver &precond)
Definition: hypre.cpp:6450
void SetTol(double tol)
Definition: hypre.cpp:6419
void SetPreconditioner(Solver &precond)
Definition: hypre.cpp:6136
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:279
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
void SetMassMatrix(Operator &M)
Definition: hypre.cpp:6192
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void SetPrintLevel(int logging)
Definition: hypre.cpp:6121
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:561
STL namespace.
(Q div u, div v) for RT elements
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:317
void SetSymmetricPattern(bool sym)
Definition: superlu.cpp:454
Subdomain representation of a topological parent in another ParMesh.
Definition: psubmesh.hpp:51
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
void SetTol(double tol)
Definition: hypre.cpp:6099
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:179
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: superlu.cpp:587
void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: superlu.cpp:475
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1590
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_iter)
Definition: hypre.cpp:6115
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
A class to handle Block diagonal preconditioners in a matrix-free implementation. ...
char vishost[]
void VectorWaveGuide(int mode, ParGridFunction &x)
Definition: ex35p.cpp:673
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:723
void SetColumnPermutation(superlu::ColPerm col_perm)
Definition: superlu.cpp:399
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:302
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1670
ParMesh * GetParMesh() const
Definition: pfespace.hpp:273
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
int main(int argc, char *argv[])
Definition: ex35p.cpp:64
void Assemble(int skip_zeros=1)
Assemble the local matrix.
ParTransferMap represents a mapping of degrees of freedom from a source ParGridFunction to a destinat...
void EliminateEssentialBCDiag(const Array< int > &bdr_attr_is_ess, double value)
Perform elimination and set the diagonal entry to the given value.
void SetPrintStatistics(bool print_stat)
Definition: superlu.cpp:385
prob_type prob
Definition: ex25.cpp:154
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:281
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:380
double omega
Definition: ex25.cpp:142
void Transfer(const ParGridFunction &src, ParGridFunction &dst) const
Transfer the source ParGridFunction to the destination ParGridFunction.
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
void SetRelTol(double rtol)
Definition: solvers.hpp:199
Scaled Operator B: x -> a A(x).
Definition: operator.hpp:726
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:5645
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:113
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 &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
void Solve()
Solve the eigenproblem.
Definition: hypre.cpp:6478
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
Definition: array.cpp:103
HYPRE_Int HYPRE_BigInt
void SetNumModes(int num_eigs)
Definition: hypre.cpp:6411
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
void SetMaxIter(int max_iter)
Definition: hypre.cpp:6435
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:909
void SetRandomSeed(int s)
Definition: hypre.hpp:2004
double a
Definition: lissajous.cpp:41
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:104
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
int dim
Definition: ex24.cpp:53
Specialization of the ComplexOperator built from a pair of HypreParMatrices.
double freq
Definition: ex24.cpp:54
void SetOperator(Operator &A)
Definition: hypre.cpp:6145
Adds new Domain Integrator. Assumes ownership of bfi.
void SetPrintLevel(int logging)
Definition: hypre.cpp:6441
Class for parallel bilinear form.
Operator * Ptr() const
Access the underlying Operator pointer.
Definition: handle.hpp:87
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
RefCoord t[3]
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:454
Vector data type.
Definition: vector.hpp:58
void SetPrecondUsageMode(int pcg_mode)
Definition: hypre.cpp:6130
void SetOperator(const HypreParMatrix &A)
Definition: hypre.cpp:6456
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
void Print(std::ostream &out=mfem::out) const override
Definition: pmesh.cpp:4825
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
Class for parallel grid function.
Definition: pgridfunc.hpp:32
static void Transfer(const ParGridFunction &src, ParGridFunction &dst)
Transfer the dofs of a ParGridFunction.
Definition: psubmesh.cpp:932
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
void SetNumModes(int num_eigs)
Definition: hypre.hpp:2002
Class for parallel meshes.
Definition: pmesh.hpp:32
void SetSingularProblem()
Set this option when solving a curl-curl problem with zero mass term.
Definition: hypre.hpp:1802
void Solve()
Solve the eigenproblem.
Definition: hypre.cpp:6259
void PseudoScalarWaveGuide(int mode, ParGridFunction &x_l2)
Definition: ex35p.cpp:731
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
const HypreParVector & GetEigenvector(unsigned int i) const
Extract a single eigenvector.
Definition: hypre.cpp:6514
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:327