MFEM  v4.2.0 Finite element discretization library
navier_mms.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, 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_Session mpi(argc, argv);
89
90  OptionsParser args(argc, argv);
92  "-rs",
93  "--refine-serial",
94  "Number of times to refine the mesh uniformly in serial.");
96  "-o",
97  "--order",
98  "Order (degree) of the finite elements.");
99  args.AddOption(&ctx.dt, "-dt", "--time-step", "Time step.");
100  args.AddOption(&ctx.t_final, "-tf", "--final-time", "Final time.");
102  "-pa",
103  "--enable-pa",
104  "-no-pa",
105  "--disable-pa",
106  "Enable partial assembly.");
108  "-ni",
109  "--enable-ni",
110  "-no-ni",
111  "--disable-ni",
112  "Enable numerical integration rules.");
114  "-vis",
115  "--visualization",
116  "-no-vis",
117  "--no-visualization",
118  "Enable or disable GLVis visualization.");
120  &ctx.checkres,
121  "-cr",
122  "--checkresult",
123  "-no-cr",
124  "--no-checkresult",
125  "Enable or disable checking of the result. Returns -1 on failure.");
126  args.Parse();
127  if (!args.Good())
128  {
129  if (mpi.Root())
130  {
131  args.PrintUsage(mfem::out);
132  }
133  return 1;
134  }
135  if (mpi.Root())
136  {
137  args.PrintOptions(mfem::out);
138  }
139
140  Mesh *mesh = new Mesh("../../data/inline-quad.mesh");
141  mesh->EnsureNodes();
142  GridFunction *nodes = mesh->GetNodes();
143  *nodes *= 2.0;
144  *nodes -= 1.0;
145
146  for (int i = 0; i < ctx.ser_ref_levels; ++i)
147  {
148  mesh->UniformRefinement();
149  }
150
151  if (mpi.Root())
152  {
153  std::cout << "Number of elements: " << mesh->GetNE() << std::endl;
154  }
155
156  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
157  delete mesh;
158
159  // Create the flow solver.
160  NavierSolver naviersolver(pmesh, ctx.order, ctx.kinvis);
161  naviersolver.EnablePA(ctx.pa);
162  naviersolver.EnableNI(ctx.ni);
163
164  // Set the initial condition.
165  ParGridFunction *u_ic = naviersolver.GetCurrentVelocity();
166  VectorFunctionCoefficient u_excoeff(pmesh->Dimension(), vel);
167  u_ic->ProjectCoefficient(u_excoeff);
168
169  FunctionCoefficient p_excoeff(p);
170
171  // Add Dirichlet boundary conditions to velocity space restricted to
172  // selected attributes on the mesh.
173  Array<int> attr(pmesh->bdr_attributes.Max());
174  attr = 1;
176
177  Array<int> domain_attr(pmesh->attributes.Max());
178  domain_attr = 1;
180
181  double t = 0.0;
182  double dt = ctx.dt;
183  double t_final = ctx.t_final;
184  bool last_step = false;
185
186  naviersolver.Setup(dt);
187
188  double err_u = 0.0;
189  double err_p = 0.0;
190  ParGridFunction *u_gf = nullptr;
191  ParGridFunction *p_gf = nullptr;
192  u_gf = naviersolver.GetCurrentVelocity();
193  p_gf = naviersolver.GetCurrentPressure();
194
195  for (int step = 0; !last_step; ++step)
196  {
197  if (t + dt >= t_final - dt / 2)
198  {
199  last_step = true;
200  }
201
202  naviersolver.Step(t, dt, step);
203
204  // Compare against exact solution of velocity and pressure.
205  u_excoeff.SetTime(t);
206  p_excoeff.SetTime(t);
207  err_u = u_gf->ComputeL2Error(u_excoeff);
208  err_p = p_gf->ComputeL2Error(p_excoeff);
209
210  if (mpi.Root())
211  {
212  printf("%11s %11s %11s %11s\n", "Time", "dt", "err_u", "err_p");
213  printf("%.5E %.5E %.5E %.5E err\n", t, dt, err_u, err_p);
214  fflush(stdout);
215  }
216  }
217
218  if (ctx.visualization)
219  {
220  char vishost[] = "localhost";
221  int visport = 19916;
222  socketstream sol_sock(vishost, visport);
223  sol_sock << "parallel " << mpi.WorldSize() << " " << mpi.WorldRank()
224  << "\n";
225  sol_sock << "solution\n" << *pmesh << *u_ic << std::flush;
226  }
227
228  naviersolver.PrintTimingData();
229
230  // Test if the result for the test run is as expected.
231  if (ctx.checkres)
232  {
233  double tol = 1e-3;
234  if (err_u > tol || err_p > tol)
235  {
236  if (mpi.Root())
237  {
238  mfem::out << "Result has a larger error than expected."
239  << std::endl;
240  }
241  return -1;
242  }
243  }
244
245  delete pmesh;
246
247  return 0;
248 }
int WorldSize() const
Return MPI_COMM_WORLD&#39;s size.
void PrintTimingData()
Print timing summary of the solving routine.
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
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...
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:150
constexpr char vishost[]
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:8382
constexpr int visport
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
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.
int Dimension() const
Definition: mesh.hpp:788
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:434
void SetTime(double t)
Set the time for time dependent coefficients.
Definition: coefficient.hpp:48
A general vector function coefficient.
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:201
int WorldRank() const
Return MPI_COMM_WORLD&#39;s rank.
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:267
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:304
A general function coefficient.
Vector data type.
Definition: vector.hpp:51
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:6603
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.
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition: mesh.hpp:199
void EnsureNodes()
Definition: mesh.cpp:4132
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:145