MFEM  v4.5.1 Finite element discretization library
navier_solver.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 #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:465
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:1579
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:1135
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:655
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Square Operator for imposing essential boundary conditions using only the action, Mult()...
Definition: operator.hpp:845
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)