MFEM  v4.4.0 Finite element discretization library
navier_mms.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
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 MMS example
13 //
14 // A manufactured solution is defined as
15 //
16 // u = [pi * sin(t) * sin(pi * x)^2 * sin(2 * pi * y),
17 // -(pi * sin(t) * sin(2 * pi * x)) * sin(pi * y)^2].
18 //
19 // p = cos(pi * x) * sin(t) * sin(pi * y)
20 //
21 // The solution is used to compute the symbolic forcing term (right hand side),
22 // of the equation. Then the numerical solution is computed and compared to the
23 // exact manufactured solution to determine the error.
24
25 #include "navier_solver.hpp"
26 #include <fstream>
27
28 using namespace mfem;
29 using namespace navier;
30
31 struct s_NavierContext
32 {
33  int ser_ref_levels = 1;
34  int order = 5;
35  double kinvis = 1.0;
36  double t_final = 10 * 0.25e-4;
37  double dt = 0.25e-4;
38  bool pa = true;
39  bool ni = false;
40  bool visualization = false;
41  bool checkres = false;
42 } ctx;
43
44 void vel(const Vector &x, double t, Vector &u)
45 {
46  double xi = x(0);
47  double yi = x(1);
48
49  u(0) = M_PI * sin(t) * pow(sin(M_PI * xi), 2.0) * sin(2.0 * M_PI * yi);
50  u(1) = -(M_PI * sin(t) * sin(2.0 * M_PI * xi) * pow(sin(M_PI * yi), 2.0));
51 }
52
53 double p(const Vector &x, double t)
54 {
55  double xi = x(0);
56  double yi = x(1);
57
58  return cos(M_PI * xi) * sin(t) * sin(M_PI * yi);
59 }
60
61 void accel(const Vector &x, double t, Vector &u)
62 {
63  double xi = x(0);
64  double yi = x(1);
65
66  u(0) = M_PI * sin(t) * sin(M_PI * xi) * sin(M_PI * yi)
67  * (-1.0
68  + 2.0 * pow(M_PI, 2.0) * sin(t) * sin(M_PI * xi)
69  * sin(2.0 * M_PI * xi) * sin(M_PI * yi))
70  + M_PI
71  * (2.0 * ctx.kinvis * pow(M_PI, 2.0)
72  * (1.0 - 2.0 * cos(2.0 * M_PI * xi)) * sin(t)
73  + cos(t) * pow(sin(M_PI * xi), 2.0))
74  * sin(2.0 * M_PI * yi);
75
76  u(1) = M_PI * cos(M_PI * yi) * sin(t)
77  * (cos(M_PI * xi)
78  + 2.0 * ctx.kinvis * pow(M_PI, 2.0) * cos(M_PI * yi)
79  * sin(2.0 * M_PI * xi))
80  - M_PI * (cos(t) + 6.0 * ctx.kinvis * pow(M_PI, 2.0) * sin(t))
81  * sin(2.0 * M_PI * xi) * pow(sin(M_PI * yi), 2.0)
82  + 4.0 * pow(M_PI, 3.0) * cos(M_PI * yi) * pow(sin(t), 2.0)
83  * pow(sin(M_PI * xi), 2.0) * pow(sin(M_PI * yi), 3.0);
84 }
85
86 int main(int argc, char *argv[])
87 {
88  Mpi::Init(argc, argv);
89  Hypre::Init();
90
91  OptionsParser args(argc, argv);
93  "-rs",
94  "--refine-serial",
95  "Number of times to refine the mesh uniformly in serial.");
97  "-o",
98  "--order",
99  "Order (degree) of the finite elements.");
100  args.AddOption(&ctx.dt, "-dt", "--time-step", "Time step.");
101  args.AddOption(&ctx.t_final, "-tf", "--final-time", "Final time.");
103  "-pa",
104  "--enable-pa",
105  "-no-pa",
106  "--disable-pa",
107  "Enable partial assembly.");
109  "-ni",
110  "--enable-ni",
111  "-no-ni",
112  "--disable-ni",
113  "Enable numerical integration rules.");
115  "-vis",
116  "--visualization",
117  "-no-vis",
118  "--no-visualization",
119  "Enable or disable GLVis visualization.");
121  &ctx.checkres,
122  "-cr",
123  "--checkresult",
124  "-no-cr",
125  "--no-checkresult",
126  "Enable or disable checking of the result. Returns -1 on failure.");
127  args.Parse();
128  if (!args.Good())
129  {
130  if (Mpi::Root())
131  {
132  args.PrintUsage(mfem::out);
133  }
134  return 1;
135  }
136  if (Mpi::Root())
137  {
138  args.PrintOptions(mfem::out);
139  }
140
141  Mesh *mesh = new Mesh("../../data/inline-quad.mesh");
142  mesh->EnsureNodes();
143  GridFunction *nodes = mesh->GetNodes();
144  *nodes *= 2.0;
145  *nodes -= 1.0;
146
147  for (int i = 0; i < ctx.ser_ref_levels; ++i)
148  {
149  mesh->UniformRefinement();
150  }
151
152  if (Mpi::Root())
153  {
154  std::cout << "Number of elements: " << mesh->GetNE() << std::endl;
155  }
156
157  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
158  delete mesh;
159
160  // Create the flow solver.
161  NavierSolver naviersolver(pmesh, ctx.order, ctx.kinvis);
162  naviersolver.EnablePA(ctx.pa);
163  naviersolver.EnableNI(ctx.ni);
164
165  // Set the initial condition.
166  ParGridFunction *u_ic = naviersolver.GetCurrentVelocity();
167  VectorFunctionCoefficient u_excoeff(pmesh->Dimension(), vel);
168  u_ic->ProjectCoefficient(u_excoeff);
169
170  FunctionCoefficient p_excoeff(p);
171
172  // Add Dirichlet boundary conditions to velocity space restricted to
173  // selected attributes on the mesh.
174  Array<int> attr(pmesh->bdr_attributes.Max());
175  attr = 1;
177
178  Array<int> domain_attr(pmesh->attributes.Max());
179  domain_attr = 1;
181
182  double t = 0.0;
183  double dt = ctx.dt;
184  double t_final = ctx.t_final;
185  bool last_step = false;
186
187  naviersolver.Setup(dt);
188
189  double err_u = 0.0;
190  double err_p = 0.0;
191  ParGridFunction *u_gf = nullptr;
192  ParGridFunction *p_gf = nullptr;
193  u_gf = naviersolver.GetCurrentVelocity();
194  p_gf = naviersolver.GetCurrentPressure();
195
196  for (int step = 0; !last_step; ++step)
197  {
198  if (t + dt >= t_final - dt / 2)
199  {
200  last_step = true;
201  }
202
203  naviersolver.Step(t, dt, step);
204
205  // Compare against exact solution of velocity and pressure.
206  u_excoeff.SetTime(t);
207  p_excoeff.SetTime(t);
208  err_u = u_gf->ComputeL2Error(u_excoeff);
209  err_p = p_gf->ComputeL2Error(p_excoeff);
210
211  if (Mpi::Root())
212  {
213  printf("%11s %11s %11s %11s\n", "Time", "dt", "err_u", "err_p");
214  printf("%.5E %.5E %.5E %.5E err\n", t, dt, err_u, err_p);
215  fflush(stdout);
216  }
217  }
218
219  if (ctx.visualization)
220  {
221  char vishost[] = "localhost";
222  int visport = 19916;
223  socketstream sol_sock(vishost, visport);
224  sol_sock << "parallel " << Mpi::WorldSize() << " "
225  << Mpi::WorldRank() << "\n";
226  sol_sock << "solution\n" << *pmesh << *u_ic << std::flush;
227  }
228
229  naviersolver.PrintTimingData();
230
231  // Test if the result for the test run is as expected.
232  if (ctx.checkres)
233  {
234  double tol = 1e-3;
235  if (err_u > tol || err_p > tol)
236  {
237  if (Mpi::Root())
238  {
239  mfem::out << "Result has a larger error than expected."
240  << std::endl;
241  }
242  return -1;
243  }
244  }
245
246  delete pmesh;
247
248  return 0;
249 }
const char vishost[]
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:67
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:928
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:525
static int WorldSize()
Return the size of MPI_COMM_WORLD.
virtual void SetTime(double t)
Set the time for time dependent coefficients.
Definition: coefficient.hpp:48
void Setup(double dt)
Initialize forms, solvers and preconditioners.
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9737
ParGridFunction * GetCurrentPressure()
Return a pointer to the current pressure ParGridFunction.
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
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.
static void Init()
Singleton creation with Mpi::Init();.
const int visport
int Dimension() const
Definition: mesh.hpp:999
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
FDualNumber< tbase > cos(const FDualNumber< tbase > &f)
cos([dual number])
Definition: fdual.hpp:471
static bool Root()
Return true if the rank in MPI_COMM_WORLD is zero.
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
Definition: fdual.hpp:543
A general vector function coefficient.
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
FDualNumber< tbase > sin(const FDualNumber< tbase > &f)
sin([dual number])
Definition: fdual.hpp:572
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
ParGridFunction * GetCurrentVelocity()
Return a pointer to the current velocity ParGridFunction.
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const
Definition: pgridfunc.hpp:281
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
RefCoord t[3]
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7841
void Step(double &time, double dt, int cur_step, bool provisional=false)
Compute solution at the next time step t+dt.
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
Class for parallel grid function.
Definition: pgridfunc.hpp:32
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
Class for parallel meshes.
Definition: pmesh.hpp:32
Transient incompressible Navier Stokes solver in a split scheme formulation.
int main()
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition: mesh.hpp:268
void EnsureNodes()
Definition: mesh.cpp:5279
void AddAccelTerm(VectorCoefficient *coeff, Array< int > &attr)
Add an acceleration term to the RHS of the equation.
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150