Example 04: Lid-driven cavity flow
Description
This section will introduce the lid-driven cavity flow example, which is a real case in Computational Fluid Dynamics (CFD). The lid-driven cavity problem tries to solve the incompressible viscous flow in a unit square domain with no-slip walls, where only the upper lid has a constant velocity from left to right. The governing equation is the Navier-Stokes equation:
with the incompressibility condition:
We assume the velocity to be zero at \(t=0\), together with the lid’s boundary condition:
We adopt the 2nd order fully implicit projection method developed by Kim et al. (2002) . This method uses Crank-Nicolson scheme both for convection and viscosity terms, leading to a fully implicit method which can use a CFL number as large as 2.0. The discretized equations can be written as:
We won’t explain in detail how this method works. Please to the original paper for better knowledge. Writing your algorithm in the discretized form is always the first thing you need to do before using OpFlow.
Implementation
Now we show how the above equations are translated into OpFlow statements. To begin with, we need to discover how many types of operators we need. Since we have already solve the heat transfer equation and Poisson equation on Cartesian mesh, the viscosity term and Poisson equation is not longer an issue. The major problem is the convection term, who’s stencil pattern is like:
The figure above shows how the term \(\partial (uv)/\partial y\) is calculated. We need to interpolate the
\(u\) field to corner in \(y\)-direction and the \(v\) field to corner in \(x\)-direction. Then,
multiply the two interpolated field and calculate its 1st order derivative in \(y\)-direction. This can be done
with the d1IntpCenterToCorner operator:
auto conv_xy = [&](auto&& _1, auto&& _2) {
return dy<D1FirstOrderCenteredUpwind>(d1IntpCenterToCorner<1>(_1) * d1IntpCenterToCorner<0>(_2));
};
The above code defines a functor doing this calculation. Defining functors for sub-expressions are very useful for
composing complex expressions. Just as intermediate expressions, they add zero cost to the actual computation.
Following the same pattern, we can define other convection functors conv_xx, conv_yx and conv_yy.
Since we have all the operators we need, we can start to composing the main algorithm. We start with defining the three implicit solvers for the \(u\), \(v\) and \(p\) equations:
StructSolverParams<StructSolverType::GMRES> params; params.tol = 1e-10; params.maxIter = 100;
StructSolverParams<StructSolverType::GMRES> poisson_params = params;
StructSolverParams<StructSolverType::PFMG> p_params {.useZeroGuess = true, .relaxType = 1, .rapType = 0,
.numPreRelax = 1, .numPostRelax = 1, .skipRelax = 0}; p_params.tol = 1e-10;
auto solver = PrecondStructSolver<StructSolverType::GMRES, StructSolverType::PFMG>(params, p_params);
auto u_handler = makeEqnSolveHandler(
[&](auto&& e) {
return e / dt + conv_xx(u, e) + 0.5 * conv_xy(e, v)
== nu * (d2x<D2SecondOrderCentered>(u) + d2y<D2SecondOrderCentered>(u))
+ 0.5 * nu * (d2x<D2SecondOrderCentered>(e) + d2y<D2SecondOrderCentered>(e))
- (conv_xx(u, u) + conv_xy(u, v)) - dx<D1FirstOrderCenteredDownwind>(p);
},
du, solver);
auto v_handler = makeEqnSolveHandler(
[&](auto&& e) {
return e / dt + conv_yy(v, e) + conv_yy(v, v) + conv_yx(u, v) + 0.5 * conv_yx(u, e)
+ 0.5 * conv_yx(du, v)
== nu * (d2x<D2SecondOrderCentered>(v) + d2y<D2SecondOrderCentered>(v))
+ 0.5 * nu * (d2x<D2SecondOrderCentered>(e) + d2y<D2SecondOrderCentered>(e))
- dy<D1FirstOrderCenteredDownwind>(p);
},
dv, solver);
poisson_params.staticMat = true; poisson_params.pinValue = true;
auto p_solver = PrecondStructSolver<StructSolverType::GMRES, StructSolverType::PFMG>(poisson_params, p_params);
auto p_handler = makeEqnSolveHandler(
[&](auto&& e) {
return d2x<D2SecondOrderCentered>(e) + d2y<D2SecondOrderCentered>(e)
== (dx<D1FirstOrderCenteredUpwind>(du) + dy<D1FirstOrderCenteredUpwind>(dv)) / dt;
},
dp, p_solver);
Recall that in example Poisson equation we mentioned that by explicitly constructing
solver handlers, we can reuse the generated solver & matrix repeatedly. Also note that since the Poisson equation
is not well defined when all boundaries use zero Neumann conditions, we pin the left corner of pressure to be zero.
This is done by specifying poisson_params.pinValue=true for the Poisson solver.
We now define some field streams for data saving:
Utils::TecplotASCIIStream uf("u.tec"), vf("v.tec"), pf("p.tec");
Finally, we start the main loop:
for (auto i = 0; i < 1000; ++i) {
u_handler.solve();
v_handler.solve();
du = du - dt * conv_xy(u, dv);
u = u + du;
v = v + dv;
p_handler.solve();
u = u - dt * dx<D1FirstOrderCenteredDownwind>(dp);
v = v - dt * dy<D1FirstOrderCenteredDownwind>(dp);
p = p + dp;
uf << Utils::TimeStamp(i) << u; vf << Utils::TimeStamp(i) << v; pf << Utils::TimeStamp(i) << p;
OP_INFO("Current step {}", i);
}
The complete code for this example can be found in examples/LidDriven/LidDriven2D.cpp. Also checkout the MPI version
& 3D version in LidDriven2D-MPI.cpp and LidDriven3D.cpp. As you can see, composing a numerical solver in OpFlow
is concise and straightforward. You can easily write a implicit solver within 99 lines of code, with intuitive mathematical
semantics. Meanwhile, all the computational details can be explicit shown in the code, without writing tedious
differential & interpolation schemes over and over again. With OpFlow’s evaluation engine, the unified field operations
are translated into element-wise operations automatically. Together with inlining, auto vectorization and other powerful
optimizations provided by modern compilers, the final compiled program can be run on both multi-thread workstations
& distributed memory clusters, giving users decent performance comparable with hard coded solvers. This is the original
aspiration of OpFlow, and a major advantage over other similar languages & software.
Visualization
Run the above example and plot the result with any compatible tools, you can get some fields plots like:
