MFEM  v4.5.2
Finite element discretization library
div_free_solver.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 #include "div_free_solver.hpp"
13 
14 using namespace std;
15 using namespace mfem;
16 using namespace blocksolvers;
17 
18 void SetOptions(IterativeSolver& solver, const IterSolveParameters& param)
19 {
20  solver.SetPrintLevel(param.print_level);
21  solver.SetMaxIter(param.max_iter);
22  solver.SetAbsTol(param.abs_tol);
23  solver.SetRelTol(param.rel_tol);
24 }
25 
27  const HypreParMatrix& P)
28 {
29  OperatorPtr R(Rt.Transpose());
30  OperatorPtr RA(ParMult(R.As<HypreParMatrix>(), &A));
31  return ParMult(RA.As<HypreParMatrix>(), &P, true);
32 }
33 
34 void GetRowColumnsRef(const SparseMatrix& A, int row, Array<int>& cols)
35 {
36  cols.MakeRef(const_cast<int*>(A.GetRowColumns(row)), A.RowSize(row));
37 }
38 
40 {
41  int* I = new int[fes.GetNE()+1];
42  copy_n(fes.GetElementToDofTable().GetI(), fes.GetNE()+1, I);
43  Array<int> J(new int[I[fes.GetNE()]], I[fes.GetNE()]);
44  copy_n(fes.GetElementToDofTable().GetJ(), J.Size(), J.begin());
45  fes.AdjustVDofs(J);
46  double* D = new double[J.Size()];
47  fill_n(D, J.Size(), 1.0);
48  return SparseMatrix(I, J, D, fes.GetNE(), fes.GetVSize());
49 }
50 
51 DFSSpaces::DFSSpaces(int order, int num_refine, ParMesh *mesh,
52  const Array<int>& ess_attr, const DFSParameters& param)
53  : hdiv_fec_(order, mesh->Dimension()), l2_fec_(order, mesh->Dimension()),
54  l2_0_fec_(0, mesh->Dimension()), ess_bdr_attr_(ess_attr), level_(0)
55 {
56  if (mesh->GetElement(0)->GetType() == Element::TETRAHEDRON && order)
57  {
58  mfem_error("DFSDataCollector: High order spaces on tetrahedra are not supported");
59  }
60 
61  data_.param = param;
62 
63  if (mesh->Dimension() == 3)
64  {
65  hcurl_fec_.reset(new ND_FECollection(order+1, mesh->Dimension()));
66  }
67  else
68  {
69  hcurl_fec_.reset(new H1_FECollection(order+1, mesh->Dimension()));
70  }
71 
72  all_bdr_attr_.SetSize(ess_attr.Size(), 1);
73  hdiv_fes_.reset(new ParFiniteElementSpace(mesh, &hdiv_fec_));
74  l2_fes_.reset(new ParFiniteElementSpace(mesh, &l2_fec_));
75  coarse_hdiv_fes_.reset(new ParFiniteElementSpace(*hdiv_fes_));
76  coarse_l2_fes_.reset(new ParFiniteElementSpace(*l2_fes_));
77  l2_0_fes_.reset(new ParFiniteElementSpace(mesh, &l2_0_fec_));
78  l2_0_fes_->SetUpdateOperatorType(Operator::MFEM_SPARSEMAT);
79  el_l2dof_.reserve(num_refine+1);
80  el_l2dof_.push_back(ElemToDof(*coarse_l2_fes_));
81 
82  data_.agg_hdivdof.resize(num_refine);
83  data_.agg_l2dof.resize(num_refine);
84  data_.P_hdiv.resize(num_refine, OperatorPtr(Operator::Hypre_ParCSR));
85  data_.P_l2.resize(num_refine, OperatorPtr(Operator::Hypre_ParCSR));
86  data_.Q_l2.resize(num_refine);
87  hdiv_fes_->GetEssentialTrueDofs(ess_attr, data_.coarsest_ess_hdivdofs);
88  data_.C.resize(num_refine+1);
89 
90  hcurl_fes_.reset(new ParFiniteElementSpace(mesh, hcurl_fec_.get()));
91  coarse_hcurl_fes_.reset(new ParFiniteElementSpace(*hcurl_fes_));
92  data_.P_hcurl.resize(num_refine, OperatorPtr(Operator::Hypre_ParCSR));
93 }
94 
96  const SparseMatrix& agg_elem,
97  const SparseMatrix& elem_dof,
98  const HypreParMatrix& dof_truedof,
99  Array<HYPRE_BigInt>& agg_starts)
100 {
101  OperatorPtr agg_dof(Mult(agg_elem, elem_dof));
102  SparseMatrix& agg_dof_ref = *agg_dof.As<SparseMatrix>();
103  OperatorPtr agg_tdof(dof_truedof.LeftDiagMult(agg_dof_ref, agg_starts));
104  OperatorPtr agg_tdof_T(agg_tdof.As<HypreParMatrix>()->Transpose());
105  SparseMatrix tdof_agg, is_shared;
106  HYPRE_BigInt* trash;
107  agg_tdof_T.As<HypreParMatrix>()->GetDiag(tdof_agg);
108  agg_tdof_T.As<HypreParMatrix>()->GetOffd(is_shared, trash);
109 
110  int * I = new int [tdof_agg.NumRows()+1]();
111  int * J = new int[tdof_agg.NumNonZeroElems()];
112 
113  Array<int> is_bdr;
114  FiniteElementSpace::ListToMarker(bdr_truedofs, tdof_agg.NumRows(), is_bdr);
115 
116  int counter = 0;
117  for (int i = 0; i < tdof_agg.NumRows(); ++i)
118  {
119  bool agg_bdr = is_bdr[i] || is_shared.RowSize(i) || tdof_agg.RowSize(i)>1;
120  if (agg_bdr) { I[i+1] = I[i]; continue; }
121  I[i+1] = I[i] + 1;
122  J[counter++] = tdof_agg.GetRowColumns(i)[0];
123  }
124 
125  double * D = new double[I[tdof_agg.NumRows()]];
126  std::fill_n(D, I[tdof_agg.NumRows()], 1.0);
127 
128  SparseMatrix intdof_agg(I, J, D, tdof_agg.NumRows(), tdof_agg.NumCols());
129  return Transpose(intdof_agg);
130 }
131 
132 void DFSSpaces::MakeDofRelationTables(int level)
133 {
134  Array<HYPRE_BigInt> agg_starts(Array<HYPRE_BigInt>(l2_0_fes_->GetDofOffsets(),
135  2));
136  auto& elem_agg = (const SparseMatrix&)*l2_0_fes_->GetUpdateOperator();
137  OperatorPtr agg_elem(Transpose(elem_agg));
138  SparseMatrix& agg_el = *agg_elem.As<SparseMatrix>();
139 
140  el_l2dof_.push_back(ElemToDof(*l2_fes_));
141  data_.agg_l2dof[level].Reset(Mult(agg_el, el_l2dof_[level+1]));
142 
143  Array<int> bdr_tdofs;
144  hdiv_fes_->GetEssentialTrueDofs(all_bdr_attr_, bdr_tdofs);
145  auto tmp = AggToInteriorDof(bdr_tdofs, agg_el, ElemToDof(*hdiv_fes_),
146  *hdiv_fes_->Dof_TrueDof_Matrix(), agg_starts);
147  data_.agg_hdivdof[level].Reset(tmp);
148 }
149 
151 {
152  auto GetP = [this](OperatorPtr& P, unique_ptr<ParFiniteElementSpace>& cfes,
153  ParFiniteElementSpace& fes, bool remove_zero)
154  {
155  fes.Update();
156  fes.GetTrueTransferOperator(*cfes, P);
157  if (remove_zero)
158  {
159  P.As<HypreParMatrix>()->DropSmallEntries(1e-16);
160  }
161  (level_ < (int)data_.P_l2.size()-1) ? cfes->Update() : cfes.reset();
162  };
163 
164  GetP(data_.P_hdiv[level_], coarse_hdiv_fes_, *hdiv_fes_, true);
165  GetP(data_.P_l2[level_], coarse_l2_fes_, *l2_fes_, false);
166  MakeDofRelationTables(level_);
167 
168  GetP(data_.P_hcurl[level_], coarse_hcurl_fes_, *hcurl_fes_, true);
169 
170  Vector trash1(hcurl_fes_->GetVSize()), trash2(hdiv_fes_->GetVSize());
171  ParDiscreteLinearOperator curl(hcurl_fes_.get(), hdiv_fes_.get());
173  curl.Assemble();
174  curl.EliminateTrialDofs(ess_bdr_attr_, trash1, trash2);
175  curl.Finalize();
176  data_.C[level_+1].Reset(curl.ParallelAssemble());
177 
178  ++level_;
179 
180  if (level_ == (int)data_.P_l2.size()) { DataFinalize(); }
181 }
182 
183 void DFSSpaces::DataFinalize()
184 {
185  ParBilinearForm mass(l2_fes_.get());
187  mass.Assemble();
188  mass.Finalize();
189  OperatorPtr W(mass.LoseMat());
190 
191  SparseMatrix P_l2;
192  for (int l = (int)data_.P_l2.size()-1; l >= 0; --l)
193  {
194  data_.P_l2[l].As<HypreParMatrix>()->GetDiag(P_l2);
195  OperatorPtr PT_l2(Transpose(P_l2));
196  auto PTW = Mult(*PT_l2.As<SparseMatrix>(), *W.As<SparseMatrix>());
197  auto cW = Mult(*PTW, P_l2);
198  auto cW_inv = new SymDirectSubBlockSolver(*cW, el_l2dof_[l]);
199  data_.Q_l2[l].Reset(new ProductOperator(cW_inv, PTW, true, true));
200  W.Reset(cW);
201  }
202 
203  l2_0_fes_.reset();
204 }
205 
207  : Solver(B.NumRows()), BBT_solver_(B.GetComm())
208 {
209  OperatorPtr BT(B.Transpose());
210  BBT_.Reset(ParMult(&B, BT.As<HypreParMatrix>()));
211  BBT_.As<HypreParMatrix>()->CopyColStarts();
212 
213  BBT_prec_.Reset(new HypreBoomerAMG(*BBT_.As<HypreParMatrix>()));
214  BBT_prec_.As<HypreBoomerAMG>()->SetPrintLevel(0);
215 
216  SetOptions(BBT_solver_, param);
217  BBT_solver_.SetOperator(*BBT_);
218  BBT_solver_.SetPreconditioner(*BBT_prec_.As<HypreBoomerAMG>());
219 }
220 
222  : Solver(M.NumRows()+B.NumRows()), local_system_(height), offset_(M.NumRows())
223 {
224  local_system_.CopyMN(M, 0, 0);
225  local_system_.CopyMN(B, offset_, 0);
226  local_system_.CopyMNt(B, 0, offset_);
227 
228  local_system_.SetRow(offset_, 0.0);
229  local_system_.SetCol(offset_, 0.0);
230  local_system_(offset_, offset_) = -1.0;
231  local_solver_.SetOperator(local_system_);
232 }
233 
234 void LocalSolver::Mult(const Vector &x, Vector &y) const
235 {
236  const double x0 = x[offset_];
237  const_cast<Vector&>(x)[offset_] = 0.0;
238 
239  y.SetSize(local_system_.NumRows());
240  local_solver_.Mult(x, y);
241 
242  const_cast<Vector&>(x)[offset_] = x0;
243 }
244 
246  const HypreParMatrix& B,
247  const SparseMatrix& agg_hdivdof,
248  const SparseMatrix& agg_l2dof,
249  const HypreParMatrix& P_l2,
250  const HypreParMatrix& Q_l2)
251  : Solver(M.NumRows() + B.NumRows()), agg_hdivdof_(agg_hdivdof),
252  agg_l2dof_(agg_l2dof), solvers_loc_(agg_l2dof.NumRows())
253 {
254  coarse_l2_projector_.Reset(new ProductOperator(&P_l2, &Q_l2, false, false));
255 
256  offsets_loc_.SetSize(3, 0);
257  offsets_.SetSize(3, 0);
258  offsets_[1] = M.NumRows();
259  offsets_[2] = M.NumRows() + B.NumRows();
260 
261  SparseMatrix M_diag, B_diag;
262  M.GetDiag(M_diag);
263  B.GetDiag(B_diag);
264 
265  DenseMatrix B_loc, M_loc;
266 
267  for (int agg = 0; agg < (int)solvers_loc_.size(); agg++)
268  {
269  GetRowColumnsRef(agg_hdivdof_, agg, hdivdofs_loc_);
270  GetRowColumnsRef(agg_l2dof_, agg, l2dofs_loc_);
271  M_loc.SetSize(hdivdofs_loc_.Size(), hdivdofs_loc_.Size());
272  B_loc.SetSize(l2dofs_loc_.Size(), hdivdofs_loc_.Size());
273  M_diag.GetSubMatrix(hdivdofs_loc_, hdivdofs_loc_, M_loc);
274  B_diag.GetSubMatrix(l2dofs_loc_, hdivdofs_loc_, B_loc);
275  solvers_loc_[agg].Reset(new LocalSolver(M_loc, B_loc));
276  }
277 }
278 
279 void SaddleSchwarzSmoother::Mult(const Vector & x, Vector & y) const
280 {
281  y.SetSize(offsets_[2]);
282  y = 0.0;
283 
284  BlockVector blk_y(y.GetData(), offsets_);
285  BlockVector Pi_x(offsets_); // aggregate-wise average free projection of x
286  static_cast<Vector&>(Pi_x) = x;
287 
288  // Right hand side: F_l = F - W_l P_l2[l] (W_{l+1})^{-1} P_l2[l]^T F
289  // This ensures the existence of solutions to the local problems
290  Vector coarse_l2_projection(Pi_x.BlockSize(1));
291  coarse_l2_projector_->MultTranspose(Pi_x.GetBlock(1), coarse_l2_projection);
292 
293  Pi_x.GetBlock(1) -= coarse_l2_projection;
294 
295  for (int agg = 0; agg < (int)solvers_loc_.size(); agg++)
296  {
297  GetRowColumnsRef(agg_hdivdof_, agg, hdivdofs_loc_);
298  GetRowColumnsRef(agg_l2dof_, agg, l2dofs_loc_);
299 
300  offsets_loc_[1] = hdivdofs_loc_.Size();
301  offsets_loc_[2] = offsets_loc_[1]+l2dofs_loc_.Size();
302 
303  BlockVector rhs_loc(offsets_loc_), sol_loc(offsets_loc_);
304  Pi_x.GetBlock(0).GetSubVector(hdivdofs_loc_, rhs_loc.GetBlock(0));
305  Pi_x.GetBlock(1).GetSubVector(l2dofs_loc_, rhs_loc.GetBlock(1));
306 
307  solvers_loc_[agg]->Mult(rhs_loc, sol_loc);
308 
309  blk_y.GetBlock(0).AddElementVector(hdivdofs_loc_, sol_loc.GetBlock(0));
310  blk_y.GetBlock(1).AddElementVector(l2dofs_loc_, sol_loc.GetBlock(1));
311  }
312 
313  coarse_l2_projector_->Mult(blk_y.GetBlock(1), coarse_l2_projection);
314  blk_y.GetBlock(1) -= coarse_l2_projection;
315 }
316 
318  IterSolveParameters param)
319  : DarcySolver(M.NumRows(), B.NumRows()), op_(offsets_), prec_(offsets_),
320  BT_(B.Transpose()), solver_(M.GetComm())
321 {
322  op_.SetBlock(0,0, &M);
323  op_.SetBlock(0,1, BT_.As<HypreParMatrix>());
324  op_.SetBlock(1,0, &B);
325 
326  Vector Md;
327  M.GetDiag(Md);
328  BT_.As<HypreParMatrix>()->InvScaleRows(Md);
329  S_.Reset(ParMult(&B, BT_.As<HypreParMatrix>()));
330  BT_.As<HypreParMatrix>()->ScaleRows(Md);
331 
332  prec_.SetDiagonalBlock(0, new HypreDiagScale(M));
333  prec_.SetDiagonalBlock(1, new HypreBoomerAMG(*S_.As<HypreParMatrix>()));
334  static_cast<HypreBoomerAMG&>(prec_.GetDiagonalBlock(1)).SetPrintLevel(0);
335  prec_.owns_blocks = true;
336 
337  SetOptions(solver_, param);
338  solver_.SetOperator(op_);
339  solver_.SetPreconditioner(prec_);
340 }
341 
342 void BDPMinresSolver::Mult(const Vector & x, Vector & y) const
343 {
344  solver_.Mult(x, y);
345  for (int dof : ess_zero_dofs_) { y[dof] = 0.0; }
346 }
347 
349  const DFSData& data)
350  : DarcySolver(M.NumRows(), B.NumRows()), data_(data), param_(data.param),
351  BT_(B.Transpose()), BBT_solver_(B, param_.BBT_solve_param),
352  ops_offsets_(data.P_l2.size()+1), ops_(ops_offsets_.size()),
353  blk_Ps_(ops_.Size()-1), smoothers_(ops_.Size())
354 {
355  ops_offsets_.back().MakeRef(DarcySolver::offsets_);
356  ops_.Last() = new BlockOperator(ops_offsets_.back());
357  ops_.Last()->SetBlock(0, 0, const_cast<HypreParMatrix*>(&M));
358  ops_.Last()->SetBlock(1, 0, const_cast<HypreParMatrix*>(&B));
359  ops_.Last()->SetBlock(0, 1, BT_.Ptr());
360 
361  for (int l = data.P_l2.size(); l >= 0; --l)
362  {
363  auto& M_f = static_cast<HypreParMatrix&>(ops_[l]->GetBlock(0, 0));
364  auto& B_f = static_cast<HypreParMatrix&>(ops_[l]->GetBlock(1, 0));
365 
366  if (l == 0)
367  {
368  SparseMatrix M_f_diag, B_f_diag;
369  M_f.GetDiag(M_f_diag);
370  B_f.GetDiag(B_f_diag);
371  for (int dof : data.coarsest_ess_hdivdofs)
372  {
373  M_f_diag.EliminateRowCol(dof);
374  B_f_diag.EliminateCol(dof);
375  }
376 
377  const IterSolveParameters& param = param_.coarse_solve_param;
378  auto coarse_solver = new BDPMinresSolver(M_f, B_f, param);
379  if (ops_.Size() > 1)
380  {
381  coarse_solver->SetEssZeroDofs(data.coarsest_ess_hdivdofs);
382  }
383  smoothers_[l] = coarse_solver;
384  continue;
385  }
386 
387  HypreParMatrix& P_hdiv_l = *data.P_hdiv[l-1].As<HypreParMatrix>();
388  HypreParMatrix& P_l2_l = *data.P_l2[l-1].As<HypreParMatrix>();
389  SparseMatrix& agg_hdivdof_l = *data.agg_hdivdof[l-1].As<SparseMatrix>();
390  SparseMatrix& agg_l2dof_l = *data.agg_l2dof[l-1].As<SparseMatrix>();
391  HypreParMatrix& Q_l2_l = *data.Q_l2[l-1].As<HypreParMatrix>();
392  HypreParMatrix* C_l = data.C[l].As<HypreParMatrix>();
393 
394  auto S0 = new SaddleSchwarzSmoother(M_f, B_f, agg_hdivdof_l,
395  agg_l2dof_l, P_l2_l, Q_l2_l);
396  if (param_.coupled_solve)
397  {
398  auto S1 = new BlockDiagonalPreconditioner(ops_offsets_[l]);
399  S1->SetDiagonalBlock(0, new AuxSpaceSmoother(M_f, C_l));
400  S1->owns_blocks = true;
401  smoothers_[l] = new ProductSolver(ops_[l], S0, S1, false, true, true);
402  }
403  else
404  {
405  smoothers_[l] = S0;
406  }
407 
408  HypreParMatrix* M_c = TwoStepsRAP(P_hdiv_l, M_f, P_hdiv_l);
409  HypreParMatrix* B_c = TwoStepsRAP(P_l2_l, B_f, P_hdiv_l);
410 
411  ops_offsets_[l-1].SetSize(3, 0);
412  ops_offsets_[l-1][1] = M_c->NumRows();
413  ops_offsets_[l-1][2] = M_c->NumRows() + B_c->NumRows();
414 
415  blk_Ps_[l-1] = new BlockOperator(ops_offsets_[l], ops_offsets_[l-1]);
416  blk_Ps_[l-1]->SetBlock(0, 0, &P_hdiv_l);
417  blk_Ps_[l-1]->SetBlock(1, 1, &P_l2_l);
418 
419  ops_[l-1] = new BlockOperator(ops_offsets_[l-1]);
420  ops_[l-1]->SetBlock(0, 0, M_c);
421  ops_[l-1]->SetBlock(1, 0, B_c);
422  ops_[l-1]->SetBlock(0, 1, B_c->Transpose());
423  ops_[l-1]->owns_blocks = true;
424  }
425 
426  Array<bool> own_ops(ops_.Size());
427  Array<bool> own_smoothers(smoothers_.Size());
428  Array<bool> own_Ps(blk_Ps_.Size());
429  own_ops = true;
430  own_smoothers = true;
431  own_Ps = true;
432 
433  if (data_.P_l2.size() == 0) { return; }
434 
435  if (param_.coupled_solve)
436  {
437  solver_.Reset(new GMRESSolver(B.GetComm()));
438  solver_.As<GMRESSolver>()->SetOperator(*(ops_.Last()));
439  prec_.Reset(new Multigrid(ops_, smoothers_, blk_Ps_,
440  own_ops, own_smoothers, own_Ps));
441  }
442  else
443  {
444  Array<HypreParMatrix*> ops(data_.P_hcurl.size()+1);
445  Array<Solver*> smoothers(ops.Size());
446  Array<HypreParMatrix*> Ps(data_.P_hcurl.size());
447  own_Ps = false;
448 
449  HypreParMatrix& C_finest = *data.C.back().As<HypreParMatrix>();
450  ops.Last() = TwoStepsRAP(C_finest, M, C_finest);
451  ops.Last()->EliminateZeroRows();
452  ops.Last()->DropSmallEntries(1e-14);
453 
454  solver_.Reset(new CGSolver(B.GetComm()));
455  solver_.As<CGSolver>()->SetOperator(*ops.Last());
456  smoothers.Last() = new HypreSmoother(*ops.Last());
457  static_cast<HypreSmoother*>(smoothers.Last())->SetOperatorSymmetry(true);
458 
459  for (int l = Ps.Size()-1; l >= 0; --l)
460  {
461  Ps[l] = data_.P_hcurl[l].As<HypreParMatrix>();
462  ops[l] = TwoStepsRAP(*Ps[l], *ops[l+1], *Ps[l]);
463  ops[l]->DropSmallEntries(1e-14);
464  smoothers[l] = new HypreSmoother(*ops[l]);
465  static_cast<HypreSmoother*>(smoothers[l])->SetOperatorSymmetry(true);
466  }
467 
468  prec_.Reset(new Multigrid(ops, smoothers, Ps, own_ops, own_smoothers, own_Ps));
469  }
470 
471  solver_.As<IterativeSolver>()->SetPreconditioner(*prec_.As<Solver>());
472  SetOptions(*solver_.As<IterativeSolver>(), param_);
473 }
474 
476 {
477  if (param_.coupled_solve) { return; }
478  for (int i = 0; i < ops_.Size(); ++i)
479  {
480  delete ops_[i];
481  delete smoothers_[i];
482  if (i == ops_.Size() - 1) { break; }
483  delete blk_Ps_[i];
484  }
485 }
486 
487 void DivFreeSolver::SolveParticular(const Vector& rhs, Vector& sol) const
488 {
489  std::vector<Vector> rhss(smoothers_.Size());
490  std::vector<Vector> sols(smoothers_.Size());
491 
492  rhss.back().SetDataAndSize(const_cast<double*>(rhs.HostRead()), rhs.Size());
493  sols.back().SetDataAndSize(sol.HostWrite(), sol.Size());
494 
495  for (int l = blk_Ps_.Size()-1; l >= 0; --l)
496  {
497  rhss[l].SetSize(blk_Ps_[l]->NumCols());
498  sols[l].SetSize(blk_Ps_[l]->NumCols());
499 
500  sols[l] = 0.0;
501  rhss[l] = 0.0;
502 
503  blk_Ps_[l]->MultTranspose(rhss[l+1], rhss[l]);
504  }
505 
506  for (int l = 0; l < smoothers_.Size(); ++l)
507  {
508  smoothers_[l]->Mult(rhss[l], sols[l]);
509  }
510 
511  for (int l = 0; l < blk_Ps_.Size(); ++l)
512  {
513  Vector P_sol(blk_Ps_[l]->NumRows());
514  blk_Ps_[l]->Mult(sols[l], P_sol);
515  sols[l+1] += P_sol;
516  }
517 }
518 
519 void DivFreeSolver::SolveDivFree(const Vector &rhs, Vector& sol) const
520 {
521  Vector rhs_divfree(data_.C.back()->NumCols());
522  data_.C.back()->MultTranspose(rhs, rhs_divfree);
523 
524  Vector potential_divfree(rhs_divfree.Size());
525  potential_divfree = 0.0;
526  solver_->Mult(rhs_divfree, potential_divfree);
527 
528  data_.C.back()->Mult(potential_divfree, sol);
529 }
530 
531 void DivFreeSolver::SolvePotential(const Vector& rhs, Vector& sol) const
532 {
533  Vector rhs_p(BT_->NumCols());
534  BT_->MultTranspose(rhs, rhs_p);
535  BBT_solver_.Mult(rhs_p, sol);
536 }
537 
538 void DivFreeSolver::Mult(const Vector & x, Vector & y) const
539 {
540  MFEM_VERIFY(x.Size() == offsets_[2], "MLDivFreeSolver: x size is invalid");
541  MFEM_VERIFY(y.Size() == offsets_[2], "MLDivFreeSolver: y size is invalid");
542 
543  if (ops_.Size() == 1) { smoothers_[0]->Mult(x, y); return; }
544 
545  BlockVector blk_y(y, offsets_);
546 
547  BlockVector resid(offsets_);
548  ops_.Last()->Mult(y, resid);
549  add(1.0, x, -1.0, resid, resid);
550 
551  BlockVector correction(offsets_);
552  correction = 0.0;
553 
554  if (param_.coupled_solve)
555  {
556  solver_->Mult(resid, correction);
557  y += correction;
558  }
559  else
560  {
561  StopWatch ch;
562  ch.Start();
563 
564  SolveParticular(resid, correction);
565  blk_y += correction;
566 
567  if (param_.verbose)
568  {
569  cout << "Particular solution found in " << ch.RealTime() << "s.\n";
570  }
571 
572  ch.Clear();
573  ch.Start();
574 
575  ops_.Last()->Mult(y, resid);
576  add(1.0, x, -1.0, resid, resid);
577 
578  SolveDivFree(resid.GetBlock(0), correction.GetBlock(0));
579  blk_y.GetBlock(0) += correction.GetBlock(0);
580 
581  if (param_.verbose)
582  {
583  cout << "Divergence free solution found in " << ch.RealTime() << "s.\n";
584  }
585 
586  ch.Clear();
587  ch.Start();
588 
589  auto M = dynamic_cast<HypreParMatrix&>(ops_.Last()->GetBlock(0, 0));
590  M.Mult(-1.0, correction.GetBlock(0), 1.0, resid.GetBlock(0));
591  SolvePotential(resid.GetBlock(0), correction.GetBlock(1));
592  blk_y.GetBlock(1) += correction.GetBlock(1);
593 
594  if (param_.verbose)
595  {
596  cout << "Scalar potential found in " << ch.RealTime() << "s.\n";
597  }
598  }
599 }
600 
602 {
603  if (ops_.Size() == 1)
604  {
605  return static_cast<BDPMinresSolver*>(smoothers_[0])->GetNumIterations();
606  }
607  return solver_.As<IterativeSolver>()->GetNumIterations();
608 }
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: operator.hpp:93
Multigrid solver class.
Definition: multigrid.hpp:135
Conjugate gradient method.
Definition: solvers.hpp:493
int * GetJ()
Definition: table.hpp:114
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
void EliminateCol(int col, DiagonalPolicy dpolicy=DIAG_ZERO)
Eliminates the column col from the matrix.
Definition: sparsemat.cpp:1741
void SetCol(int c, const double *col)
Definition: densemat.cpp:2193
void GetDiag(Vector &d) const
Returns the Diagonal of A.
Definition: sparsemat.cpp:669
A class to handle Vectors in a block fashion.
Definition: blockvector.hpp:30
void SetRow(int r, const double *row)
Definition: densemat.cpp:2178
DivFreeSolver(const HypreParMatrix &M, const HypreParMatrix &B, const DFSData &data)
int Dimension() const
Definition: mesh.hpp:1047
void Mult(const double *x, double *y) const
Matrix vector multiplication with the inverse of dense matrix.
Definition: densemat.cpp:3963
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:512
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:452
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
Definition: sparsemat.cpp:391
int RowSize(const int i) const
Returns the number of elements in row i.
Definition: sparsemat.cpp:344
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
std::vector< OperatorPtr > Q_l2
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1612
Abstract solver class for Darcy&#39;s flow.
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual double * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:460
int NumRows() const
Get the number of rows (size of output) of the Operator. Synonym with Height().
Definition: operator.hpp:69
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
const Element * GetElement(int i) const
Definition: mesh.hpp:1081
double RealTime()
Definition: tic_toc.cpp:426
Abstract parallel finite element space.
Definition: pfespace.hpp:28
SaddleSchwarzSmoother(const HypreParMatrix &M, const HypreParMatrix &B, const SparseMatrix &agg_hdivdof, const SparseMatrix &agg_l2dof, const HypreParMatrix &P_l2, const HypreParMatrix &Q_l2)
STL namespace.
std::vector< OperatorPtr > P_hcurl
void GetSubMatrix(const Array< int > &rows, const Array< int > &cols, DenseMatrix &subm) const
Definition: sparsemat.cpp:2886
HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B, bool own_matrix)
Definition: hypre.cpp:2850
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:314
void AddDomainInterpolator(DiscreteInterpolator *di)
Adds a domain interpolator. Assumes ownership of di.
IterSolveParameters coarse_solve_param
LocalSolver(const DenseMatrix &M, const DenseMatrix &B)
Operator & GetDiagonalBlock(int iblock)
Return a reference to block i,i.
ID for class SparseMatrix.
Definition: operator.hpp:286
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1590
Data for the divergence free solver.
General product operator: x -> (A*B)(x) = A(B(x)).
Definition: operator.hpp:774
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.hpp:616
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
A class to handle Block diagonal preconditioners in a matrix-free implementation. ...
Data type sparse matrix.
Definition: sparsemat.hpp:50
Wrapper for the block-diagonal-preconditioned MINRES defined in ex5p.cpp.
Jacobi preconditioner in hypre.
Definition: hypre.hpp:1421
SparseMatrix ElemToDof(const ParFiniteElementSpace &fes)
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
Timing object.
Definition: tic_toc.hpp:34
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
std::vector< OperatorPtr > P_hdiv
HypreParMatrix * LeftDiagMult(const SparseMatrix &D, HYPRE_BigInt *row_starts=NULL) const
Multiply the HypreParMatrix on the left by a block-diagonal parallel matrix D and return the result a...
Definition: hypre.cpp:1912
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Parallel smoothers in hypre.
Definition: hypre.hpp:971
virtual void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: densemat.cpp:3956
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1591
void GetDiag(Vector &diag) const
Get the local diagonal of the matrix.
Definition: hypre.cpp:1481
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:614
void Start()
Start the stopwatch. The elapsed time is not cleared.
Definition: tic_toc.cpp:411
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:208
BDPMinresSolver(HypreParMatrix &M, HypreParMatrix &B, IterSolveParameters param)
MPI_Comm GetComm() const
MPI communicator.
Definition: hypre.hpp:534
void SetAbsTol(double atol)
Definition: solvers.hpp:200
void SetRelTol(double rtol)
Definition: solvers.hpp:199
Abstract base class for iterative solver.
Definition: solvers.hpp:66
OperatorHandle OperatorPtr
Add an alternative name for OperatorHandle – OperatorPtr.
Definition: handle.hpp:212
SparseMatrix * AggToInteriorDof(const Array< int > &bdr_truedofs, const SparseMatrix &agg_elem, const SparseMatrix &elem_dof, const HypreParMatrix &dof_truedof, Array< HYPRE_BigInt > &agg_starts)
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:684
void EliminateRowCol(int rc, const double sol, Vector &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate row rc and column rc and modify the rhs using sol.
Definition: sparsemat.cpp:1854
HYPRE_Int HYPRE_BigInt
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition: table.cpp:413
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
std::vector< OperatorPtr > agg_hdivdof
std::vector< OperatorPtr > agg_l2dof
GMRES method.
Definition: solvers.hpp:525
Parameters for the divergence free solver.
int NumCols() const
Get the number of columns (size of input) of the Operator. Synonym with Width().
Definition: operator.hpp:75
HypreParMatrix * Transpose() const
Returns the transpose of *this.
Definition: hypre.cpp:1608
std::vector< OperatorPtr > C
void CopyMNt(const DenseMatrix &A, int row_offset, int col_offset)
Copy matrix A^t to the location in *this at row_offset, col_offset.
Definition: densemat.cpp:1618
void SetOptions(IterativeSolver &solver, const IterSolveParameters &param)
HypreParMatrix * TwoStepsRAP(const HypreParMatrix &Rt, const HypreParMatrix &A, const HypreParMatrix &P)
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:104
Block diagonal solver for symmetric A, each block is inverted by direct solver.
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:587
void GetRowColumnsRef(const SparseMatrix &A, int row, Array< int > &cols)
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
static void ListToMarker(const Array< int > &list, int marker_size, Array< int > &marker, int mark_val=-1)
Convert an array of indices (list) to a Boolean marker array where all indices in the list are marked...
Definition: fespace.cpp:617
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:869
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:507
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:447
Vector data type.
Definition: vector.hpp:60
ID for class HypreParMatrix.
Definition: operator.hpp:287
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
Parameters for iterative solver.
void DropSmallEntries(double tol)
Wrapper for hypre_ParCSRMatrixDropSmallEntries in different versions of hypre. Drop off-diagonal entr...
Definition: hypre.cpp:2227
int * GetI()
Definition: table.hpp:113
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:252
void CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco)
Copy the m x n submatrix of A at row/col offsets Aro/Aco to *this.
Definition: densemat.cpp:1592
Base class for solvers.
Definition: operator.hpp:682
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
Solver for local problems in SaddleSchwarzSmoother.
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
A class to handle Block systems in a matrix-free implementation.
BBTSolver(const HypreParMatrix &B, IterSolveParameters param)
Class for parallel meshes.
Definition: pmesh.hpp:32
const Table & GetElementToDofTable() const
Return a reference to the internal Table that stores the lists of scalar dofs, for each mesh element...
Definition: fespace.hpp:750
virtual Type GetType() const =0
Returns element&#39;s type.
std::vector< OperatorPtr > P_l2
static void AdjustVDofs(Array< int > &vdofs)
Definition: fespace.cpp:267
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0) const
Computes y = alpha * A * x + beta * y.
Definition: hypre.cpp:1733
void Clear()
Clear the elapsed time on the stopwatch and restart it if it&#39;s running.
Definition: tic_toc.cpp:406
void SetBlock(int iRow, int iCol, Operator *op, double c=1.0)
Add a block op in the block-entry (iblock, jblock).
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Definition: handle.hpp:145
void SetDiagonalBlock(int iblock, Operator *op)
Add a square block op in the block-entry (iblock, iblock).
Vector & GetBlock(int i)
Get the i-th vector in the block.
Definition: blockvector.hpp:90