MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
navier_turbchan.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 // DNS of a channel flow at Re_tau = 180 (variable). A detailed description of
13 // the test case can be found at [1]. Like described in the reference, the
14 // initial condition is based on the Reichardt function.
15 //
16 // [1] https://how5.cenaero.be/content/ws2-les-plane-channel-ret550
17 
18 #include "navier_solver.hpp"
19 #include <fstream>
20 #include <cmath>
21 
22 using namespace mfem;
23 using namespace navier;
24 
25 struct s_NavierContext
26 {
27  int order = 5;
28  double Re_tau = 180.0;
29  double kin_vis = 1.0 / Re_tau;
30  double t_final = 50.0;
31  double dt = -1.0;
32 } ctx;
33 
34 double mesh_stretching_func(const double y)
35 {
36  double C = 1.8;
37  double delta = 1.0;
38 
39  return delta * tanh(C * (2.0 * y - 1.0)) / tanh(C);
40 }
41 
42 void accel(const Vector &x, double t, Vector &f)
43 {
44  f(0) = 1.0;
45  f(1) = 0.0;
46  f(2) = 0.0;
47 }
48 
49 void vel_ic_reichardt(const Vector &coords, double t, Vector &u)
50 {
51  double yp;
52  double x = coords(0);
53  double y = coords(1);
54  double z = coords(2);
55 
56  double C = 5.17;
57  double k = 0.4;
58  double eps = 1e-2;
59 
60  if (y < 0)
61  {
62  yp = (1.0 + y) * ctx.Re_tau;
63  }
64  else
65  {
66  yp = (1.0 - y) * ctx.Re_tau;
67  }
68 
69  u(0) = 1.0 / k * log(1.0 + k * yp) + (C - (1.0 / k) * log(k)) * (1 - exp(
70  -yp / 11.0) - yp / 11.0 * exp(-yp / 3.0));
71 
72  double kx = 23.0;
73  double kz = 13.0;
74 
75  double alpha = kx * 2.0 * M_PI / 2.0 * M_PI;
76  double beta = kz * 2.0 * M_PI / M_PI;
77 
78  u(0) += eps * beta * sin(alpha * x) * cos(beta * z);
79  u(1) = eps * sin(alpha * x) * sin(beta * z);
80  u(2) = -eps * alpha * cos(alpha * x) * sin(beta * z);
81 }
82 
83 void vel_wall(const Vector &x, double t, Vector &u)
84 {
85  u(0) = 0.0;
86  u(1) = 0.0;
87  u(2) = 0.0;
88 }
89 
90 int main(int argc, char *argv[])
91 {
92  Mpi::Init();
93  Hypre::Init();
94 
95  double Lx = 2.0 * M_PI;
96  double Ly = 1.0;
97  double Lz = M_PI;
98 
99  int N = ctx.order + 1;
100  int NL = static_cast<int>(std::round(64.0 / N)); // Coarse
101  // int NL = std::round(96.0 / N); // Baseline
102  // int NL = std::round(128.0 / N); // Fine
103  double LC = M_PI / NL;
104  int NX = 2 * NL;
105  int NY = 2 * static_cast<int>(std::round(48.0 / N));
106  int NZ = NL;
107 
108  Mesh mesh = Mesh::MakeCartesian3D(NX, NY, NZ, Element::HEXAHEDRON, Lx, Ly, Lz);
109 
110  for (int i = 0; i < mesh.GetNV(); ++i)
111  {
112  double *v = mesh.GetVertex(i);
113  v[1] = mesh_stretching_func(v[1]);
114  }
115 
116  // Create translation vectors defining the periodicity
117  Vector x_translation({Lx, 0.0, 0.0});
118  Vector z_translation({0.0, 0.0, Lz});
119  std::vector<Vector> translations = {x_translation, z_translation};
120 
121  // Create the periodic mesh using the vertex mapping defined by the translation vectors
122  Mesh periodic_mesh = Mesh::MakePeriodic(mesh,
123  mesh.CreatePeriodicVertexMapping(translations));
124 
125  if (Mpi::Root())
126  {
127  printf("NL=%d NX=%d NY=%d NZ=%d dx+=%f\n", NL, NX, NY, NZ, LC * ctx.Re_tau);
128  std::cout << "Number of elements: " << mesh.GetNE() << std::endl;
129  }
130 
131  double hmin, hmax, kappa_min, kappa_max;
132  periodic_mesh.GetCharacteristics(hmin, hmax, kappa_min, kappa_max);
133 
134  double umax = 22.0;
135  ctx.dt = 1.0 / pow(ctx.order, 1.5) * hmin / umax;
136 
137  auto *pmesh = new ParMesh(MPI_COMM_WORLD, periodic_mesh);
138 
139  // Create the flow solver.
140  NavierSolver flowsolver(pmesh, ctx.order, ctx.kin_vis);
141  flowsolver.EnablePA(true);
142 
143  // Set the initial condition.
144  ParGridFunction *u_gf = flowsolver.GetCurrentVelocity();
145  ParGridFunction *p_gf = flowsolver.GetCurrentPressure();
146 
147  VectorFunctionCoefficient u_ic_coef(pmesh->Dimension(), vel_ic_reichardt);
148  u_gf->ProjectCoefficient(u_ic_coef);
149 
150  Array<int> domain_attr(pmesh->attributes);
151  domain_attr = 1;
152  flowsolver.AddAccelTerm(accel, domain_attr);
153 
154  Array<int> attr(pmesh->bdr_attributes.Max());
155  attr[1] = 1;
156  attr[3] = 1;
157  flowsolver.AddVelDirichletBC(vel_wall, attr);
158 
159  double t = 0.0;
160  double dt = ctx.dt;
161  double t_final = ctx.t_final;
162  bool last_step = false;
163 
164  flowsolver.Setup(dt);
165 
166  ParaViewDataCollection pvdc("turbchan", pmesh);
168  pvdc.SetHighOrderOutput(true);
169  pvdc.SetLevelsOfDetail(ctx.order);
170  pvdc.SetCycle(0);
171  pvdc.SetTime(t);
172  pvdc.RegisterField("velocity", u_gf);
173  pvdc.RegisterField("pressure", p_gf);
174  pvdc.Save();
175 
176  for (int step = 0; !last_step; ++step)
177  {
178  if (t + dt >= t_final - dt / 2)
179  {
180  last_step = true;
181  }
182 
183  flowsolver.Step(t, dt, step);
184 
185  if (step % 1000 == 0)
186  {
187  pvdc.SetCycle(step);
188  pvdc.SetTime(t);
189  pvdc.Save();
190  }
191 
192  if (t > 5.0)
193  {
194  dt = 1e-2;
195  }
196 
197  if (Mpi::Root())
198  {
199  printf("%11s %11s\n", "Time", "dt");
200  printf("%.5E %.5E\n", t, dt);
201  fflush(stdout);
202  }
203  }
204 
205  flowsolver.PrintTimingData();
206 
207  delete pmesh;
208 
209  return 0;
210 }
static Mesh MakePeriodic(const Mesh &orig_mesh, const std::vector< int > &v2v)
Create a periodic mesh by identifying vertices of orig_mesh.
Definition: mesh.cpp:4891
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
const double * GetVertex(int i) const
Return pointer to vertex i&#39;s coordinates.
Definition: mesh.hpp:1012
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void SetDataFormat(VTKFormat fmt)
void GetCharacteristics(double &h_min, double &h_max, double &kappa_min, double &kappa_max, Vector *Vh=NULL, Vector *Vk=NULL)
Definition: mesh.cpp:193
Helper class for ParaView visualization data.
static Mesh MakeCartesian3D(int nx, int ny, int nz, Element::Type type, double sx=1.0, double sy=1.0, double sz=1.0, bool sfc_ordering=true)
Definition: mesh.cpp:3723
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
void Setup(double dt)
Initialize forms, solvers and preconditioners.
std::vector< int > CreatePeriodicVertexMapping(const std::vector< Vector > &translations, double tol=1e-8) const
Creates a mapping v2v from the vertex indices of the mesh such that coincident vertices under the giv...
Definition: mesh.cpp:4925
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
ParGridFunction * GetCurrentPressure()
Return a pointer to the current pressure ParGridFunction.
double delta
Definition: lissajous.cpp:43
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 SetHighOrderOutput(bool high_order_output_)
static void Init()
Singleton creation with Mpi::Init();.
void SetTime(double t)
Set physical time (for time-dependent simulations)
static bool Root()
Return true if the rank in MPI_COMM_WORLD is zero.
A general vector function coefficient.
ParGridFunction * GetCurrentVelocity()
Return a pointer to the current velocity ParGridFunction.
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:920
RefCoord t[3]
void SetLevelsOfDetail(int levels_of_detail_)
const double alpha
Definition: ex15.cpp:369
Vector data type.
Definition: vector.hpp:60
void Step(double &time, double dt, int cur_step, bool provisional=false)
Compute solution at the next time step t+dt.
virtual void Save() override
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Class for parallel meshes.
Definition: pmesh.hpp:32
Transient incompressible Navier Stokes solver in a split scheme formulation.
int main()
void AddAccelTerm(VectorCoefficient *coeff, Array< int > &attr)
Add an acceleration term to the RHS of the equation.