MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
navier_solver.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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
19namespace mfem
20{
21namespace navier
22{
23using VecFuncT = void(const Vector &x, real_t t, Vector &u);
24using ScalarFuncT = real_t(const Vector &x, real_t t);
25
26/// Container for a Dirichlet boundary condition of the velocity field.
28{
29public:
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
45
48};
49
50/// Container for a Dirichlet boundary condition of the pressure field.
52{
53public:
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
69
72};
73
74/// Container for an acceleration term.
76{
77public:
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 $M_v^{-1}$. On
108 * the other hand this introduces a CFL stability condition on the maximum
109 * timestep.
110 *
111 * 2. A Poisson solve $S_p^{-1}$.
112 *
113 * 3. A Helmholtz like solve $(M_v - \partial t K_v)^{-1}$.
114 *
115 * The numerical solver setup for each step are as follows.
116 *
117 * $M_v^{-1}$ is solved using CG with Jacobi as preconditioner.
118 *
119 * $S_p^{-1}$ 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 * $(M_v - \partial t K_v)^{-1}$ 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{
143public:
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 * $(P_N)^d P_N$ 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 */
157
158 /// Initialize forms, solvers and preconditioners.
159 void Setup(real_t 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(real_t &time, real_t 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.
196
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
241
242 /// Compute $\nabla \times \nabla \times u$ for $u \in (H^1)^2$.
244 ParGridFunction &cu,
245 bool assume_scalar = false);
246
247 /// Compute $\nabla \times \nabla \times u$ for $u \in (H^1)^3$.
249
250 /// Remove mean from a Vector.
251 /**
252 * Modify the Vector @a v by subtracting its mean using
253 * $v = v - \frac{\sum_i^N v_i}{N}$
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 * $v = v - \int_\Omega \frac{v}{vol(\Omega)} dx$.
261 */
262 void MeanZero(ParGridFunction &v);
263
264 /// Rotate entries in the time step and solution history arrays.
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
272
273 /// Set the number of modes to cut off in the interpolation filter
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 */
285
286protected:
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).
330
332
333 /// Velocity $H^1$ finite element collection.
335
336 /// Pressure $H^1$ finite element collection.
338
339 /// Velocity $(H^1)^d$ finite element space.
341
342 /// Pressure $H^1$ finite element space.
344
345 ParNonlinearForm *N = nullptr;
346
348
350
352
354
356
358
360
362
364
365 /// Linear form to compute the mass matrix in various subroutines.
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
394
396
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
420 int cur_step = 0;
421 std::vector<real_t> dthist = {0.0, 0.0, 0.0};
422
423 // BDFk/EXTk coefficients.
424 real_t bd0 = 0.0;
425 real_t bd1 = 0.0;
426 real_t bd2 = 0.0;
427 real_t bd3 = 0.0;
428 real_t ab1 = 0.0;
429 real_t ab2 = 0.0;
430 real_t 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#if defined(MFEM_USE_DOUBLE)
446#elif defined(MFEM_USE_SINGLE)
447 real_t rtol_mvsolve = 1e-9;
448 real_t rtol_spsolve = 1e-5;
449 real_t rtol_hsolve = 1e-7;
450#else
451#error "Only single and double precision are supported!"
452 real_t rtol_mvsolve = 1e-12;
453 real_t rtol_spsolve = 1e-6;
454 real_t rtol_hsolve = 1e-8;
455#endif
456
457 // Iteration counts.
459
460 // Residuals.
462
463 // LOR related.
465
466 // Filter-based stabilization
473};
474
475} // namespace navier
476
477} // namespace mfem
478
479#endif
Definition solvers.hpp:513
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
A coefficient that is constant across space and time.
Square Operator for imposing essential boundary conditions using only the action, Mult(),...
Definition operator.hpp:895
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
Container class for integration rules.
Definition intrules.hpp:416
Pointer to an Operator of a specified type.
Definition handle.hpp:34
Abstract operator.
Definition operator.hpp:25
Solver wrapper which orthogonalizes the input and output vector.
Definition solvers.hpp:1218
Class for parallel bilinear form.
Abstract parallel finite element space.
Definition pfespace.hpp:29
Class for parallel grid function.
Definition pgridfunc.hpp:33
Create and assemble a low-order refined version of a ParBilinearForm.
Definition lor.hpp:166
Class for parallel linear form.
Class for parallel meshes.
Definition pmesh.hpp:34
Class for parallel bilinear form using different test and trial FE spaces.
Parallel non-linear operator on the true dofs.
Base class for solvers.
Definition operator.hpp:683
Timing object.
Definition tic_toc.hpp:36
Base class for vector Coefficients that optionally depend on time and space.
Vector coefficient defined by a vector GridFunction.
Vector data type.
Definition vector.hpp:80
Container for an acceleration term.
AccelTerm_T(AccelTerm_T &&obj)
AccelTerm_T(Array< int > attr, VectorCoefficient *coeff)
VectorCoefficient * coeff
Transient incompressible Navier Stokes solver in a split scheme formulation.
ParGridFunction * GetProvisionalVelocity()
Return a pointer to the provisional velocity ParGridFunction.
void Setup(real_t dt)
Initialize forms, solvers and preconditioners.
ParFiniteElementSpace * vfes_filter
real_t kin_vis
Kinematic viscosity (dimensionless).
real_t ComputeCFL(ParGridFunction &u, real_t dt)
Compute CFL.
ParLORDiscretization * lor
bool verbose
Enable/disable verbose output.
void AddVelDirichletBC(VectorCoefficient *coeff, Array< int > &attr)
Add a Dirichlet boundary condition to the velocity field.
ConstantCoefficient Sp_coeff
FiniteElementCollection * vfec
Velocity finite element collection.
std::vector< PresDirichletBC_T > pres_dbcs
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.
VectorGridFunctionCoefficient * FText_gfcoeff
void SetMaxBDFOrder(int maxbdforder)
Set the maximum order to use for the BDF method.
int order
The order of the velocity and pressure space.
bool partial_assembly
Enable/disable partial assembly of forms.
ParGridFunction * GetCurrentPressure()
Return a pointer to the current pressure ParGridFunction.
void PrintTimingData()
Print timing summary of the solving routine.
void AddPresDirichletBC(Coefficient *coeff, Array< int > &attr)
Add a Dirichlet boundary condition to the pressure field.
FiniteElementCollection * vfec_filter
ConstantCoefficient onecoeff
void EnablePA(bool pa)
Enable partial assembly for every operator.
ParFiniteElementSpace * pfes
Pressure finite element space.
ConstantCoefficient nlcoeff
void PrintInfo()
Print information about the Navier version.
void SetFilterAlpha(real_t a)
Set the interpolation filter parameter a.
std::vector< AccelTerm_T > accel_terms
void SetCutoffModes(int c)
Set the number of modes to cut off in the interpolation filter.
void MeanZero(ParGridFunction &v)
Remove the mean from a ParGridFunction.
void ComputeCurl2D(ParGridFunction &u, ParGridFunction &cu, bool assume_scalar=false)
Compute for .
ParGridFunction * GetCurrentVelocity()
Return a pointer to the current velocity ParGridFunction.
bool numerical_integ
Enable/disable numerical integration rules of forms.
ParMixedBilinearForm * D_form
void UpdateTimestepHistory(real_t dt)
Rotate entries in the time step and solution history arrays.
ParLinearForm * mass_lf
Linear form to compute the mass matrix in various subroutines.
bool debug
Enable/disable debug output.
ConstantCoefficient H_lincoeff
ParMesh * pmesh
The parallel mesh.
ConstantCoefficient H_bdfcoeff
void AddAccelTerm(VectorCoefficient *coeff, Array< int > &attr)
Add an acceleration term to the RHS of the equation.
std::vector< VelDirichletBC_T > vel_dbcs
void Step(real_t &time, real_t dt, int cur_step, bool provisional=false)
Compute solution at the next time step t+dt.
ParMixedBilinearForm * G_form
FiniteElementCollection * pfec
Pressure finite element collection.
void SetTimeIntegrationCoefficients(int step)
Update the EXTk/BDF time integration coefficient.
void Orthogonalize(Vector &v)
Remove mean from a Vector.
ParFiniteElementSpace * vfes
Velocity finite element space.
void ComputeCurl3D(ParGridFunction &u, ParGridFunction &cu)
Compute for .
NavierSolver(ParMesh *mesh, int order, real_t kin_vis)
Initialize data structures, set FE space order and kinematic viscosity.
std::vector< real_t > dthist
Container for a Dirichlet boundary condition of the pressure field.
PresDirichletBC_T(Array< int > attr, Coefficient *coeff)
PresDirichletBC_T(PresDirichletBC_T &&obj)
Container for a Dirichlet boundary condition of the velocity field.
VelDirichletBC_T(VelDirichletBC_T &&obj)
VelDirichletBC_T(Array< int > attr, VectorCoefficient *coeff)
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
real_t(const Vector &x, real_t t) ScalarFuncT
void(const Vector &x, real_t t, Vector &u) VecFuncT
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
RefCoord t[3]