MFEM  v4.6.0
Finite element discretization library
navier_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 "navier_solver.hpp"
13 #include "../../general/forall.hpp"
14 #include <fstream>
15 #include <iomanip>
16 
17 using namespace mfem;
18 using namespace navier;
19 
21 {
23  for (int i = 0; i < bffis->Size(); ++i)
24  {
25  dst->AddDomainIntegrator((*bffis)[i]);
26  }
27 }
28 
29 NavierSolver::NavierSolver(ParMesh *mesh, int order, double kin_vis)
30  : pmesh(mesh), order(order), kin_vis(kin_vis),
31  gll_rules(0, Quadrature1D::GaussLobatto)
32 {
34  pfec = new H1_FECollection(order);
37 
38  // Check if fully periodic mesh
39  if (!(pmesh->bdr_attributes.Size() == 0))
40  {
42  vel_ess_attr = 0;
43 
45  pres_ess_attr = 0;
46  }
47 
48  int vfes_truevsize = vfes->GetTrueVSize();
49  int pfes_truevsize = pfes->GetTrueVSize();
50 
51  un.SetSize(vfes_truevsize);
52  un = 0.0;
53  un_next.SetSize(vfes_truevsize);
54  un_next = 0.0;
55  unm1.SetSize(vfes_truevsize);
56  unm1 = 0.0;
57  unm2.SetSize(vfes_truevsize);
58  unm2 = 0.0;
59  fn.SetSize(vfes_truevsize);
60  Nun.SetSize(vfes_truevsize);
61  Nun = 0.0;
62  Nunm1.SetSize(vfes_truevsize);
63  Nunm1 = 0.0;
64  Nunm2.SetSize(vfes_truevsize);
65  Nunm2 = 0.0;
66  Fext.SetSize(vfes_truevsize);
67  FText.SetSize(vfes_truevsize);
68  Lext.SetSize(vfes_truevsize);
69  resu.SetSize(vfes_truevsize);
70 
71  tmp1.SetSize(vfes_truevsize);
72 
73  pn.SetSize(pfes_truevsize);
74  pn = 0.0;
75  resp.SetSize(pfes_truevsize);
76  resp = 0.0;
77  FText_bdr.SetSize(pfes_truevsize);
78  g_bdr.SetSize(pfes_truevsize);
79 
81  un_gf = 0.0;
83  un_next_gf = 0.0;
84 
90 
92  pn_gf = 0.0;
94 
95  cur_step = 0;
96 
97  PrintInfo();
98 }
99 
100 void NavierSolver::Setup(double dt)
101 {
102  if (verbose && pmesh->GetMyRank() == 0)
103  {
104  mfem::out << "Setup" << std::endl;
105  if (partial_assembly)
106  {
107  mfem::out << "Using Partial Assembly" << std::endl;
108  }
109  else
110  {
111  mfem::out << "Using Full Assembly" << std::endl;
112  }
113  }
114 
115  sw_setup.Start();
116 
119 
120  Array<int> empty;
121 
122  // GLL integration rule (Numerical Integration)
123  const IntegrationRule &ir_ni = gll_rules.Get(vfes->GetFE(0)->GetGeomType(),
124  2 * order - 1);
125 
126  nlcoeff.constant = -1.0;
127  N = new ParNonlinearForm(vfes);
128  auto *nlc_nlfi = new VectorConvectionNLFIntegrator(nlcoeff);
129  if (numerical_integ)
130  {
131  nlc_nlfi->SetIntRule(&ir_ni);
132  }
133  N->AddDomainIntegrator(nlc_nlfi);
134  if (partial_assembly)
135  {
137  N->Setup();
138  }
139 
141  auto *mv_blfi = new VectorMassIntegrator;
142  if (numerical_integ)
143  {
144  mv_blfi->SetIntRule(&ir_ni);
145  }
146  Mv_form->AddDomainIntegrator(mv_blfi);
147  if (partial_assembly)
148  {
150  }
151  Mv_form->Assemble();
152  Mv_form->FormSystemMatrix(empty, Mv);
153 
155  auto *sp_blfi = new DiffusionIntegrator;
156  if (numerical_integ)
157  {
158  sp_blfi->SetIntRule(&ir_ni);
159  }
160  Sp_form->AddDomainIntegrator(sp_blfi);
161  if (partial_assembly)
162  {
164  }
165  Sp_form->Assemble();
167 
169  auto *vd_mblfi = new VectorDivergenceIntegrator();
170  if (numerical_integ)
171  {
172  vd_mblfi->SetIntRule(&ir_ni);
173  }
174  D_form->AddDomainIntegrator(vd_mblfi);
175  if (partial_assembly)
176  {
178  }
179  D_form->Assemble();
180  D_form->FormRectangularSystemMatrix(empty, empty, D);
181 
183  auto *g_mblfi = new GradientIntegrator();
184  if (numerical_integ)
185  {
186  g_mblfi->SetIntRule(&ir_ni);
187  }
188  G_form->AddDomainIntegrator(g_mblfi);
189  if (partial_assembly)
190  {
192  }
193  G_form->Assemble();
194  G_form->FormRectangularSystemMatrix(empty, empty, G);
195 
197  H_bdfcoeff.constant = 1.0 / dt;
198  H_form = new ParBilinearForm(vfes);
199  auto *hmv_blfi = new VectorMassIntegrator(H_bdfcoeff);
200  auto *hdv_blfi = new VectorDiffusionIntegrator(H_lincoeff);
201  if (numerical_integ)
202  {
203  hmv_blfi->SetIntRule(&ir_ni);
204  hdv_blfi->SetIntRule(&ir_ni);
205  }
206  H_form->AddDomainIntegrator(hmv_blfi);
207  H_form->AddDomainIntegrator(hdv_blfi);
208  if (partial_assembly)
209  {
211  }
212  H_form->Assemble();
214 
217  auto *ftext_bnlfi = new BoundaryNormalLFIntegrator(*FText_gfcoeff);
218  if (numerical_integ)
219  {
220  ftext_bnlfi->SetIntRule(&ir_ni);
221  }
223 
225  for (auto &vel_dbc : vel_dbcs)
226  {
227  auto *gbdr_bnlfi = new BoundaryNormalLFIntegrator(*vel_dbc.coeff);
228  if (numerical_integ)
229  {
230  gbdr_bnlfi->SetIntRule(&ir_ni);
231  }
232  g_bdr_form->AddBoundaryIntegrator(gbdr_bnlfi, vel_dbc.attr);
233  }
234 
235  f_form = new ParLinearForm(vfes);
236  for (auto &accel_term : accel_terms)
237  {
238  auto *vdlfi = new VectorDomainLFIntegrator(*accel_term.coeff);
239  // @TODO: This order should always be the same as the nonlinear forms one!
240  // const IntegrationRule &ir = IntRules.Get(vfes->GetFE(0)->GetGeomType(),
241  // 4 * order);
242  // vdlfi->SetIntRule(&ir);
243  if (numerical_integ)
244  {
245  vdlfi->SetIntRule(&ir_ni);
246  }
247  f_form->AddDomainIntegrator(vdlfi);
248  }
249 
250  if (partial_assembly)
251  {
252  Vector diag_pa(vfes->GetTrueVSize());
253  Mv_form->AssembleDiagonal(diag_pa);
254  MvInvPC = new OperatorJacobiSmoother(diag_pa, empty);
255  }
256  else
257  {
259  dynamic_cast<HypreSmoother *>(MvInvPC)->SetType(HypreSmoother::Jacobi, 1);
260  }
261  MvInv = new CGSolver(vfes->GetComm());
262  MvInv->iterative_mode = false;
263  MvInv->SetOperator(*Mv);
266  MvInv->SetRelTol(1e-12);
267  MvInv->SetMaxIter(200);
268 
269  if (partial_assembly)
270  {
274  SpInvPC->Mult(resp, pn);
277  }
278  else
279  {
284  }
285  SpInv = new CGSolver(vfes->GetComm());
286  SpInv->iterative_mode = true;
287  SpInv->SetOperator(*Sp);
288  if (pres_dbcs.empty())
289  {
291  }
292  else
293  {
295  }
298  SpInv->SetMaxIter(200);
299 
300  if (partial_assembly)
301  {
302  Vector diag_pa(vfes->GetTrueVSize());
303  H_form->AssembleDiagonal(diag_pa);
305  }
306  else
307  {
309  dynamic_cast<HypreSmoother *>(HInvPC)->SetType(HypreSmoother::Jacobi, 1);
310  }
311  HInv = new CGSolver(vfes->GetComm());
312  HInv->iterative_mode = true;
313  HInv->SetOperator(*H);
317  HInv->SetMaxIter(200);
318 
319  // If the initial condition was set, it has to be aligned with dependent
320  // Vectors and GridFunctions
322  un_next = un;
324 
325  // Set initial time step in the history array
326  dthist[0] = dt;
327 
328  if (filter_alpha != 0.0)
329  {
331  pmesh->Dimension());
333  vfec_filter,
334  pmesh->Dimension());
335 
337  un_NM1_gf = 0.0;
338 
340  un_filtered_gf = 0.0;
341  }
342 
343  sw_setup.Stop();
344 }
345 
347 {
348  // Rotate values in time step history
349  dthist[2] = dthist[1];
350  dthist[1] = dthist[0];
351  dthist[0] = dt;
352 
353  // Rotate values in nonlinear extrapolation history
354  Nunm2 = Nunm1;
355  Nunm1 = Nun;
356 
357  // Rotate values in solution history
358  unm2 = unm1;
359  unm1 = un;
360 
361  // Update the current solution and corresponding GridFunction
363  un = un_next;
365 }
366 
367 void NavierSolver::Step(double &time, double dt, int current_step,
368  bool provisional)
369 {
370  sw_step.Start();
371 
372  SetTimeIntegrationCoefficients(current_step);
373 
374  // Set current time for velocity Dirichlet boundary conditions.
375  for (auto &vel_dbc : vel_dbcs)
376  {
377  vel_dbc.coeff->SetTime(time + dt);
378  }
379 
380  // Set current time for pressure Dirichlet boundary conditions.
381  for (auto &pres_dbc : pres_dbcs)
382  {
383  pres_dbc.coeff->SetTime(time + dt);
384  }
385 
386  H_bdfcoeff.constant = bd0 / dt;
387  H_form->Update();
388  H_form->Assemble();
390 
391  HInv->SetOperator(*H);
392  if (partial_assembly)
393  {
394  delete HInvPC;
395  Vector diag_pa(vfes->GetTrueVSize());
396  H_form->AssembleDiagonal(diag_pa);
399  }
400 
401  // Extrapolated f^{n+1}.
402  for (auto &accel_term : accel_terms)
403  {
404  accel_term.coeff->SetTime(time + dt);
405  }
406 
407  f_form->Assemble();
409 
410  // Nonlinear extrapolated terms.
411  sw_extrap.Start();
412 
413  N->Mult(un, Nun);
414  N->Mult(unm1, Nunm1);
415  N->Mult(unm2, Nunm2);
416 
417  {
418  const auto d_Nun = Nun.Read();
419  const auto d_Nunm1 = Nunm1.Read();
420  const auto d_Nunm2 = Nunm2.Read();
421  auto d_Fext = Fext.Write();
422  const auto ab1_ = ab1;
423  const auto ab2_ = ab2;
424  const auto ab3_ = ab3;
425  mfem::forall(Fext.Size(), [=] MFEM_HOST_DEVICE (int i)
426  {
427  d_Fext[i] = ab1_ * d_Nun[i] +
428  ab2_ * d_Nunm1[i] +
429  ab3_ * d_Nunm2[i];
430  });
431  }
432 
433  Fext.Add(1.0, fn);
434 
435  // Fext = M^{-1} (F(u^{n}) + f^{n+1})
436  MvInv->Mult(Fext, tmp1);
439  Fext.Set(1.0, tmp1);
440 
441  // Compute BDF terms.
442  {
443  const double bd1idt = -bd1 / dt;
444  const double bd2idt = -bd2 / dt;
445  const double bd3idt = -bd3 / dt;
446  const auto d_un = un.Read();
447  const auto d_unm1 = unm1.Read();
448  const auto d_unm2 = unm2.Read();
449  auto d_Fext = Fext.ReadWrite();
450  mfem::forall(Fext.Size(), [=] MFEM_HOST_DEVICE (int i)
451  {
452  d_Fext[i] += bd1idt * d_un[i] +
453  bd2idt * d_unm1[i] +
454  bd3idt * d_unm2[i];
455  });
456  }
457 
458  sw_extrap.Stop();
459 
460  // Pressure Poisson.
461  sw_curlcurl.Start();
462  {
463  const auto d_un = un.Read();
464  const auto d_unm1 = unm1.Read();
465  const auto d_unm2 = unm2.Read();
466  auto d_Lext = Lext.Write();
467  const auto ab1_ = ab1;
468  const auto ab2_ = ab2;
469  const auto ab3_ = ab3;
470  mfem::forall(Lext.Size(), [=] MFEM_HOST_DEVICE (int i)
471  {
472  d_Lext[i] = ab1_ * d_un[i] +
473  ab2_ * d_unm1[i] +
474  ab3_ * d_unm2[i];
475  });
476  }
477 
479  if (pmesh->Dimension() == 2)
480  {
483  }
484  else
485  {
488  }
489 
491  Lext *= kin_vis;
492 
493  sw_curlcurl.Stop();
494 
495  // \tilde{F} = F - \nu CurlCurl(u)
496  FText.Set(-1.0, Lext);
497  FText.Add(1.0, Fext);
498 
499  // p_r = \nabla \cdot FText
500  D->Mult(FText, resp);
501  resp.Neg();
502 
503  // Add boundary terms.
507 
508  g_bdr_form->Assemble();
510  resp.Add(1.0, FText_bdr);
511  resp.Add(-bd0 / dt, g_bdr);
512 
513  if (pres_dbcs.empty())
514  {
516  }
517 
518  for (auto &pres_dbc : pres_dbcs)
519  {
520  pn_gf.ProjectBdrCoefficient(*pres_dbc.coeff, pres_dbc.attr);
521  }
522 
524 
525  Vector X1, B1;
526  if (partial_assembly)
527  {
528  auto *SpC = Sp.As<ConstrainedOperator>();
529  EliminateRHS(*Sp_form, *SpC, pres_ess_tdof, pn_gf, resp_gf, X1, B1, 1);
530  }
531  else
532  {
534  }
535  sw_spsolve.Start();
536  SpInv->Mult(B1, X1);
537  sw_spsolve.Stop();
541 
542  // If the boundary conditions on the pressure are pure Neumann remove the
543  // nullspace by removing the mean of the pressure solution. This is also
544  // ensured by the OrthoSolver wrapper for the preconditioner which removes
545  // the nullspace after every application.
546  if (pres_dbcs.empty())
547  {
548  MeanZero(pn_gf);
549  }
550 
552 
553  // Project velocity.
554  G->Mult(pn, resu);
555  resu.Neg();
556  Mv->Mult(Fext, tmp1);
557  resu.Add(1.0, tmp1);
558 
559  // un_next_gf = un_gf;
560 
561  for (auto &vel_dbc : vel_dbcs)
562  {
564  }
565 
567 
568  Vector X2, B2;
569  if (partial_assembly)
570  {
571  auto *HC = H.As<ConstrainedOperator>();
572  EliminateRHS(*H_form, *HC, vel_ess_tdof, un_next_gf, resu_gf, X2, B2, 1);
573  }
574  else
575  {
577  }
578  sw_hsolve.Start();
579  HInv->Mult(B2, X2);
580  sw_hsolve.Stop();
584 
586 
587  // If the current time step is not provisional, accept the computed solution
588  // and update the time step history by default.
589  if (!provisional)
590  {
592  time += dt;
593  }
594 
595  if (filter_alpha != 0.0)
596  {
599  const auto d_un_filtered_gf = un_filtered_gf.Read();
600  auto d_un_gf = un_gf.ReadWrite();
601  const auto filter_alpha_ = filter_alpha;
602  mfem::forall(un_gf.Size(), [=] MFEM_HOST_DEVICE (int i)
603  {
604  d_un_gf[i] = (1.0 - filter_alpha_) * d_un_gf[i]
605  + filter_alpha_ * d_un_filtered_gf[i];
606  });
607  }
608 
609  sw_step.Stop();
610 
611  if (verbose && pmesh->GetMyRank() == 0)
612  {
613  mfem::out << std::setw(7) << "" << std::setw(3) << "It" << std::setw(8)
614  << "Resid" << std::setw(12) << "Reltol"
615  << "\n";
616  // If numerical integration is active, there is no solve (thus no
617  // iterations), on the inverse velocity mass application.
618  if (!numerical_integ)
619  {
620  mfem::out << std::setw(5) << "MVIN " << std::setw(5) << std::fixed
621  << iter_mvsolve << " " << std::setw(3)
622  << std::setprecision(2) << std::scientific << res_mvsolve
623  << " " << 1e-12 << "\n";
624  }
625  mfem::out << std::setw(5) << "PRES " << std::setw(5) << std::fixed
626  << iter_spsolve << " " << std::setw(3) << std::setprecision(2)
627  << std::scientific << res_spsolve << " " << rtol_spsolve
628  << "\n";
629  mfem::out << std::setw(5) << "HELM " << std::setw(5) << std::fixed
630  << iter_hsolve << " " << std::setw(3) << std::setprecision(2)
631  << std::scientific << res_hsolve << " " << rtol_hsolve
632  << "\n";
633  mfem::out << std::setprecision(8);
634  mfem::out << std::fixed;
635  }
636 }
637 
639 {
640  // Make sure not to recompute the inner product linear form every
641  // application.
642  if (mass_lf == nullptr)
643  {
644  onecoeff.constant = 1.0;
645  mass_lf = new ParLinearForm(v.ParFESpace());
646  auto *dlfi = new DomainLFIntegrator(onecoeff);
647  if (numerical_integ)
648  {
649  const IntegrationRule &ir_ni = gll_rules.Get(vfes->GetFE(0)->GetGeomType(),
650  2 * order - 1);
651  dlfi->SetIntRule(&ir_ni);
652  }
654  mass_lf->Assemble();
655 
656  ParGridFunction one_gf(v.ParFESpace());
658 
659  volume = mass_lf->operator()(one_gf);
660  }
661 
662  double integ = mass_lf->operator()(v);
663 
664  v -= integ / volume;
665 }
666 
668  ConstrainedOperator &constrainedA,
669  const Array<int> &ess_tdof_list,
670  Vector &x,
671  Vector &b,
672  Vector &X,
673  Vector &B,
674  int copy_interior)
675 {
676  const Operator *Po = A.GetOutputProlongation();
677  const Operator *Pi = A.GetProlongation();
678  const Operator *Ri = A.GetRestriction();
679  A.InitTVectors(Po, Ri, Pi, x, b, X, B);
680  if (!copy_interior)
681  {
682  X.SetSubVectorComplement(ess_tdof_list, 0.0);
683  }
684  constrainedA.EliminateRHS(X, B);
685 }
686 
688 {
689  double loc_sum = v.Sum();
690  double global_sum = 0.0;
691  int loc_size = v.Size();
692  int global_size = 0;
693 
694  MPI_Allreduce(&loc_sum, &global_sum, 1, MPI_DOUBLE, MPI_SUM, pfes->GetComm());
695  MPI_Allreduce(&loc_size, &global_size, 1, MPI_INT, MPI_SUM, pfes->GetComm());
696 
697  v -= global_sum / static_cast<double>(global_size);
698 }
699 
701 {
702  FiniteElementSpace *fes = u.FESpace();
703 
704  // AccumulateAndCountZones.
705  Array<int> zones_per_vdof;
706  zones_per_vdof.SetSize(fes->GetVSize());
707  zones_per_vdof = 0;
708 
709  cu = 0.0;
710 
711  // Local interpolation.
712  int elndofs;
713  Array<int> vdofs;
714  Vector vals;
715  Vector loc_data;
716  int vdim = fes->GetVDim();
717  DenseMatrix grad_hat;
718  DenseMatrix dshape;
719  DenseMatrix grad;
720  Vector curl;
721 
722  for (int e = 0; e < fes->GetNE(); ++e)
723  {
724  fes->GetElementVDofs(e, vdofs);
725  u.GetSubVector(vdofs, loc_data);
726  vals.SetSize(vdofs.Size());
728  const FiniteElement *el = fes->GetFE(e);
729  elndofs = el->GetDof();
730  int dim = el->GetDim();
731  dshape.SetSize(elndofs, dim);
732 
733  for (int dof = 0; dof < elndofs; ++dof)
734  {
735  // Project.
736  const IntegrationPoint &ip = el->GetNodes().IntPoint(dof);
737  tr->SetIntPoint(&ip);
738 
739  // Eval and GetVectorGradientHat.
740  el->CalcDShape(tr->GetIntPoint(), dshape);
741  grad_hat.SetSize(vdim, dim);
742  DenseMatrix loc_data_mat(loc_data.GetData(), elndofs, vdim);
743  MultAtB(loc_data_mat, dshape, grad_hat);
744 
745  const DenseMatrix &Jinv = tr->InverseJacobian();
746  grad.SetSize(grad_hat.Height(), Jinv.Width());
747  Mult(grad_hat, Jinv, grad);
748 
749  curl.SetSize(3);
750  curl(0) = grad(2, 1) - grad(1, 2);
751  curl(1) = grad(0, 2) - grad(2, 0);
752  curl(2) = grad(1, 0) - grad(0, 1);
753 
754  for (int j = 0; j < curl.Size(); ++j)
755  {
756  vals(elndofs * j + dof) = curl(j);
757  }
758  }
759 
760  // Accumulate values in all dofs, count the zones.
761  for (int j = 0; j < vdofs.Size(); j++)
762  {
763  int ldof = vdofs[j];
764  cu(ldof) += vals[j];
765  zones_per_vdof[ldof]++;
766  }
767  }
768 
769  // Communication
770 
771  // Count the zones globally.
772  GroupCommunicator &gcomm = u.ParFESpace()->GroupComm();
773  gcomm.Reduce<int>(zones_per_vdof, GroupCommunicator::Sum);
774  gcomm.Bcast(zones_per_vdof);
775 
776  // Accumulate for all vdofs.
777  gcomm.Reduce<double>(cu.GetData(), GroupCommunicator::Sum);
778  gcomm.Bcast<double>(cu.GetData());
779 
780  // Compute means.
781  for (int i = 0; i < cu.Size(); i++)
782  {
783  const int nz = zones_per_vdof[i];
784  if (nz)
785  {
786  cu(i) /= nz;
787  }
788  }
789 }
790 
792  ParGridFunction &cu,
793  bool assume_scalar)
794 {
795  FiniteElementSpace *fes = u.FESpace();
796 
797  // AccumulateAndCountZones.
798  Array<int> zones_per_vdof;
799  zones_per_vdof.SetSize(fes->GetVSize());
800  zones_per_vdof = 0;
801 
802  cu = 0.0;
803 
804  // Local interpolation.
805  int elndofs;
806  Array<int> vdofs;
807  Vector vals;
808  Vector loc_data;
809  int vdim = fes->GetVDim();
810  DenseMatrix grad_hat;
811  DenseMatrix dshape;
812  DenseMatrix grad;
813  Vector curl;
814 
815  for (int e = 0; e < fes->GetNE(); ++e)
816  {
817  fes->GetElementVDofs(e, vdofs);
818  u.GetSubVector(vdofs, loc_data);
819  vals.SetSize(vdofs.Size());
821  const FiniteElement *el = fes->GetFE(e);
822  elndofs = el->GetDof();
823  int dim = el->GetDim();
824  dshape.SetSize(elndofs, dim);
825 
826  for (int dof = 0; dof < elndofs; ++dof)
827  {
828  // Project.
829  const IntegrationPoint &ip = el->GetNodes().IntPoint(dof);
830  tr->SetIntPoint(&ip);
831 
832  // Eval and GetVectorGradientHat.
833  el->CalcDShape(tr->GetIntPoint(), dshape);
834  grad_hat.SetSize(vdim, dim);
835  DenseMatrix loc_data_mat(loc_data.GetData(), elndofs, vdim);
836  MultAtB(loc_data_mat, dshape, grad_hat);
837 
838  const DenseMatrix &Jinv = tr->InverseJacobian();
839  grad.SetSize(grad_hat.Height(), Jinv.Width());
840  Mult(grad_hat, Jinv, grad);
841 
842  if (assume_scalar)
843  {
844  curl.SetSize(2);
845  curl(0) = grad(0, 1);
846  curl(1) = -grad(0, 0);
847  }
848  else
849  {
850  curl.SetSize(2);
851  curl(0) = grad(1, 0) - grad(0, 1);
852  curl(1) = 0.0;
853  }
854 
855  for (int j = 0; j < curl.Size(); ++j)
856  {
857  vals(elndofs * j + dof) = curl(j);
858  }
859  }
860 
861  // Accumulate values in all dofs, count the zones.
862  for (int j = 0; j < vdofs.Size(); j++)
863  {
864  int ldof = vdofs[j];
865  cu(ldof) += vals[j];
866  zones_per_vdof[ldof]++;
867  }
868  }
869 
870  // Communication.
871 
872  // Count the zones globally.
873  GroupCommunicator &gcomm = u.ParFESpace()->GroupComm();
874  gcomm.Reduce<int>(zones_per_vdof, GroupCommunicator::Sum);
875  gcomm.Bcast(zones_per_vdof);
876 
877  // Accumulate for all vdofs.
878  gcomm.Reduce<double>(cu.GetData(), GroupCommunicator::Sum);
879  gcomm.Bcast<double>(cu.GetData());
880 
881  // Compute means.
882  for (int i = 0; i < cu.Size(); i++)
883  {
884  const int nz = zones_per_vdof[i];
885  if (nz)
886  {
887  cu(i) /= nz;
888  }
889  }
890 }
891 
893 {
894  ParMesh *pmesh_u = u.ParFESpace()->GetParMesh();
895  FiniteElementSpace *fes = u.FESpace();
896  int vdim = fes->GetVDim();
897 
898  Vector ux, uy, uz;
899  Vector ur, us, ut;
900  double cflx = 0.0;
901  double cfly = 0.0;
902  double cflz = 0.0;
903  double cflm = 0.0;
904  double cflmax = 0.0;
905 
906  for (int e = 0; e < fes->GetNE(); ++e)
907  {
908  const FiniteElement *fe = fes->GetFE(e);
909  const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
910  fe->GetOrder());
912 
913  u.GetValues(e, ir, ux, 1);
914  ur.SetSize(ux.Size());
915  u.GetValues(e, ir, uy, 2);
916  us.SetSize(uy.Size());
917  if (vdim == 3)
918  {
919  u.GetValues(e, ir, uz, 3);
920  ut.SetSize(uz.Size());
921  }
922 
923  double hmin = pmesh_u->GetElementSize(e, 1) /
924  (double) fes->GetElementOrder(0);
925 
926  for (int i = 0; i < ir.GetNPoints(); ++i)
927  {
928  const IntegrationPoint &ip = ir.IntPoint(i);
929  tr->SetIntPoint(&ip);
930  const DenseMatrix &invJ = tr->InverseJacobian();
931  const double detJinv = 1.0 / tr->Jacobian().Det();
932 
933  if (vdim == 2)
934  {
935  ur(i) = (ux(i) * invJ(0, 0) + uy(i) * invJ(1, 0)) * detJinv;
936  us(i) = (ux(i) * invJ(0, 1) + uy(i) * invJ(1, 1)) * detJinv;
937  }
938  else if (vdim == 3)
939  {
940  ur(i) = (ux(i) * invJ(0, 0) + uy(i) * invJ(1, 0)
941  + uz(i) * invJ(2, 0))
942  * detJinv;
943  us(i) = (ux(i) * invJ(0, 1) + uy(i) * invJ(1, 1)
944  + uz(i) * invJ(2, 1))
945  * detJinv;
946  ut(i) = (ux(i) * invJ(0, 2) + uy(i) * invJ(1, 2)
947  + uz(i) * invJ(2, 2))
948  * detJinv;
949  }
950 
951  cflx = fabs(dt * ux(i) / hmin);
952  cfly = fabs(dt * uy(i) / hmin);
953  if (vdim == 3)
954  {
955  cflz = fabs(dt * uz(i) / hmin);
956  }
957  cflm = cflx + cfly + cflz;
958  cflmax = fmax(cflmax, cflm);
959  }
960  }
961 
962  double cflmax_global = 0.0;
963  MPI_Allreduce(&cflmax,
964  &cflmax_global,
965  1,
966  MPI_DOUBLE,
967  MPI_MAX,
968  pmesh_u->GetComm());
969 
970  return cflmax_global;
971 }
972 
974 {
975  vel_dbcs.emplace_back(attr, coeff);
976 
977  if (verbose && pmesh->GetMyRank() == 0)
978  {
979  mfem::out << "Adding Velocity Dirichlet BC to attributes ";
980  for (int i = 0; i < attr.Size(); ++i)
981  {
982  if (attr[i] == 1)
983  {
984  mfem::out << i << " ";
985  }
986  }
987  mfem::out << std::endl;
988  }
989 
990  for (int i = 0; i < attr.Size(); ++i)
991  {
992  MFEM_ASSERT((vel_ess_attr[i] && attr[i]) == 0,
993  "Duplicate boundary definition deteceted.");
994  if (attr[i] == 1)
995  {
996  vel_ess_attr[i] = 1;
997  }
998  }
999 }
1000 
1002 {
1004 }
1005 
1007 {
1008  pres_dbcs.emplace_back(attr, coeff);
1009 
1010  if (verbose && pmesh->GetMyRank() == 0)
1011  {
1012  mfem::out << "Adding Pressure Dirichlet BC to attributes ";
1013  for (int i = 0; i < attr.Size(); ++i)
1014  {
1015  if (attr[i] == 1)
1016  {
1017  mfem::out << i << " ";
1018  }
1019  }
1020  mfem::out << std::endl;
1021  }
1022 
1023  for (int i = 0; i < attr.Size(); ++i)
1024  {
1025  MFEM_ASSERT((pres_ess_attr[i] && attr[i]) == 0,
1026  "Duplicate boundary definition deteceted.");
1027  if (attr[i] == 1)
1028  {
1029  pres_ess_attr[i] = 1;
1030  }
1031  }
1032 }
1033 
1035 {
1037 }
1038 
1040 {
1041  accel_terms.emplace_back(attr, coeff);
1042 
1043  if (verbose && pmesh->GetMyRank() == 0)
1044  {
1045  mfem::out << "Adding Acceleration term to attributes ";
1046  for (int i = 0; i < attr.Size(); ++i)
1047  {
1048  if (attr[i] == 1)
1049  {
1050  mfem::out << i << " ";
1051  }
1052  }
1053  mfem::out << std::endl;
1054  }
1055 }
1056 
1058 {
1060 }
1061 
1063 {
1064  // Maximum BDF order to use at current time step
1065  // step + 1 <= order <= max_bdf_order
1066  int bdf_order = std::min(step + 1, max_bdf_order);
1067 
1068  // Ratio of time step history at dt(t_{n}) - dt(t_{n-1})
1069  double rho1 = 0.0;
1070 
1071  // Ratio of time step history at dt(t_{n-1}) - dt(t_{n-2})
1072  double rho2 = 0.0;
1073 
1074  rho1 = dthist[0] / dthist[1];
1075 
1076  if (bdf_order == 3)
1077  {
1078  rho2 = dthist[1] / dthist[2];
1079  }
1080 
1081  if (step == 0 && bdf_order == 1)
1082  {
1083  bd0 = 1.0;
1084  bd1 = -1.0;
1085  bd2 = 0.0;
1086  bd3 = 0.0;
1087  ab1 = 1.0;
1088  ab2 = 0.0;
1089  ab3 = 0.0;
1090  }
1091  else if (step >= 1 && bdf_order == 2)
1092  {
1093  bd0 = (1.0 + 2.0 * rho1) / (1.0 + rho1);
1094  bd1 = -(1.0 + rho1);
1095  bd2 = pow(rho1, 2.0) / (1.0 + rho1);
1096  bd3 = 0.0;
1097  ab1 = 1.0 + rho1;
1098  ab2 = -rho1;
1099  ab3 = 0.0;
1100  }
1101  else if (step >= 2 && bdf_order == 3)
1102  {
1103  bd0 = 1.0 + rho1 / (1.0 + rho1)
1104  + (rho2 * rho1) / (1.0 + rho2 * (1 + rho1));
1105  bd1 = -1.0 - rho1 - (rho2 * rho1 * (1.0 + rho1)) / (1.0 + rho2);
1106  bd2 = pow(rho1, 2.0) * (rho2 + 1.0 / (1.0 + rho1));
1107  bd3 = -(pow(rho2, 3.0) * pow(rho1, 2.0) * (1.0 + rho1))
1108  / ((1.0 + rho2) * (1.0 + rho2 + rho2 * rho1));
1109  ab1 = ((1.0 + rho1) * (1.0 + rho2 * (1.0 + rho1))) / (1.0 + rho2);
1110  ab2 = -rho1 * (1.0 + rho2 * (1.0 + rho1));
1111  ab3 = (pow(rho2, 2.0) * rho1 * (1.0 + rho1)) / (1.0 + rho2);
1112  }
1113 }
1114 
1116 {
1117  double my_rt[6], rt_max[6];
1118 
1119  my_rt[0] = sw_setup.RealTime();
1120  my_rt[1] = sw_step.RealTime();
1121  my_rt[2] = sw_extrap.RealTime();
1122  my_rt[3] = sw_curlcurl.RealTime();
1123  my_rt[4] = sw_spsolve.RealTime();
1124  my_rt[5] = sw_hsolve.RealTime();
1125 
1126  MPI_Reduce(my_rt, rt_max, 6, MPI_DOUBLE, MPI_MAX, 0, pmesh->GetComm());
1127 
1128  if (pmesh->GetMyRank() == 0)
1129  {
1130  mfem::out << std::setw(10) << "SETUP" << std::setw(10) << "STEP"
1131  << std::setw(10) << "EXTRAP" << std::setw(10) << "CURLCURL"
1132  << std::setw(10) << "PSOLVE" << std::setw(10) << "HSOLVE"
1133  << "\n";
1134 
1135  mfem::out << std::setprecision(3) << std::setw(10) << my_rt[0]
1136  << std::setw(10) << my_rt[1] << std::setw(10) << my_rt[2]
1137  << std::setw(10) << my_rt[3] << std::setw(10) << my_rt[4]
1138  << std::setw(10) << my_rt[5] << "\n";
1139 
1140  mfem::out << std::setprecision(3) << std::setw(10) << " " << std::setw(10)
1141  << my_rt[1] / my_rt[1] << std::setw(10) << my_rt[2] / my_rt[1]
1142  << std::setw(10) << my_rt[3] / my_rt[1] << std::setw(10)
1143  << my_rt[4] / my_rt[1] << std::setw(10) << my_rt[5] / my_rt[1]
1144  << "\n";
1145 
1146  mfem::out << std::setprecision(8);
1147  }
1148 }
1149 
1151 {
1152  int fes_size0 = vfes->GlobalVSize();
1153  int fes_size1 = pfes->GlobalVSize();
1154 
1155  if (pmesh->GetMyRank() == 0)
1156  {
1157  mfem::out << "NAVIER version: " << NAVIER_VERSION << std::endl
1158  << "MFEM version: " << MFEM_VERSION << std::endl
1159  << "MFEM GIT: " << MFEM_GIT_STRING << std::endl
1160  << "Velocity #DOFs: " << fes_size0 << std::endl
1161  << "Pressure #DOFs: " << fes_size1 << std::endl;
1162  }
1163 }
1164 
1166 {
1167  delete FText_gfcoeff;
1168  delete g_bdr_form;
1169  delete FText_bdr_form;
1170  delete mass_lf;
1171  delete Mv_form;
1172  delete N;
1173  delete Sp_form;
1174  delete D_form;
1175  delete G_form;
1176  delete HInvPC;
1177  delete HInv;
1178  delete H_form;
1179  delete SpInv;
1180  delete MvInvPC;
1181  delete SpInvOrthoPC;
1182  delete SpInvPC;
1183  delete lor;
1184  delete f_form;
1185  delete MvInv;
1186  delete vfec;
1187  delete pfec;
1188  delete vfes;
1189  delete pfes;
1190  delete vfec_filter;
1191  delete vfes_filter;
1192 }
Abstract class for all finite elements.
Definition: fe_base.hpp:233
FiniteElementCollection * vfec_filter
void PrintTimingData()
Print timing summary of the solving routine.
void Orthogonalize(Vector &v)
Remove mean from a Vector.
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:1031
Conjugate gradient method.
Definition: solvers.hpp:493
FiniteElementCollection * pfec
Pressure finite element collection.
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:253
Create and assemble a low-order refined version of a ParBilinearForm.
Definition: lor.hpp:165
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
void EliminateRHS(Operator &A, ConstrainedOperator &constrainedA, const Array< int > &ess_tdof_list, Vector &x, Vector &b, Vector &X, Vector &B, int copy_interior=0)
Eliminate essential BCs in an Operator and apply to RHS.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:980
Base class for vector Coefficients that optionally depend on time and space.
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(...
void PrintInfo()
Print information about the Navier version.
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
Array< BilinearFormIntegrator * > * GetDBFI()
Access all the integrators added with AddDomainIntegrator().
void Bcast(T *ldata, int layout) const
Broadcast within each group where the master is the root.
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
virtual void Setup()
Setup the NonlinearForm: based on the current AssemblyLevel and the current mesh, optionally...
ConstantCoefficient nlcoeff
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
void AddDomainIntegrator(NonlinearFormIntegrator *nlfi)
Adds new Domain Integrator.
virtual const Operator * GetRestriction() const
Restriction operator from input vectors for the operator to linear algebra (linear system) vectors...
Definition: operator.hpp:143
void(const Vector &x, double t, Vector &u) VecFuncT
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
virtual void Mult(const Vector &x, Vector &y) const
Evaluate the action of the NonlinearForm.
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:718
double Det() const
Definition: densemat.cpp:487
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
void InitTVectors(const Operator *Po, const Operator *Ri, const Operator *Pi, Vector &x, Vector &b, Vector &X, Vector &B) const
Initializes memory for true vectors of linear system.
Definition: operator.cpp:22
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
ParGridFunction un_filtered_gf
Parallel non-linear operator on the true dofs.
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:453
virtual void SetIntRule(const IntegrationRule *ir)
Prescribe a fixed IntegrationRule to use (when ir != NULL) or let the integrator choose (when ir == N...
Definition: nonlininteg.hpp:57
double(const Vector &x, double t) ScalarFuncT
ParLinearForm * FText_bdr_form
std::vector< PresDirichletBC_T > pres_dbcs
const DenseMatrix & InverseJacobian()
Return the inverse of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:145
bool numerical_integ
Enable/disable numerical integration rules of forms.
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
double RealTime()
Return the number of real seconds elapsed since the stopwatch was started.
Definition: tic_toc.cpp:429
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the ParGridFunction.
Definition: pgridfunc.cpp:96
Abstract parallel finite element space.
Definition: pfespace.hpp:28
ConstantCoefficient H_lincoeff
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:561
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:686
A class container for 1D quadrature type constants.
Definition: intrules.hpp:390
const IntegrationPoint & GetIntPoint()
Get a const reference to the currently set integration point. This will return NULL if no integration...
Definition: eltrans.hpp:98
NavierSolver(ParMesh *mesh, int order, double kin_vis)
Initialize data structures, set FE space order and kinematic viscosity.
virtual const SparseMatrix * GetRestrictionMatrix() const
Get the R matrix which restricts a local dof vector to true dof vector.
Definition: pfespace.hpp:379
ParMixedBilinearForm * G_form
ParLORDiscretization * lor
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2841
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:480
double GetFinalNorm() const
Returns the final residual norm after termination of the solver during the last call to Mult()...
Definition: solvers.hpp:265
virtual void Update(FiniteElementSpace *nfes=NULL)
Update the FiniteElementSpace and delete all data associated with the old one.
void Setup(double dt)
Initialize forms, solvers and preconditioners.
void MeanZero(ParGridFunction &v)
Remove the mean from a ParGridFunction.
std::function< double(const Vector &)> f(double mass_coeff)
Definition: lor_mms.hpp:30
void Stop()
Stop the stopwatch.
Definition: tic_toc.cpp:419
bool partial_assembly
Enable/disable partial assembly of forms.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds a domain integrator. Assumes ownership of bfi.
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:320
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1590
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:256
void ComputeCurl3D(ParGridFunction &u, ParGridFunction &cu)
Compute for .
Class for parallel linear form.
Definition: plinearform.hpp:26
Communicator performing operations within groups defined by a GroupTopology with arbitrary-size data ...
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
FiniteElementCollection * vfec
Velocity finite element collection.
double GetElementSize(ElementTransformation *T, int type=0)
Definition: mesh.cpp:83
static void Sum(OpData< T >)
Reduce operation Sum, instantiated for int and double.
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: pgridfunc.hpp:169
virtual void AssembleDiagonal(Vector &diag) const
Assemble the diagonal of the bilinear form into diag. Note that diag is a true-dof Vector...
double ComputeCFL(ParGridFunction &u, double dt)
Compute CFL.
Class for boundary integration .
Definition: lininteg.hpp:210
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:461
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:302
int GetNumIterations() const
Returns the number of iterations taken during the last call to Mult()
Definition: solvers.hpp:250
void Assemble()
Assembles the ParLinearForm i.e. sums over all domain/bdr integrators.
Definition: plinearform.cpp:46
double b
Definition: lissajous.cpp:42
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1670
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:119
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
ConstantCoefficient onecoeff
void Assemble(int skip_zeros=1)
double Sum() const
Return the sum of the vector entries.
Definition: vector.cpp:1286
void Assemble(int skip_zeros=1)
Assemble the local matrix.
std::vector< double > dthist
void SetTimeIntegrationCoefficients(int step)
Update the EXTk/BDF time integration coefficient.
virtual void FormRectangularSystemMatrix(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, OperatorHandle &A)
Return in A a parallel (on truedofs) version of this operator.
void AddVelDirichletBC(VectorCoefficient *coeff, Array< int > &attr)
Add a Dirichlet boundary condition to the velocity field.
Parallel smoothers in hypre.
Definition: hypre.hpp:971
void AddBoundaryIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:72
void ComputeCurl2D(ParGridFunction &u, ParGridFunction &cu, bool assume_scalar=false)
Compute for .
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::LEGACY.
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:736
virtual const Operator * GetProlongation() const
Prolongation operator from linear algebra (linear system) vectors, to input vectors for the operator...
Definition: operator.hpp:139
MPI_Comm GetComm() const
Definition: pmesh.hpp:351
void Start()
Start the stopwatch. The elapsed time is not cleared.
Definition: tic_toc.cpp:408
void SetSubVectorComplement(const Array< int > &dofs, const double val)
Set all vector entries NOT in the dofs Array to the given val.
Definition: vector.cpp:740
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:206
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:702
A general vector function coefficient.
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::LEGACY.
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i&#39;th element. The returned indices are offsets into an ...
Definition: fespace.cpp:281
int GetMyRank() const
Definition: pmesh.hpp:353
HypreParMatrix & GetAssembledMatrix() const
Return the assembled LOR operator as a HypreParMatrix.
Definition: lor.cpp:522
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:41
virtual void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:908
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
void SetRelTol(double rtol)
Definition: solvers.hpp:199
void EliminateRHS(const Vector &x, Vector &b) const
Eliminate "essential boundary condition" values specified in x from the given right-hand side b...
Definition: operator.cpp:520
ParFiniteElementSpace * vfes
Velocity finite element space.
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe_base.hpp:311
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
virtual const FiniteElement * GetFE(int i) const
Definition: pfespace.cpp:534
void forall(int N, lambda &&body)
Definition: forall.hpp:742
void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff, Array< int > &attr)
Definition: pgridfunc.cpp:669
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
HYPRE_BigInt GlobalVSize() const
Definition: pfespace.hpp:279
ParLinearForm * mass_lf
Linear form to compute the mass matrix in various subroutines.
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
void ProjectGridFunction(const GridFunction &src)
Project the src GridFunction to this GridFunction, both of which must be on the same mesh...
Definition: gridfunc.cpp:1849
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:113
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
Vector & Set(const double a, const Vector &x)
(*this) = a * x
Definition: vector.cpp:263
std::vector< AccelTerm_T > accel_terms
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
virtual int GetTrueVSize() const
Return the number of local vector true dofs.
Definition: pfespace.hpp:285
bool verbose
Enable/disable verbose output.
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
Solver wrapper which orthogonalizes the input and output vector.
Definition: solvers.hpp:1180
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:469
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:248
double kin_vis
Kinematic viscosity (dimensionless).
Class for integration point with weight.
Definition: intrules.hpp:31
Class for parallel bilinear form using different test and trial FE spaces.
int GetElementOrder(int i) const
Returns the order of the i&#39;th finite element.
Definition: fespace.cpp:178
std::vector< VelDirichletBC_T > vel_dbcs
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:104
VectorGridFunctionCoefficient * FText_gfcoeff
void UpdateTimestepHistory(double dt)
Rotate entries in the time step and solution history arrays.
ParFiniteElementSpace * pfes
Pressure finite element space.
int dim
Definition: ex24.cpp:53
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void Reduce(T *ldata, void(*Op)(OpData< T >)) const
Reduce within each group where the master is the root.
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:709
Class for parallel bilinear form.
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
ParMixedBilinearForm * D_form
MPI_Comm GetComm() const
Definition: pfespace.hpp:269
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:507
A general function coefficient.
Vector data type.
Definition: vector.hpp:58
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
Definition: pgridfunc.cpp:152
void SetSolver(Solver &s)
Set the solver used by the OrthoSolver.
Definition: solvers.cpp:3433
Vector coefficient defined by a vector GridFunction.
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
virtual const Operator * GetOutputProlongation() const
Prolongation operator from linear algebra (linear system) vectors, to output vectors for the operator...
Definition: operator.hpp:147
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve the linear system Ax=b.
Definition: hypre.cpp:3901
void Step(double &time, double dt, int cur_step, bool provisional=false)
Compute solution at the next time step t+dt.
ParFiniteElementSpace * vfes_filter
void AddPresDirichletBC(Coefficient *coeff, Array< int > &attr)
Add a Dirichlet boundary condition to the pressure field.
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
Class for parallel grid function.
Definition: pgridfunc.hpp:32
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
Square Operator for imposing essential boundary conditions using only the action, Mult()...
Definition: operator.hpp:872
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:389
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Class for parallel meshes.
Definition: pmesh.hpp:32
void ParallelAssemble(Vector &tv)
Assemble the vector on the true dofs, i.e. P^t v.
Definition: plinearform.cpp:95
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:327
ParMesh * pmesh
The parallel mesh.
void AddAccelTerm(VectorCoefficient *coeff, Array< int > &attr)
Add an acceleration term to the RHS of the equation.
int order
The order of the velocity and pressure space.
ConstantCoefficient H_bdfcoeff
void Neg()
(*this) = -(*this)
Definition: vector.cpp:301
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:769