MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
navier_kovasznay_vs.cpp
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 // Navier Kovasznay example with variable time step
13 //
14 // Solve for the steady Kovasznay flow at Re = 40 defined by
15 //
16 // u = [1 - exp(L * x) * cos(2 * pi * y),
17 // L / (2 * pi) * exp(L * x) * sin(2 * pi * y)],
18 //
19 // p = 1/2 * (1 - exp(2 * L * x)),
20 //
21 // with L = Re/2 - sqrt(Re^2/4 + 4 * pi^2).
22 //
23 // The problem domain is set up like this
24 //
25 // +-------------+
26 // | |
27 // | |
28 // | |
29 // | |
30 // Inflow -> | | -> Outflow
31 // | |
32 // | |
33 // | |
34 // | |
35 // | |
36 // +-------------+
37 //
38 // and Dirichlet boundary conditions are applied for the velocity on every
39 // boundary. The problem, although steady state, is time integrated up to the
40 // final time and the solution is compared with the known exact solution.
41 //
42 // Additionally, this example shows the usage of variable time steps with
43 // Navier. A basic sample algorithm for determining the next time step is
44 // provided based on a CFL restriction on the velocity.
45 
46 #include "navier_solver.hpp"
47 #include <fstream>
48 
49 using namespace mfem;
50 using namespace navier;
51 
52 struct s_NavierContext
53 {
54  int ser_ref_levels = 1;
55  int order = 6;
56  double kinvis = 1.0 / 40.0;
57  double t_final = 10 * 0.001;
58  double dt = 0.001;
59  double reference_pressure = 0.0;
60  double reynolds = 1.0 / kinvis;
61  double lam = 0.5 * reynolds
62  - sqrt(0.25 * reynolds * reynolds + 4.0 * M_PI * M_PI);
63  bool pa = true;
64  bool ni = false;
65  bool visualization = false;
66  bool checkres = false;
67 } ctx;
68 
69 void vel_kovasznay(const Vector &x, double t, Vector &u)
70 {
71  double xi = x(0);
72  double yi = x(1);
73 
74  u(0) = 1.0 - exp(ctx.lam * xi) * cos(2.0 * M_PI * yi);
75  u(1) = ctx.lam / (2.0 * M_PI) * exp(ctx.lam * xi) * sin(2.0 * M_PI * yi);
76 }
77 
78 double pres_kovasznay(const Vector &x, double t)
79 {
80  double xi = x(0);
81 
82  return 0.5 * (1.0 - exp(2.0 * ctx.lam * xi)) + ctx.reference_pressure;
83 }
84 
85 int main(int argc, char *argv[])
86 {
87  Mpi::Init(argc, argv);
88  Hypre::Init();
89 
90  OptionsParser args(argc, argv);
91  args.AddOption(&ctx.ser_ref_levels,
92  "-rs",
93  "--refine-serial",
94  "Number of times to refine the mesh uniformly in serial.");
95  args.AddOption(&ctx.order,
96  "-o",
97  "--order",
98  "Order (degree) of the finite elements.");
99  args.AddOption(&ctx.dt, "-dt", "--time-step", "Time step.");
100  args.AddOption(&ctx.t_final, "-tf", "--final-time", "Final time.");
101  args.AddOption(&ctx.pa,
102  "-pa",
103  "--enable-pa",
104  "-no-pa",
105  "--disable-pa",
106  "Enable partial assembly.");
107  args.AddOption(&ctx.ni,
108  "-ni",
109  "--enable-ni",
110  "-no-ni",
111  "--disable-ni",
112  "Enable numerical integration rules.");
113  args.AddOption(&ctx.visualization,
114  "-vis",
115  "--visualization",
116  "-no-vis",
117  "--no-visualization",
118  "Enable or disable GLVis visualization.");
119  args.AddOption(
120  &ctx.checkres,
121  "-cr",
122  "--checkresult",
123  "-no-cr",
124  "--no-checkresult",
125  "Enable or disable checking of the result. Returns -1 on failure.");
126  args.Parse();
127  if (!args.Good())
128  {
129  if (Mpi::Root())
130  {
131  args.PrintUsage(mfem::out);
132  }
133  return 1;
134  }
135  if (Mpi::Root())
136  {
137  args.PrintOptions(mfem::out);
138  }
139 
140  Mesh mesh = Mesh::MakeCartesian2D(2, 4, Element::QUADRILATERAL, false, 1.5,
141  2.0);
142 
143  mesh.EnsureNodes();
144  GridFunction *nodes = mesh.GetNodes();
145  *nodes -= 0.5;
146 
147  for (int i = 0; i < ctx.ser_ref_levels; ++i)
148  {
149  mesh.UniformRefinement();
150  }
151 
152  if (Mpi::Root())
153  {
154  std::cout << "Number of elements: " << mesh.GetNE() << std::endl;
155  }
156 
157  auto *pmesh = new ParMesh(MPI_COMM_WORLD, mesh);
158  mesh.Clear();
159 
160  // Create the flow solver.
161  NavierSolver flowsolver(pmesh, ctx.order, ctx.kinvis);
162  flowsolver.EnablePA(ctx.pa);
163  flowsolver.EnableNI(ctx.ni);
164 
165  // Set the initial condition.
166  ParGridFunction *u_ic = flowsolver.GetCurrentVelocity();
167  VectorFunctionCoefficient u_excoeff(pmesh->Dimension(), vel_kovasznay);
168  u_ic->ProjectCoefficient(u_excoeff);
169 
171 
172  // Add Dirichlet boundary conditions to velocity space restricted to
173  // selected attributes on the mesh.
174  Array<int> attr(pmesh->bdr_attributes.Max());
175  attr = 1;
176  flowsolver.AddVelDirichletBC(vel_kovasznay, attr);
177 
178  double t = 0.0;
179  double dt = ctx.dt;
180  double t_final = ctx.t_final;
181  bool last_step = false;
182 
183  flowsolver.Setup(dt);
184 
185  double err_u = 0.0;
186  double err_p = 0.0;
187  ParGridFunction *u_next_gf = nullptr;
188  ParGridFunction *u_gf = nullptr;
189  ParGridFunction *p_gf = nullptr;
190 
191  ParGridFunction p_ex_gf(flowsolver.GetCurrentPressure()->ParFESpace());
192  GridFunctionCoefficient p_ex_gf_coeff(&p_ex_gf);
193 
194  double cfl_max = 0.8;
195  double cfl_tol = 1e-4;
196 
197  for (int step = 0; !last_step; ++step)
198  {
199  if (t + dt >= t_final - dt / 2)
200  {
201  last_step = true;
202  }
203 
204  // Take a provisional step
205  flowsolver.Step(t, dt, step, true);
206 
207  // Retrieve the computed provisional velocity
208  u_next_gf = flowsolver.GetProvisionalVelocity();
209 
210  // Compute the CFL based on the provisional velocity
211  double cfl = flowsolver.ComputeCFL(*u_next_gf, dt);
212 
213  double error_est = cfl / (cfl_max + cfl_tol);
214  if (error_est >= 1.0)
215  {
216  // Reject the time step
217  if (Mpi::Root())
218  {
219  std::cout
220  << "Step reached maximum CFL, retrying with smaller step size..."
221  << std::endl;
222  }
223  dt *= 0.5;
224  step -= 1;
225  }
226  else
227  {
228  // Accept the time step
229  t += dt;
230 
231  // Predict new step size
232  double fac_safety = 2.0;
233  double eta = pow(1.0 / (fac_safety * error_est), 1.0 / (1.0 + 3.0));
234  double fac_min = 0.1;
235  double fac_max = 1.4;
236  dt = dt * std::min(fac_max, std::max(fac_min, eta));
237 
238  // Queue new time step in the history array
239  flowsolver.UpdateTimestepHistory(dt);
240  }
241 
242  u_gf = flowsolver.GetCurrentVelocity();
243  p_gf = flowsolver.GetCurrentPressure();
244 
245  u_excoeff.SetTime(t);
246  p_excoeff.SetTime(t);
247 
248  // Remove mean value from exact pressure solution.
249  p_ex_gf.ProjectCoefficient(p_excoeff);
250  flowsolver.MeanZero(p_ex_gf);
251 
252  err_u = u_gf->ComputeL2Error(u_excoeff);
253  err_p = p_gf->ComputeL2Error(p_ex_gf_coeff);
254 
255  if (Mpi::Root())
256  {
257  printf("%5s %8s %8s %8s %11s %11s\n",
258  "Order",
259  "CFL",
260  "Time",
261  "dt",
262  "err_u",
263  "err_p");
264  printf("%5.2d %8.2E %.2E %.2E %.5E %.5E err\n",
265  ctx.order,
266  cfl,
267  t,
268  dt,
269  err_u,
270  err_p);
271  fflush(stdout);
272  }
273  }
274 
275  if (ctx.visualization)
276  {
277  char vishost[] = "localhost";
278  int visport = 19916;
279  socketstream sol_sock(vishost, visport);
280  sol_sock.precision(8);
281  sol_sock << "parallel " << Mpi::WorldSize() << " "
282  << Mpi::WorldRank() << "\n";
283  sol_sock << "solution\n" << *pmesh << *u_ic << std::flush;
284  }
285 
286  flowsolver.PrintTimingData();
287 
288  // Test if the result for the test run is as expected.
289  if (ctx.checkres)
290  {
291  double tol_u = 1e-6;
292  double tol_p = 1e-5;
293  if (err_u > tol_u || err_p > tol_p)
294  {
295  if (Mpi::Root())
296  {
297  mfem::out << "Result has a larger error than expected."
298  << std::endl;
299  }
300  return -1;
301  }
302  }
303 
304  delete pmesh;
305 
306  return 0;
307 }
const char vishost[]
void PrintTimingData()
Print timing summary of the solving routine.
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(), hypre will be finalized automatically at program exit.
Definition: hypre.hpp:67
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:928
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:525
FDualNumber< tbase > exp(const FDualNumber< tbase > &f)
exp([dual number])
Definition: fdual.hpp:499
static int WorldSize()
Return the size of MPI_COMM_WORLD.
virtual void SetTime(double t)
Set the time for time dependent coefficients.
Definition: coefficient.hpp:48
void Setup(double dt)
Initialize forms, solvers and preconditioners.
void MeanZero(ParGridFunction &v)
Remove the mean from a ParGridFunction.
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
double ComputeCFL(ParGridFunction &u, double dt)
Compute CFL.
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9737
ParGridFunction * GetCurrentPressure()
Return a pointer to the current pressure ParGridFunction.
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.
static void Init()
Singleton creation with Mpi::Init();.
const int visport
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
FDualNumber< tbase > cos(const FDualNumber< tbase > &f)
cos([dual number])
Definition: fdual.hpp:471
static bool Root()
Return true if the rank in MPI_COMM_WORLD is zero.
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
Definition: fdual.hpp:543
A general vector function coefficient.
FDualNumber< tbase > sin(const FDualNumber< tbase > &f)
sin([dual number])
Definition: fdual.hpp:572
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
static Mesh MakeCartesian2D(int nx, int ny, Element::Type type, bool generate_edges=false, double sx=1.0, double sy=1.0, bool sfc_ordering=true)
Definition: mesh.cpp:3701
ParGridFunction * GetCurrentVelocity()
Return a pointer to the current velocity ParGridFunction.
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const
Definition: pgridfunc.hpp:281
void UpdateTimestepHistory(double dt)
Rotate entries in the time step and solution history arrays.
FDualNumber< tbase > sqrt(const FDualNumber< tbase > &f)
sqrt([dual number])
Definition: fdual.hpp:600
ParGridFunction * GetProvisionalVelocity()
Return a pointer to the provisional velocity ParGridFunction.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
RefCoord t[3]
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:904
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7841
void Step(double &time, double dt, int cur_step, bool provisional=false)
Compute solution at the next time step t+dt.
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
Class for parallel grid function.
Definition: pgridfunc.hpp:32
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
Class for parallel meshes.
Definition: pmesh.hpp:32
Transient incompressible Navier Stokes solver in a split scheme formulation.
int main()
void EnsureNodes()
Definition: mesh.cpp:5279
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:113
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150