MFEM  v4.5.2
Finite element discretization library
navier_solver.hpp
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 #ifndef MFEM_NAVIER_SOLVER_HPP
13 #define MFEM_NAVIER_SOLVER_HPP
14 
15 #define NAVIER_VERSION 0.1
16 
17 #include "mfem.hpp"
18 
19 namespace mfem
20 {
21 namespace navier
22 {
23 using VecFuncT = void(const Vector &x, double t, Vector &u);
24 using ScalarFuncT = double(const Vector &x, double t);
25 
26 /// Container for a Dirichlet boundary condition of the velocity field.
28 {
29 public:
31  : attr(attr), coeff(coeff)
32  {}
33 
35  {
36  // Deep copy the attribute array
37  this->attr = obj.attr;
38 
39  // Move the coefficient pointer
40  this->coeff = obj.coeff;
41  obj.coeff = nullptr;
42  }
43 
44  ~VelDirichletBC_T() { delete coeff; }
45 
48 };
49 
50 /// Container for a Dirichlet boundary condition of the pressure field.
52 {
53 public:
55  : attr(attr), coeff(coeff)
56  {}
57 
59  {
60  // Deep copy the attribute array
61  this->attr = obj.attr;
62 
63  // Move the coefficient pointer
64  this->coeff = obj.coeff;
65  obj.coeff = nullptr;
66  }
67 
68  ~PresDirichletBC_T() { delete coeff; }
69 
72 };
73 
74 /// Container for an acceleration term.
76 {
77 public:
79  : attr(attr), coeff(coeff)
80  {}
81 
83  {
84  // Deep copy the attribute array
85  this->attr = obj.attr;
86 
87  // Move the coefficient pointer
88  this->coeff = obj.coeff;
89  obj.coeff = nullptr;
90  }
91 
92  ~AccelTerm_T() { delete coeff; }
93 
96 };
97 
98 /// Transient incompressible Navier Stokes solver in a split scheme formulation.
99 /**
100  * This implementation of a transient incompressible Navier Stokes solver uses
101  * the non-dimensionalized formulation. The coupled momentum and
102  * incompressibility equations are decoupled using the split scheme described in
103  * [1]. This leads to three solving steps.
104  *
105  * 1. An extrapolation step for all nonlinear terms which are treated
106  * explicitly. This step avoids a fully coupled nonlinear solve and only
107  * requires a solve of the mass matrix in velocity space \f$M_v^{-1}\f$. On
108  * the other hand this introduces a CFL stability condition on the maximum
109  * timestep.
110  *
111  * 2. A Poisson solve \f$S_p^{-1}\f$.
112  *
113  * 3. A Helmholtz like solve \f$(M_v - \partial t K_v)^{-1}\f$.
114  *
115  * The numerical solver setup for each step are as follows.
116  *
117  * \f$M_v^{-1}\f$ is solved using CG with Jacobi as preconditioner.
118  *
119  * \f$S_p^{-1}\f$ is solved using CG with AMG applied to the low order refined
120  * (LOR) assembled pressure Poisson matrix. To avoid assembling a matrix for
121  * preconditioning, one can use p-MG as an alternative (NYI).
122  *
123  * \f$(M_v - \partial t K_v)^{-1}\f$ due to the CFL condition we expect the time
124  * step to be small. Therefore this is solved using CG with Jacobi as
125  * preconditioner. For large time steps a preconditioner like AMG or p-MG should
126  * be used (NYI).
127  *
128  * Statements marked with NYI mean this feature is planned but Not Yet
129  * Implemented.
130  *
131  * A detailed description is available in [1] section 4.2. The algorithm is
132  * originated from [2].
133  *
134  * [1] Michael Franco, Jean-Sylvain Camier, Julian Andrej, Will Pazner (2020)
135  * High-order matrix-free incompressible flow solvers with GPU acceleration and
136  * low-order refined preconditioners (https://arxiv.org/abs/1910.03032)
137  *
138  * [2] A. G. Tomboulides, J. C. Y. Lee & S. A. Orszag (1997) Numerical
139  * Simulation of Low Mach Number Reactive Flows
140  */
142 {
143 public:
144  /// Initialize data structures, set FE space order and kinematic viscosity.
145  /**
146  * The ParMesh @a mesh can be a linear or curved parallel mesh. The @a order
147  * of the finite element spaces is this algorithm is of equal order
148  * \f$(P_N)^d P_N\f$ for velocity and pressure respectively. This means the
149  * pressure is in discretized in the same space (just scalar instead of a
150  * vector space) as the velocity.
151  *
152  * Kinematic viscosity (dimensionless) is set using @a kin_vis and
153  * automatically converted to the Reynolds number. If you want to set the
154  * Reynolds number directly, you can provide the inverse.
155  */
156  NavierSolver(ParMesh *mesh, int order, double kin_vis);
157 
158  /// Initialize forms, solvers and preconditioners.
159  void Setup(double dt);
160 
161  /// Compute solution at the next time step t+dt.
162  /**
163  * This method can be called with the default value @a provisional which
164  * always accepts the computed time step by automatically calling
165  * UpdateTimestepHistory.
166  *
167  * If @a provisional is set to true, the solution at t+dt is not accepted
168  * automatically and the application code has to call UpdateTimestepHistory
169  * and update the @a time variable accordingly.
170  *
171  * The application code can check the provisional step by retrieving the
172  * GridFunction with the method GetProvisionalVelocity. If the check fails,
173  * it is possible to retry the step with a different time step by not
174  * calling UpdateTimestepHistory and calling this method with the previous
175  * @a time and @a cur_step.
176  *
177  * The method and parameter choices are based on [1].
178  *
179  * [1] D. Wang, S.J. Ruuth (2008) Variable step-size implicit-explicit
180  * linear multistep methods for time-dependent partial differential
181  * equations
182  */
183  void Step(double &time, double dt, int cur_step, bool provisional = false);
184 
185  /// Return a pointer to the provisional velocity ParGridFunction.
187 
188  /// Return a pointer to the current velocity ParGridFunction.
190 
191  /// Return a pointer to the current pressure ParGridFunction.
193 
194  /// Add a Dirichlet boundary condition to the velocity field.
195  void AddVelDirichletBC(VectorCoefficient *coeff, Array<int> &attr);
196 
197  void AddVelDirichletBC(VecFuncT *f, Array<int> &attr);
198 
199  /// Add a Dirichlet boundary condition to the pressure field.
200  void AddPresDirichletBC(Coefficient *coeff, Array<int> &attr);
201 
203 
204  /// Add an acceleration term to the RHS of the equation.
205  /**
206  * The VecFuncT @a f is evaluated at the current time t and extrapolated
207  * together with the nonlinear parts of the Navier Stokes equation.
208  */
209  void AddAccelTerm(VectorCoefficient *coeff, Array<int> &attr);
210 
211  void AddAccelTerm(VecFuncT *f, Array<int> &attr);
212 
213  /// Enable partial assembly for every operator.
214  void EnablePA(bool pa) { partial_assembly = pa; }
215 
216  /// Enable numerical integration rules. This means collocated quadrature at
217  /// the nodal points.
218  void EnableNI(bool ni) { numerical_integ = ni; }
219 
220  /// Print timing summary of the solving routine.
221  /**
222  * The summary shows the timing in seconds in the first row of
223  *
224  * 1. SETUP: Time spent for the setup of all forms, solvers and
225  * preconditioners.
226  * 2. STEP: Time spent computing a full time step. It includes all three
227  * solves.
228  * 3. EXTRAP: Time spent for extrapolation of all forcing and nonlinear
229  * terms.
230  * 4. CURLCURL: Time spent for computing the curl curl term in the pressure
231  * Poisson equation (see references for detailed explanation).
232  * 5. PSOLVE: Time spent in the pressure Poisson solve.
233  * 6. HSOLVE: Time spent in the Helmholtz solve.
234  *
235  * The second row shows a proportion of a column relative to the whole
236  * time step.
237  */
238  void PrintTimingData();
239 
240  ~NavierSolver();
241 
242  /// Compute \f$\nabla \times \nabla \times u\f$ for \f$u \in (H^1)^2\f$.
244  ParGridFunction &cu,
245  bool assume_scalar = false);
246 
247  /// Compute \f$\nabla \times \nabla \times u\f$ for \f$u \in (H^1)^3\f$.
249 
250  /// Remove mean from a Vector.
251  /**
252  * Modify the Vector @a v by subtracting its mean using
253  * \f$v = v - \frac{\sum_i^N v_i}{N} \f$
254  */
255  void Orthogonalize(Vector &v);
256 
257  /// Remove the mean from a ParGridFunction.
258  /**
259  * Modify the ParGridFunction @a v by subtracting its mean using
260  * \f$ v = v - \int_\Omega \frac{v}{vol(\Omega)} dx \f$.
261  */
262  void MeanZero(ParGridFunction &v);
263 
264  /// Rotate entries in the time step and solution history arrays.
265  void UpdateTimestepHistory(double dt);
266 
267  /// Set the maximum order to use for the BDF method.
268  void SetMaxBDFOrder(int maxbdforder) { max_bdf_order = maxbdforder; };
269 
270  /// Compute CFL
271  double ComputeCFL(ParGridFunction &u, double dt);
272 
273  /// Set the number of modes to cut off in the interpolation filter
274  void SetCutoffModes(int c) { filter_cutoff_modes = c; }
275 
276  /// Set the interpolation filter parameter @a a
277  /**
278  * If @a a is > 0, the filtering algorithm for the velocity field after every
279  * time step from [1] is used. The parameter should be 0 > @a >= 1.
280  *
281  * [1] Paul Fischer, Julia Mullen (2001) Filter-based stabilization of
282  * spectral element methods
283  */
284  void SetFilterAlpha(double a) { filter_alpha = a; }
285 
286 protected:
287  /// Print information about the Navier version.
288  void PrintInfo();
289 
290  /// Update the EXTk/BDF time integration coefficient.
291  /**
292  * Depending on which time step the computation is in, the EXTk/BDF time
293  * integration coefficients have to be set accordingly. This allows
294  * bootstrapping with a BDF scheme of order 1 and increasing the order each
295  * following time step, up to order 3 (or whichever order is set in
296  * SetMaxBDFOrder).
297  */
298  void SetTimeIntegrationCoefficients(int step);
299 
300  /// Eliminate essential BCs in an Operator and apply to RHS.
301  void EliminateRHS(Operator &A,
302  ConstrainedOperator &constrainedA,
303  const Array<int> &ess_tdof_list,
304  Vector &x,
305  Vector &b,
306  Vector &X,
307  Vector &B,
308  int copy_interior = 0);
309 
310  /// Enable/disable debug output.
311  bool debug = false;
312 
313  /// Enable/disable verbose output.
314  bool verbose = true;
315 
316  /// Enable/disable partial assembly of forms.
317  bool partial_assembly = false;
318 
319  /// Enable/disable numerical integration rules of forms.
320  bool numerical_integ = false;
321 
322  /// The parallel mesh.
323  ParMesh *pmesh = nullptr;
324 
325  /// The order of the velocity and pressure space.
326  int order;
327 
328  /// Kinematic viscosity (dimensionless).
329  double kin_vis;
330 
332 
333  /// Velocity \f$H^1\f$ finite element collection.
335 
336  /// Pressure \f$H^1\f$ finite element collection.
338 
339  /// Velocity \f$(H^1)^d\f$ finite element space.
341 
342  /// Pressure \f$H^1\f$ finite element space.
344 
345  ParNonlinearForm *N = nullptr;
346 
348 
350 
352 
354 
356 
358 
360 
361  ParLinearForm *f_form = nullptr;
362 
364 
365  /// Linear form to compute the mass matrix in various subroutines.
366  ParLinearForm *mass_lf = nullptr;
368  double volume = 0.0;
369 
374 
380 
381  Solver *MvInvPC = nullptr;
382  CGSolver *MvInv = nullptr;
383 
386  CGSolver *SpInv = nullptr;
387 
388  Solver *HInvPC = nullptr;
389  CGSolver *HInv = nullptr;
390 
392  resu;
394 
396 
398  resu_gf;
399 
401 
402  // All essential attributes.
405 
406  // All essential true dofs.
409 
410  // Bookkeeping for velocity dirichlet bcs.
411  std::vector<VelDirichletBC_T> vel_dbcs;
412 
413  // Bookkeeping for pressure dirichlet bcs.
414  std::vector<PresDirichletBC_T> pres_dbcs;
415 
416  // Bookkeeping for acceleration (forcing) terms.
417  std::vector<AccelTerm_T> accel_terms;
418 
419  int max_bdf_order = 3;
420  int cur_step = 0;
421  std::vector<double> dthist = {0.0, 0.0, 0.0};
422 
423  // BDFk/EXTk coefficients.
424  double bd0 = 0.0;
425  double bd1 = 0.0;
426  double bd2 = 0.0;
427  double bd3 = 0.0;
428  double ab1 = 0.0;
429  double ab2 = 0.0;
430  double ab3 = 0.0;
431 
432  // Timers.
434 
435  // Print levels.
436  int pl_mvsolve = 0;
437  int pl_spsolve = 0;
438  int pl_hsolve = 0;
439  int pl_amg = 0;
440 
441  // Relative tolerances.
442  double rtol_spsolve = 1e-6;
443  double rtol_hsolve = 1e-8;
444 
445  // Iteration counts.
447 
448  // Residuals.
449  double res_mvsolve = 0.0, res_spsolve = 0.0, res_hsolve = 0.0;
450 
451  // LOR related.
453 
454  // Filter-based stabilization
456  double filter_alpha = 0.0;
461 };
462 
463 } // namespace navier
464 
465 } // namespace mfem
466 
467 #endif
FiniteElementCollection * vfec_filter
VectorCoefficient * coeff
Container for a Dirichlet boundary condition of the velocity field.
void PrintTimingData()
Print timing summary of the solving routine.
void Orthogonalize(Vector &v)
Remove mean from a Vector.
Conjugate gradient method.
Definition: solvers.hpp:493
FiniteElementCollection * pfec
Pressure finite element collection.
Create and assemble a low-order refined version of a ParBilinearForm.
Definition: lor.hpp:165
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.
Container for a Dirichlet boundary condition of the pressure field.
Base class for vector Coefficients that optionally depend on time and space.
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
void PrintInfo()
Print information about the Navier version.
void SetFilterAlpha(double a)
Set the interpolation filter parameter a.
ConstantCoefficient nlcoeff
void SetMaxBDFOrder(int maxbdforder)
Set the maximum order to use for the BDF method.
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
void(const Vector &x, double t, Vector &u) VecFuncT
Container class for integration rules.
Definition: intrules.hpp:311
ParGridFunction un_filtered_gf
Parallel non-linear operator on the true dofs.
double(const Vector &x, double t) ScalarFuncT
ParLinearForm * FText_bdr_form
std::vector< PresDirichletBC_T > pres_dbcs
bool numerical_integ
Enable/disable numerical integration rules of forms.
Abstract parallel finite element space.
Definition: pfespace.hpp:28
ConstantCoefficient H_lincoeff
Container for an acceleration term.
NavierSolver(ParMesh *mesh, int order, double kin_vis)
Initialize data structures, set FE space order and kinematic viscosity.
AccelTerm_T(Array< int > attr, VectorCoefficient *coeff)
PresDirichletBC_T(PresDirichletBC_T &&obj)
ParMixedBilinearForm * G_form
ParLORDiscretization * lor
void Setup(double dt)
Initialize forms, solvers and preconditioners.
void MeanZero(ParGridFunction &v)
Remove the mean from a ParGridFunction.
VelDirichletBC_T(VelDirichletBC_T &&obj)
bool partial_assembly
Enable/disable partial assembly of forms.
ConstantCoefficient Sp_coeff
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1590
void ComputeCurl3D(ParGridFunction &u, ParGridFunction &cu)
Compute for .
Class for parallel linear form.
Definition: plinearform.hpp:26
FiniteElementCollection * vfec
Velocity finite element collection.
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
double ComputeCFL(ParGridFunction &u, double dt)
Compute CFL.
double b
Definition: lissajous.cpp:42
ParGridFunction * GetCurrentPressure()
Return a pointer to the current pressure ParGridFunction.
Timing object.
Definition: tic_toc.hpp:34
ConstantCoefficient onecoeff
void SetCutoffModes(int c)
Set the number of modes to cut off in the interpolation filter.
std::vector< double > dthist
void SetTimeIntegrationCoefficients(int step)
Update the EXTk/BDF time integration coefficient.
void EnablePA(bool pa)
Enable partial assembly for every operator.
void AddVelDirichletBC(VectorCoefficient *coeff, Array< int > &attr)
Add a Dirichlet boundary condition to the velocity field.
void ComputeCurl2D(ParGridFunction &u, ParGridFunction &cu, bool assume_scalar=false)
Compute for .
ParFiniteElementSpace * vfes
Velocity finite element space.
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
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
ParGridFunction * GetCurrentVelocity()
Return a pointer to the current velocity ParGridFunction.
std::vector< AccelTerm_T > accel_terms
bool verbose
Enable/disable verbose output.
double a
Definition: lissajous.cpp:41
Solver wrapper which orthogonalizes the input and output vector.
Definition: solvers.hpp:1163
double kin_vis
Kinematic viscosity (dimensionless).
Class for parallel bilinear form using different test and trial FE spaces.
std::vector< VelDirichletBC_T > vel_dbcs
VectorGridFunctionCoefficient * FText_gfcoeff
void UpdateTimestepHistory(double dt)
Rotate entries in the time step and solution history arrays.
ParFiniteElementSpace * pfes
Pressure finite element space.
ParGridFunction * GetProvisionalVelocity()
Return a pointer to the provisional velocity ParGridFunction.
PresDirichletBC_T(Array< int > attr, Coefficient *coeff)
Class for parallel bilinear form.
AccelTerm_T(AccelTerm_T &&obj)
RefCoord t[3]
ParMixedBilinearForm * D_form
Vector data type.
Definition: vector.hpp:60
Vector coefficient defined by a vector GridFunction.
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:24
Base class for solvers.
Definition: operator.hpp:682
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Square Operator for imposing essential boundary conditions using only the action, Mult()...
Definition: operator.hpp:872
Abstract operator.
Definition: operator.hpp:24
bool debug
Enable/disable debug output.
Class for parallel meshes.
Definition: pmesh.hpp:32
Transient incompressible Navier Stokes solver in a split scheme formulation.
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
VelDirichletBC_T(Array< int > attr, VectorCoefficient *coeff)