MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
navier_kovasznay.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
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 #include "navier_solver.hpp"
43 #include <fstream>
44 
45 using namespace mfem;
46 using namespace navier;
47 
48 struct s_NavierContext
49 {
50  int ser_ref_levels = 1;
51  int order = 6;
52  double kinvis = 1.0 / 40.0;
53  double t_final = 10 * 0.001;
54  double dt = 0.001;
55  double reference_pressure = 0.0;
56  double reynolds = 1.0 / kinvis;
57  double lam = 0.5 * reynolds
58  - sqrt(0.25 * reynolds * reynolds + 4.0 * M_PI * M_PI);
59  bool pa = true;
60  bool ni = false;
61  bool visualization = false;
62  bool checkres = false;
63 } ctx;
64 
65 void vel_kovasznay(const Vector &x, double t, Vector &u)
66 {
67  double xi = x(0);
68  double yi = x(1);
69 
70  u(0) = 1.0 - exp(ctx.lam * xi) * cos(2.0 * M_PI * yi);
71  u(1) = ctx.lam / (2.0 * M_PI) * exp(ctx.lam * xi) * sin(2.0 * M_PI * yi);
72 }
73 
74 double pres_kovasznay(const Vector &x, double t)
75 {
76  double xi = x(0);
77 
78  return 0.5 * (1.0 - exp(2.0 * ctx.lam * xi)) + ctx.reference_pressure;
79 }
80 
81 int main(int argc, char *argv[])
82 {
83  Mpi::Init(argc, argv);
84  Hypre::Init();
85 
86  OptionsParser args(argc, argv);
87  args.AddOption(&ctx.ser_ref_levels,
88  "-rs",
89  "--refine-serial",
90  "Number of times to refine the mesh uniformly in serial.");
91  args.AddOption(&ctx.order,
92  "-o",
93  "--order",
94  "Order (degree) of the finite elements.");
95  args.AddOption(&ctx.dt, "-dt", "--time-step", "Time step.");
96  args.AddOption(&ctx.t_final, "-tf", "--final-time", "Final time.");
97  args.AddOption(&ctx.pa,
98  "-pa",
99  "--enable-pa",
100  "-no-pa",
101  "--disable-pa",
102  "Enable partial assembly.");
103  args.AddOption(&ctx.ni,
104  "-ni",
105  "--enable-ni",
106  "-no-ni",
107  "--disable-ni",
108  "Enable numerical integration rules.");
109  args.AddOption(&ctx.visualization,
110  "-vis",
111  "--visualization",
112  "-no-vis",
113  "--no-visualization",
114  "Enable or disable GLVis visualization.");
115  args.AddOption(
116  &ctx.checkres,
117  "-cr",
118  "--checkresult",
119  "-no-cr",
120  "--no-checkresult",
121  "Enable or disable checking of the result. Returns -1 on failure.");
122  args.Parse();
123  if (!args.Good())
124  {
125  if (Mpi::Root())
126  {
127  args.PrintUsage(mfem::out);
128  }
129  return 1;
130  }
131  if (Mpi::Root())
132  {
133  args.PrintOptions(mfem::out);
134  }
135 
136  Mesh mesh = Mesh::MakeCartesian2D(2, 4, Element::QUADRILATERAL, false, 1.5,
137  2.0);
138 
139  mesh.EnsureNodes();
140  GridFunction *nodes = mesh.GetNodes();
141  *nodes -= 0.5;
142 
143  for (int i = 0; i < ctx.ser_ref_levels; ++i)
144  {
145  mesh.UniformRefinement();
146  }
147 
148  if (Mpi::Root())
149  {
150  std::cout << "Number of elements: " << mesh.GetNE() << std::endl;
151  }
152 
153  auto *pmesh = new ParMesh(MPI_COMM_WORLD, mesh);
154  mesh.Clear();
155 
156  // Create the flow solver.
157  NavierSolver flowsolver(pmesh, ctx.order, ctx.kinvis);
158  flowsolver.EnablePA(ctx.pa);
159  flowsolver.EnableNI(ctx.ni);
160 
161  // Set the initial condition.
162  ParGridFunction *u_ic = flowsolver.GetCurrentVelocity();
163  VectorFunctionCoefficient u_excoeff(pmesh->Dimension(), vel_kovasznay);
164  u_ic->ProjectCoefficient(u_excoeff);
165 
167 
168  // Add Dirichlet boundary conditions to velocity space restricted to
169  // selected attributes on the mesh.
170  Array<int> attr(pmesh->bdr_attributes.Max());
171  attr = 1;
172  flowsolver.AddVelDirichletBC(vel_kovasznay, attr);
173 
174  double t = 0.0;
175  double dt = ctx.dt;
176  double t_final = ctx.t_final;
177  bool last_step = false;
178 
179  flowsolver.Setup(dt);
180 
181  double err_u = 0.0;
182  double err_p = 0.0;
183  ParGridFunction *u_gf = nullptr;
184  ParGridFunction *p_gf = nullptr;
185 
186  ParGridFunction p_ex_gf(flowsolver.GetCurrentPressure()->ParFESpace());
187  GridFunctionCoefficient p_ex_gf_coeff(&p_ex_gf);
188 
189  for (int step = 0; !last_step; ++step)
190  {
191  if (t + dt >= t_final - dt / 2)
192  {
193  last_step = true;
194  }
195 
196  flowsolver.Step(t, dt, step);
197 
198  // Compare against exact solution of velocity and pressure.
199  u_gf = flowsolver.GetCurrentVelocity();
200  p_gf = flowsolver.GetCurrentPressure();
201 
202  u_excoeff.SetTime(t);
203  p_excoeff.SetTime(t);
204 
205  // Remove mean value from exact pressure solution.
206  p_ex_gf.ProjectCoefficient(p_excoeff);
207  flowsolver.MeanZero(p_ex_gf);
208 
209  err_u = u_gf->ComputeL2Error(u_excoeff);
210  err_p = p_gf->ComputeL2Error(p_ex_gf_coeff);
211 
212  double cfl = flowsolver.ComputeCFL(*u_gf, dt);
213 
214  if (Mpi::Root())
215  {
216  printf("%5s %8s %8s %8s %11s %11s\n",
217  "Order",
218  "CFL",
219  "Time",
220  "dt",
221  "err_u",
222  "err_p");
223  printf("%5.2d %8.2E %.2E %.2E %.5E %.5E err\n",
224  ctx.order,
225  cfl,
226  t,
227  dt,
228  err_u,
229  err_p);
230  fflush(stdout);
231  }
232  }
233 
234  if (ctx.visualization)
235  {
236  char vishost[] = "localhost";
237  int visport = 19916;
238  socketstream sol_sock(vishost, visport);
239  sol_sock.precision(8);
240  sol_sock << "parallel " << Mpi::WorldSize() << " "
241  << Mpi::WorldRank() << "\n";
242  sol_sock << "solution\n" << *pmesh << *u_ic << std::flush;
243  }
244 
245  flowsolver.PrintTimingData();
246 
247  // Test if the result for the test run is as expected.
248  if (ctx.checkres)
249  {
250  double tol_u = 1e-6;
251  double tol_p = 1e-5;
252  if (err_u > tol_u || err_p > tol_p)
253  {
254  if (Mpi::Root())
255  {
256  mfem::out << "Result has a larger error than expected."
257  << std::endl;
258  }
259  return -1;
260  }
261  }
262 
263  delete pmesh;
264 
265  return 0;
266 }
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:80
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:923
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:525
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:50
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:9816
ParGridFunction * GetCurrentPressure()
Return a pointer to the current pressure ParGridFunction.
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
Definition: pgridfunc.hpp:281
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
static bool Root()
Return true if the rank in MPI_COMM_WORLD is zero.
A general vector function coefficient.
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:3713
ParGridFunction * GetCurrentVelocity()
Return a pointer to the current 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:899
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7908
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:5291
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:113
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150