MFEM  v4.6.0
Finite element discretization library
navier_shear.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 // Navier double shear layer example
13 //
14 // Solve the double shear problem in the following configuration.
15 //
16 // +-------------------+
17 // | |
18 // | u0 = ua |
19 // | |
20 // -------------------------------- y = 0.5
21 // | |
22 // | u0 = ub |
23 // | |
24 // +-------------------+
25 //
26 // The initial condition u0 is chosen to be a varying velocity in the y
27 // direction. It includes a perturbation at x = 0.5 which leads to an
28 // instability and the dynamics of the flow. The boundary conditions are fully
29 // periodic.
30 
31 #include "navier_solver.hpp"
32 #include <fstream>
33 
34 using namespace mfem;
35 using namespace navier;
36 
37 struct s_NavierContext
38 {
39  int order = 6;
40  double kinvis = 1.0 / 100000.0;
41  double t_final = 10 * 1e-3;
42  double dt = 1e-3;
43 } ctx;
44 
45 void vel_shear_ic(const Vector &x, double t, Vector &u)
46 {
47  double xi = x(0);
48  double yi = x(1);
49 
50  double rho = 30.0;
51  double delta = 0.05;
52 
53  if (yi <= 0.5)
54  {
55  u(0) = tanh(rho * (yi - 0.25));
56  }
57  else
58  {
59  u(0) = tanh(rho * (0.75 - yi));
60  }
61 
62  u(1) = delta * sin(2.0 * M_PI * xi);
63 }
64 
65 int main(int argc, char *argv[])
66 {
67  Mpi::Init(argc, argv);
68  Hypre::Init();
69 
70  int serial_refinements = 2;
71 
72  Mesh *mesh = new Mesh("../../data/periodic-square.mesh");
73  mesh->EnsureNodes();
74  GridFunction *nodes = mesh->GetNodes();
75  *nodes -= -1.0;
76  *nodes /= 2.0;
77 
78  for (int i = 0; i < serial_refinements; ++i)
79  {
80  mesh->UniformRefinement();
81  }
82 
83  if (Mpi::Root())
84  {
85  std::cout << "Number of elements: " << mesh->GetNE() << std::endl;
86  }
87 
88  auto *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
89  delete mesh;
90 
91  // Create the flow solver.
92  NavierSolver flowsolver(pmesh, ctx.order, ctx.kinvis);
93  flowsolver.EnablePA(true);
94 
95  // Set the initial condition.
96  ParGridFunction *u_ic = flowsolver.GetCurrentVelocity();
97  VectorFunctionCoefficient u_excoeff(pmesh->Dimension(), vel_shear_ic);
98  u_ic->ProjectCoefficient(u_excoeff);
99 
100  double t = 0.0;
101  double dt = ctx.dt;
102  double t_final = ctx.t_final;
103  bool last_step = false;
104 
105  flowsolver.Setup(dt);
106 
107  ParGridFunction *u_gf = flowsolver.GetCurrentVelocity();
108  ParGridFunction *p_gf = flowsolver.GetCurrentPressure();
109 
110  ParGridFunction w_gf(*u_gf);
111  flowsolver.ComputeCurl2D(*u_gf, w_gf);
112 
113  ParaViewDataCollection pvdc("shear_output", pmesh);
115  pvdc.SetHighOrderOutput(true);
116  pvdc.SetLevelsOfDetail(ctx.order);
117  pvdc.SetCycle(0);
118  pvdc.SetTime(t);
119  pvdc.RegisterField("velocity", u_gf);
120  pvdc.RegisterField("pressure", p_gf);
121  pvdc.RegisterField("vorticity", &w_gf);
122  pvdc.Save();
123 
124  for (int step = 0; !last_step; ++step)
125  {
126  if (t + dt >= t_final - dt / 2)
127  {
128  last_step = true;
129  }
130 
131  flowsolver.Step(t, dt, step);
132 
133  if (step % 10 == 0)
134  {
135  flowsolver.ComputeCurl2D(*u_gf, w_gf);
136  pvdc.SetCycle(step);
137  pvdc.SetTime(t);
138  pvdc.Save();
139  }
140 
141  if (Mpi::Root())
142  {
143  printf("%11s %11s\n", "Time", "dt");
144  printf("%.5E %.5E\n", t, dt);
145  fflush(stdout);
146  }
147  }
148 
149  flowsolver.PrintTimingData();
150 
151  delete pmesh;
152 
153  return 0;
154 }
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
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void SetDataFormat(VTKFormat fmt)
Helper class for ParaView visualization data.
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:561
void Setup(double dt)
Initialize forms, solvers and preconditioners.
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
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 SetHighOrderOutput(bool high_order_output_)
static void Init()
Singleton creation with Mpi::Init();.
void ComputeCurl2D(ParGridFunction &u, ParGridFunction &cu, bool assume_scalar=false)
Compute for .
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 GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
RefCoord t[3]
void SetLevelsOfDetail(int levels_of_detail_)
Vector data type.
Definition: vector.hpp:58
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:8302
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:22
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.
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:5583