MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
MultilevelHcurlHdivSolver.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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
45using namespace mfem;
46using namespace parelag;
47using namespace std;
48
49void bdrfunc(const Vector &, Vector &);
50void rhsfunc(const Vector &, Vector &);
51
52int 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();
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 real_t 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<real_t> alpha_vals = prob_list.Get("alpha values",
172 vector<real_t> {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<real_t> beta_vals = prob_list.Get("beta values",
178 vector<real_t> {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 real_t local_norm = tmp.Norml2();
717 local_norm *= local_norm;
718 real_t 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 real_t local_norm = tmp.Norml2();
745 local_norm *= local_norm;
746 real_t 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.
783void bdrfunc(const Vector &p, Vector &F)
784{
785 F = 0.0;
786}
787
788// The right hand side.
789void rhsfunc(const Vector &p, Vector &f)
790{
791 f = 1.0;
792}
void bdrfunc(const Vector &, Vector &)
void rhsfunc(const Vector &, Vector &)
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
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).
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.
A piecewise constant coefficient with the constants keyed off the element attribute numbers.
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Definition lininteg.hpp:347
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
real_t Norml2() const
Returns the l2 norm of the vector.
Definition vector.cpp:832
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
Definition vector.hpp:227
Vector beta
void visualize(ostream &os, Mesh *mesh, GridFunction *deformed_nodes, GridFunction *field, const char *field_name=NULL, bool init_vis=false)
Definition ex10.cpp:381
const real_t alpha
Definition ex15.cpp:369
int main()
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
void Add(const DenseMatrix &A, const DenseMatrix &B, real_t alpha, DenseMatrix &C)
C = A + alpha*B.
real_t p(const Vector &x, real_t t)