MFEM  v4.6.0
Finite element discretization library
navier_3dfoc.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 flow over a cylinder benchmark example
13 
14 #include "navier_solver.hpp"
15 #include <fstream>
16 
17 using namespace mfem;
18 using namespace navier;
19 
20 struct s_NavierContext
21 {
22  int order = 4;
23  double kin_vis = 0.001;
24  double t_final = 8.0;
25  double dt = 1e-3;
26 } ctx;
27 
28 void vel(const Vector &x, double t, Vector &u)
29 {
30  double xi = x(0);
31  double yi = x(1);
32  double zi = x(2);
33 
34  double U = 2.25;
35 
36  if (xi <= 1e-8)
37  {
38  u(0) = 16.0 * U * yi * zi * sin(M_PI * t / 8.0) * (0.41 - yi)
39  * (0.41 - zi) / pow(0.41, 4.0);
40  }
41  else
42  {
43  u(0) = 0.0;
44  }
45  u(1) = 0.0;
46  u(2) = 0.0;
47 }
48 
49 int main(int argc, char *argv[])
50 {
51  Mpi::Init(argc, argv);
52  Hypre::Init();
53 
54  int serial_refinements = 0;
55 
56  Mesh *mesh = new Mesh("box-cylinder.mesh");
57 
58  for (int i = 0; i < serial_refinements; ++i)
59  {
60  mesh->UniformRefinement();
61  }
62 
63  if (Mpi::Root())
64  {
65  std::cout << "Number of elements: " << mesh->GetNE() << std::endl;
66  }
67 
68  auto *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
69  delete mesh;
70 
71  // Create the flow solver.
72  NavierSolver flowsolver(pmesh, ctx.order, ctx.kin_vis);
73  flowsolver.EnablePA(true);
74 
75  // Set the initial condition.
76  ParGridFunction *u_ic = flowsolver.GetCurrentVelocity();
77  VectorFunctionCoefficient u_excoeff(pmesh->Dimension(), vel);
78  u_ic->ProjectCoefficient(u_excoeff);
79 
80  // Add Dirichlet boundary conditions to velocity space restricted to
81  // selected attributes on the mesh.
82  Array<int> attr(pmesh->bdr_attributes.Max());
83  // Inlet is attribute 1.
84  attr[0] = 1;
85  // Walls is attribute 3.
86  attr[2] = 1;
87  flowsolver.AddVelDirichletBC(vel, attr);
88 
89  double t = 0.0;
90  double dt = ctx.dt;
91  double t_final = ctx.t_final;
92  bool last_step = false;
93 
94  flowsolver.Setup(dt);
95 
96  ParGridFunction *u_gf = flowsolver.GetCurrentVelocity();
97  ParGridFunction *p_gf = flowsolver.GetCurrentPressure();
98 
99  ParaViewDataCollection pvdc("3dfoc", pmesh);
101  pvdc.SetHighOrderOutput(true);
102  pvdc.SetLevelsOfDetail(ctx.order);
103  pvdc.SetCycle(0);
104  pvdc.SetTime(t);
105  pvdc.RegisterField("velocity", u_gf);
106  pvdc.RegisterField("pressure", p_gf);
107  pvdc.Save();
108 
109  for (int step = 0; !last_step; ++step)
110  {
111  if (t + dt >= t_final - dt / 2)
112  {
113  last_step = true;
114  }
115 
116  flowsolver.Step(t, dt, step);
117 
118  if (step % 10 == 0)
119  {
120  pvdc.SetCycle(step);
121  pvdc.SetTime(t);
122  pvdc.Save();
123  }
124 
125  if (Mpi::Root())
126  {
127  printf("%11s %11s\n", "Time", "dt");
128  printf("%.5E %.5E\n", t, dt);
129  fflush(stdout);
130  }
131  }
132 
133  flowsolver.PrintTimingData();
134 
135  delete pmesh;
136 
137  return 0;
138 }
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
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.
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 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 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.