(sec:tut_cavity)= # Lid-driven Cavity This tutorial describes how to set up and run the first non-trivial flow problem. The lid-driven cavity flow is a standard test case for numerical schemes, and a number of results have been published in literature, see, e.g., {cite}`ghia1982high`, {cite}`gao2016efficient`. This tutorial assumes that you have completed the previous tutorial, know how to edit files and post-process the solution with your favorite visualization tool, e.g., **ParaView**. Also, the later parts of the tutorial assume that you have access to a computer with an MPI-based parallelization with at least $4$ computing cores - otherwise, it will just take a lot longer :). This tutorial is divided into two sections. The [Basic](#sec:tut_cavity_basic) section introduces the setup process and guides you through running simulations, providing a solid foundation for using the code. The [Advanced](#sec:tut_cavity_advanced) section builds on this, offering insights into code modifications that enable more complex simulations and the addition of custom features. If you’re mainly interested in running the code as provided, feel free to skip the Advanced section or only explore the parts that interest you. ## Flow Description The flow under consideration is essentially incompressible and two-dimensional, but we will use the three-dimensional code for the compressible Navier-Stokes equations to solve it here. This is not the most efficient way to compute this flow, but it works well as an example how to set up and run a simulation in **FLEXI**. The computation is conducted in a three-dimensional, square domain with periodic boundary conditions in the "third" direction. The walls of the cavity are modeled as isothermal walls, and a fixed flow is prescribed at the upper boundary, i.e., the lid of the domain. For the Reynolds numbers investigated here, this generates a steady, vortical flow field in the cavity. {numref}`fig:cavity_re400_velmag` shows the resulting velocity field and streamlines for $Re = 400$. ```{figure} ./figures/cavity_Re400_streamlines_cropped.jpg :name: fig:cavity_re400_velmag :align: center :width: 45% :alt: Velocity magnitude contours for $Re=400$ Contours of velocity magnitude for the Re = 400 lid-driven cavity case. ``` ### Compiler Options **FLEXI** should be compiled with the `cavity` preset using the following commands. ```bash cmake -B build --preset cavity cmake --build build ``` (sec:tut_cavity_basic)= ## Basic Tutorial | Flow at Re=100 The basic tutorial is contained in the `Basic_Re100` subfolder of the `tutorials\cavity` directory. ### Mesh Generation The domain of interest consists of a square 2D geometry. Although the flow field is two- dimensional, we will create a three-dimensional domain here and apply periodic boundary conditions in the z−direction. Also, we will only use one element in that direction to save computational costs. In the tutorial directory, we provide the necessary mesh files, along with a parameter files for **HOPR** to generate these meshes. You can recreate any mesh by running the following command. A full tutorial on how to run **HOPR** is available at the [HOPR documentation][hopr]. ```bash hopr parameter_hopr.ini ``` For this basic tutorial, the simple meshes shown in ({numref}`fig:cavity_meshes`) will be used. ```{figure} ./figures/Cavity_Re100_Meshes.jpg :name: fig:cavity_meshes :align: center :width: 45% :alt: Meshes for the basic lid-driven cavity tutorial. Meshes for the basic lid-driven cavity tutorial. ``` ### Simulation Parameters The parameter file to run the simulation is supplied as `parameter_flexi.ini`. The parameters specific to each topic can be found in the labeled section of the file. ###### Output ```ini ! ============================================================ ! ! OUTPUT ! ============================================================ ! ProjectName = Tutorial_Cavity_Re100 OutputFormat = 0 ``` The `ProjectName` parameter defines the prefix for all files generated by the simulation. For example, if the simulation saves the state at time $0.3$, it will produce a file named `Tutorial_Cavity_Re100_State_0000000.300000000.h5` which contains the conserved variable values at each node. In this tutorial, the `OutputFormat` parameter is set to $0$, which disables on-the-fly output. Although turning on visualization output provides real-time insight, it’s generally not recommended because it can significantly slow down the code, especially in parallel executions. Regardless of this setting, HDF5 state files are always generated. For post-simulation visualization, the recommended approach is to use the **posti_visu** tool to create **ParaView**-compatible files. Currently, only the VTK output format is supported (by setting `OutputFormat=3`), as Tecplot output is unavailable due to GPL licensing restrictions. ###### Interpolation ```ini ! ============================================================ ! ! INTERPOLATION ! ============================================================ ! N = 3 NAnalyze = 10 ``` The parameter `N` sets the degree of the solution polynomial. In this example, the solution is approximated by a polynomial of degree $3$ in each spatial direction. This results in $(N+1)^3 = 64$ degrees of freedom for each (3D) element. In general, *N* can be chosen to be any integer greater or equal to $1$, however, the discretization and the timestep calculation has not extensively been tested beyond $N\approx 23$. Usually, for a good compromise of performance and accuracy is found for $N\in[3,..,9]$. `NAnalyze` determines the polynomial degree for the analysis routines, i.e., the accuracy of the calculation of error norms or test case specific integrals during the computation. A good rule of thumb is to set $NAnalyze=2\times N$. ###### Numerical Mesh ```ini ! ============================================================ ! ! MESH ! ============================================================ ! MeshFile = cavity2x2_mesh.h5 useCurveds = F BoundaryName = BC_free BoundaryType = (/2,1/) BoundaryName = BC_wall_left BoundaryType = (/4,1/) BoundaryName = BC_wall_right BoundaryType = (/4,1/) ! BoundaryName = BC_wall_lower ! BoundaryType = (/4,1/) ``` The parameter `MeshFile` contains the name of the **HOPR** mesh file in HDF5 format (and/or the full path to it). `UseCurveds` indicates whether the mesh is considered to be curved, i.e., if high-order mesh information should be used. Setting this to `F` can be used to discard high-order information in the mesh file and treat it as a linear mesh. For the current tutorial, the meshes are linear by design. The boundary conditions are set via the `BoundaryName` identifier, which must be present in the mesh file and thus match the `BoundaryName` identified used during mesh creation. Each line containing the boundary name must be followed by a line containing the `BoundaryType` identified that what kind of boundary is to be applied to the face. A list of types available for the Navier-Stokes equations can be found in table {numref}`tab:boundaryconditions`. For types that require additional information (like Dirichlet boundaries), the second index in `BoundaryType` refers to the `RefState` (short for reference state) which is used to determine the unknowns / quantities from the outside for this boundary condition. For example, for a Dirichlet inflow boundary (Type 2), the full reference state is set at the boundary. Here, `BoundaryType= (/2,1/)` indicates that the *first* reference state vector listed below is set at this boundary (the lid part of the cavity). Note that the reference vectors are always in primitive variables, i.e., $(\rho, u, v, w, p)^T$ *unless* specified otherwise. The same reference state is also chosen for the boundaries `BC_wall_left` and `BC_wall_right`. These boundaries selected as isothermal walls (Type 4), so a wall temperature needs to be specified - this is computed from the primitive variables in the associated reference state. For the $z-$oriented faces omitted here, periodic boundary conditions are selected which must be specified in **HOPR** to pre-compute the correct connectivity information. Note that the lines for the lower wall boundary are commented out. In this case, the boundary conditions set in **HOPR** will be retained and not overwritten here. Later, when running **FLEXI** with these settings, it is good practice to inspect the boundary condition information as understood by **FLEXI**. In this case, the output of **FLEXI** to the console should look like the following. ``` ---------------------------------------------------------------- | BoundaryName | BC_wall_left | *CUSTOM | | BoundaryType | (/ 4, 1 /) | *CUSTOM | | BoundaryName | BC_wall_right | *CUSTOM | | BoundaryType | (/ 4, 1 /) | *CUSTOM | | BoundaryName | BC_free | *CUSTOM | | BoundaryType | (/ 2, 1 /) | *CUSTOM | | Boundary in HDF file found | BC_free | was | 2 0 | is set to | 2 1 | Boundary in HDF file found | BC_wall_left | was | 4 1 | is set to | 4 1 | Boundary in HDF file found | BC_wall_right | was | 4 1 | is set to | 4 1 ............................................................. BOUNDARY CONDITIONS| Name Type State Alpha | | BC_zminus 1 0 1 | | BC_zplus 1 0 -1 | | BC_wall_lower 4 1 0 | | BC_free 2 1 0 | | BC_wall_left 4 1 0 | | BC_wall_right 4 1 0 ............................................................. ``` ###### Equation System ```ini ! ============================================================ ! ! EQUATION ! ============================================================ ! IniExactFunc = 1 IniRefState = 2 RefState = (/1.0,1.,0.,0.,71.4285714286/) RefState = (/1.0,0.,0.,0.,71.4285714286/) mu0 = 0.01 R = 1 Pr = 0.72 kappa = 1.4 ``` The equation system in use is set during compilation. However, these equations are unclosed without initial conditions and reference data for the boundary conditions. The parameter `IniExactFunc` specifies which solution or function should be used to fill the initial solution vector, i.e., it specifies what the starting flow field looks like. Setting this to $1$ selects a uniform initial state in the whole domain. Note that this solution is also used to compute the errors norms and can also be time-dependent. The reference state itself used for the initial function is defined by the parameter `IniRefstate`, in this case, the second one is used. As described above, each `RefState` is given in primitive variables. Here, the second reference state describes a fluid at rest and is used to initialize a resting fluid in the cavity. The first state is used to determine the driving flow at the top of the cavity (and to compute the wall temperatures for the boundaries). Constant flow properties like the gas constant correspond to the values given in {numref}`tab:cavity_flow_prop`. These define the gas dynamics in combination with the ideal gas law $p=\rho R T$. ```{important} **FLEXI** does not distinguish between dimensional and non-dimensional quantities. It is the user's responsibility to set all data consistently. For anything other than an ideal gas with constant viscosity and heat conductivity, physically meaningful quantities should be set. ``` ```{list-table} Material properties for the lid-driven cavity tutorial. :header-rows: 1 :name: tab:cavity_flow_prop :align: center :width: 100% :widths: 30 30 40 * - Variable - Value - Description * - mu0 - 0.1 - Dynamic viscosity $\mu$ * - R - 1 - Ideal gas constant $R$ * - kappa - 1.4 - Isentropic coefficient $\kappa$ * - Pr - 0.72 - Prandtl number $\Pr$ ``` From these settings, the Mach and Reynolds number can be computed as follows, taking into account a reference cavity length of $1$ and the magnitude of the driving velocity. Since we are comparing against an incompressible reference solution, setting the Mach number to $0.1$ is a good compromise between accuracy and efficiency of the explicit time integration. ```{math} Mach = u/c=1.0/\sqrt{\kappa \frac{p}{\rho}}=1.0/10.0=0.1 ``` ```{math} Re = \frac{u L \rho}{\mu_0}=\frac{1.0}{0.01}=100 ``` ### Temporal Discretization ```ini ! ============================================================ ! ! TIMEDISC ! ============================================================ ! TEnd = 5.0 Analyze_dt = 0.1 nWriteData = 1 CFLscale = 0.9 DFLscale = 0.4 ``` The parameter `TEnd` determines the end time of the solution, `Analyze_dt` the interval at which the analysis routines (like error computation, checking of wall velocities etc. see below) are called. The multiplier `nWriteData` determines the interval at which full solution state files in HDF5 format are written to the file system, e.g., in this case `nWriteData` $\times$ `Analyze_dt`$ = 0.1$ is the interval for writing to disc. The CFL and DFL numbers determine the explicit time step restriction for the advective and viscous parts. Note that these values should always be chosen to be $<1$. However, since the determination of the timestep includes some heuristics, both values might require to be chosen even more conservatively. (sec:tut_cavity_running_and_results)= ### Simulation and Results We proceed by running the code with the following command. ```bash flexi parameter_flexi.ini ``` If **FLEXI** was compiled with MPI support, it can also be run in parallel with the following command. Here, `` is an integer denoting the number of processes to be used in parallel. ```bash mpirun -np flexi parameter_flexi.ini ``` ```{important} **FLEXI** uses an element-based domain decomposition approach for parallelization. Consequently, the minimum load per process is *one* grid element, i.e. do not use more processes than grid elements! ``` Running the code prints all output to `STDOUT`. If the run completes successfully, the last lines should appear similar to the following (condensed) output. After a successful run, the directory should contain additional generated files named `_State_.h5`. Each file stores the solution vector of the conserved variables at each interpolation node at a specific time, corresponding to multiples of `Analyze_dt`. ``` -------------------------------------------------------------- Sys date : 07.11.2024 13:02:23 CALCULATION TIME PER STAGE/DOF: [ 4.64195E-07 sec ] EFFICIENCY: CALCULATION TIME [s]/[Core-h]: [ 1.23651E+04 sec/h ] Timestep : 4.1665427E-03 #Timesteps : 1.2000000E+03 WRITE STATE TO HDF5 FILE... DONE! [ 0.00 sec ] Sim time : 5.000E+00 L_2 : 1.269E-03 1.810E-01 1.127E-01 5.762E-14 1.568E-01 L_inf : 8.090E-03 9.497E-01 3.924E-01 2.020E-13 1.713E+00 BodyForces (Pressure, Friction) : BC_wall_lower 0.00E+00 -7.14E+01 0.00E+00 -2.73E-03 1.80E-05 3.63E-15 BC_wall_left -7.13E+01 0.00E+00 0.00E+00 2.77E-04 1.91E-02 3.34E-15 BC_wall_right 7.14E+01 -8.74E-15 -1.61E-15 9.38E-04 -4.15E-02 -1.78E-15 Wall Velocities (mean/min/max) : BC_wall_lower 1.138E-02 3.164E-03 2.183E-02 BC_wall_left 3.604E-02 2.169E-03 1.488E-01 BC_wall_right 1.114E-01 2.742E-03 3.930E-01 MeanFlux through boundaries : BC_zminus 3.539E-14 -1.047E-14 -5.820E-15 7.139E+01 8.846E-12 BC_wall_lower 0.000E+00 -2.737E-03 -7.140E+01 3.633E-15 8.844E-04 BC_free 3.913E-05 -1.048E-01 7.142E+01 2.564E-15 -4.560E-02 BC_wall_left 0.000E+00 -7.137E+01 1.914E-02 3.345E-15 1.024E-02 BC_wall_right 0.000E+00 7.148E+01 -4.159E-02 -3.400E-15 3.555E-02 -------------------------------------------------------------- Time = 0.5000E+01 dt = 0.4169E-02 ETA [d:h:m]:<1 min |==>| [100.00%] ============================================================= FLEXI RUNNING Tutorial_Cavity_Re100_mesh2x2...! [ 0.38 sec ] [ 0:00:00:00 ] ============================================================ ============================================================ FLEXI FINISHED!! [ 0.38 sec ] [ 0:00:00:00 ] ============================================================ ``` Since we start the simulation from a fluid at rest, it will take some iterations / time steps to achieve a steady state solution. One way to check if the solution has converged to a steady state is to check some characteristic quantities. Note that since the boundary conditions are applied *weakly* in a DG setting, a velocity slip at walls can occur, with its magnitude depending on the local wall resolution. Thus, the velocities at the walls offer a measure of the solution convergence. In {numref}`fig:cavity_re100_wallvel`, the temporal evolution of the velocities at the lower wall are plotted over time for $4$ different simulations. For all cases, a sufficiently stationary solution has been achieved at $t_{end}=5$. ```{figure} ./figures/Cavity_Re100_WallVel_BC_wall_lower.jpg :name: fig:cavity_re100_wallvel :align: center :width: 70% :alt: Time evolution of wall velocity at lower wall for the $Re=100$ lid-driven cavity case. Time evolution of wall velocity at lower wall for the $Re=100$ lid-driven cavity case. ``` A contour plot of the velocity magnitude at the end time is given in {numref}`fig:cavity_re100_velcontours`. To generate this plot, convert the *State* files to a **ParaView** format using the **posti_visu** tool. ```{figure} ./figures/Re100_streamlines_cropped.jpg :name: fig:cavity_re100_velcontours :align: center :width: 45% :alt: Contours of velocity magnitude for the $Re=100$ lid-driven cavity case. Contours of velocity magnitude for the $Re=100$ lid-driven cavity case. ``` For a more quantitative comparison with published data, you can generate a plot of the $u$-velocity on the centerline ($x=0$) of the cavity. {numref}`fig:cavity_re100_u_over_y` shows the results for the 4 simulations run here, along with published data available in {cite}`ghia1982high`, {cite}`gao2016efficient`. ```{figure} ./figures/cavity_u_over_y_Re100_plot.jpg :name: fig:cavity_re100_u_over_y :align: center :width: 45% :alt: Comparison of centerline velocities for the $Re=100$ lid-driven cavity case with literature. Comparison of centerline velocities for the $Re=100$ lid-driven cavity case with literature. ``` For simulation $1$, the agreement with literature results is fair. This is due to the coarse resolution with $2\times (3+1) = 8$ degrees of freedom in $x$- and $y$-direction. Doubling the grid elements results in visible improvement. Doubling the grid resolution again (simulation $3$), the agreement with the published data is excellent. It should be noted that the same accuracy can be achieved by increasing $N$ and keeping the coarse grid. Simulation $3$ and $4$ have nearly identical results, although the number of degrees of freedom differs by a factor of $2$. This is an indication of the excellent convergence properties of high-order schemes for smooth problems. (sec:tut_cavity_advanced)= ## Advanced Tutorial | Flow at Re=400 In this section, we build on the concepts covered in the basic tutorial. While the general setup of the simulation remains the same, we increase the Reynolds number, which requires a new, higher-resolution mesh to capture the finer flow details. Additionally, this part introduces basic code customization by showing how to add a new function for custom initial or boundary conditions. Before diving in, it is recommended that you have completed the basic tutorial, have access to at least four computational cores (or be prepared for longer run times), and be comfortable with the modern **Fortran** syntax. The basic tutorial is contained in the `Advanced_Re400` subfolder of the `tutorials\cavity` directory. (sec:tut_cav_meshgen_adv)= ### Mesh Generation To account for the increased Reynolds number, the number of elements in the $x-y$-plane is increased to $12\times 12$. Also, a stretching in the $y-$direction is introduced, as depicted in {numref}`fig:cavity_re400_stretchmesh`. You can generate your own mesh or re-use the provided one, labeled `cavity12x12_stretch_mesh.h5`. In `parameter_hopr.ini`, the following line introduces an exponential stretching. ```ini factor = (/1.0,-1.2,1./) ! stretching with constant growth factor ! (+/- changes direction) ``` ```{figure} ./figures/Cavity_Re400_Mesh.jpg :name: fig:cavity_re400_stretchmesh :align: center :width: 45% :alt: Stretched mesh for the $Re=400$ lid-driven cavity case. Stretched mesh for the $Re=400$ lid-driven cavity case. ``` ### Custom Initial / Boundary Function To set up a custom boundary condition for the top boundary, which drives the cavity flow, we define a function for the velocity profile. In this, we follow the suggestions from {cite}`gao2016efficient`, where the $u(x)$ velocity at the lid is given as ```{math} :label: eq:cavity_code u(x)= \begin{cases} c_1 x^4 + c_2 x^3 + c_3 x^2 +c_4 x,& \text{if } 0\le x <0.2\\ d_1 x^4 + d_2 x^3 + d_3 x^2 +d_4 x + d_5,& \text{if } 0.8 < x \le 1.0\\ 1, & \text{otherwise} \end{cases} ``` with ```{math} [c_1,c_2,c_3,c_4]=& 1000\times[4.9333,-1.4267,0.1297,-0.0033]\\ [d_1,d_2,d_3,d_4,d_5]=& 10000\times[0.4933,-1.8307,2.5450,1-.5709,0.3633] ``` This assumes the top boundary to be from $x \in [0, 1]$, as is the case in our domain. To add this new function to **FLEXI**, locate the file `src/equations/navierstokes/idealgas/exactfunc.f90` and open it in the text editor of your choice. Locate the `SUBROUTINE ExactFunc`, which provides functions with analytic solutions for the boundary condition, initial condition, and analysis routines of the code. The header of the routine you are looking for is shown below. ```{code-block} fortran :caption: Header of the ExactFunc subroutine. !============================================================ !> Specifies all the initial conditions. The state in conservative !> variables is returned. t is the actual time. dt is only needed !> to compute the time dependent boundary values for the RK scheme !> for each function resu and the first and second time derivative !> resu_t and resu_tt have to be defined (is trivial for constants) !============================================================ SUBROUTINE ExactFunc(ExactFunction,tIn,x,resu,RefStateOpt) ! MODULES USE MOD_Preproc ,ONLY: PP_PI USE MOD_Globals ,ONLY: Abort USE MOD_Mathtools ,ONLY: CROSS USE MOD_Eos_Vars ,ONLY: Kappa,sKappaM1,KappaM1,KappaP1,R USE MOD_Exactfunc_Vars ,ONLY: IniCenter,IniHalfwidth,IniAmplitude,IniFrequency,IniAxis,AdvVel USE MOD_Exactfunc_Vars ,ONLY: MachShock,PreShockDens USE MOD_Exactfunc_Vars ,ONLY: P_Parameter,U_Parameter USE MOD_Exactfunc_Vars ,ONLY: JetRadius,JetEnd,JetAmplitude USE MOD_Equation_Vars ,ONLY: IniRefState,RefStateCons,RefStatePrim USE MOD_Timedisc_Vars ,ONLY: fullBoundaryOrder,CurrentStage,dt,RKb,RKc,t USE MOD_TestCase ,ONLY: ExactFuncTestcase USE MOD_EOS ,ONLY: PrimToCons,ConsToPrim #if PARABOLIC USE MOD_Eos_Vars ,ONLY: mu0 USE MOD_Exactfunc_Vars ,ONLY: delta99_in,x_in #endif IMPLICIT NONE !------------------------------------------------------------- ! INPUT/OUTPUT VARIABLES INTEGER,INTENT(IN) :: ExactFunction REAL,INTENT(IN) :: x(3) REAL,INTENT(IN) :: tIn REAL,INTENT(OUT) :: Resu(PP_nVar) INTEGER,INTENT(IN),OPTIONAL :: RefStateOpt !------------------------------------------------------------- ! LOCAL VARIABLES INTEGER :: RefState REAL :: tEval REAL :: Resu_t(PP_nVar),Resu_tt(PP_nVar),ov REAL :: Frequency,Amplitude REAL :: Omega REAL :: Vel(3),Cent(3),a REAL :: Prim(PP_nVarPrim) REAL :: r_len REAL :: Ms,xs REAL :: Resul(PP_nVar),Resur(PP_nVar) REAL :: random REAL :: du, dTemp, RT, r2 REAL :: pi_loc,phi,radius REAL :: h,sRT,pexit,pentry #if PARABOLIC ! needed for blasius BL INTEGER :: nSteps,i REAL :: eta,deta,deta2,f,fp,fpp,fppp,fbar,fpbar,fppbar,fpppbar REAL :: x_eff(3),x_offset(3) #endif !============================================================= ``` To add equation {eq}`eq:cavity_code` to the code, add a new `CASE` to the routine. In this `CASE`, define the state vector for the primitive variables, and then convert them to conservative ones (see e.g., `CASE(8)` for how this is done). You might also need to introduce some new local variables for this routine. To check if your changes are syntactically correct, compile the code with your changes. If the compilation process was not successful, check the compiler output for any error messages to diagnose the issue. To benefit from this tutorial, it is recommended that you do try to complete this programming task. For reference, the listing below shows an example code for a correct implementation as `CASE(9)`. Note that in the example code, we do not specify the full primitive state vector within the routine, but re-use the `IniRefState` and just overwrite the $u$-velocity. This is a matter of choice, but it allows imposing the Mach number by setting the `IniRefState` accordingly. ```{code-block} fortran :caption: Example code for a correct implementation from {cite}`gao2016efficient`. CASE(9) ! Lid driven cavity flow from Gao, Hesthaven, Warburton ! "Absorbing layers for weakly compressible flows", JSC, 2016 ! Special "regularized" driven cavity BC to prevent singularities ! at corners. Top BC assumed to be in x-direction from 0..1 Prim = RefStatePrim(:,RefState) IF (x(1).LT.0.2) THEN prim(VEL1)=1000*4.9333*x(1)**4-1.4267*1000*x(1)**3+0.1297*1000*x(1)**2-0.0033*1000*x(1) ELSEIF (x(1).LE.0.8) THEN prim(VEL1)=1.0 ELSE prim(VEL1)=1000*4.9333*x(1)**4-1.8307*10000*x(1)**3+2.5450*10000*x(1)**2-1.5709*10000*x(1)+10000*0.3633 ENDIF CALL PrimToCons(prim,Resu) ``` ```{note} From now on, we refer to the above `CASE(9)` as the case number in question while yours might differ. ``` (sec:tut_cav_prepflowsim_adv)= ### Simulation Parameters To set up the simulation, you can either modify the `parameter_flexi.ini` files from the basic tutorial or use the ones provided in the `Advanced_Re400` directory. We will conduct two simulations for the advanced configuration. The first one uses the constant driving flow boundary conditions as before while the second one utilizes the new custom equation specified as equation `CASE(9)`. In both cases, we need to modify the mesh file, the fluid viscosity (in order to set the Reynolds number), and the end time of the simulation. Here, we chose `TEnd=100` to account for the longer accommodation period required for the simulation to reach a quasi-steady state. As will be seen later, $100$ is very conservative, so feel free to lower this value as you see fit. ```ini MeshFile = cavity12x12_stretch_mesh.h5 [...] mu0 = 0.0025 [...] TEnd = 100.0 ``` For the custom boundary condition case, the parameter file needs to be adjusted to include the new boundary settings. According to {numref}`tab:boundaryconditions`, Dirichlet boundary conditions with a specified reference equation (instead of a reference _state_) are of type $22$. The second index in the entry thus refers to the equation (`CASE(9)`) which we have programmed above. ```ini BoundaryName = BC_free BoundaryType = (/22,9/) ``` (sec:tut_cav_simres_adv)= ### Simulation and Results We proceed by running the code with the following command. If you chose to use the provided .ini files, change the parameter file to `parameter_flexi_custombc.ini` for the second run. ```bash flexi parameter_flexi.ini ``` If **FLEXI** was compiled with MPI support, it can also be run in parallel with the following command. Here, `` is an integer denoting the number of processes to be used in parallel. ```bash mpirun -np flexi parameter_flexi.ini ``` The evolution of the mean velocity of the lower wall is given in {numref}`fig:cavity_re400_wallvel`. Note the difference in the axis scales compared to {numref}`fig:cavity_re100_wallvel`. It is evident that the flow reaches steady state after about $t=35$. ```{figure} ./figures/Cavity_Re400_WallVel_BC_wall_lower.jpg :name: fig:cavity_re400_wallvel :align: center :width: 70% :alt: Time evolution of wall velocity at lower wall for the $Re=400$ lid-driven cavity case. Time evolution of wall velocity at lower wall for the $Re=400$ lid-driven cavity case. ``` In the following, {numref}`fig:cavity_re400_velcomp` and {numref}`fig:cavity_re400_u_over_y` show the flow field and comparison of the centerline velocities with published results {cite}`ghia1982high`, {cite}`gao2016efficient`. ```{figure} ./figures/cavity_Re400_constant_costum.jpg :name: fig:cavity_re400_velcomp :align: center :width: 70% :alt: Steady state solution of velocity magnitude of $Re=400$ lid driven cavity. Left: constant boundary condition, right: custom boundary condition. Steady state solution of velocity magnitude of $Re=400$ lid driven cavity. Left: constant boundary condition, right: custom boundary condition. ``` ```{figure} ./figures/cavity_u_over_y_Re400_plot.jpg :name: fig:cavity_re400_u_over_y :align: center :width: 45% :alt: Evolution of wall velocities at the lower wall for $Re=400$ lid driven cavity simulations. Evolution of wall velocities at the lower wall for $Re=400$ lid driven cavity simulations. ``` [hopr]: https://hopr.readthedocs.io/en/latest/