MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
navier_mms.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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
28using namespace mfem;
29using namespace navier;
30
31struct s_NavierContext
32{
33 int ser_ref_levels = 1;
34 int order = 5;
35 real_t kinvis = 1.0;
36 real_t t_final = 10 * 0.25e-4;
37 real_t dt = 0.25e-4;
38 bool pa = true;
39 bool ni = false;
40 bool visualization = false;
41 bool checkres = false;
43
44void vel(const Vector &x, real_t t, Vector &u)
45{
46 real_t xi = x(0);
47 real_t 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
54{
55 real_t xi = x(0);
56 real_t yi = x(1);
57
58 return cos(M_PI * xi) * sin(t) * sin(M_PI * yi);
59}
60
61void accel(const Vector &x, real_t t, Vector &u)
62{
63 real_t xi = x(0);
64 real_t 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
86int main(int argc, char *argv[])
87{
88 Mpi::Init(argc, argv);
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 {
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 real_t t = 0.0;
183 real_t dt = ctx.dt;
184 real_t t_final = ctx.t_final;
185 bool last_step = false;
186
187 naviersolver.Setup(dt);
188
189 real_t err_u = 0.0;
190 real_t 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 real_t 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}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
virtual void SetTime(real_t t)
Set the time for time dependent coefficients.
A general function coefficient.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
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:6159
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void GetNodes(Vector &node_coord) const
Definition mesh.cpp:8973
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition mesh.hpp:280
static bool Root()
Return true if the rank in MPI_COMM_WORLD is zero.
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static int WorldSize()
Return the size of MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Class for parallel grid function.
Definition pgridfunc.hpp:33
real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const override
void ProjectCoefficient(Coefficient &coeff) override
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Class for parallel meshes.
Definition pmesh.hpp:34
virtual void SetTime(real_t t)
Set the time for time dependent coefficients.
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
Transient incompressible Navier Stokes solver in a split scheme formulation.
void Setup(real_t dt)
Initialize forms, solvers and preconditioners.
void AddVelDirichletBC(VectorCoefficient *coeff, Array< int > &attr)
Add a Dirichlet boundary condition to the velocity field.
ParGridFunction * GetCurrentPressure()
Return a pointer to the current pressure ParGridFunction.
void PrintTimingData()
Print timing summary of the solving routine.
void EnablePA(bool pa)
Enable partial assembly for every operator.
ParGridFunction * GetCurrentVelocity()
Return a pointer to the current velocity ParGridFunction.
void AddAccelTerm(VectorCoefficient *coeff, Array< int > &attr)
Add an acceleration term to the RHS of the equation.
void Step(real_t &time, real_t dt, int cur_step, bool provisional=false)
Compute solution at the next time step t+dt.
int main()
const int visport
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
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
float real_t
Definition config.hpp:43
const char vishost[]
struct s_NavierContext ctx
void vel(const Vector &x, real_t t, Vector &u)
real_t p(const Vector &x, real_t t)
void accel(const Vector &x, real_t t, Vector &u)
RefCoord t[3]