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