MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
block-solvers.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 // Block Solvers Miniapp: Compare Saddle Point System Solvers
14 // ----------------------------------------------------------
15 //
16 // This miniapp compares various linear solvers for the saddle point system
17 // obtained from mixed finite element discretization of the simple mixed Darcy
18 // problem in ex5p
19 //
20 // k*u + grad p = f
21 // - div u = g
22 //
23 // with natural boundary condition -p = <given pressure>. We use a given exact
24 // solution (u,p) and compute the corresponding r.h.s. (f,g). We discretize
25 // with Raviart-Thomas finite elements (velocity u) and piecewise discontinuous
26 // polynomials (pressure p).
27 //
28 // The solvers being compared include:
29 // 1. The divergence free solver (couple and decoupled modes)
30 // 2. MINRES preconditioned by a block diagonal preconditioner
31 //
32 // We recommend viewing example 5 before viewing this miniapp.
33 //
34 // Sample runs:
35 //
36 // mpirun -np 8 block-solvers -r 2 -o 0
37 // mpirun -np 8 block-solvers -m anisotropic.mesh -c anisotropic.coeff -be anisotropic.bdr
38 //
39 //
40 // NOTE: The coefficient file (provided through -c) defines a piecewise constant
41 // scalar coefficient k. The number of entries in this file should equal
42 // to the number of "element attributes" in the mesh file. The value of
43 // the coefficient in elements with the i-th attribute is given by the
44 // i-th entry of the coefficient file.
45 //
46 //
47 // NOTE: The essential boundary attribute file (provided through -eb) defines
48 // which attributes to impose essential boundary condition (on u). The
49 // number of entries in this file should equal to the number of "boundary
50 // attributes" in the mesh file. If the i-th entry of the file is nonzero
51 // (respectively 0), essential (respectively natural) boundary condition
52 // will be imposed on boundary with the i-th attribute.
53 
54 #include "mfem.hpp"
55 #include "div_free_solver.hpp"
56 #include <fstream>
57 #include <iostream>
58 #include <memory>
59 
60 using namespace std;
61 using namespace mfem;
62 using namespace blocksolvers;
63 
64 // Exact solution, u and p, and r.h.s., f and g.
65 void u_exact(const Vector & x, Vector & u);
66 double p_exact(const Vector & x);
67 void f_exact(const Vector & x, Vector & f);
68 double g_exact(const Vector & x);
69 double natural_bc(const Vector & x);
70 
71 /** Wrapper for assembling the discrete Darcy problem (ex5p)
72  [ M B^T ] [u] = [f]
73  [ B 0 ] [p] = [g]
74  where:
75  M = \int_\Omega (k u_h) \cdot v_h dx,
76  B = -\int_\Omega (div_h u_h) q_h dx,
77  f = \int_\Omega f_exact v_h dx + \int_D natural_bc v_h dS,
78  g = \int_\Omega g_exact q_h dx,
79  u_h, v_h \in R_h (Raviart-Thomas finite element space),
80  q_h \in W_h (piecewise discontinuous polynomials),
81  D: subset of the boundary where natural boundary condition is imposed. */
82 class DarcyProblem
83 {
84  OperatorPtr M_;
85  OperatorPtr B_;
86  Vector rhs_;
87  Vector ess_data_;
88  ParGridFunction u_;
89  ParGridFunction p_;
90  ParMesh mesh_;
92  FunctionCoefficient pcoeff_;
93  DFSSpaces dfs_spaces_;
94  const IntegrationRule *irs_[Geometry::NumGeom];
95 public:
96  DarcyProblem(Mesh &mesh, int num_refines, int order, const char *coef_file,
97  Array<int> &ess_bdr, DFSParameters param);
98 
99  HypreParMatrix& GetM() { return *M_.As<HypreParMatrix>(); }
100  HypreParMatrix& GetB() { return *B_.As<HypreParMatrix>(); }
101  const Vector& GetRHS() { return rhs_; }
102  const Vector& GetEssentialBC() { return ess_data_; }
103  const DFSData& GetDFSData() const { return dfs_spaces_.GetDFSData(); }
104  void ShowError(const Vector &sol, bool verbose);
105  void VisualizeSolution(const Vector &sol, string tag);
106 };
107 
108 DarcyProblem::DarcyProblem(Mesh &mesh, int num_refs, int order,
109  const char *coef_file, Array<int> &ess_bdr,
110  DFSParameters dfs_param)
111  : mesh_(MPI_COMM_WORLD, mesh), ucoeff_(mesh.Dimension(), u_exact),
112  pcoeff_(p_exact), dfs_spaces_(order, num_refs, &mesh_, ess_bdr, dfs_param)
113 {
114  for (int l = 0; l < num_refs; l++)
115  {
116  mesh_.UniformRefinement();
117  dfs_spaces_.CollectDFSData();
118  }
119 
120  Vector coef_vector(mesh.GetNE());
121  coef_vector = 1.0;
122  if (std::strcmp(coef_file, ""))
123  {
124  ifstream coef_str(coef_file);
125  coef_vector.Load(coef_str, mesh.GetNE());
126  }
127  PWConstCoefficient mass_coeff(coef_vector);
128  VectorFunctionCoefficient fcoeff(mesh_.Dimension(), f_exact);
131 
132  u_.SetSpace(dfs_spaces_.GetHdivFES());
133  p_.SetSpace(dfs_spaces_.GetL2FES());
134  p_ = 0.0;
135  u_ = 0.0;
136  u_.ProjectBdrCoefficientNormal(ucoeff_, ess_bdr);
137 
138  ParLinearForm fform(dfs_spaces_.GetHdivFES());
140  fform.AddBoundaryIntegrator(new VectorFEBoundaryFluxLFIntegrator(natcoeff));
141  fform.Assemble();
142 
143  ParLinearForm gform(dfs_spaces_.GetL2FES());
144  gform.AddDomainIntegrator(new DomainLFIntegrator(gcoeff));
145  gform.Assemble();
146 
147  ParBilinearForm mVarf(dfs_spaces_.GetHdivFES());
148  ParMixedBilinearForm bVarf(dfs_spaces_.GetHdivFES(), dfs_spaces_.GetL2FES());
149 
150  mVarf.AddDomainIntegrator(new VectorFEMassIntegrator(mass_coeff));
151  mVarf.Assemble();
152  mVarf.EliminateEssentialBC(ess_bdr, u_, fform);
153  mVarf.Finalize();
154  M_.Reset(mVarf.ParallelAssemble());
155 
156  bVarf.AddDomainIntegrator(new VectorFEDivergenceIntegrator);
157  bVarf.Assemble();
158  bVarf.SpMat() *= -1.0;
159  bVarf.EliminateTrialDofs(ess_bdr, u_, gform);
160  bVarf.Finalize();
161  B_.Reset(bVarf.ParallelAssemble());
162 
163  rhs_.SetSize(M_->NumRows() + B_->NumRows());
164  Vector rhs_block0(rhs_.GetData(), M_->NumRows());
165  Vector rhs_block1(rhs_.GetData()+M_->NumRows(), B_->NumRows());
166  fform.ParallelAssemble(rhs_block0);
167  gform.ParallelAssemble(rhs_block1);
168 
169  ess_data_.SetSize(M_->NumRows() + B_->NumRows());
170  ess_data_ = 0.0;
171  Vector ess_data_block0(ess_data_.GetData(), M_->NumRows());
172  u_.ParallelProject(ess_data_block0);
173 
174  int order_quad = max(2, 2*order+1);
175  for (int i=0; i < Geometry::NumGeom; ++i)
176  {
177  irs_[i] = &(IntRules.Get(i, order_quad));
178  }
179 }
180 
181 void DarcyProblem::ShowError(const Vector& sol, bool verbose)
182 {
183  u_.Distribute(Vector(sol.GetData(), M_->NumRows()));
184  p_.Distribute(Vector(sol.GetData()+M_->NumRows(), B_->NumRows()));
185 
186  double err_u = u_.ComputeL2Error(ucoeff_, irs_);
187  double norm_u = ComputeGlobalLpNorm(2, ucoeff_, mesh_, irs_);
188  double err_p = p_.ComputeL2Error(pcoeff_, irs_);
189  double norm_p = ComputeGlobalLpNorm(2, pcoeff_, mesh_, irs_);
190 
191  if (!verbose) { return; }
192  cout << "|| u_h - u_ex || / || u_ex || = " << err_u / norm_u << "\n";
193  cout << "|| p_h - p_ex || / || p_ex || = " << err_p / norm_p << "\n";
194 }
195 
196 void DarcyProblem::VisualizeSolution(const Vector& sol, string tag)
197 {
198  int num_procs, myid;
199  MPI_Comm_size(mesh_.GetComm(), &num_procs);
200  MPI_Comm_rank(mesh_.GetComm(), &myid);
201 
202  u_.Distribute(Vector(sol.GetData(), M_->NumRows()));
203  p_.Distribute(Vector(sol.GetData()+M_->NumRows(), B_->NumRows()));
204 
205  const char vishost[] = "localhost";
206  const int visport = 19916;
207  socketstream u_sock(vishost, visport);
208  u_sock << "parallel " << num_procs << " " << myid << "\n";
209  u_sock.precision(8);
210  u_sock << "solution\n" << mesh_ << u_ << "window_title 'Velocity ("
211  << tag << " solver)'" << endl;
212  MPI_Barrier(mesh_.GetComm());
213  socketstream p_sock(vishost, visport);
214  p_sock << "parallel " << num_procs << " " << myid << "\n";
215  p_sock.precision(8);
216  p_sock << "solution\n" << mesh_ << p_ << "window_title 'Pressure ("
217  << tag << " solver)'" << endl;
218 }
219 
220 bool IsAllNeumannBoundary(const Array<int>& ess_bdr_attr)
221 {
222  for (int attr : ess_bdr_attr) { if (attr == 0) { return false; } }
223  return true;
224 }
225 
226 int main(int argc, char *argv[])
227 {
228 #ifdef HYPRE_USING_GPU
229  cout << "\nAs of mfem-4.3 and hypre-2.22.0 (July 2021) this miniapp\n"
230  << "is NOT supported with the GPU version of hypre.\n\n";
231  return 242;
232 #endif
233 
234  // Initialize MPI and HYPRE.
235  Mpi::Init(argc, argv);
236  Hypre::Init();
237 
238  StopWatch chrono;
239  auto ResetTimer = [&chrono]() { chrono.Clear(); chrono.Start(); };
240 
241  // Parse command-line options.
242  const char *mesh_file = "../../data/beam-hex.mesh";
243  const char *coef_file = "";
244  const char *ess_bdr_attr_file = "";
245  int order = 0;
246  int par_ref_levels = 2;
247  bool show_error = false;
248  bool visualization = false;
249  DFSParameters param;
250  OptionsParser args(argc, argv);
251  args.AddOption(&mesh_file, "-m", "--mesh",
252  "Mesh file to use.");
253  args.AddOption(&order, "-o", "--order",
254  "Finite element order (polynomial degree).");
255  args.AddOption(&par_ref_levels, "-r", "--ref",
256  "Number of parallel refinement steps.");
257  args.AddOption(&coef_file, "-c", "--coef",
258  "Coefficient file to use.");
259  args.AddOption(&ess_bdr_attr_file, "-eb", "--ess-bdr",
260  "Essential boundary attribute file to use.");
261  args.AddOption(&show_error, "-se", "--show-error", "-no-se",
262  "--no-show-error",
263  "Show or not show approximation error.");
264  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
265  "--no-visualization",
266  "Enable or disable GLVis visualization.");
267  args.Parse();
268  if (!args.Good())
269  {
270  if (Mpi::Root()) { args.PrintUsage(cout); }
271  return 1;
272  }
273  if (Mpi::Root()) { args.PrintOptions(cout); }
274 
275  if (Mpi::Root() && par_ref_levels == 0)
276  {
277  std::cout << "WARNING: DivFree solver is equivalent to BDPMinresSolver "
278  << "when par_ref_levels == 0.\n";
279  }
280 
281  // Initialize the mesh, boundary attributes, and solver parameters
282  Mesh *mesh = new Mesh(mesh_file, 1, 1);
283  int dim = mesh->Dimension();
284  int ser_ref_lvls =
285  (int)ceil(log(Mpi::WorldSize()/mesh->GetNE())/log(2.)/dim);
286  for (int i = 0; i < ser_ref_lvls; ++i)
287  {
288  mesh->UniformRefinement();
289  }
290 
291  Array<int> ess_bdr(mesh->bdr_attributes.Max());
292  ess_bdr = 0;
293  if (std::strcmp(ess_bdr_attr_file, ""))
294  {
295  ifstream ess_bdr_attr_str(ess_bdr_attr_file);
296  ess_bdr.Load(mesh->bdr_attributes.Max(), ess_bdr_attr_str);
297  }
298  if (IsAllNeumannBoundary(ess_bdr))
299  {
300  if (Mpi::Root())
301  {
302  cout << "\nSolution is not unique when Neumann boundary condition is "
303  << "imposed on the entire boundary. \nPlease provide a different "
304  << "boundary condition.\n";
305  }
306  delete mesh;
307  return 0;
308  }
309 
310  string line = "**********************************************************\n";
311 
312  ResetTimer();
313 
314  // Generate components of the saddle point problem
315  DarcyProblem darcy(*mesh, par_ref_levels, order, coef_file, ess_bdr, param);
316  HypreParMatrix& M = darcy.GetM();
317  HypreParMatrix& B = darcy.GetB();
318  const DFSData& DFS_data = darcy.GetDFSData();
319  delete mesh;
320 
321  if (Mpi::Root())
322  {
323  cout << line << "System assembled in " << chrono.RealTime() << "s.\n";
324  cout << "Dimension of the physical space: " << dim << "\n";
325  cout << "Size of the discrete Darcy system: " << M.M() + B.M() << "\n";
326  if (par_ref_levels > 0)
327  {
328  cout << "Dimension of the divergence free subspace: "
329  << DFS_data.C.back().Ptr()->NumCols() << "\n\n";
330  }
331  }
332 
333  // Setup various solvers for the discrete problem
334  std::map<const DarcySolver*, double> setup_time;
335  ResetTimer();
336  BDPMinresSolver bdp(M, B, param);
337  setup_time[&bdp] = chrono.RealTime();
338 
339  ResetTimer();
340  DivFreeSolver dfs_dm(M, B, DFS_data);
341  setup_time[&dfs_dm] = chrono.RealTime();
342 
343  ResetTimer();
344  const_cast<bool&>(DFS_data.param.coupled_solve) = true;
345  DivFreeSolver dfs_cm(M, B, DFS_data);
346  setup_time[&dfs_cm] = chrono.RealTime();
347 
348  std::map<const DarcySolver*, std::string> solver_to_name;
349  solver_to_name[&bdp] = "Block-diagonal-preconditioned MINRES";
350  solver_to_name[&dfs_dm] = "Divergence free (decoupled mode)";
351  solver_to_name[&dfs_cm] = "Divergence free (coupled mode)";
352 
353  // Solve the problem using all solvers
354  for (const auto& solver_pair : solver_to_name)
355  {
356  auto& solver = solver_pair.first;
357  auto& name = solver_pair.second;
358 
359  Vector sol = darcy.GetEssentialBC();
360 
361  ResetTimer();
362  solver->Mult(darcy.GetRHS(), sol);
363  chrono.Stop();
364 
365  if (Mpi::Root())
366  {
367  cout << line << name << " solver:\n Setup time: "
368  << setup_time[solver] << "s.\n Solve time: "
369  << chrono.RealTime() << "s.\n Total time: "
370  << setup_time[solver] + chrono.RealTime() << "s.\n"
371  << " Iteration count: " << solver->GetNumIterations() <<"\n\n";
372  }
373  if (show_error && std::strcmp(coef_file, "") == 0)
374  {
375  darcy.ShowError(sol, Mpi::Root());
376  }
377  else if (show_error && Mpi::Root())
378  {
379  cout << "Exact solution is unknown for coefficient '" << coef_file
380  << "'.\nApproximation error is computed in this case!\n\n";
381  }
382 
383  if (visualization) { darcy.VisualizeSolution(sol, name); }
384 
385  }
386 
387  return 0;
388 }
389 
390 void u_exact(const Vector & x, Vector & u)
391 {
392  double xi(x(0));
393  double yi(x(1));
394  double zi(x.Size() == 3 ? x(2) : 0.0);
395 
396  u(0) = - exp(xi)*sin(yi)*cos(zi);
397  u(1) = - exp(xi)*cos(yi)*cos(zi);
398  if (x.Size() == 3)
399  {
400  u(2) = exp(xi)*sin(yi)*sin(zi);
401  }
402 }
403 
404 double p_exact(const Vector & x)
405 {
406  double xi(x(0));
407  double yi(x(1));
408  double zi(x.Size() == 3 ? x(2) : 0.0);
409  return exp(xi)*sin(yi)*cos(zi);
410 }
411 
412 void f_exact(const Vector & x, Vector & f)
413 {
414  f = 0.0;
415 }
416 
417 double g_exact(const Vector & x)
418 {
419  if (x.Size() == 3) { return -p_exact(x); }
420  return 0;
421 }
422 
423 double natural_bc(const Vector & x)
424 {
425  return (-p_exact(x));
426 }
const char vishost[]
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:97
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:67
void Load(std::istream &in, int fmt=0)
Read an Array from the stream in using format fmt. The format fmt can be:
Definition: array.cpp:53
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
static const int NumGeom
Definition: geom.hpp:42
bool IsAllNeumannBoundary(const Array< int > &ess_bdr_attr)
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:521
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:928
double RealTime()
Definition: tic_toc.cpp:426
FDualNumber< tbase > exp(const FDualNumber< tbase > &f)
exp([dual number])
Definition: fdual.hpp:499
static int WorldSize()
Return the size of MPI_COMM_WORLD.
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:208
HYPRE_BigInt M() const
Returns the global number of rows.
Definition: hypre.hpp:577
void Stop()
Stop the stopwatch.
Definition: tic_toc.cpp:416
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds a domain integrator. Assumes ownership of bfi.
Data for the divergence free solver.
Class for parallel linear form.
Definition: plinearform.hpp:26
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
Wrapper for the block-diagonal-preconditioned MINRES defined in ex5p.cpp.
double natural_bc(const Vector &x)
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9737
Timing object.
Definition: tic_toc.hpp:34
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
static void Init()
Singleton creation with Mpi::Init();.
const int visport
int Dimension() const
Definition: mesh.hpp:999
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
void Start()
Start the stopwatch. The elapsed time is not cleared.
Definition: tic_toc.cpp:411
FDualNumber< tbase > cos(const FDualNumber< tbase > &f)
cos([dual number])
Definition: fdual.hpp:471
static bool Root()
Return true if the rank in MPI_COMM_WORLD is zero.
double p_exact(const Vector &x)
Definition: ex24.cpp:379
A general vector function coefficient.
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:39
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
void u_exact(const Vector &x, Vector &u)
FDualNumber< tbase > sin(const FDualNumber< tbase > &f)
sin([dual number])
Definition: fdual.hpp:572
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 f_exact(const Vector &, Vector &)
Definition: ex3.cpp:258
Parameters for the divergence free solver.
double ComputeGlobalLpNorm(double p, Coefficient &coeff, ParMesh &pmesh, const IntegrationRule *irs[])
Compute the global Lp norm of a function f. .
FDualNumber< tbase > log(const FDualNumber< tbase > &f)
log([dual number])
Definition: fdual.hpp:515
A piecewise constant coefficient with the constants keyed off the element attribute numbers...
Definition: coefficient.hpp:94
std::vector< OperatorPtr > C
Class for parallel bilinear form using different test and trial FE spaces.
int dim
Definition: ex24.cpp:53
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
Class for parallel bilinear form.
double g_exact(const Vector &x)
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Definition: lininteg.hpp:266
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:337
Class for parallel meshes.
Definition: pmesh.hpp:32
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
int main()
void Clear()
Clear the elapsed time on the stopwatch and restart it if it&#39;s running.
Definition: tic_toc.cpp:406
double f(const Vector &p)
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150