MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
MultilevelHcurlHdivSolver.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 //
12 // ---------------------------------------------------
13 // ParELAG Miniapp: AMGe solver for H(curl) and H(div)
14 // ---------------------------------------------------
15 //
16 // This miniapp employs MFEM and ParELAG to solve H(curl)- and H(div)-elliptic
17 // forms by an element based algebraic multigrid (AMGe). It uses: a multilevel
18 // hierarchy of de Rham complexes of finite element spaces, built by ParELAG;
19 // Hiptmair-type smoothers, implemented in ParELAG; and AMS (Auxiliary-space
20 // Maxwell Solver) or ADS (Auxiliary-space Divergence Solver), from HYPRE, for
21 // preconditioning or solving on the coarsest levels. See the README file in
22 // this directory for more details.
23 //
24 // Compile with: see README
25 //
26 // Sample runs: MultilevelHcurlHdivSolver -curl -f \
27 // MultilevelHcurlSolver_pipe_example_parameters.xml
28 // MultilevelHcurlHdivSolver -div -f \
29 // MultilevelHdivSolver_pipe_example_parameters.xml
30 //
31 // We recommend viewing MFEM's examples 3 and 4 before viewing this miniapp.
32 
33 #include <fstream>
34 #include <sstream>
35 #include <ostream>
36 #include <string>
37 #include <vector>
38 #include <memory>
39 
40 #include <mpi.h>
41 
42 #include "elag.hpp"
43 #include "utilities/MPIDataTypes.hpp"
44 
45 using namespace mfem;
46 using namespace parelag;
47 using namespace std;
48 
49 void bdrfunc(const Vector &, Vector &);
50 void rhsfunc(const Vector &, Vector &);
51 
52 int main(int argc, char *argv[])
53 {
54  // Initialize MPI and HYPRE.
55  Mpi::Init();
56  int num_ranks = Mpi::WorldSize();
57  int myid = Mpi::WorldRank();
58  Hypre::Init();
59 
60  Timer total_timer = TimeManager::AddTimer("Program Execution -- Total");
61  Timer init_timer = TimeManager::AddTimer("Initial Setup");
62 
63  if (!myid)
64  cout << "-- This is an example of using a geometric-like multilevel "
65  "hierarchy, constructed by ParELAG,\n"
66  "-- to solve respective finite element H(curl) and H(div) forms: \n"
67  "(alpha curl u, curl v) + (beta u, v);\n"
68  "(alpha div u, div v) + (beta u, v).\n\n";
69 
70  // Get basic parameters from command line.
71  const char *xml_file_c = NULL;
72  bool hcurl = true;
73  bool visualize = false;
74  double tolSVD = 1e-3;
75  OptionsParser args(argc, argv);
76  args.AddOption(&xml_file_c, "-f", "--xml-file",
77  "XML parameter list (an XML file with detailed parameters).");
78  args.AddOption(&hcurl, "-curl", "--hcurl", "-div", "--hdiv",
79  "Whether the H(curl) or H(div) form is being solved.");
80  args.AddOption(&visualize, "-v", "--visualize", "-nv", "--no-visualize",
81  "Use GLVis to visualize the final solution and the "
82  "agglomerates.");
83  args.AddOption(&tolSVD, "-s", "--svd-tol",
84  "SVD tolerance. It is used for filtering out local linear "
85  "dependencies in the basis construction and extension "
86  "process in ParELAG. Namely, right singular vectors with "
87  "singular values smaller than this tolerance are removed.");
88  args.Parse();
89  if (!args.Good())
90  {
91  if (!myid)
92  {
93  args.PrintUsage(cout);
94  }
95  return EXIT_FAILURE;
96  }
97  if (!xml_file_c)
98  {
99  if (hcurl)
100  {
101  xml_file_c = "MultilevelHcurlSolver_cube_example_parameters.xml";
102  }
103  else
104  {
105  xml_file_c = "MultilevelHdivSolver_cube_example_parameters.xml";
106  }
107  if (!myid)
108  {
109  cout << "No XML parameter list provided! Using default "
110  << xml_file_c << "." << endl;
111  }
112  }
113  if (!myid)
114  {
115  args.PrintOptions(cout);
116  }
117  string xml_file(xml_file_c);
118 
119  // Read and parse the detailed parameter list from file.
120  unique_ptr<ParameterList> master_list;
121  ifstream xml_in(xml_file);
122  if (!xml_in.good())
123  {
124  if (!myid)
125  {
126  cerr << "ERROR: Cannot read from input file: " << xml_file << ".\n";
127  }
128  return EXIT_FAILURE;
129  }
130  SimpleXMLParameterListReader reader;
131  master_list = reader.GetParameterList(xml_in);
132  xml_in.close();
133 
134  // General parameters for the problem.
135  ParameterList& prob_list = master_list->Sublist("Problem parameters", true);
136 
137  // The file from which to read the mesh.
138  const string meshfile = prob_list.Get("Mesh file", "");
139 
140  // The number of times to refine in serial.
141  // Negative means refine until mesh is big enough to distribute, i.e.,
142  // until the number of elements is 6 times the number of processes.
143  int ser_ref_levels = prob_list.Get("Serial refinement levels", -1);
144 
145  // The number of times to refine in parallel. This determines the
146  // number of levels in the AMGe hierarchy, if amge_levels is set to 0.
147  const int par_ref_levels = prob_list.Get("Parallel refinement levels", 2);
148 
149  // Number of levels in the AMGe hierarchy. Should not be larger than
150  // par_ref_levels + 1. If set to 0, it will be interpreted as equal to
151  // par_ref_levels + 1.
152  const int amge_levels = prob_list.Get("AMGe levels", 0);
153 
154  // The order of the finite elements on the finest level.
155  const int feorder = prob_list.Get("Finite element order", 0);
156 
157  // The order of the polynomials to include in the coarse spaces
158  // (after interpolating them onto the fine space).
159  const int upscalingOrder = prob_list.Get("Upscaling order", 0);
160 
161  // A list of 1s and 0s stating which boundary attribute is appointed as
162  // essential. If only a single entry is given, it is applied to the whole
163  // boundary. That is, if a single 0 is given the whole boundary is "natural",
164  // while a single 1 means that the whole boundary is essential.
165  vector<int> par_ess_attr = prob_list.Get("Essential attributes",
166  vector<int> {1});
167 
168  // A list of (piecewise) constant values for the coefficient 'alpha', in
169  // accordance with the mesh attributes. If only a single entry is given, it
170  // is applied to the whole mesh/domain.
171  vector<double> alpha_vals = prob_list.Get("alpha values",
172  vector<double> {1.0});
173 
174  // A list of (piecewise) constant values for the coefficient 'beta', in
175  // accordance with the mesh attributes. If only a single entry is given, it
176  // is applied to the whole mesh/domain.
177  vector<double> beta_vals = prob_list.Get("beta values",
178  vector<double> {1.0});
179 
180  // The list of solvers to invoke.
181  auto list_of_solvers = prob_list.Get<list<string>>("List of linear solvers");
182 
183  ostringstream mesh_msg;
184  if (!myid)
185  {
186  mesh_msg << '\n' << string(50, '*') << '\n'
187  << "* Mesh: " << meshfile << "\n*\n"
188  << "* FE order: " << feorder << '\n'
189  << "* Upscaling order: " << upscalingOrder << "\n*\n";
190  }
191 
192  // Read the (serial) mesh from the given mesh file and uniformly refine it.
193  shared_ptr<ParMesh> pmesh;
194  {
195  if (!myid)
196  {
197  cout << "\nReading and refining serial mesh...\n";
198  cout << "Times to refine mesh in serial: " << ser_ref_levels << ".\n";
199  }
200 
201  ifstream imesh(meshfile);
202  if (!imesh)
203  {
204  if (!myid)
205  {
206  cerr << "ERROR: Cannot open mesh file: " << meshfile << ".\n";
207  }
208  return EXIT_FAILURE;
209  }
210 
211  auto mesh = make_unique<Mesh>(imesh, true, true);
212  imesh.close();
213 
214  for (int l = 0; l < ser_ref_levels; ++l)
215  {
216  if (!myid)
217  {
218  cout << "Refining mesh in serial: " << l + 1 << "...\n";
219  }
220  mesh->UniformRefinement();
221  }
222 
223  if (ser_ref_levels < 0)
224  {
225  ser_ref_levels = 0;
226  for (; mesh->GetNE() < 6 * num_ranks; ++ser_ref_levels)
227  {
228  if (!myid)
229  {
230  cout << "Refining mesh in serial: " << ser_ref_levels + 1
231  << "...\n";
232  }
233  mesh->UniformRefinement();
234  }
235  }
236 
237  if (!myid)
238  {
239  cout << "Times refined mesh in serial: " << ser_ref_levels << ".\n";
240  cout << "Building and refining parallel mesh...\n";
241  cout << "Times to refine mesh in parallel: " << par_ref_levels
242  << ".\n";
243  mesh_msg << "* Serial refinements: " << ser_ref_levels << '\n'
244  << "* Coarse mesh size: " << mesh->GetNE() << "\n*\n";
245  }
246 
247  pmesh = make_shared<ParMesh>(MPI_COMM_WORLD, *mesh);
248  }
249 
250  // Mark essential boundary attributes.
251  MFEM_VERIFY(par_ess_attr.size() <= 1 ||
252  par_ess_attr.size() == (unsigned) pmesh->bdr_attributes.Max(),
253  "Incorrect size of the essential attributes vector in parameters"
254  " input.");
255  vector<Array<int>> ess_attr(1);
256  ess_attr[0].SetSize(pmesh->bdr_attributes.Max());
257  if (par_ess_attr.size() == 0)
258  {
259  ess_attr[0] = 1;
260  }
261  else if (par_ess_attr.size() == 1)
262  {
263  ess_attr[0] = par_ess_attr[0];
264  }
265  else
266  {
267  for (unsigned i = 0; i < par_ess_attr.size(); ++i)
268  {
269  ess_attr[0][i] = par_ess_attr[i];
270  }
271  }
272 
273  // Initialize piecewise constant coefficients in the form.
274  MFEM_VERIFY(alpha_vals.size() <= 1 ||
275  alpha_vals.size() == (unsigned) pmesh->attributes.Max(),
276  "Incorrect size of the 'alpha' local values vector in parameters"
277  " input.");
278  MFEM_VERIFY(alpha_vals.size() <= 1 ||
279  alpha_vals.size() == (unsigned) pmesh->attributes.Max(),
280  "Incorrect size of the 'alpha' local values vector in parameters"
281  " input.");
282  PWConstCoefficient alpha(pmesh->attributes.Max());
283  PWConstCoefficient beta(pmesh->attributes.Max());
284 
285  if (alpha_vals.size() == 0)
286  {
287  alpha = 1.0;
288  }
289  else if (alpha_vals.size() == 1)
290  {
291  alpha = alpha_vals[0];
292  }
293  else
294  {
295  for (unsigned i = 0; i < alpha_vals.size(); ++i)
296  {
297  alpha(i+1) = alpha_vals[i];
298  }
299  }
300 
301  if (beta_vals.size() == 0)
302  {
303  beta = 1.0;
304  }
305  else if (beta_vals.size() == 1)
306  {
307  beta = beta_vals[0];
308  }
309  else
310  {
311  for (unsigned i = 0; i < beta_vals.size(); ++i)
312  {
313  beta(i+1) = beta_vals[i];
314  }
315  }
316 
317  // Refine the mesh in parallel.
318  const int nDimensions = pmesh->Dimension();
319 
320  // This is mainly because AMS and ADS (at least the way ParELAG uses them)
321  // are bound to be used in 3D. Note that, for the purpose of demonstration,
322  // some of the code below is still constructed in a way that is applicable in
323  // 2D as well, taking into account that case as well. Also, in 2D, ParELAG
324  // defaults to H(div) interpretation of form 1.
325  MFEM_VERIFY(nDimensions == 3, "Only 3D problems are currently supported.");
326 
327  const int nLevels = amge_levels <= 0 ? par_ref_levels + 1 : amge_levels;
328  MFEM_VERIFY(nLevels <= par_ref_levels + 1,
329  "Number of AMGe levels too high relative to parallel"
330  " refinements.");
331  vector<int> level_nElements(nLevels);
332  for (int l = 0; l < par_ref_levels; ++l)
333  {
334  if (!myid)
335  {
336  cout << "Refining mesh in parallel: " << l + 1
337  << (par_ref_levels - l > nLevels ? " (not in hierarchy)"
338  : " (in hierarchy)")
339  << "...\n";
340  }
341  if (par_ref_levels - l < nLevels)
342  {
343  level_nElements[par_ref_levels - l] = pmesh->GetNE();
344  }
345  pmesh->UniformRefinement();
346  }
347  level_nElements[0] = pmesh->GetNE();
348 
349  if (!myid)
350  {
351  cout << "Times refined mesh in parallel: " << par_ref_levels << ".\n";
352  }
353 
354  {
355  size_t local_num_elmts = pmesh->GetNE(), global_num_elmts;
356  MPI_Reduce(&local_num_elmts, &global_num_elmts, 1, GetMPIType<size_t>(0),
357  MPI_SUM, 0, MPI_COMM_WORLD);
358  if (!myid)
359  {
360  mesh_msg << "* Parallel refinements: " << par_ref_levels << '\n'
361  << "* Fine mesh size: " << global_num_elmts << '\n'
362  << "* Total levels: " << nLevels << '\n'
363  << string(50, '*') << "\n\n";
364  }
365  }
366 
367  if (!myid)
368  {
369  cout << mesh_msg.str();
370  }
371  init_timer.Stop();
372 
373  // Obtain the hierarchy of agglomerate topologies.
374  Timer agg_timer = TimeManager::AddTimer("Mesh Agglomeration -- Total");
375  Timer agg0_timer = TimeManager::AddTimer("Mesh Agglomeration -- Level 0");
376  if (!myid)
377  {
378  cout << "Agglomerating topology for " << nLevels - 1
379  << " coarse levels...\n";
380  }
381 
382  constexpr auto AT_elem = AgglomeratedTopology::ELEMENT;
383  // This partitioner simply geometrically coarsens the mesh by recovering the
384  // geometric coarse elements as agglomerate elements. That is, it reverts the
385  // MFEM uniform refinement procedure to provide agglomeration.
386  MFEMRefinedMeshPartitioner partitioner(nDimensions);
387  vector<shared_ptr<AgglomeratedTopology>> topology(nLevels);
388 
389  if (!myid)
390  {
391  cout << "Agglomerating level: 0...\n";
392  }
393 
394  topology[0] = make_shared<AgglomeratedTopology>(pmesh, nDimensions);
395 
396  if (!myid)
397  {
398  cout << "Level 0 global number of mesh entities: "
399  << topology[0]->
400  GetNumberGlobalTrueEntities((AgglomeratedTopology::Entity)0);
401  for (int j = 1; j <= nDimensions; ++j)
402  cout << ", " << topology[0]->
403  GetNumberGlobalTrueEntities((AgglomeratedTopology::Entity)j);
404  cout << endl;
405  }
406 
407  agg0_timer.Stop();
408 
409  for (int l = 0; l < nLevels - 1; ++l)
410  {
411  Timer aggl_timer = TimeManager::AddTimer(std::string("Mesh "
412  "Agglomeration -- Level ").
413  append(std::to_string(l+1)));
414  Array<int> partitioning(topology[l]->GetNumberLocalEntities(AT_elem));
415  partitioner.Partition(topology[l]->GetNumberLocalEntities(AT_elem),
416  level_nElements[l + 1], partitioning);
417 
418  if (!myid)
419  {
420  cout << "Agglomerating level: " << l + 1 << "...\n";
421  }
422 
423  topology[l + 1] = topology[l]->CoarsenLocalPartitioning(partitioning,
424  false, false, 2);
425  if (!myid)
426  {
427  cout << "Level " << l + 1 << " global number of mesh entities: "
428  << topology[l + 1]->
429  GetNumberGlobalTrueEntities((AgglomeratedTopology::Entity)0);
430  for (int j = 1; j <= nDimensions; ++j)
431  {
432  cout << ", " << topology[l + 1]->
433  GetNumberGlobalTrueEntities((AgglomeratedTopology::Entity)j);
434  }
435  cout << endl;
436  }
437  }
438 
439  agg_timer.Stop();
440 
441  if (visualize && nDimensions <= 3)
442  {
443  for (int l = 1; l < nLevels; ++l)
444  {
445  ShowTopologyAgglomeratedElements(topology[l].get(), pmesh.get());
446  }
447  }
448 
449  // Construct the hierarchy of spaces, thus forming a hierarchy of (partial)
450  // de Rham sequences.
451  Timer derham_timer = TimeManager::AddTimer("DeRhamSequence Construction -- "
452  "Total");
453  Timer derham0_timer = TimeManager::AddTimer("DeRhamSequence Construction -- "
454  "Level 0");
455  if (!myid)
456  {
457  cout << "Building the fine-level de Rham sequence...\n";
458  }
459 
460  vector<shared_ptr<DeRhamSequence>> sequence(topology.size());
461 
462  const int jform = DeRhamSequence::GetForm(nDimensions,
463  hcurl ? DeRhamSequence::HCURL :
464  DeRhamSequence::HDIV);
465  if (nDimensions == 3)
466  {
467  sequence[0] = make_shared<DeRhamSequence3D_FE>(topology[0], pmesh.get(),
468  feorder, true, false);
469  }
470  else
471  {
472  MFEM_VERIFY(nDimensions == 2, "Only 2D or 3D problems are supported "
473  "by the utilized ParELAG.");
474  if (hcurl)
475  {
476  MFEM_ABORT("No H(curl) 2D interpretation of form 1 is implemented.");
477  }
478  sequence[0] = make_shared<DeRhamSequence2D_Hdiv_FE>(topology[0],
479  pmesh.get(), feorder,
480  true, false);
481  }
482 
483  // To build H(curl) (form 1 in 3D), it is needed to obtain all forms and
484  // spaces with larger indices. To use the so called "Hiptmair smoothers", a
485  // one form lower is needed (H1, form 0). Anyway, to use AMS all forms and
486  // spaces to H1 (0 form) are needed. Therefore, the entire de Rham complex is
487  // constructed.
488  // To build H(div) (form 2 in 3D), it is needed to obtain all forms and
489  // spaces with larger indices. To use the so called "Hiptmair smoothers", a
490  // one form lower is needed (H(curl), form 1, in 3D). To use AMS and ADS, all
491  // forms and spaces to H1 (0 form) are needed. Therefore, the entire de Rham
492  // complex is constructed.
493  sequence[0]->SetjformStart(0);
494 
495  DeRhamSequenceFE *DRSequence_FE = sequence[0]->FemSequence();
496  MFEM_VERIFY(DRSequence_FE,
497  "Failed to obtain the fine-level de Rham sequence.");
498 
499  if (!myid)
500  {
501  cout << "Level 0 global number of dofs: "
502  << DRSequence_FE->GetDofHandler(0)->GetDofTrueDof().
503  GetTrueGlobalSize();
504  for (int j = 1; j <= nDimensions; ++j)
505  {
506  cout << ", " << DRSequence_FE->GetDofHandler(j)->GetDofTrueDof().
507  GetTrueGlobalSize();
508  }
509  cout << endl;
510  }
511 
512  if (!myid)
513  {
514  cout << "Setting coefficients and computing fine-level local "
515  << "matrices...\n";
516  }
517 
518  DRSequence_FE->ReplaceMassIntegrator(AT_elem, jform,
519  make_unique<VectorFEMassIntegrator>(beta), false);
520  if (hcurl && nDimensions == 3)
521  {
522  DRSequence_FE->ReplaceMassIntegrator(AT_elem, jform + 1,
523  make_unique<VectorFEMassIntegrator>(alpha), true);
524  }
525  else
526  {
527  DRSequence_FE->ReplaceMassIntegrator(AT_elem, jform + 1,
528  make_unique<MassIntegrator>(alpha), true);
529  }
530 
531  if (!myid)
532  {
533  cout << "Interpolating and setting polynomial targets...\n";
534  }
535 
536  DRSequence_FE->SetUpscalingTargets(nDimensions, upscalingOrder);
537  derham0_timer.Stop();
538 
539  if (!myid)
540  {
541  cout << "Building the coarse-level de Rham sequences...\n";
542  }
543 
544  for (int l = 0; l < nLevels - 1; ++l)
545  {
546  Timer derhaml_timer = TimeManager::AddTimer(std::string("DeRhamSequence "
547  "Construction -- Level ").
548  append(std::to_string(l+1)));
549  if (!myid)
550  {
551  cout << "Building the level " << l + 1 << " de Rham sequences...\n";
552  }
553 
554  sequence[l]->SetSVDTol(tolSVD);
555  sequence[l + 1] = sequence[l]->Coarsen();
556 
557  if (!myid)
558  {
559  auto DRSequence = sequence[l + 1];
560  cout << "Level " << l + 1 << " global number of dofs: "
561  << DRSequence->GetDofHandler(0)->GetDofTrueDof().
562  GetTrueGlobalSize();
563  for (int j = 1; j <= nDimensions; ++j)
564  {
565  cout << ", " << DRSequence->GetDofHandler(j)->GetDofTrueDof().
566  GetTrueGlobalSize();
567  }
568  cout << endl;
569  }
570  }
571  derham_timer.Stop();
572 
573  Timer assemble_timer = TimeManager::AddTimer("Fine Matrix Assembly");
574  if (!myid)
575  {
576  cout << "Assembling the fine-level system...\n";
577  }
578 
579  VectorFunctionCoefficient rhscoeff(nDimensions, rhsfunc);
580  VectorFunctionCoefficient solcoeff(nDimensions, bdrfunc);
581 
582  // Take the vector FE space and construct a RHS linear form on it. Then, move
583  // the linear form to a vector. This is local, i.e. on all known dofs for the
584  // process.
585  FiniteElementSpace *fespace = DRSequence_FE->GetFeSpace(jform);
586  auto rhsform = make_unique<LinearForm>(fespace);
587  rhsform->AddDomainIntegrator(new VectorFEDomainLFIntegrator(rhscoeff));
588  rhsform->Assemble();
589  unique_ptr<Vector> rhs = move(rhsform);
590 
591  // Obtain the boundary data. This is local, i.e. on all known dofs for the
592  // process.
593  auto solgf = make_unique<GridFunction>(fespace);
594  solgf->ProjectCoefficient(solcoeff);
595  unique_ptr<Vector> sol = move(solgf);
596 
597  // Create the parallel linear system.
598  const SharingMap& hcurlhdiv_dofTrueDof =
599  sequence[0]->GetDofHandler(jform)->GetDofTrueDof();
600 
601  // System RHS, B. It is defined on the true dofs owned by the process.
602  Vector B(hcurlhdiv_dofTrueDof.GetTrueLocalSize());
603 
604  // System matrix, A.
605  shared_ptr<HypreParMatrix> A;
606  {
607  // Get the mass and derivative operators.
608  // For H(curl):
609  // M1 represents the form (beta u, v) on H(curl) vector fields.
610  // M2 represents the form (alpha u, v) on H(div) vector fields, in 3D.
611  // D1 is the curl operator from H(curl) vector fields to H(div) vector
612  // fields, in 3D.
613  // In 2D, instead of considering H(div) vector fields, L2 scalar fields
614  // are to be considered.
615  // Thus, D1^T * M2 * D1 represents the form (alpha curl u, curl v) on
616  // H(curl) vector fields.
617  // For H(div):
618  // M1 represents the form (beta u, v) on H(div) vector fields.
619  // M2 represents the form (alpha u, v) on L2 scalar fields.
620  // D1 is the divergence operator from H(div) vector fields to L2 scalar
621  // fields.
622  // Thus, D1^T * M2 * D1 represents the form (alpha div u, div v) on H(div)
623  // vector fields.
624  auto M1 = sequence[0]->ComputeMassOperator(jform),
625  M2 = sequence[0]->ComputeMassOperator(jform + 1);
626  auto D1 = sequence[0]->GetDerivativeOperator(jform);
627 
628  // spA = D1^T * M2 * D1 + M1 represents the respective H(curl) or H(div)
629  // form:
630  // (alpha curl u, curl v) + (beta u, v), on H(curl) vector fields;
631  // (alpha div u, div v) + (beta u, v), on H(div) vector fields.
632  // This is local, i.e. on all known dofs for the process.
633  auto spA = ToUnique(Add(*M1, *ToUnique(RAP(*D1, *M2, *D1))));
634 
635  // Eliminate the boundary conditions
636  Array<int> marker(spA->Height());
637  marker = 0;
638  sequence[0]->GetDofHandler(jform)->MarkDofsOnSelectedBndr(ess_attr[0],
639  marker);
640 
641  for (int i = 0; i < spA->Height(); ++i)
642  {
643  if (marker[i])
644  {
645  spA->EliminateRowCol(i, sol->Elem(i), *rhs);
646  }
647  }
648 
649  A = Assemble(hcurlhdiv_dofTrueDof, *spA, hcurlhdiv_dofTrueDof);
650  hcurlhdiv_dofTrueDof.Assemble(*rhs, B);
651  }
652  if (!myid)
653  {
654  cout << "A size: " << A->GetGlobalNumRows() << 'x'
655  << A->GetGlobalNumCols() << '\n' << " A NNZ: " << A->NNZ() << '\n';
656  }
657  MFEM_VERIFY(B.Size() == A->Height(),
658  "Matrix and vector size are incompatible.");
659  assemble_timer.Stop();
660 
661  // Perform the solves.
662  Timer solvers_timer = TimeManager::AddTimer("Solvers -- Total");
663  if (!myid)
664  {
665  cout << "\nRunning fine-level solvers...\n\n";
666  }
667 
668  // Create the solver library.
669  auto lib = SolverLibrary::CreateLibrary(
670  master_list->Sublist("Preconditioner Library"));
671 
672  // Loop through the solvers.
673  for (const auto& solver_name : list_of_solvers)
674  {
675  Timer solver_timer = TimeManager::AddTimer(std::string("Solver \"").
676  append(solver_name).
677  append("\" -- Total"));
678  // Get the solver factory.
679  auto solver_factory = lib->GetSolverFactory(solver_name);
680  auto solver_state = solver_factory->GetDefaultState();
681  solver_state->SetDeRhamSequence(sequence[0]);
682  solver_state->SetBoundaryLabels(ess_attr);
683  solver_state->SetForms({jform});
684 
685  // Build the solver.
686  Timer build_timer = TimeManager::AddTimer(std::string("Solver \"").
687  append(solver_name).
688  append("\" -- Build"));
689  if (!myid)
690  {
691  cout << "Building solver \"" << solver_name << "\"...\n";
692  }
693  unique_ptr<Solver> solver = solver_factory->BuildSolver(A, *solver_state);
694  build_timer.Stop();
695 
696  // Run the solver.
697  Timer pre_timer = TimeManager::AddTimer(std::string("Solver \"").
698  append(solver_name).
699  append("\" -- Pre-solve"));
700  if (!myid)
701  {
702  cout << "Solving system with \"" << solver_name << "\"...\n";
703  }
704 
705  // Note that X is on true dofs owned by the process, while x is on local
706  // dofs that are known to the process.
707  Vector X(A->Width()), x(sequence[0]->GetNumberOfDofs(jform));
708  X=0.0;
709 
710  {
711  Vector tmp(A->Height());
712  A->Mult(X, tmp);
713  tmp *= -1.0;
714  tmp += B;
715 
716  double local_norm = tmp.Norml2();
717  local_norm *= local_norm;
718  double global_norm;
719  MPI_Reduce(&local_norm, &global_norm, 1, GetMPIType(local_norm),
720  MPI_SUM, 0, MPI_COMM_WORLD);
721  if (!myid)
722  {
723  cout << "Initial residual l2 norm: " << sqrt(global_norm) << '\n';
724  }
725  }
726  pre_timer.Stop();
727 
728  // Perform the solve.
729  Timer solve_timer = TimeManager::AddTimer(std::string("Solver \"").
730  append(solver_name).
731  append("\" -- Solve"));
732  solver->Mult(B, X);
733  solve_timer.Stop();
734 
735  Timer post_timer = TimeManager::AddTimer(std::string("Solver \"").
736  append(solver_name).
737  append("\" -- Post-solve"));
738  {
739  Vector tmp(A->Height());
740  A->Mult(X, tmp);
741  tmp *= -1.0;
742  tmp += B;
743 
744  double local_norm = tmp.Norml2();
745  local_norm *= local_norm;
746  double global_norm;
747  MPI_Reduce(&local_norm, &global_norm, 1, GetMPIType(local_norm),
748  MPI_SUM, 0, MPI_COMM_WORLD);
749  if (!myid)
750  {
751  cout << "Final residual l2 norm: " << sqrt(global_norm) << '\n';
752  }
753  }
754 
755  if (!myid)
756  {
757  cout << "Solver \"" << solver_name << "\" finished.\n";
758  }
759 
760  // Visualize the solution.
761  if (visualize)
762  {
763  hcurlhdiv_dofTrueDof.Distribute(X, x);
764  MultiVector tmp(x.GetData(), 1, x.Size());
765  sequence[0]->show(jform, tmp);
766  }
767  post_timer.Stop();
768  }
769  solvers_timer.Stop();
770 
771  total_timer.Stop();
772  TimeManager::Print();
773 
774  if (!myid)
775  {
776  cout << "\nFinished.\n";
777  }
778 
779  return EXIT_SUCCESS;
780 }
781 
782 // A vector field, used for setting boundary conditions.
783 void bdrfunc(const Vector &p, Vector &F)
784 {
785  F = 0.0;
786 }
787 
788 // The right hand side.
789 void rhsfunc(const Vector &p, Vector &f)
790 {
791  f = 1.0;
792 }
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(), hypre will be finalized automatically at program exit.
Definition: hypre.hpp:80
double & Elem(int i)
Access Vector entries. Index i = 0 .. size-1.
Definition: vector.cpp:101
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:812
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
static int WorldSize()
Return the size of MPI_COMM_WORLD.
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2282
void rhsfunc(const Vector &, Vector &)
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
void visualize(ostream &os, Mesh *mesh, GridFunction *deformed_nodes, GridFunction *field, const char *field_name=NULL, bool init_vis=false)
Definition: ex10.cpp:381
static void Init()
Singleton creation with Mpi::Init();.
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
Definition: densemat.cpp:3213
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
A general vector function coefficient.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
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
A piecewise constant coefficient with the constants keyed off the element attribute numbers...
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
void bdrfunc(const Vector &, Vector &)
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Definition: lininteg.hpp:346
const double alpha
Definition: ex15.cpp:369
Vector data type.
Definition: vector.hpp:60
int main()
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150