29 using namespace navier;
31 struct s_NavierContext
33 int ser_ref_levels = 1;
36 double t_final = 10 * 0.25e-4;
40 bool visualization =
false;
41 bool checkres =
false;
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));
58 return cos(M_PI * xi) * sin(t) * sin(M_PI * yi);
66 u(0) = M_PI * sin(t) * sin(M_PI * xi) * sin(M_PI * yi)
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))
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);
76 u(1) = M_PI * cos(M_PI * yi) * sin(t)
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);
86 int main(
int argc,
char *argv[])
94 "Number of times to refine the mesh uniformly in serial.");
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.");
106 "Enable partial assembly.");
112 "Enable numerical integration rules.");
117 "--no-visualization",
118 "Enable or disable GLVis visualization.");
125 "Enable or disable checking of the result. Returns -1 on failure.");
140 Mesh *mesh =
new Mesh(
"../../data/inline-quad.mesh");
146 for (
int i = 0; i <
ctx.ser_ref_levels; ++i)
153 std::cout <<
"Number of elements: " << mesh->
GetNE() << std::endl;
183 double t_final =
ctx.t_final;
184 bool last_step =
false;
186 naviersolver.
Setup(dt);
195 for (
int step = 0; !last_step; ++step)
197 if (t + dt >= t_final - dt / 2)
202 naviersolver.
Step(t, dt, step);
205 u_excoeff.SetTime(t);
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);
218 if (
ctx.visualization)
225 sol_sock <<
"solution\n" << *pmesh << *u_ic << std::flush;
234 if (err_u > tol || err_p > tol)
238 mfem::out <<
"Result has a larger error than expected."
int WorldSize() const
Return MPI_COMM_WORLD's size.
void PrintTimingData()
Print timing summary of the solving routine.
Class for grid function - Vector with associated FE space.
int GetNE() const
Returns number of elements.
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
void Setup(double dt)
Initialize forms, solvers and preconditioners.
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 ...
struct s_NavierContext ctx
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
ParGridFunction * GetCurrentPressure()
Return a pointer to the current pressure ParGridFunction.
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
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 PrintUsage(std::ostream &out) const
Print the usage message.
void SetTime(double t)
Set the time for time dependent coefficients.
A general vector function coefficient.
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
double p(const Vector &x, double t)
int WorldRank() const
Return MPI_COMM_WORLD'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 'var' to receive the value. Enable/disable tags are used to set the bool...
ParGridFunction * GetCurrentVelocity()
Return a pointer to the current velocity ParGridFunction.
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const
void vel(const Vector &x, double t, Vector &u)
void PrintOptions(std::ostream &out) const
Print the options.
void accel(const Vector &x, double t, Vector &u)
A general function coefficient.
void GetNodes(Vector &node_coord) const
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)
Class for parallel grid function.
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Class for parallel meshes.
Transient incompressible Navier Stokes solver in a split scheme formulation.
Array< int > attributes
A list of all unique element attributes used by the Mesh.
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.