MFEM  v4.2.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
navier_3dfoc.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, 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_Session mpi(argc, argv);
52 
53  int serial_refinements = 0;
54 
55  Mesh *mesh = new Mesh("box-cylinder.mesh");
56 
57  for (int i = 0; i < serial_refinements; ++i)
58  {
59  mesh->UniformRefinement();
60  }
61 
62  if (mpi.Root())
63  {
64  std::cout << "Number of elements: " << mesh->GetNE() << std::endl;
65  }
66 
67  auto *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
68  delete mesh;
69 
70  // Create the flow solver.
71  NavierSolver flowsolver(pmesh, ctx.order, ctx.kin_vis);
72  flowsolver.EnablePA(true);
73 
74  // Set the initial condition.
75  ParGridFunction *u_ic = flowsolver.GetCurrentVelocity();
76  VectorFunctionCoefficient u_excoeff(pmesh->Dimension(), vel);
77  u_ic->ProjectCoefficient(u_excoeff);
78 
79  // Add Dirichlet boundary conditions to velocity space restricted to
80  // selected attributes on the mesh.
81  Array<int> attr(pmesh->bdr_attributes.Max());
82  // Inlet is attribute 1.
83  attr[0] = 1;
84  // Walls is attribute 3.
85  attr[2] = 1;
86  flowsolver.AddVelDirichletBC(vel, attr);
87 
88  double t = 0.0;
89  double dt = ctx.dt;
90  double t_final = ctx.t_final;
91  bool last_step = false;
92 
93  flowsolver.Setup(dt);
94 
95  ParGridFunction *u_gf = flowsolver.GetCurrentVelocity();
96  ParGridFunction *p_gf = flowsolver.GetCurrentPressure();
97 
98  ParaViewDataCollection pvdc("3dfoc", pmesh);
100  pvdc.SetHighOrderOutput(true);
101  pvdc.SetLevelsOfDetail(ctx.order);
102  pvdc.SetCycle(0);
103  pvdc.SetTime(t);
104  pvdc.RegisterField("velocity", u_gf);
105  pvdc.RegisterField("pressure", p_gf);
106  pvdc.Save();
107 
108  for (int step = 0; !last_step; ++step)
109  {
110  if (t + dt >= t_final - dt / 2)
111  {
112  last_step = true;
113  }
114 
115  flowsolver.Step(t, dt, step);
116 
117  if (step % 10 == 0)
118  {
119  pvdc.SetCycle(step);
120  pvdc.SetTime(t);
121  pvdc.Save();
122  }
123 
124  if (mpi.Root())
125  {
126  printf("%11s %11s\n", "Time", "dt");
127  printf("%.5E %.5E\n", t, dt);
128  fflush(stdout);
129  }
130  }
131 
132  flowsolver.PrintTimingData();
133 
134  delete pmesh;
135 
136  return 0;
137 }
void PrintTimingData()
Print timing summary of the solving routine.
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void SetDataFormat(VTKFormat fmt)
Helper class for ParaView visualization data.
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:737
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: pgridfunc.cpp:474
int main(int argc, char *argv[])
Definition: ex1.cpp:66
void Setup(double dt)
Initialize forms, solvers and preconditioners.
void Step(double &time, double dt, int cur_step)
Compute solution at the next time step t+dt.
A simple convenience class that calls MPI_Init() at construction and MPI_Finalize() at destruction...
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:8382
ParGridFunction * GetCurrentPressure()
Return a pointer to the current pressure ParGridFunction.
bool Root() const
Return true if WorldRank() == 0.
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_)
void SetTime(double t)
Set physical time (for time-dependent simulations)
A general vector function coefficient.
ParGridFunction * GetCurrentVelocity()
Return a pointer to the current velocity ParGridFunction.
void SetLevelsOfDetail(int levels_of_detail_)
Vector data type.
Definition: vector.hpp:51
virtual void Save() override
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.