MFEM  v4.5.2
Finite element discretization library
navier_tgv.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 // 3D Taylor-Green vortex benchmark example at Re=1600
13 // Unsteady flow of a decaying vortex is computed and compared against a known,
14 // analytical solution.
15 
16 #include "navier_solver.hpp"
17 #include <fstream>
18 
19 using namespace mfem;
20 using namespace navier;
21 
22 struct s_NavierContext
23 {
24  int element_subdivisions = 1;
25  int order = 4;
26  double kinvis = 1.0 / 1600.0;
27  double t_final = 10 * 1e-3;
28  double dt = 1e-3;
29  bool pa = true;
30  bool ni = false;
31  bool visualization = false;
32  bool checkres = false;
33 } ctx;
34 
35 void vel_tgv(const Vector &x, double t, Vector &u)
36 {
37  double xi = x(0);
38  double yi = x(1);
39  double zi = x(2);
40 
41  u(0) = sin(xi) * cos(yi) * cos(zi);
42  u(1) = -cos(xi) * sin(yi) * cos(zi);
43  u(2) = 0.0;
44 }
45 
46 class QuantitiesOfInterest
47 {
48 public:
49  QuantitiesOfInterest(ParMesh *pmesh)
50  {
51  H1_FECollection h1fec(1);
52  ParFiniteElementSpace h1fes(pmesh, &h1fec);
53 
54  onecoeff.constant = 1.0;
55  mass_lf = new ParLinearForm(&h1fes);
56  mass_lf->AddDomainIntegrator(new DomainLFIntegrator(onecoeff));
57  mass_lf->Assemble();
58 
59  ParGridFunction one_gf(&h1fes);
60  one_gf.ProjectCoefficient(onecoeff);
61 
62  volume = mass_lf->operator()(one_gf);
63  };
64 
65  double ComputeKineticEnergy(ParGridFunction &v)
66  {
67  Vector velx, vely, velz;
68  double integ = 0.0;
69  const FiniteElement *fe;
71  FiniteElementSpace *fes = v.FESpace();
72 
73  for (int i = 0; i < fes->GetNE(); i++)
74  {
75  fe = fes->GetFE(i);
76  int intorder = 2 * fe->GetOrder();
77  const IntegrationRule *ir = &IntRules.Get(fe->GetGeomType(), intorder);
78 
79  v.GetValues(i, *ir, velx, 1);
80  v.GetValues(i, *ir, vely, 2);
81  v.GetValues(i, *ir, velz, 3);
82 
83  T = fes->GetElementTransformation(i);
84  for (int j = 0; j < ir->GetNPoints(); j++)
85  {
86  const IntegrationPoint &ip = ir->IntPoint(j);
87  T->SetIntPoint(&ip);
88 
89  double vel2 = velx(j) * velx(j) + vely(j) * vely(j)
90  + velz(j) * velz(j);
91 
92  integ += ip.weight * T->Weight() * vel2;
93  }
94  }
95 
96  double global_integral = 0.0;
97  MPI_Allreduce(&integ,
98  &global_integral,
99  1,
100  MPI_DOUBLE,
101  MPI_SUM,
102  MPI_COMM_WORLD);
103 
104  return 0.5 * global_integral / volume;
105  };
106 
107  ~QuantitiesOfInterest() { delete mass_lf; };
108 
109 private:
110  ConstantCoefficient onecoeff;
111  ParLinearForm *mass_lf;
112  double volume;
113 };
114 
115 template<typename T>
116 T sq(T x)
117 {
118  return x * x;
119 }
120 
121 // Computes Q = 0.5*(tr(\nabla u)^2 - tr(\nabla u \cdot \nabla u))
123 {
124  FiniteElementSpace *v_fes = u.FESpace();
125  FiniteElementSpace *fes = q.FESpace();
126 
127  // AccumulateAndCountZones
128  Array<int> zones_per_vdof;
129  zones_per_vdof.SetSize(fes->GetVSize());
130  zones_per_vdof = 0;
131 
132  q = 0.0;
133 
134  // Local interpolation
135  int elndofs;
136  Array<int> v_dofs, dofs;
137  Vector vals;
138  Vector loc_data;
139  int vdim = v_fes->GetVDim();
140  DenseMatrix grad_hat;
141  DenseMatrix dshape;
142  DenseMatrix grad;
143 
144  for (int e = 0; e < fes->GetNE(); ++e)
145  {
146  fes->GetElementVDofs(e, dofs);
147  v_fes->GetElementVDofs(e, v_dofs);
148  u.GetSubVector(v_dofs, loc_data);
149  vals.SetSize(dofs.Size());
151  const FiniteElement *el = fes->GetFE(e);
152  elndofs = el->GetDof();
153  int dim = el->GetDim();
154  dshape.SetSize(elndofs, dim);
155 
156  for (int dof = 0; dof < elndofs; ++dof)
157  {
158  // Project
159  const IntegrationPoint &ip = el->GetNodes().IntPoint(dof);
160  tr->SetIntPoint(&ip);
161 
162  // Eval
163  // GetVectorGradientHat
164  el->CalcDShape(tr->GetIntPoint(), dshape);
165  grad_hat.SetSize(vdim, dim);
166  DenseMatrix loc_data_mat(loc_data.GetData(), elndofs, vdim);
167  MultAtB(loc_data_mat, dshape, grad_hat);
168 
169  const DenseMatrix &Jinv = tr->InverseJacobian();
170  grad.SetSize(grad_hat.Height(), Jinv.Width());
171  Mult(grad_hat, Jinv, grad);
172 
173  double q_val = 0.5 * (sq(grad(0, 0)) + sq(grad(1, 1)) + sq(grad(2, 2)))
174  + grad(0, 1) * grad(1, 0) + grad(0, 2) * grad(2, 0)
175  + grad(1, 2) * grad(2, 1);
176 
177  vals(dof) = q_val;
178  }
179 
180  // Accumulate values in all dofs, count the zones.
181  for (int j = 0; j < dofs.Size(); j++)
182  {
183  int ldof = dofs[j];
184  q(ldof) += vals[j];
185  zones_per_vdof[ldof]++;
186  }
187  }
188 
189  // Communication
190 
191  // Count the zones globally.
192  GroupCommunicator &gcomm = q.ParFESpace()->GroupComm();
193  gcomm.Reduce<int>(zones_per_vdof, GroupCommunicator::Sum);
194  gcomm.Bcast(zones_per_vdof);
195 
196  // Accumulate for all vdofs.
197  gcomm.Reduce<double>(q.GetData(), GroupCommunicator::Sum);
198  gcomm.Bcast<double>(q.GetData());
199 
200  // Compute means
201  for (int i = 0; i < q.Size(); i++)
202  {
203  const int nz = zones_per_vdof[i];
204  if (nz)
205  {
206  q(i) /= nz;
207  }
208  }
209 }
210 
211 int main(int argc, char *argv[])
212 {
213  Mpi::Init(argc, argv);
214  Hypre::Init();
215 
216  OptionsParser args(argc, argv);
217  args.AddOption(&ctx.element_subdivisions,
218  "-es",
219  "--element-subdivisions",
220  "Number of 1d uniform subdivisions for each element.");
221  args.AddOption(&ctx.order,
222  "-o",
223  "--order",
224  "Order (degree) of the finite elements.");
225  args.AddOption(&ctx.dt, "-dt", "--time-step", "Time step.");
226  args.AddOption(&ctx.t_final, "-tf", "--final-time", "Final time.");
227  args.AddOption(&ctx.pa,
228  "-pa",
229  "--enable-pa",
230  "-no-pa",
231  "--disable-pa",
232  "Enable partial assembly.");
233  args.AddOption(&ctx.ni,
234  "-ni",
235  "--enable-ni",
236  "-no-ni",
237  "--disable-ni",
238  "Enable numerical integration rules.");
239  args.AddOption(&ctx.visualization,
240  "-vis",
241  "--visualization",
242  "-no-vis",
243  "--no-visualization",
244  "Enable or disable GLVis visualization.");
245  args.AddOption(
246  &ctx.checkres,
247  "-cr",
248  "--checkresult",
249  "-no-cr",
250  "--no-checkresult",
251  "Enable or disable checking of the result. Returns -1 on failure.");
252  args.Parse();
253  if (!args.Good())
254  {
255  if (Mpi::Root())
256  {
257  args.PrintUsage(mfem::out);
258  }
259  return 1;
260  }
261  if (Mpi::Root())
262  {
263  args.PrintOptions(mfem::out);
264  }
265 
266  Mesh orig_mesh("../../data/periodic-cube.mesh");
267  Mesh mesh = Mesh::MakeRefined(orig_mesh, ctx.element_subdivisions,
269  orig_mesh.Clear();
270 
271  mesh.EnsureNodes();
272  GridFunction *nodes = mesh.GetNodes();
273  *nodes *= M_PI;
274 
275  int nel = mesh.GetNE();
276  if (Mpi::Root())
277  {
278  mfem::out << "Number of elements: " << nel << std::endl;
279  }
280 
281  auto *pmesh = new ParMesh(MPI_COMM_WORLD, mesh);
282  mesh.Clear();
283 
284  // Create the flow solver.
285  NavierSolver flowsolver(pmesh, ctx.order, ctx.kinvis);
286  flowsolver.EnablePA(ctx.pa);
287  flowsolver.EnableNI(ctx.ni);
288 
289  // Set the initial condition.
290  ParGridFunction *u_ic = flowsolver.GetCurrentVelocity();
291  VectorFunctionCoefficient u_excoeff(pmesh->Dimension(), vel_tgv);
292  u_ic->ProjectCoefficient(u_excoeff);
293 
294  double t = 0.0;
295  double dt = ctx.dt;
296  double t_final = ctx.t_final;
297  bool last_step = false;
298 
299  flowsolver.Setup(dt);
300 
301  ParGridFunction *u_gf = flowsolver.GetCurrentVelocity();
302  ParGridFunction *p_gf = flowsolver.GetCurrentPressure();
303 
304  ParGridFunction w_gf(*u_gf);
305  ParGridFunction q_gf(*p_gf);
306  flowsolver.ComputeCurl3D(*u_gf, w_gf);
307  ComputeQCriterion(*u_gf, q_gf);
308 
309  QuantitiesOfInterest kin_energy(pmesh);
310 
311  ParaViewDataCollection pvdc("tgv_output", pmesh);
313  pvdc.SetHighOrderOutput(true);
314  pvdc.SetLevelsOfDetail(ctx.order);
315  pvdc.SetCycle(0);
316  pvdc.SetTime(t);
317  pvdc.RegisterField("velocity", u_gf);
318  pvdc.RegisterField("pressure", p_gf);
319  pvdc.RegisterField("vorticity", &w_gf);
320  pvdc.RegisterField("qcriterion", &q_gf);
321  pvdc.Save();
322 
323  double u_inf_loc = u_gf->Normlinf();
324  double p_inf_loc = p_gf->Normlinf();
325  double u_inf = GlobalLpNorm(infinity(), u_inf_loc, MPI_COMM_WORLD);
326  double p_inf = GlobalLpNorm(infinity(), p_inf_loc, MPI_COMM_WORLD);
327  double ke = kin_energy.ComputeKineticEnergy(*u_gf);
328 
329  std::string fname = "tgv_out_p_" + std::to_string(ctx.order) + ".txt";
330  FILE *f = NULL;
331 
332  if (Mpi::Root())
333  {
334  int nel1d = static_cast<int>(std::round(pow(nel, 1.0 / 3.0)));
335  int ngridpts = p_gf->ParFESpace()->GlobalVSize();
336  printf("%11s %11s %11s %11s %11s\n", "Time", "dt", "u_inf", "p_inf", "ke");
337  printf("%.5E %.5E %.5E %.5E %.5E\n", t, dt, u_inf, p_inf, ke);
338 
339  f = fopen(fname.c_str(), "w");
340  fprintf(f, "3D Taylor Green Vortex\n");
341  fprintf(f, "order = %d\n", ctx.order);
342  fprintf(f, "grid = %d x %d x %d\n", nel1d, nel1d, nel1d);
343  fprintf(f, "dofs per component = %d\n", ngridpts);
344  fprintf(f, "=================================================\n");
345  fprintf(f, " time kinetic energy\n");
346  fprintf(f, "%20.16e %20.16e\n", t, ke);
347  fflush(f);
348  fflush(stdout);
349  }
350 
351  for (int step = 0; !last_step; ++step)
352  {
353  if (t + dt >= t_final - dt / 2)
354  {
355  last_step = true;
356  }
357 
358  flowsolver.Step(t, dt, step);
359 
360  if ((step + 1) % 100 == 0 || last_step)
361  {
362  flowsolver.ComputeCurl3D(*u_gf, w_gf);
363  ComputeQCriterion(*u_gf, q_gf);
364  pvdc.SetCycle(step);
365  pvdc.SetTime(t);
366  pvdc.Save();
367  }
368 
369  u_inf_loc = u_gf->Normlinf();
370  p_inf_loc = p_gf->Normlinf();
371  u_inf = GlobalLpNorm(infinity(), u_inf_loc, MPI_COMM_WORLD);
372  p_inf = GlobalLpNorm(infinity(), p_inf_loc, MPI_COMM_WORLD);
373  ke = kin_energy.ComputeKineticEnergy(*u_gf);
374  if (Mpi::Root())
375  {
376  printf("%.5E %.5E %.5E %.5E %.5E\n", t, dt, u_inf, p_inf, ke);
377  fprintf(f, "%20.16e %20.16e\n", t, ke);
378  fflush(f);
379  fflush(stdout);
380  }
381  }
382 
383  flowsolver.PrintTimingData();
384 
385  // Test if the result for the test run is as expected.
386  if (ctx.checkres)
387  {
388  double tol = 2e-5;
389  double ke_expected = 1.25e-1;
390  if (fabs(ke - ke_expected) > tol)
391  {
392  if (Mpi::Root())
393  {
394  mfem::out << "Result has a larger error than expected."
395  << std::endl;
396  }
397  return -1;
398  }
399  }
400 
401  delete pmesh;
402 
403  return 0;
404 }
Abstract class for all finite elements.
Definition: fe_base.hpp:232
void PrintTimingData()
Print timing summary of the solving routine.
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
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
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
void SetDataFormat(VTKFormat fmt)
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
int Dimension() const
Definition: mesh.hpp:1047
void Bcast(T *ldata, int layout) const
Broadcast within each group where the master is the root.
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:512
Helper class for ParaView visualization data.
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
static Mesh MakeRefined(Mesh &orig_mesh, int ref_factor, int ref_type)
Create a refined (by any factor) version of orig_mesh.
Definition: mesh.cpp:3746
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
const DenseMatrix & InverseJacobian()
Return the inverse of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:145
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:525
const IntegrationPoint & GetIntPoint()
Get a const reference to the currently set integration point. This will return NULL if no integration...
Definition: eltrans.hpp:98
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2783
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
void Setup(double dt)
Initialize forms, solvers and preconditioners.
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:319
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
void ComputeCurl3D(ParGridFunction &u, ParGridFunction &cu)
Compute for .
Class for parallel linear form.
Definition: plinearform.hpp:26
Communicator performing operations within groups defined by a GroupTopology with arbitrary-size data ...
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void GetValues(int i, const IntegrationRule &ir, Vector &vals, int vdim=1) const
Definition: gridfunc.cpp:524
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
static void Sum(OpData< T >)
Reduce operation Sum, instantiated for int and double.
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
ParGridFunction * GetCurrentPressure()
Return a pointer to the current pressure ParGridFunction.
void EnablePA(bool pa)
Enable partial assembly for every operator.
void SetHighOrderOutput(bool high_order_output_)
static void Init()
Singleton creation with Mpi::Init();.
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:683
void SetTime(double t)
Set physical time (for time-dependent simulations)
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:614
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:208
static bool Root()
Return true if the rank in MPI_COMM_WORLD is zero.
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:581
A general vector function coefficient.
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:281
double GlobalLpNorm(const double p, double loc_norm, MPI_Comm comm)
Compute a global Lp norm from the local Lp norms computed by each processor.
Definition: pgridfunc.cpp:1065
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe_base.hpp:310
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
HYPRE_BigInt GlobalVSize() const
Definition: pfespace.hpp:283
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
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:113
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
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:684
GroupCommunicator & GroupComm()
Return a reference to the internal GroupCommunicator (on VDofs)
Definition: pfespace.hpp:338
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:322
ParGridFunction * GetCurrentVelocity()
Return a pointer to the current velocity ParGridFunction.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:936
Class for integration point with weight.
Definition: intrules.hpp:25
int dim
Definition: ex24.cpp:53
void Reduce(T *ldata, void(*Op)(OpData< T >)) const
Reduce within each group where the master is the root.
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:587
double infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition: vector.hpp:46
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
RefCoord t[3]
void SetLevelsOfDetail(int levels_of_detail_)
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:912
Nodes: x_i = i/(n-1), i=0,...,n-1.
Definition: fe_base.hpp:34
Vector data type.
Definition: vector.hpp:60
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:252
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7949
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
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:388
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Class for parallel meshes.
Definition: pmesh.hpp:32
Transient incompressible Navier Stokes solver in a split scheme formulation.
double Normlinf() const
Returns the l_infinity norm of the vector.
Definition: vector.cpp:822
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:326
void EnsureNodes()
Make sure that the mesh has valid nodes, i.e. its geometry is described by a vector finite element gr...
Definition: mesh.cpp:5304
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:647