Skip to content
Snippets Groups Projects
Commit ef376aaf authored by Mohsen Lahooti's avatar Mohsen Lahooti
Browse files

an xml tag is corrected

parent d1bb8476
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
02-incns-solver/incns-solver.ipynb
+ 72
− 72
View file @ ef376aaf
%% Cell type:markdown id:214d1e5a tags: %% Cell type:markdown id:0b018ab6 tags:
# **Incompressible Flow Simulation** # **Incompressible Flow Simulation**
%% Cell type:markdown id:fb7db30b tags: %% Cell type:markdown id:98e5cc81 tags:
## Learning objectives <a id='learningObjectives'></a> ## Learning objectives <a id='learningObjectives'></a>
Welcome to the tutorial on the incompressible flow simulation using the Nektar++ framework. This tutorial aims to show how to setup and run a simulation for incompressible flows using the Nektar++ `IncNavierStokesSolver` for the well-known problem of external flow over a cylinder. Welcome to the tutorial on the incompressible flow simulation using the Nektar++ framework. This tutorial aims to show how to setup and run a simulation for incompressible flows using the Nektar++ `IncNavierStokesSolver` for the well-known problem of external flow over a cylinder.
After completetion of this tutorial you will be familiar with: After completetion of this tutorial you will be familiar with:
- Configuring a XML file suitable for incompressible Navier-Stokes in Nektar++. - Configuring a XML file suitable for incompressible Navier-Stokes in Nektar++.
- Simulating 2D and quasi-3D incompressible flows. - Simulating 2D and quasi-3D incompressible flows.
- Running the simulation in serial and parallel. - Running the simulation in serial and parallel.
- Computing aerodynamic forces and modal energy as a function of time. - Computing aerodynamic forces and modal energy as a function of time.
- Post-processing the results and converting the Nektar++ field outputs for visualization. - Post-processing the results and converting the Nektar++ field outputs for visualization.
%% Cell type:markdown id:d4e7482f tags: %% Cell type:markdown id:348ebe82 tags:
# Introduction <a id='introduction'></a> # Introduction <a id='introduction'></a>
In this tutorial we will walk through how to configure a session file to be able to solve incompressible Navier-Stokes equations using a simple flow example. We consider flow over a circular cylinder. At the beginning we consider a two-dimensional case. Upon completion of the 2D case, we extend the problem to a quasi-3D case. By quasi-3D, we mean that the 2D mesh will be extruded in the spanwise direction, with periodic boundary conditions applied in the spanwise direction. In this tutorial we will walk through how to configure a session file to be able to solve incompressible Navier-Stokes equations using a simple flow example. We consider flow over a circular cylinder. At the beginning we consider a two-dimensional case. Upon completion of the 2D case, we extend the problem to a quasi-3D case. By quasi-3D, we mean that the 2D mesh will be extruded in the spanwise direction, with periodic boundary conditions applied in the spanwise direction.
%% Cell type:markdown id:cb5727a3 tags: %% Cell type:markdown id:5ee722c7 tags:
## Governing equation <a id='geqn' ></a> ## Governing equation <a id='geqn' ></a>
The incompressible solver implemented in Nektar++ uses the convective form of the Navier-Stokes equations for viscous incompressible and Newtonian fluids, which can be written as: The incompressible solver implemented in Nektar++ uses the convective form of the Navier-Stokes equations for viscous incompressible and Newtonian fluids, which can be written as:
\begin{equation} \begin{equation}
\frac{\partial \mathbf{V}}{\partial t}+ \mathbf{V} \cdot \nabla \mathbf{V}=-\nabla p + \nu \nabla ^{2} \mathbf{V} + \mathbf{f} \frac{\partial \mathbf{V}}{\partial t}+ \mathbf{V} \cdot \nabla \mathbf{V}=-\nabla p + \nu \nabla ^{2} \mathbf{V} + \mathbf{f}
\end{equation} \end{equation}
\begin{equation} \begin{equation}
\nabla \cdot \mathbf{V} =0 \nabla \cdot \mathbf{V} =0
\end{equation} \end{equation}
where $\mathbf{V} = (u,v,w)$ is the velocity of the fluid in the $x$, $y$, and $z$ directions, respectively, $p$ is the specific pressure (including density), $\nu$ is the kinematic viscosity, and $\mathbf{f}$ is the forcing term. where $\mathbf{V} = (u,v,w)$ is the velocity of the fluid in the $x$, $y$, and $z$ directions, respectively, $p$ is the specific pressure (including density), $\nu$ is the kinematic viscosity, and $\mathbf{f}$ is the forcing term.
%% Cell type:markdown id:beb687a7 tags: %% Cell type:markdown id:eebfd817 tags:
# Problem description: 2D case <a id='pdes' ></a> # Problem description: 2D case <a id='pdes' ></a>
In this tutorial we use air to simulate the flow past a 2D circular cylinder with diameter $D=1$. For our study, we will use a Reynolds number $Re = 40$ with $Re = UD/\nu$ where $U$ is the velocity scale, $D$ is the length scale and $\nu$ is the kinematic viscosity. In this tutorial we use air to simulate the flow past a 2D circular cylinder with diameter $D=1$. For our study, we will use a Reynolds number $Re = 40$ with $Re = UD/\nu$ where $U$ is the velocity scale, $D$ is the length scale and $\nu$ is the kinematic viscosity.
The flow domain is a rectangle of sizes $[-5,10] \times [-5,5]$ as shown in the figure below. The flow domain is a rectangle of sizes $[-5,10] \times [-5,5]$ as shown in the figure below.
<table class="center"> <table class="center">
<tr> <tr>
<td> <img src="CylinderGeometry.png" /> </td> <td> <img src="CylinderGeometry.png" /> </td>
<td> <img src="CylinderGrid.png" /> </td> <td> <img src="CylinderGrid.png" /> </td>
</tr> </tr>
<tr> <tr>
<td style="text-align:center"> Computational Domain </td> <td style="text-align:center"> Computational Domain </td>
<td style="text-align:center"> Computational Mesh </td> <td style="text-align:center"> Computational Mesh </td>
</tr> </tr>
</table> </table>
We apply the following boundary conditions (BCs): We apply the following boundary conditions (BCs):
- no−slip conditions on the cylinder surface (i.e., homogeneous Dirichlet BC $\mathbf{V}=\mathbf{0}$); - no−slip conditions on the cylinder surface (i.e., homogeneous Dirichlet BC $\mathbf{V}=\mathbf{0}$);
- far−field with Dirichlet BC: $(u, v) = (1, 0)$ at the bottom and top boundaries; - far−field with Dirichlet BC: $(u, v) = (1, 0)$ at the bottom and top boundaries;
- a corresonding inflow with Dirichlet BC: $(u, v) = (1, 0)$ at the left boundary; - a corresonding inflow with Dirichlet BC: $(u, v) = (1, 0)$ at the left boundary;
- and an outflow at the right boundary with static pressure $p=0$. - and an outflow at the right boundary with static pressure $p=0$.
<div class="alert alert-warning"> <div class="alert alert-warning">
**Note:** Mathematically we do not require a boundary condition for the pressure $p$, asides from at the outflow. However, due to the velocity-correction scheme that Nektar++ uses, a numerical Neumann BC is required at any boundary where a Dirichlet condition is given for $\mathbf{V}$, in order to preserve the time integration accuracy of the scheme. This is referred to as the 'high-order pressure BC' in the below. **Note:** Mathematically we do not require a boundary condition for the pressure $p$, asides from at the outflow. However, due to the velocity-correction scheme that Nektar++ uses, a numerical Neumann BC is required at any boundary where a Dirichlet condition is given for $\mathbf{V}$, in order to preserve the time integration accuracy of the scheme. This is referred to as the 'high-order pressure BC' in the below.
</div> </div>
%% Cell type:markdown id:ada34927 tags: %% Cell type:markdown id:b81b944e tags:
# Pre-processing <a id='preprocessing'></a> # Pre-processing <a id='preprocessing'></a>
In this section we discuss how to setup the `.xml` files suitable for Nektar++ simulation. Here we use two separate XML files: In this section we discuss how to setup the `.xml` files suitable for Nektar++ simulation. Here we use two separate XML files:
1. A **geometry XML file**, which contains the mesh used to represent the computational domain; 1. A **geometry XML file**, which contains the mesh used to represent the computational domain;
2. A **session** or **conditions** XML file, which contains simulation conditions, problem parameters, filters, etc. 2. A **session** or **conditions** XML file, which contains simulation conditions, problem parameters, filters, etc.
<div class="alert alert-info"> <div class="alert alert-info">
**Note:** It is also possible to use a single xml file instead of two. In that case, one needs to copy the `<GEOMETRY>` tag from the geometry XML file into the session (condition) XML file to make a single XML file. **Note:** It is also possible to use a single xml file instead of two. In that case, one needs to copy the <code>Geometry</code> tag from the geometry XML file into the session (condition) XML file to make a single XML file.
</div> </div>
Please open the condition file [**session2d.xml**](session2d.xml). You will see **three** main tags (`<EXPANSIONS>`,`<CONDITIONS>` and `<FILTERS>`) within the `<NEKTAR>` tag as shown below: Please open the condition file [**session2d.xml**](session2d.xml). You will see **three** main tags (`<EXPANSIONS>`,`<CONDITIONS>` and `<FILTERS>`) within the `<NEKTAR>` tag as shown below:
```xml ```xml
<NEKTAR> <NEKTAR>
<EXPANSIONS> <EXPANSIONS>
<!-- The basis to be used in each element. --> <!-- The basis to be used in each element. -->
</EXPANSIONS> </EXPANSIONS>
<CONDITIONS> <CONDITIONS>
<!-- Solver-specific information, including boundary and initial conditions. --> <!-- Solver-specific information, including boundary and initial conditions. -->
</CONDITIONS> </CONDITIONS>
<FILTERS> <FILTERS>
<!-- Calculate a variety of quantities as the solution evolves in time. --> <!-- Calculate a variety of quantities as the solution evolves in time. -->
</FILTERS> </FILTERS>
</NEKTAR> </NEKTAR>
``` ```
%% Cell type:markdown id:4f4354ed tags: %% Cell type:markdown id:f9bc0470 tags:
## Expansion bases <a id='expansionBases'></a> ## Expansion bases <a id='expansionBases'></a>
First, let us complete the `<EXPANSIONS>` tags. First, let us complete the `<EXPANSIONS>` tags.
<span style="color:Chocolate"><b> Task 1: </b> </span> <span style="color:Chocolate"><b> Task 1: </b> </span>
Please write the missing parts in `<EXPANSIONS>` tag. Please write the missing parts in `<EXPANSIONS>` tag.
A completed `<EXPANSIONS>` tag has the following standard form: A completed `<EXPANSIONS>` tag has the following standard form:
```xml ```xml
<EXPANSIONS> <EXPANSIONS>
<E COMPOSITE="C[100,101]" NUMMODES="3" FIELDS="u,v,p" TYPE="MODIFIED" /> <E COMPOSITE="C[100,101]" NUMMODES="3" FIELDS="u,v,p" TYPE="MODIFIED" />
</EXPANSIONS>` </EXPANSIONS>`
``` ```
Here, Here,
- `COMPOSITE` refers to the ID assigned to the geometric domain (or sub-domain) defined in the geometry xml file. - `COMPOSITE` refers to the ID assigned to the geometric domain (or sub-domain) defined in the geometry xml file.
- `NUMMODES` is used to indentify the polynomial order of the expansion basis functions as $P = \texttt{NUMMODES}-1$. - `NUMMODES` is used to indentify the polynomial order of the expansion basis functions as $P = \texttt{NUMMODES}-1$.
- `TYPE` denotes the basis functions. Possible choices are `MODIFIED`, `GLL_LAGRANGE`,`GLL_LAGRANGE_SEM`. One may note that the first basis type the _modified_ modal expansion of Karniadakis & Sherwin, whereas the remaining two consider nodal expansion using Lagrange basis functions. - `TYPE` denotes the basis functions. Possible choices are `MODIFIED`, `GLL_LAGRANGE`,`GLL_LAGRANGE_SEM`. One may note that the first basis type the _modified_ modal expansion of Karniadakis & Sherwin, whereas the remaining two consider nodal expansion using Lagrange basis functions.
%% Cell type:markdown id:7226e6fd tags: %% Cell type:markdown id:2d35a8c7 tags:
## Conditions setup <a id='conditionsSetup'></a> ## Conditions setup <a id='conditionsSetup'></a>
Now, we will complete the `<CONDITIONS>` tag. There are **five** sections in this tag: Now, we will complete the `<CONDITIONS>` tag. There are **five** sections in this tag:
- `SOLVERINFO` defines solver setup (e.g. the equations to solve); - `SOLVERINFO` defines solver setup (e.g. the equations to solve);
- `PARAMETERS` allow us to define solver parameters, such as $Re$ or the timestep size; - `PARAMETERS` allow us to define solver parameters, such as $Re$ or the timestep size;
- `VARIABLES` define the scalar variables that will be solved for (i.e. $u$, $v$ and $p$); - `VARIABLES` define the scalar variables that will be solved for (i.e. $u$, $v$ and $p$);
- `BOUNDARYREGIONS` and `BOUNDARYCONDITIONS` define the boundary conditions; - `BOUNDARYREGIONS` and `BOUNDARYCONDITIONS` define the boundary conditions;
- finally we need a function to define the initial conditions for the problem. - finally we need a function to define the initial conditions for the problem.
%% Cell type:markdown id:9d8567dd tags: %% Cell type:markdown id:41931951 tags:
### Solver info <a id='solverInfo'></a> ### Solver info <a id='solverInfo'></a>
<span style="color:Chocolate"><b> Task 2: </b> </span> Please complete the properties `SolverType`, `EQTYPE`, `EvolutionOperator` as shown below, <span style="color:Chocolate"><b> Task 2: </b> </span> Please complete the properties `SolverType`, `EQTYPE`, `EvolutionOperator` as shown below,
```xml ```xml
<I PROPERTY="SolverType" VALUE="VelocityCorrectionScheme"/> <I PROPERTY="SolverType" VALUE="VelocityCorrectionScheme"/>
<I PROPERTY="EQTYPE" VALUE="UnsteadyNavierStokes" /> <I PROPERTY="EQTYPE" VALUE="UnsteadyNavierStokes" />
<I PROPERTY="EvolutionOperator" VALUE="Nonlinear" /> <I PROPERTY="EvolutionOperator" VALUE="Nonlinear" />
``` ```
Here, Here,
- `SolverType` sets the algorithm we want to use to solve the governing equations. - `SolverType` sets the algorithm we want to use to solve the governing equations.
- `EQTYPE` sets the kind of equations we want to solve. Here we use `UnsteadyNavierStokes` since we want to solve an unsteady-in-time Navier-Stokes solver; other options might be e.g. `UnsteadyStokes`, but this is not suitable for our problem. - `EQTYPE` sets the kind of equations we want to solve. Here we use `UnsteadyNavierStokes` since we want to solve an unsteady-in-time Navier-Stokes solver; other options might be e.g. `UnsteadyStokes`, but this is not suitable for our problem.
- `EvolutionOperator` sets the choice of the convection operator (i.e. the term $\mathbf{V}\cdot\nabla\mathbf{V}$): - `EvolutionOperator` sets the choice of the convection operator (i.e. the term $\mathbf{V}\cdot\nabla\mathbf{V}$):
- `Nonlinear` for standard non-linear Navier-Stokes equations - `Nonlinear` for standard non-linear Navier-Stokes equations
- `Direct` for linearised Navier-Stokes equations - `Direct` for linearised Navier-Stokes equations
- `Adjoint ` for adjoint Navier-Stokes equations - `Adjoint ` for adjoint Navier-Stokes equations
- `TransientGrowth` for transient growth analysis - `TransientGrowth` for transient growth analysis
- `SkewSymmetric` for skew symmetric form of the Navier-Stokes equations. - `SkewSymmetric` for skew symmetric form of the Navier-Stokes equations.
%% Cell type:markdown id:0d9a2a43 tags: %% Cell type:markdown id:29d8bcc8 tags:
### Parameters <a id='parameters'></a> ### Parameters <a id='parameters'></a>
In `<PARAMETERS>` tag, we can define the parameters present in the flow problem as well as the parameters required by various solvers of Nektar++. Hence a parameter can either be defined in the context of the problem we are solving (a 'dummy' or utility variable), or be predefined within Nektar++ (solver parameter). One should note that the names of the solver parameters can't be altered, and some solvers will expect certain parameters to be defined. In `<PARAMETERS>` tag, we can define the parameters present in the flow problem as well as the parameters required by various solvers of Nektar++. Hence a parameter can either be defined in the context of the problem we are solving (a 'dummy' or utility variable), or be predefined within Nektar++ (solver parameter). One should note that the names of the solver parameters can't be altered, and some solvers will expect certain parameters to be defined.
<span style="color:Chocolate"><b> Task 3: </b> </span> Assign values to the following parameters as <span style="color:Chocolate"><b> Task 3: </b> </span> Assign values to the following parameters as
```xml ```xml
<P> TimeStep = 0.025 </P> <P> TimeStep = 0.025 </P>
<P> FinTime = 20.0 </P> <P> FinTime = 20.0 </P>
<P> IO_CheckSteps = NumSteps/20 </P> <P> IO_CheckSteps = NumSteps/20 </P>
<P> Re = 40 </P> <P> Re = 40 </P>
``` ```
Here, Here,
- `TimeStep` sets the time-step size for the time integration. One should select magnitude of `TimeStep` strategically to avoid the situation of CFL (Courant-Friedrichs-Lewy) number to be too large! - `TimeStep` sets the time-step size for the time integration. One should select magnitude of `TimeStep` strategically to avoid the situation of CFL (Courant-Friedrichs-Lewy) number to be too large!
- `FinTime` is used for define the final time of the simulation. Equivalently, one could specify `NumSteps` to set the number of time-steps for the simulation. Given two of `TimeStep`, `FinTime` and `NumSteps`, Nektar++ will calculate the remaining quantity automatically. - `FinTime` is used for define the final time of the simulation. Equivalently, one could specify `NumSteps` to set the number of time-steps for the simulation. Given two of `TimeStep`, `FinTime` and `NumSteps`, Nektar++ will calculate the remaining quantity automatically.
- `IO_CheckSteps` sets the number of steps between successive checkpoint files (which are saved as files ending in the extension `.chk`) to save the solution at equal-interval intermediate time levels. For example, since `TimeStep` is `0.025` and `FinTime` is `20.0`, we have that `NumSteps` is `800`. Since `IO_CheckSteps` is then set as `NumSteps/20` which will generate twenty successive checkpoint files. - `IO_CheckSteps` sets the number of steps between successive checkpoint files (which are saved as files ending in the extension `.chk`) to save the solution at equal-interval intermediate time levels. For example, since `TimeStep` is `0.025` and `FinTime` is `20.0`, we have that `NumSteps` is `800`. Since `IO_CheckSteps` is then set as `NumSteps/20` which will generate twenty successive checkpoint files.
- `Re` denotes the Reynolds number used for the simulation. This is not explicitly required by the solver, i.e. it is a 'dummy' parameter. - `Re` denotes the Reynolds number used for the simulation. This is not explicitly required by the solver, i.e. it is a 'dummy' parameter.
- `Kinvis` denotes kinematic viscosity of the fluid used and is required as a solver parameter. One may directly supply the value of `Kinvis` without defining `Re` in the XML file, but generally it is more convenient to adjust $Re$ as the quantity of interest. - `Kinvis` denotes kinematic viscosity of the fluid used and is required as a solver parameter. One may directly supply the value of `Kinvis` without defining `Re` in the XML file, but generally it is more convenient to adjust $Re$ as the quantity of interest.
%% Cell type:markdown id:c21edca9 tags: %% Cell type:markdown id:ba2a0976 tags:
### Variables <a id='variables'></a> ### Variables <a id='variables'></a>
In `<VARIABLES>` tag, we define the number (and name) of solution variables. We prescribed a unique ID to each variable. In `<VARIABLES>` tag, we define the number (and name) of solution variables. We prescribed a unique ID to each variable.
<span style="color:Chocolate"><b> Task 4: </b> </span> Provide a ID for the v-velocity component as <span style="color:Chocolate"><b> Task 4: </b> </span> Provide a ID for the v-velocity component as
```xml ```xml
<V ID="1">v</V> <V ID="1">v</V>
``` ```
%% Cell type:markdown id:8e2cbe01 tags: %% Cell type:markdown id:c50dbd5a tags:
### Boundary conditions <a id='bc'></a> ### Boundary conditions <a id='bc'></a>
Boundary conditions are defined by two tags. First we define the boundary regions in the `<BOUNDARYREGIONS>` tag in terms of composite `C[ ]` entities from the `<GEOMETRY>` tag section of the [mesh xml](cyl2DMesh.xml) file. Each boundary region has a unique ID. For the current computational domain there are eight composites (see towards the bottom section of [mesh xml](cyl2DMesh.xml) file) marked as `C[1]` to `C[8]`. Boundary conditions are defined by two tags. First we define the boundary regions in the `<BOUNDARYREGIONS>` tag in terms of composite `C[ ]` entities from the `<GEOMETRY>` tag section of the [mesh xml](cyl2DMesh.xml) file. Each boundary region has a unique ID. For the current computational domain there are eight composites (see towards the bottom section of [mesh xml](cyl2DMesh.xml) file) marked as `C[1]` to `C[8]`.
<span style="color:Chocolate"><b> Task 5: </b> </span> Define the farfield region using `C[2]` and `C[4]` as shown below <span style="color:Chocolate"><b> Task 5: </b> </span> Define the farfield region using `C[2]` and `C[4]` as shown below
```xml ```xml
<B ID="3"> C[2,4] </B> <!-- farfield --> <B ID="3"> C[2,4] </B> <!-- farfield -->
``` ```
The next tag we define is `<BOUNDARYCONDITIONS>` tag, in which specify the boundary conditions for each boundary ID specified in the above `<BOUNDARYREGIONS>` tag. The next tag we define is `<BOUNDARYCONDITIONS>` tag, in which specify the boundary conditions for each boundary ID specified in the above `<BOUNDARYREGIONS>` tag.
<span style="color:Chocolate"><b> Task 6: </b> </span> Define $u=1$, $v=0$ and higher-order pressure boundary condition for the inflow section as <span style="color:Chocolate"><b> Task 6: </b> </span> Define $u=1$, $v=0$ and higher-order pressure boundary condition for the inflow section as
```xml ```xml
<REGION REF="1"> <!-- inflow --> <REGION REF="1"> <!-- inflow -->
<D VAR="u" VALUE="1" /> <D VAR="u" VALUE="1" />
<D VAR="v" VALUE="0" /> <D VAR="v" VALUE="0" />
<N VAR="p" VALUE="0" USERDEFINEDTYPE="H" /> <N VAR="p" VALUE="0" USERDEFINEDTYPE="H" />
</REGION> </REGION>
``` ```
Note the use of `USERDEFINEDTYPE="H"`, which tells the `IncNavierStokesSolver` to set the high-order pressure boundary condition for $p$. In this case the `VALUE` attribute is ignored. Note the use of `USERDEFINEDTYPE="H"`, which tells the `IncNavierStokesSolver` to set the high-order pressure boundary condition for $p$. In this case the `VALUE` attribute is ignored.
%% Cell type:markdown id:34a4461e tags: %% Cell type:markdown id:2f397f7d tags:
### Initial conditions <a id='ic'></a> ### Initial conditions <a id='ic'></a>
<span style="color:Chocolate"><b> Task 7: </b> </span> Provide the initial conditions $u(t=0)=1$, $v(t=0)=0$ and $p(t=0)=0$ in the "InitialConditions" `FUNCTION` as, <span style="color:Chocolate"><b> Task 7: </b> </span> Provide the initial conditions $u(t=0)=1$, $v(t=0)=0$ and $p(t=0)=0$ in the "InitialConditions" `FUNCTION` as,
```xml ```xml
<E VAR="u" VALUE="1" /> <E VAR="u" VALUE="1" />
<E VAR="v" VALUE="0" /> <E VAR="v" VALUE="0" />
<E VAR="p" VALUE="0" /> <E VAR="p" VALUE="0" />
``` ```
%% Cell type:markdown id:1e9e49b0 tags: %% Cell type:markdown id:510769ee tags:
## Special parameters for solution stabilisation ## Special parameters for solution stabilisation
One of the main challenges of running fluid simulations at high-order is numerical stability, since the low dissipation and diffusion at high orders can lead to accumulation of numerical errors in the presence of only mild under-resolution. The mesh we will use today is deliberably under-resolved to make it cheap to simulate. For this reason, we need to now enable some parameters in the `<SOLVERINFO>` tag to allow for the stabilisation of the incompressible flow solver. One of the main challenges of running fluid simulations at high-order is numerical stability, since the low dissipation and diffusion at high orders can lead to accumulation of numerical errors in the presence of only mild under-resolution. The mesh we will use today is deliberably under-resolved to make it cheap to simulate. For this reason, we need to now enable some parameters in the `<SOLVERINFO>` tag to allow for the stabilisation of the incompressible flow solver.
%% Cell type:markdown id:65e2fb26 tags: %% Cell type:markdown id:4e91dd3b tags:
### Spectral vanishing viscosity ### Spectral vanishing viscosity
```xml ```xml
<I PROPERTY="SpectralVanishingViscosity" VALUE="DGKernel" /> <I PROPERTY="SpectralVanishingViscosity" VALUE="DGKernel" />
``` ```
- `SpectralVanishingViscosity` activates a stabilization technique that increases the viscosity on the modes with the highest frequencies. Possible options availabe are `True`, `ExpKernel`, `DGKernel` for exponential, power and DG Kernel, respectively. - `SpectralVanishingViscosity` activates a stabilization technique that increases the viscosity on the modes with the highest frequencies. Possible options availabe are `True`, `ExpKernel`, `DGKernel` for exponential, power and DG Kernel, respectively.
- To override the default values, we define the parameters `SVVCutoffRatio` in the `<PARAMETERS>` tag as - To override the default values, we define the parameters `SVVCutoffRatio` in the `<PARAMETERS>` tag as
```xml ```xml
<P> SVVCutoffRatio = 0.5 </P> <P> SVVCutoffRatio = 0.5 </P>
``` ```
To explore more about spectral vanishing viscoity one may consult with Ref. [2] and Section 11.3.1 of Ref. [1]. To explore more about spectral vanishing viscoity one may consult with Ref. [2] and Section 11.3.1 of Ref. [1].
%% Cell type:markdown id:df160715 tags: %% Cell type:markdown id:8bccf39f tags:
### Spectral dealiasing ### Spectral dealiasing
```xml ```xml
<I PROPERTY="SPECTRALHPDEALIASING" VALUE="True" /> <I PROPERTY="SPECTRALHPDEALIASING" VALUE="True" />
``` ```
- `SPECTRALHPDEALIASING` activates the dealiasing (or over-integration) to stabilize the simulation. This dealiasing technique helps to reduce the numerical instability associated with the integration of the nonlinear terms [3]. One may note that Nektar++ also uses another dealiasing which is associated with the Fourier expansion considered in quasi-3D simulations and activates using the property `DEALIASING`. - `SPECTRALHPDEALIASING` activates the dealiasing (or over-integration) to stabilize the simulation. This dealiasing technique helps to reduce the numerical instability associated with the integration of the nonlinear terms [3]. One may note that Nektar++ also uses another dealiasing which is associated with the Fourier expansion considered in quasi-3D simulations and activates using the property `DEALIASING`.
%% Cell type:markdown id:917620df tags: %% Cell type:markdown id:eee0421c tags:
## Filter setup<a id='filtersSetup'></a> ## Filter setup<a id='filtersSetup'></a>
`<FILTER>` tags are used for calculating a variety of useful quantities from the field variables as the solution evolves in **time**. There are several filters predefined in Nektar++, such as, aerodynamic forces, time-averaged fields and extracting the field variables at certain points inside the domain, to name a few. `<FILTER>` tags are used for calculating a variety of useful quantities from the field variables as the solution evolves in **time**. There are several filters predefined in Nektar++, such as, aerodynamic forces, time-averaged fields and extracting the field variables at certain points inside the domain, to name a few.
Each filter is defined in a `<FILTER>` tag inside a `<FILTERS>` block which lies in the main `<NEKTAR>` tag. Each filter is defined in a `<FILTER>` tag inside a `<FILTERS>` block which lies in the main `<NEKTAR>` tag.
In this tutorial we will discuss how to set up a`AeroForces` filter to calculate aerodynamic forces (such as drag and lift forces) with time. In this tutorial we will discuss how to set up a`AeroForces` filter to calculate aerodynamic forces (such as drag and lift forces) with time.
%% Cell type:markdown id:f6b3f655 tags: %% Cell type:markdown id:f18949ad tags:
### AeroForces filter <a id='aeroForces'></a> ### AeroForces filter <a id='aeroForces'></a>
<span style="color:Chocolate"><b> Task 8: </b> </span> <span style="color:Chocolate"><b> Task 8: </b> </span>
Provide the name of the `OutputFile` as `DragLift`, `OutputFrequency` = 5 and `Boundary` = B[0] since we are interested to calculate aeroforces along the cylinder wall which has a boundary ID=0. Provide the name of the `OutputFile` as `DragLift`, `OutputFrequency` = 5 and `Boundary` = B[0] since we are interested to calculate aeroforces along the cylinder wall which has a boundary ID=0.
```xml ```xml
<PARAM NAME="OutputFile">DragLift</PARAM> <PARAM NAME="OutputFile">DragLift</PARAM>
<PARAM NAME="OutputFrequency">5</PARAM> <PARAM NAME="OutputFrequency">5</PARAM>
<PARAM NAME="Boundary">B[0]</PARAM> <PARAM NAME="Boundary">B[0]</PARAM>
``` ```
Here, Here,
- Parameter `"OutputFile"` is used to name a output file in which data will be stored during the execution. File name can be supplied by the user, here we named it as `DragLift`, hence a file `DragLift.fce` will be created. - Parameter `"OutputFile"` is used to name a output file in which data will be stored during the execution. File name can be supplied by the user, here we named it as `DragLift`, hence a file `DragLift.fce` will be created.
- Parameter `"OutputFrequency"` specifies the number of timesteps after which output is written. - Parameter `"OutputFrequency"` specifies the number of timesteps after which output is written.
- Parameter `"Boundary"` specifies the Boundary surface(s) on which the forces are to be evaluated. For example, `B[0]` means the forces will be evaluated at Boundary ID=0 (cylinder wall) defined in the `<BOUNDARYREGIONS>` tag. - Parameter `"Boundary"` specifies the Boundary surface(s) on which the forces are to be evaluated. For example, `B[0]` means the forces will be evaluated at Boundary ID=0 (cylinder wall) defined in the `<BOUNDARYREGIONS>` tag.
%% Cell type:markdown id:1adb2f72 tags: %% Cell type:markdown id:cdca0c97 tags:
# Running the solver <a id='basic_runningSolver'></a> # Running the solver <a id='basic_runningSolver'></a>
The two xml files provided in this tutorial are: The two xml files provided in this tutorial are:
- Geometry file: [cyl2DMesh.xml](cyl2DMesh.xml) - Geometry file: [cyl2DMesh.xml](cyl2DMesh.xml)
- Condition file: [session2d.xml](session2d.xml) - Condition file: [session2d.xml](session2d.xml)
The flow solvers in Nektar++ can be run either in serial or parallel. First we run the solver in serial. Here we use `!` (the jupyter notebook syntax for bash command) to run the solver in this online interactive tutorial which used Docker imange and Binder to run Nektar++ in the background. The flow solvers in Nektar++ can be run either in serial or parallel. First we run the solver in serial. Here we use `!` (the jupyter notebook syntax for bash command) to run the solver in this online interactive tutorial which used Docker imange and Binder to run Nektar++ in the background.
<span style="color:Chocolate"><b> Task 9: </b> </span> Go to the next cell and click the `Run` icon to run the solver. <span style="color:Chocolate"><b> Task 9: </b> </span> Go to the next cell and click the `Run` icon to run the solver.
<div class="alert alert-warning"> <div class="alert alert-warning">
**Note:** Since the number of steps are 800, it might take a couple of minutes for the simulation to finish. **Note:** Since the number of steps are 800, it might take a couple of minutes for the simulation to finish.
</div> </div>
<div class="alert alert-success"> <div class="alert alert-success">
**Note:** Once you have run this simulation, you should see that simulation produces output files as `cyl2DMesh_N.chk` where $N=1,2,3,\dots$ and `cyl2DMesh.fld`, which is the solution at the final time $t=20$. **Note:** Once you have run this simulation, you should see that simulation produces output files as `cyl2DMesh_N.chk` where $N=1,2,3,\dots$ and `cyl2DMesh.fld`, which is the solution at the final time $t=20$.
</div> </div>
%% Cell type:code id:122b67e9 tags: %% Cell type:code id:86733406 tags:
``` python ``` python
!IncNavierStokesSolver cyl2DMesh.xml session2d.xml !IncNavierStokesSolver cyl2DMesh.xml session2d.xml
``` ```
%% Cell type:markdown id:92cc59b2 tags: %% Cell type:markdown id:1d82a467 tags:
# Post-processing # Post-processing
%% Cell type:markdown id:31e7f047 tags: %% Cell type:markdown id:426e6b7b tags:
## Post-processing the flow field<a id='post_filters'></a> ## Post-processing the flow field<a id='post_filters'></a>
To visualise the flowfield first we convert the **.fld** file (generated by the solver) to a **.vtu** format file by using the `FieldConvert` ulitity available in the Nektar++. To visualise the flowfield first we convert the **.fld** file (generated by the solver) to a **.vtu** format file by using the `FieldConvert` ulitity available in the Nektar++.
<span style="color:Chocolate"><b> Task 10: </b> </span> `Run` the following code cell to convet the **.fld** into a **.vtu** file. <span style="color:Chocolate"><b> Task 10: </b> </span> `Run` the following code cell to convet the **.fld** into a **.vtu** file.
%% Cell type:code id:9af97b18 tags: %% Cell type:code id:9819a13a tags:
``` python ``` python
!FieldConvert -f cyl2DMesh.xml session2d.xml cyl2DMesh.fld flowfield.vtu !FieldConvert -f cyl2DMesh.xml session2d.xml cyl2DMesh.fld flowfield.vtu
``` ```
%% Cell type:markdown id:2540fadd tags: %% Cell type:markdown id:99fd2bf5 tags:
<div class="alert alert-info"> <div class="alert alert-info">
**Note:** `-f` is supplied to force output, otherwise the user would be prompted whether to overwrite an existing file. **Note:** `-f` is supplied to force output, otherwise the user would be prompted whether to overwrite an existing file.
</div> </div>
<div class="alert alert-warning"> <div class="alert alert-warning">
**Note:** It is strongly recommended to use `-f` option when one is running the commands in notebook cells. **Note:** It is strongly recommended to use `-f` option when one is running the commands in notebook cells.
</div> </div>
%% Cell type:markdown id:8525713d tags: %% Cell type:markdown id:a1f6966a tags:
<span style="color:Chocolate"><b> Task 11: </b> </span> `Run` the following `python` script which uses `pyvista` combined with `itkwidgets` to visualise the flowfield in-browser. <span style="color:Chocolate"><b> Task 11: </b> </span> `Run` the following `python` script which uses `pyvista` combined with `itkwidgets` to visualise the flowfield in-browser.
%% Cell type:code id:cd89c851 tags: %% Cell type:code id:91674843 tags:
``` python ``` python
import pyvista as pv import pyvista as pv
# First read the VTK file # First read the VTK file
mesh = pv.read('flowfield.vtu') mesh = pv.read('flowfield.vtu')
# Now create an itkwidget to visualise this in-browser. # Now create an itkwidget to visualise this in-browser.
pl = pv.PlotterITK() pl = pv.PlotterITK()
pl.add_mesh(mesh, scalars=mesh['u'], smooth_shading=True) pl.add_mesh(mesh, scalars=mesh['u'], smooth_shading=True)
pl.show(True) pl.show(True)
``` ```
%% Cell type:markdown id:27b61c94 tags: %% Cell type:markdown id:1ce8863d tags:
## Post-processing the Filter <a id='post_filters'></a> ## Post-processing the Filter <a id='post_filters'></a>
<span style="color:Chocolate"><b> Task 12: </b> </span> `Run` the following `python` script to plot `drag force vs time`, <span style="color:Chocolate"><b> Task 12: </b> </span> `Run` the following `python` script to plot `drag force vs time`,
%% Cell type:code id:78ae382b tags: %% Cell type:code id:4066e97a tags:
``` python ``` python
%matplotlib inline %matplotlib inline
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
# Read the aerodynamic force file (DragLift.fce). We set skiprows=5, since first 5 lines in the data file contain notes. # Read the aerodynamic force file (DragLift.fce). We set skiprows=5, since first 5 lines in the data file contain notes.
fdata = np.loadtxt('DragLift.fce', skiprows=5, usecols=(0,3)) fdata = np.loadtxt('DragLift.fce', skiprows=5, usecols=(0,3))
# Plot time (in the first column) and force (in the second column) # Plot time (in the first column) and force (in the second column)
plt.plot(fdata[:,0], fdata[:,1]) plt.plot(fdata[:,0], fdata[:,1])
plt.xlabel("Time") plt.xlabel("Time")
plt.ylabel("Force") plt.ylabel("Force")
plt.show() plt.show()
``` ```
%% Cell type:markdown id:529d617e tags: %% Cell type:markdown id:8a7bf34b tags:
## **Calculating the pressure coefficient** ## **Calculating the pressure coefficient**
To calculate the pressure coefficient $C_p$ along the cylinder wall first we recall that we had defined the cylinder wall in the `<BOUNDARYREGIONS>` as, To calculate the pressure coefficient $C_p$ along the cylinder wall first we recall that we had defined the cylinder wall in the `<BOUNDARYREGIONS>` as,
```xml ```xml
<BOUNDARYREGIONS> <BOUNDARYREGIONS>
<B ID="0"> C[5-8] </B> <!-- wall --> <B ID="0"> C[5-8] </B> <!-- wall -->
... ...
</BOUNDARYREGIONS> </BOUNDARYREGIONS>
``` ```
Next, we follow the following steps: Next, we follow the following steps:
**Step 1:** Extract the wall boundary surface in a xml file using the `NekMesh` utility as **Step 1:** Extract the wall boundary surface in a xml file using the `NekMesh` utility as
%% Cell type:code id:8f2716fe tags: %% Cell type:code id:dfcf299d tags:
``` python ``` python
!NekMesh -v -m extract:surf=5,6,7,8 cyl2DMesh.xml WallBndSurf.xml !NekMesh -v -m extract:surf=5,6,7,8 cyl2DMesh.xml WallBndSurf.xml
``` ```
%% Cell type:markdown id:c81b935f tags: %% Cell type:markdown id:e2e89f0f tags:
Here the option `surf` uses the wall surface composite number defined as 5,6,7,8 in the `<BOUNDARYREGIONS>`. Here the option `surf` uses the wall surface composite number defined as 5,6,7,8 in the `<BOUNDARYREGIONS>`.
**Step 2:** Extract the field data over the cylinder wall using the `FieldConvert` utility as **Step 2:** Extract the field data over the cylinder wall using the `FieldConvert` utility as
%% Cell type:code id:72fbee70 tags: %% Cell type:code id:7676adab tags:
``` python ``` python
!FieldConvert -f -m extract:bnd=0 cyl2DMesh.xml session2d.xml cyl2DMesh.fld outFieldData.fld !FieldConvert -f -m extract:bnd=0 cyl2DMesh.xml session2d.xml cyl2DMesh.fld outFieldData.fld
``` ```
%% Cell type:markdown id:afb2fea1 tags: %% Cell type:markdown id:cc9c5d5a tags:
which will produce a `outFieldData_b0.fld` file containing the field data over the cylinder wall. The option `bnd` specifies the wall boundary ID defined as 0 in the `<BOUNDARYREGIONS>`. which will produce a `outFieldData_b0.fld` file containing the field data over the cylinder wall. The option `bnd` specifies the wall boundary ID defined as 0 in the `<BOUNDARYREGIONS>`.
**Step 3** To process extracted data file into `vtu` format one can use **Step 3** To process extracted data file into `vtu` format one can use
%% Cell type:code id:620a4dab tags: %% Cell type:code id:f9c725e3 tags:
``` python ``` python
!FieldConvert -f WallBndSurf.xml outFieldData_b0.fld outFieldData_b0.vtu !FieldConvert -f WallBndSurf.xml outFieldData_b0.fld outFieldData_b0.vtu
``` ```
%% Cell type:markdown id:30f06489 tags: %% Cell type:markdown id:ed2699ed tags:
Further processing is required in visualising software to calculate the pressure coefficient. On the other hand, one can use the following simple procedure and utilise a simple Python script as shown below. Further processing is required in visualising software to calculate the pressure coefficient. On the other hand, one can use the following simple procedure and utilise a simple Python script as shown below.
**Step 1**: Prepare a point file either as `pts` in XML format or `csv` as comma delimited containing the $(x,y)$ co-ordinates of the cylinder wall. Here we have supplied the [points_cyl.pts](points_cyl.pts) file with the tutorial. As described in [Section 3](#pdes), the cylinder has a diameter $D=1$ (radius $R=0.5$), so the co-ordinates on the cylinder wall can easily be constructed using $x=R\cos\theta$ and $y=R\sin\theta$ with $0\leq\theta <2\pi$. **Step 1**: Prepare a point file either as `pts` in XML format or `csv` as comma delimited containing the $(x,y)$ co-ordinates of the cylinder wall. Here we have supplied the [points_cyl.pts](points_cyl.pts) file with the tutorial. As described in [Section 3](#pdes), the cylinder has a diameter $D=1$ (radius $R=0.5$), so the co-ordinates on the cylinder wall can easily be constructed using $x=R\cos\theta$ and $y=R\sin\theta$ with $0\leq\theta <2\pi$.
**Step 2**: Interpolate the field data at the points specified in [points_cyl.pts](points_cyl.pts) using the `interppoints` module defied in `FieldConvert` utility. For this, we have to use geometry (cyl2DMesh.xml) and conditions (session2d.xml) in a single file [SessionAndGeo.xml](SessionAndGeo.xml) **Step 2**: Interpolate the field data at the points specified in [points_cyl.pts](points_cyl.pts) using the `interppoints` module defied in `FieldConvert` utility. For this, we have to use geometry (cyl2DMesh.xml) and conditions (session2d.xml) in a single file [SessionAndGeo.xml](SessionAndGeo.xml)
%% Cell type:code id:4b351df8 tags: %% Cell type:code id:25b98464 tags:
``` python ``` python
!FieldConvert -f -m interppoints:fromxml=SessionAndGeo.xml:fromfld=cyl2DMesh.fld:topts=points_cyl.pts cylWall.dat !FieldConvert -f -m interppoints:fromxml=SessionAndGeo.xml:fromfld=cyl2DMesh.fld:topts=points_cyl.pts cylWall.dat
``` ```
%% Cell type:markdown id:98d97059 tags: %% Cell type:markdown id:8584dc89 tags:
The output file `cylWall.dat` contains the variables $x,y,u,v,p$ in each column. The following Python script can be used to plot the pressure over the boundary and calculate the pressure coefficient $C_p$. The output file `cylWall.dat` contains the variables $x,y,u,v,p$ in each column. The following Python script can be used to plot the pressure over the boundary and calculate the pressure coefficient $C_p$.
%% Cell type:code id:92a1def3 tags: %% Cell type:code id:1d1c1986 tags:
``` python ``` python
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import numpy as np import numpy as np
# Open up the output file: skip over the first 3 rows. # Open up the output file: skip over the first 3 rows.
fdata = np.loadtxt('cylWall.dat', skiprows=3) fdata = np.loadtxt('cylWall.dat', skiprows=3)
x = fdata[:,0] # first column x = fdata[:,0] # first column
y = fdata[:,1] # second column y = fdata[:,1] # second column
p = fdata[:,4] # pressure in the final column p = fdata[:,4] # pressure in the final column
# plot pressure along the cylinder wall # plot pressure along the cylinder wall
plt.plot(x, p, 'r') plt.plot(x, p, 'r')
plt.xlabel("$x$") plt.xlabel("$x$")
plt.ylabel("$p$") plt.ylabel("$p$")
plt.show() plt.show()
# calculate Cp using pressure data # calculate Cp using pressure data
U_inf = 1 # reference velocity U_inf = 1 # reference velocity
rho_inf = 1 # reference density rho_inf = 1 # reference density
q2 = 0.5 * rho_inf * U_inf * U_inf q2 = 0.5 * rho_inf * U_inf * U_inf
Cp = p / q2 Cp = p / q2
# plot Cp # plot Cp
plt.plot(x, Cp, 'b') plt.plot(x, Cp, 'b')
plt.xlabel("$x$") plt.xlabel("$x$")
plt.ylabel("$C_p$") plt.ylabel("$C_p$")
plt.show() plt.show()
``` ```
%% Cell type:markdown id:549a174f tags: %% Cell type:markdown id:22239b53 tags:
## Animate the flow field ## Animate the flow field
One can animate the flow field from the `.chk` files generated by the solver. `chk` files capture the flow field at prescribed time interval defined by IO_CheckSteps. In this example we have defined `IO_CheckSteps` as `NumSteps/20`, which generates 20+1=21 checkpoint files (including the initial conditions). These `chk` files assume the names as `cyl2DMesh_0.chk`, `cyl2DMesh_1.chk`, ..., `cyl2DMesh_20.chk`, as per the name of the mesh XML file. Note that `cyl2DMesh_0.chk` contains the flow field prescribed by the initial conditions. We can convert these `chk` files to `vtu` format using `FieldConvert` utility as follows. This will generate output files `out_t_0.vtu` , `out_t_1.vtu` ... One can animate the flow field from the `.chk` files generated by the solver. `chk` files capture the flow field at prescribed time interval defined by IO_CheckSteps. In this example we have defined `IO_CheckSteps` as `NumSteps/20`, which generates 20+1=21 checkpoint files (including the initial conditions). These `chk` files assume the names as `cyl2DMesh_0.chk`, `cyl2DMesh_1.chk`, ..., `cyl2DMesh_20.chk`, as per the name of the mesh XML file. Note that `cyl2DMesh_0.chk` contains the flow field prescribed by the initial conditions. We can convert these `chk` files to `vtu` format using `FieldConvert` utility as follows. This will generate output files `out_t_0.vtu` , `out_t_1.vtu` ...
%% Cell type:code id:6db8dd37 tags: %% Cell type:code id:0564990e tags:
``` python ``` python
!for i in {0..20}; do FieldConvert -f cyl2DMesh.xml session2d.xml cyl2DMesh_${i}.chk out_t_${i}.vtu; done !for i in {0..20}; do FieldConvert -f cyl2DMesh.xml session2d.xml cyl2DMesh_${i}.chk out_t_${i}.vtu; done
``` ```
%% Output %% Output
/bin/bash: FieldConvert: command not found /bin/bash: FieldConvert: command not found
/bin/bash: FieldConvert: command not found /bin/bash: FieldConvert: command not found
/bin/bash: FieldConvert: command not found /bin/bash: FieldConvert: command not found
/bin/bash: FieldConvert: command not found /bin/bash: FieldConvert: command not found
/bin/bash: FieldConvert: command not found /bin/bash: FieldConvert: command not found
/bin/bash: FieldConvert: command not found /bin/bash: FieldConvert: command not found
/bin/bash: FieldConvert: command not found /bin/bash: FieldConvert: command not found
/bin/bash: FieldConvert: command not found /bin/bash: FieldConvert: command not found
/bin/bash: FieldConvert: command not found /bin/bash: FieldConvert: command not found
/bin/bash: FieldConvert: command not found /bin/bash: FieldConvert: command not found
/bin/bash: FieldConvert: command not found /bin/bash: FieldConvert: command not found
/bin/bash: FieldConvert: command not found /bin/bash: FieldConvert: command not found
/bin/bash: FieldConvert: command not found /bin/bash: FieldConvert: command not found
/bin/bash: FieldConvert: command not found /bin/bash: FieldConvert: command not found
/bin/bash: FieldConvert: command not found /bin/bash: FieldConvert: command not found
/bin/bash: FieldConvert: command not found /bin/bash: FieldConvert: command not found
/bin/bash: FieldConvert: command not found /bin/bash: FieldConvert: command not found
/bin/bash: FieldConvert: command not found /bin/bash: FieldConvert: command not found
/bin/bash: FieldConvert: command not found /bin/bash: FieldConvert: command not found
/bin/bash: FieldConvert: command not found /bin/bash: FieldConvert: command not found
/bin/bash: FieldConvert: command not found /bin/bash: FieldConvert: command not found
%% Cell type:markdown id:c5ed52c1 tags: %% Cell type:markdown id:fc3f7d45 tags:
Now that we have obtained the converted vtu files, one can visulaize the evolution of `u`-velocity field with respect to time using the following script. Now that we have obtained the converted vtu files, one can visulaize the evolution of `u`-velocity field with respect to time using the following script.
<div class="alert alert-warning"> <div class="alert alert-warning">
**Note:** The cell below will continue to execute indefinitely until you hit the 'stop' button at the top of the JupyterLab interface. **Note:** The cell below will continue to execute indefinitely until you hit the 'stop' button at the top of the JupyterLab interface.
</div> </div>
%% Cell type:code id:6b1ca730 tags: %% Cell type:code id:e64f3d43 tags:
``` python ``` python
import pyvista as pv import pyvista as pv
import time import time
import numpy as np import numpy as np
import ipywidgets as widgets import ipywidgets as widgets
nfiles = 21 # number of vtu files to be used for animation nfiles = 21 # number of vtu files to be used for animation
Fvtu = [ pv.read('out_t_%d.vtu' % i) for i in range(0, nfiles-1) ] Fvtu = [ pv.read('out_t_%d.vtu' % i) for i in range(0, nfiles-1) ]
plotter = pv.PlotterITK(notebook=True) plotter = pv.PlotterITK(notebook=True)
plotter.add_mesh(Fvtu[0], scalars='u') plotter.add_mesh(Fvtu[0], scalars='u')
viewer = plotter.show(False) viewer = plotter.show(False)
i = 1 i = 1
while True: while True:
viewer.geometries = Fvtu[i-1] viewer.geometries = Fvtu[i-1]
i = (i+1) % nfiles i = (i+1) % nfiles
time.sleep(0.25) time.sleep(0.25)
``` ```
%% Cell type:markdown id:a87f3de3 tags: %% Cell type:markdown id:92487626 tags:
## Optional additional tasks <a id='selfpracticeTasks'> </a> ## Optional additional tasks <a id='selfpracticeTasks'> </a>
<span style="color:Chocolate"><b> Exercise 1 </b> </span> <span style="color:Chocolate"><b> Exercise 1 </b> </span>
Run the solver by changing the flow Reynolds number at $Re = 60$. One might expect to observe Von Karman vortex shedding. Run the solver by changing the flow Reynolds number at $Re = 60$. One might expect to observe Von Karman vortex shedding.
1. Please observe the CFL condition while changing Re as the user may need to adjust timestep size. One can get the information about the CFL number by defining the IO_CFLSteps parameter the <PARAMETER> tag, and we suggest that you set `IO_CFLStep` to `10`. This means that the solver will print the maximum (approximate) CFL number every 10 steps. One may note that increase of Reynolds number might require to decrease the timestep. 1. Please observe the CFL condition while changing Re as the user may need to adjust timestep size. One can get the information about the CFL number by defining the IO_CFLSteps parameter the <PARAMETER> tag, and we suggest that you set `IO_CFLStep` to `10`. This means that the solver will print the maximum (approximate) CFL number every 10 steps. One may note that increase of Reynolds number might require to decrease the timestep.
2. Instead to using the initial conditions at $t=0$ defined in `Task 7`, one might start from a *previous* solution and progress in time, or initialise a *previous* solution for a different Re, as 2. Instead to using the initial conditions at $t=0$ defined in `Task 7`, one might start from a *previous* solution and progress in time, or initialise a *previous* solution for a different Re, as
```xml ```xml
<FUNCTION NAME="InitialConditions"> <FUNCTION NAME="InitialConditions">
<F VAR="u,v,p" FILE="cyl2DMesh.fld" /> <F VAR="u,v,p" FILE="cyl2DMesh.fld" />
</FUNCTION> </FUNCTION>
``` ```
As a good practice, one might rename `cyl2DMesh.fld` to `cyl2DMesh.rst` and use `.rst` file in the above initial conditions instead of `.fld`. Since each time we run the solver, the `fld` file stores the solution field of the most recent run. Here, `rst` stands for restart. As a good practice, one might rename `cyl2DMesh.fld` to `cyl2DMesh.rst` and use `.rst` file in the above initial conditions instead of `.fld`. Since each time we run the solver, the `fld` file stores the solution field of the most recent run. Here, `rst` stands for restart.
<span style="color:Chocolate"><b> Exercise 2 </b> </span> <span style="color:Chocolate"><b> Exercise 2 </b> </span>
Instead of using no-slip boundary condition on the cylinder wall, prescribe slip boundary condition on the cylinder wall, and run the solver. Recall that Navier slip condition is given as $$\frac{\partial u}{\partial y}=-\frac{u_s}{b}=-S,$$ where $u_s$ is the slip velocity and $b$ is the slip length. Here the slip condition constitutes a Neumann boundary condition for the $u$-velocity. Assuming $S=0.001$ (as an arbitrary value for this tutorial) and keeping the boundary conditions for $v$ and $p$ as before, one can assign the boundary conditions on the cylinder wall as Instead of using no-slip boundary condition on the cylinder wall, prescribe slip boundary condition on the cylinder wall, and run the solver. Recall that Navier slip condition is given as $$\frac{\partial u}{\partial y}=-\frac{u_s}{b}=-S,$$ where $u_s$ is the slip velocity and $b$ is the slip length. Here the slip condition constitutes a Neumann boundary condition for the $u$-velocity. Assuming $S=0.001$ (as an arbitrary value for this tutorial) and keeping the boundary conditions for $v$ and $p$ as before, one can assign the boundary conditions on the cylinder wall as
```xml ```xml
<REGION REF="0"> <!-- Wall --> <REGION REF="0"> <!-- Wall -->
<N VAR="u" VALUE="-0.001" /> <!-- Note: N for Neumann BC and D for Dirichlet BC --> <N VAR="u" VALUE="-0.001" /> <!-- Note: N for Neumann BC and D for Dirichlet BC -->
<D VAR="v" VALUE="0" /> <D VAR="v" VALUE="0" />
<N VAR="p" VALUE="0" USERDEFINEDTYPE="H" /> <N VAR="p" VALUE="0" USERDEFINEDTYPE="H" />
</REGION> </REGION>
``` ```
<span style="color:Chocolate"><b> Exercise 3 </b> </span> <span style="color:Chocolate"><b> Exercise 3 </b> </span>
In this tutorial we have kept the flow domain artificially small (and coarse grid resoluation) so that solution could be obtained in few seconds without overwhelming the online system. Interested users can use a more appropriate flow domain and grid resolution to obtain the desired accuracy. Here we have supplemented the [Gmsh](https://gmsh.info) file [cylinder_small.geo](cylinder_small.geo). After using the desired geometric quantities, one can convert `.geo` to `.msh` file in Gmsh. Then using `NekMesh`, perform the boundary layer splitting, and obtain a final XML file. Details are available at Section 4.4.7 of the User Guide. In this tutorial we have kept the flow domain artificially small (and coarse grid resoluation) so that solution could be obtained in few seconds without overwhelming the online system. Interested users can use a more appropriate flow domain and grid resolution to obtain the desired accuracy. Here we have supplemented the [Gmsh](https://gmsh.info) file [cylinder_small.geo](cylinder_small.geo). After using the desired geometric quantities, one can convert `.geo` to `.msh` file in Gmsh. Then using `NekMesh`, perform the boundary layer splitting, and obtain a final XML file. Details are available at Section 4.4.7 of the User Guide.
%% Cell type:markdown id:f15de96a tags: %% Cell type:markdown id:aa51c8a1 tags:
# Quasi-3D case # Quasi-3D case
Now that we have completed the two-dimensional simulation, next we discuss a three-dimensional flow problem by extruding the two-dimensional flow domain in the third ($z$) direction. Instead of using finite elements in this direction, we use a Fourier expansion; this imposes periodic boundary conditions in the spanwise direction of the cylinder and thus makes the problem a quasi-3D flow problem. The boundary conditions remain unchanged from the 2D simulation, asides from the inclusion of a third velocity component $w$: Now that we have completed the two-dimensional simulation, next we discuss a three-dimensional flow problem by extruding the two-dimensional flow domain in the third ($z$) direction. Instead of using finite elements in this direction, we use a Fourier expansion; this imposes periodic boundary conditions in the spanwise direction of the cylinder and thus makes the problem a quasi-3D flow problem. The boundary conditions remain unchanged from the 2D simulation, asides from the inclusion of a third velocity component $w$:
- No−slip condition on the cylinder surface (i.e., homogeneous Dirichlet BC on the velocity and high-order Neumann BC on the pressure), - No−slip condition on the cylinder surface (i.e., homogeneous Dirichlet BC on the velocity and high-order Neumann BC on the pressure),
- Far−field with Dirichlet BC: $(u, v, w) = (1, 0, 0)$ at the bottom and top boundaries, - Far−field with Dirichlet BC: $(u, v, w) = (1, 0, 0)$ at the bottom and top boundaries,
- Inflow with Dirichlet BC: $(u, v, w) = (1, 0, 0)$ at the left boundary, - Inflow with Dirichlet BC: $(u, v, w) = (1, 0, 0)$ at the left boundary,
- Outflow (higher order) at the right boundary. - Outflow (higher order) at the right boundary.
Since we are extruding the geometry, we can use the same 2D geometry XML file for a quasi-3D simulation. Since the output files `.chk` and `.fld` generated by the simulation are named as per the name of the geometry XML file, here we just renamed the [cyl2DMesh.xml](cyl2DMesh.xml) to [cyl3DMesh.xml](cyl3DMesh.xml). Since we are extruding the geometry, we can use the same 2D geometry XML file for a quasi-3D simulation. Since the output files `.chk` and `.fld` generated by the simulation are named as per the name of the geometry XML file, here we just renamed the [cyl2DMesh.xml](cyl2DMesh.xml) to [cyl3DMesh.xml](cyl3DMesh.xml).
%% Cell type:markdown id:5b4a42df tags: %% Cell type:markdown id:0244b007 tags:
<span style="color:Chocolate"><b> Task 13: </b> </span> Let us open the condition file [**session3d.xml**](session3d.xml), and complete the `<EXPANSIONS>` tag to run a simulation at order $P=2$. A completed version should have the following form: <span style="color:Chocolate"><b> Task 13: </b> </span> Let us open the condition file [**session3d.xml**](session3d.xml), and complete the `<EXPANSIONS>` tag to run a simulation at order $P=2$. A completed version should have the following form:
```xml ```xml
<EXPANSIONS> <EXPANSIONS>
<E COMPOSITE="C[100,101]" NUMMODES="3" FIELDS="u,v,w,p" TYPE="MODIFIED" /> <E COMPOSITE="C[100,101]" NUMMODES="3" FIELDS="u,v,w,p" TYPE="MODIFIED" />
</EXPANSIONS>` </EXPANSIONS>`
``` ```
%% Cell type:markdown id:1a58e92c tags: %% Cell type:markdown id:d801ab9b tags:
<span style="color:Chocolate"><b> Task 14: </b> </span> In the `<SOLVERINFO>` tag, we use the same solver properties as discussed in the 2D case tutorial. Now, add the following two extra lines for the simulation of a quasi-3D case, <span style="color:Chocolate"><b> Task 14: </b> </span> In the `<SOLVERINFO>` tag, we use the same solver properties as discussed in the 2D case tutorial. Now, add the following two extra lines for the simulation of a quasi-3D case,
```xml ```xml
<I PROPERTY="Homogeneous" VALUE="1D" /> <I PROPERTY="Homogeneous" VALUE="1D" />
<I PROPERTY="UseFFT" VALUE="FFTW" /> <I PROPERTY="UseFFT" VALUE="FFTW" />
``` ```
Here, Here,
- `Homogeneous` dictates how many dimensions of the geomety to be defied as homogeneous. For `1D` homogeneous case the $z$-direction is assumed to be homogeneous. - `Homogeneous` dictates how many dimensions of the geomety to be defied as homogeneous. For `1D` homogeneous case the $z$-direction is assumed to be homogeneous.
- `UseFFT` activates the use of fast Fourier transform. (By default Nektar++ will use a standard discrete Fourier transform, which reduces the requirement for an additional dependency on FFTW, but can be significantly slower.) - `UseFFT` activates the use of fast Fourier transform. (By default Nektar++ will use a standard discrete Fourier transform, which reduces the requirement for an additional dependency on FFTW, but can be significantly slower.)
%% Cell type:markdown id:63064e64 tags: %% Cell type:markdown id:5a534086 tags:
<span style="color:Chocolate"><b> Task 15: </b> </span> In `<PARAMETERS>` tag, add the following extra parameters needed for a quasi-3D case <span style="color:Chocolate"><b> Task 15: </b> </span> In `<PARAMETERS>` tag, add the following extra parameters needed for a quasi-3D case
```xml ```xml
<P> HomModesZ = 4 </P> <P> HomModesZ = 4 </P>
<P> LZ = 2*PI </P> <P> LZ = 2*PI </P>
``` ```
Here, Here,
- `HomModesZ` specifies the number of modes used in the harmonic (Fourier) expansion along homogeneous z-direction, and must be an **even** number. - `HomModesZ` specifies the number of modes used in the harmonic (Fourier) expansion along homogeneous z-direction, and must be an **even** number.
- `LZ` denotes the homogeneous length along the $z$-direction. - `LZ` denotes the homogeneous length along the $z$-direction.
%% Cell type:markdown id:d4e4880a tags: %% Cell type:markdown id:020ddbe1 tags:
<span style="color:Chocolate"><b> Task 16: </b> </span> <span style="color:Chocolate"><b> Task 16: </b> </span>
In `<VARIABLES>` tag, add the field variable $w$ and change the `ID` of the pressure variable $p$ to `3`. The completed section should look like: In `<VARIABLES>` tag, add the field variable $w$ and change the `ID` of the pressure variable $p$ to `3`. The completed section should look like:
```xml ```xml
<VARIABLES> <VARIABLES>
<V ID="0">u</V> <V ID="0">u</V>
<V ID="1">v</V> <V ID="1">v</V>
<V ID="2">w</V> <V ID="2">w</V>
<V ID="3">p</V> <V ID="3">p</V>
</VARIABLES> </VARIABLES>
``` ```
%% Cell type:markdown id:a673800f tags: %% Cell type:markdown id:5a528b62 tags:
<span style="color:Chocolate"><b> Task 17: </b> </span> In `<BOUNDARYCONDITIONS>` tag, add boundary conditions for $w$-velocity for each boundary region as $w$=0. <span style="color:Chocolate"><b> Task 17: </b> </span> In `<BOUNDARYCONDITIONS>` tag, add boundary conditions for $w$-velocity for each boundary region as $w$=0.
%% Cell type:markdown id:a8af74b4 tags: %% Cell type:markdown id:d8b74712 tags:
<span style="color:Chocolate"><b> Task 18: </b> </span> In `InitialConditions` function, provide the initial condition for the $w$ component of velocity as <span style="color:Chocolate"><b> Task 18: </b> </span> In `InitialConditions` function, provide the initial condition for the $w$ component of velocity as
```xml ```xml
<E VAR="w" VALUE="0.02*awgn(1.0)" /> <E VAR="w" VALUE="0.02*awgn(1.0)" />
``` ```
Here we provide the initial condition for $w$-velocity using a Gaussian noise generator defined as `A*awgn(sigma)`, where `A` is the amplitude of the noise and *sigma* is the variance of the normal distribution with zero-mean. Here we provide the initial condition for $w$-velocity using a Gaussian noise generator defined as `A*awgn(sigma)`, where `A` is the amplitude of the noise and *sigma* is the variance of the normal distribution with zero-mean.
For $u$, $v$ and $p$, initial conditions can be given as the simulation results obtained from 2D case as For $u$, $v$ and $p$, initial conditions can be given as the simulation results obtained from 2D case as
```xml ```xml
<F VAR="u,v,p" FILE="cyl2DMesh.fld" /> <F VAR="u,v,p" FILE="cyl2DMesh.fld" />
``` ```
<div class="alert alert-warning"> <div class="alert alert-warning">
**Note:** For this example the `NumSteps=200` is set with `TimeStep=0.025` which results in the total time of 5 units just for demonstrations. Feel free to increas `NumSteps` but it obviously take longer time for the simulaiton to finish. **Note:** For this example the `NumSteps=200` is set with `TimeStep=0.025` which results in the total time of 5 units just for demonstrations. Feel free to increas `NumSteps` but it obviously take longer time for the simulaiton to finish.
</div> </div>
%% Cell type:markdown id:6157e3e1 tags: %% Cell type:markdown id:c3c3cb9f tags:
## Solution stabilisation parameters ## Solution stabilisation parameters
To stabilise the solution, we define the following paramters in `SOLVERINFO`: To stabilise the solution, we define the following paramters in `SOLVERINFO`:
```xml ```xml
<I PROPERTY="SpectralVanishingViscosity" VALUE="DGKernel" /> <I PROPERTY="SpectralVanishingViscosity" VALUE="DGKernel" />
<I PROPERTY="SpectralVanishingViscosityHomo1D" VALUE="ExpKernel" /> <I PROPERTY="SpectralVanishingViscosityHomo1D" VALUE="ExpKernel" />
<I PROPERTY="SPECTRALHPDEALIASING" VALUE="True" /> <I PROPERTY="SPECTRALHPDEALIASING" VALUE="True" />
<I PROPERTY="DEALIASING" VALUE="True" /> <I PROPERTY="DEALIASING" VALUE="True" />
``` ```
In a quasi-3D simulation, `SpectralVanishingViscosity` affects both the Fourier and the spectral/$hp$ expansions. We also use a seperate parameter `SpectralVanishingViscosityHomo1D` for the Fourier or $z$-directionm with a corresponding parameter value in `PARAMETERS`. For dealiasing, we use both spectral/$hp$ dealiasing and Fourier dealiasing. In a quasi-3D simulation, `SpectralVanishingViscosity` affects both the Fourier and the spectral/$hp$ expansions. We also use a seperate parameter `SpectralVanishingViscosityHomo1D` for the Fourier or $z$-directionm with a corresponding parameter value in `PARAMETERS`. For dealiasing, we use both spectral/$hp$ dealiasing and Fourier dealiasing.
%% Cell type:markdown id:172d090c tags: %% Cell type:markdown id:eccdfe19 tags:
## ModalEnergy filter setup ## ModalEnergy filter setup
In order to calculate the time evolution of the kinetic energy of each Fourier modes, we use `ModalEnergy` filter. In order to calculate the time evolution of the kinetic energy of each Fourier modes, we use `ModalEnergy` filter.
<span style="color:Chocolate"><b> Task 19: </b> </span> <span style="color:Chocolate"><b> Task 19: </b> </span>
Provide the name of the `OutputFile` as `EnergyFile`, `OutputFrequency` = 10 as shown below Provide the name of the `OutputFile` as `EnergyFile`, `OutputFrequency` = 10 as shown below
```xml ```xml
<FILTER TYPE="ModalEnergy"> <FILTER TYPE="ModalEnergy">
<PARAM NAME="OutputFile">EnergyFile</PARAM> <PARAM NAME="OutputFile">EnergyFile</PARAM>
<PARAM NAME="OutputFrequency">10</PARAM> <PARAM NAME="OutputFrequency">10</PARAM>
</FILTER> </FILTER>
``` ```
<div class="alert alert-info"> <div class="alert alert-info">
**Note:** In the case of quasi-3D simulation with one homogeneous direction, the modal kinetic energy is defined as $ E _{m} (t) = \frac{1}{2} \int_{\Omega} ||\mathbf{\hat {V}} _{m} ||^{2}~\mathrm{d} \mathbf{x}$, where $ \mathbf{\hat {V}} _{m}$ is a complex Fourier co-efficient obtained from the velocity field defined as $\mathbf {V} (\mathbf{x},t)=\sum _{m=0} ^{N}\mathbf{\hat {V}} _{m} (t) e ^{i 2 \pi m \mathbf{x}} $. With $m=0,1,2,... **Note:** In the case of quasi-3D simulation with one homogeneous direction, the modal kinetic energy is defined as $ E _{m} (t) = \frac{1}{2} \int_{\Omega} ||\mathbf{\hat {V}} _{m} ||^{2}~\mathrm{d} \mathbf{x}$, where $ \mathbf{\hat {V}} _{m}$ is a complex Fourier co-efficient obtained from the velocity field defined as $\mathbf {V} (\mathbf{x},t)=\sum _{m=0} ^{N}\mathbf{\hat {V}} _{m} (t) e ^{i 2 \pi m \mathbf{x}} $. With $m=0,1,2,...
</div> </div>
%% Cell type:markdown id:3cca6892 tags: %% Cell type:markdown id:e4e9b3ef tags:
## Running the solver in serial ## Running the solver in serial
The two XML files that we use to run the solver are: The two XML files that we use to run the solver are:
- Geometry file: [cyl3DMesh.xml](cyl3DMesh.xml) - Geometry file: [cyl3DMesh.xml](cyl3DMesh.xml)
- Condition file: [session3d.xml](session3d.xml) - Condition file: [session3d.xml](session3d.xml)
<span style="color:Chocolate"><b> Task 20: </b> </span> `Run` the incompressible flow solver in Serial using the following code cell <span style="color:Chocolate"><b> Task 20: </b> </span> `Run` the incompressible flow solver in Serial using the following code cell
<div class="alert alert-success"> <div class="alert alert-success">
**Note:** You should see that simulation produces output files as `cyl3DMesh_N.chk` where $N=1,2,3,\dots$ and a final field file `cyl3DMesh.fld`. **Note:** You should see that simulation produces output files as `cyl3DMesh_N.chk` where $N=1,2,3,\dots$ and a final field file `cyl3DMesh.fld`.
</div> </div>
%% Cell type:code id:ee7855ac tags: %% Cell type:code id:f4801fc8 tags:
``` python ``` python
!IncNavierStokesSolver cyl3DMesh.xml session3d.xml !IncNavierStokesSolver cyl3DMesh.xml session3d.xml
``` ```
%% Cell type:markdown id:1bb2ac3d tags: %% Cell type:markdown id:edb8c19f tags:
## Post-processing the flow field<a id='post_filters'></a> ## Post-processing the flow field<a id='post_filters'></a>
<span style="color:Chocolate"><b> Task 21: </b> </span> `Run` the following code cell to convert the `.fld` into a `.vtu` file for visualisation, and then execute the subsequent cell to visualise the output. <span style="color:Chocolate"><b> Task 21: </b> </span> `Run` the following code cell to convert the `.fld` into a `.vtu` file for visualisation, and then execute the subsequent cell to visualise the output.
%% Cell type:code id:2b78faf2 tags: %% Cell type:code id:09deeff9 tags:
``` python ``` python
!FieldConvert -f -m vorticity cyl3DMesh.xml session3d.xml cyl3DMesh.fld flowfield3DH.vtu !FieldConvert -f -m vorticity cyl3DMesh.xml session3d.xml cyl3DMesh.fld flowfield3DH.vtu
``` ```
%% Cell type:code id:a2a4a6c1 tags: %% Cell type:code id:1d58152d tags:
``` python ``` python
import pyvista as pv import pyvista as pv
# First read the VTK file # First read the VTK file
mesh = pv.read('flowfield3DH.vtu') mesh = pv.read('flowfield3DH.vtu')
# Now create an itkwidget to visualise this in-browser. # Now create an itkwidget to visualise this in-browser.
pl = pv.PlotterITK() pl = pv.PlotterITK()
pl.add_mesh(mesh, scalars=mesh['u'], smooth_shading=True) pl.add_mesh(mesh, scalars=mesh['u'], smooth_shading=True)
pl.show(True) pl.show(True)
``` ```
%% Cell type:markdown id:86c70f79 tags: %% Cell type:markdown id:51c53571 tags:
### Q-criterion ### Q-criterion
The $Q$-criterion is a common post-processing technique to reveal flow features such as vortices. We can calculate the $Q$-criterion for a field output from the Navier-Stokes solver through the use of the `QCriterion` module in FieldConvert as follows: The $Q$-criterion is a common post-processing technique to reveal flow features such as vortices. We can calculate the $Q$-criterion for a field output from the Navier-Stokes solver through the use of the `QCriterion` module in FieldConvert as follows:
%% Cell type:code id:6c10e1ec tags: %% Cell type:code id:c5c3d530 tags:
``` python ``` python
!FieldConvert -f -m QCriterion cyl3DMesh.xml session3d.xml cyl3DMesh.fld cyl3DMesh-Qcrit.fld !FieldConvert -f -m QCriterion cyl3DMesh.xml session3d.xml cyl3DMesh.fld cyl3DMesh-Qcrit.fld
``` ```
%% Cell type:markdown id:f6d7c7a4 tags: %% Cell type:markdown id:480afc0f tags:
The generated`cyl3DMesh-Qcrit.fld` file can easily be visualised using `Task 20`. The generated`cyl3DMesh-Qcrit.fld` file can easily be visualised using `Task 20`.
%% Cell type:markdown id:f964b265 tags: %% Cell type:markdown id:03c595df tags:
### Extract the mean-mode in Fourier space ### Extract the mean-mode in Fourier space
A 2D expansion containing the mean mode (plane zero in Fourier space) of a quasi-3D field file can be obtained as A 2D expansion containing the mean mode (plane zero in Fourier space) of a quasi-3D field file can be obtained as
%% Cell type:code id:d7770bbc tags: %% Cell type:code id:ef5fe801 tags:
``` python ``` python
!FieldConvert -f -m meanmode cyl3DMesh.xml session3d.xml cyl3DMesh.fld cyl3DMesh-meanmode.fld !FieldConvert -f -m meanmode cyl3DMesh.xml session3d.xml cyl3DMesh.fld cyl3DMesh-meanmode.fld
``` ```
%% Cell type:markdown id:c531e19f tags: %% Cell type:markdown id:3c37d264 tags:
The generated`cyl3DMesh-meanmode.fld` file can easily be visualised using `Task 20`. The generated`cyl3DMesh-meanmode.fld` file can easily be visualised using `Task 20`.
%% Cell type:markdown id:c131c960 tags: %% Cell type:markdown id:cf1dc168 tags:
## Post-processing ModalEnergy filter ## Post-processing ModalEnergy filter
Finally, we can plot the kinetic energy of each Fourier mode as a function of time using the file `EnergyFile.mdl`, which is generated from the modal energy filter. The following cell will visualise the output on a semi-log scale. Finally, we can plot the kinetic energy of each Fourier mode as a function of time using the file `EnergyFile.mdl`, which is generated from the modal energy filter. The following cell will visualise the output on a semi-log scale.
%% Cell type:code id:1d04d276 tags: %% Cell type:code id:a240ff94 tags:
``` python ``` python
import numpy as np import numpy as np
import matplotlib import matplotlib
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
#------- Please input --------------- #------- Please input ---------------
HomModesZ = 4 HomModesZ = 4
#------------------------------------- #-------------------------------------
n_modes=int(HomModesZ/2) n_modes=int(HomModesZ/2)
print(n_modes) print(n_modes)
#------------------------------------- #-------------------------------------
# Read the Energy file data (EnergyFile.mdl) # Read the Energy file data (EnergyFile.mdl)
fileName = 'EnergyFile.mdl' fileName = 'EnergyFile.mdl'
fdata = np.loadtxt(fileName,skiprows=1) fdata = np.loadtxt(fileName,skiprows=1)
# skiprows=1, since first 1 line in the data file contain notes. # skiprows=1, since first 1 line in the data file contain notes.
# the data file contains 3 columns: Time, Mode, Energy # the data file contains 3 columns: Time, Mode, Energy
time=fdata[:,0] time=fdata[:,0]
Modes=fdata[:,1] Modes=fdata[:,1]
E=fdata[:,2] E=fdata[:,2]
nRow=time.size nRow=time.size
counter= np.arange(0,nRow,n_modes) counter= np.arange(0,nRow,n_modes)
time=time[counter] time=time[counter]
mode0e=E[counter] mode0e=E[counter]
mode1e=E[counter+1] mode1e=E[counter+1]
plt.plot(time,mode0e,label="mode 0") plt.plot(time,mode0e,label="mode 0")
plt.plot(time,mode1e,label="mode 1") plt.plot(time,mode1e,label="mode 1")
plt.yscale("log") plt.yscale("log")
plt.xlabel("Time") plt.xlabel("Time")
plt.ylabel("E") plt.ylabel("E")
plt.legend() plt.legend()
plt.show() plt.show()
``` ```
%% Cell type:markdown id:f059909f tags: %% Cell type:markdown id:bbf9c019 tags:
## Parallel run ## Parallel run
Nektar++ solvers can also be run in parallel. For a quasi-3D case, the parallel running of the solver is slightly tricky, because parallelisation can be applied in the spectral/$hp$ element plane, in the Fourier direction, or in a combination of both. The general syntax for parallel running of a quasi-3D case is Nektar++ solvers can also be run in parallel. For a quasi-3D case, the parallel running of the solver is slightly tricky, because parallelisation can be applied in the spectral/$hp$ element plane, in the Fourier direction, or in a combination of both. The general syntax for parallel running of a quasi-3D case is
`mpirun -np N IncNavierStokesSolver --npz NZ cyl3DMesh.xml session3d.xml` `mpirun -np N IncNavierStokesSolver --npz NZ cyl3DMesh.xml session3d.xml`
In the above command: In the above command:
- `N` is the total number of cores used for the simulation; - `N` is the total number of cores used for the simulation;
- `NZ` is the number of cores (out of the total `N` cores) that we want to assign for the Fourier modes used in the $z$-direction. - `NZ` is the number of cores (out of the total `N` cores) that we want to assign for the Fourier modes used in the $z$-direction.
For example if we have a mesh with $256$ elements and $32$ Fourier planes, and we run on `N = 16` processors then: For example if we have a mesh with $256$ elements and $32$ Fourier planes, and we run on `N = 16` processors then:
- `NZ = 1` means that parallelisation occurs only in the spectral element plane: every processor will be given a partition that has $16$ elements and $32$ planes; - `NZ = 1` means that parallelisation occurs only in the spectral element plane: every processor will be given a partition that has $16$ elements and $32$ planes;
- `NZ = 16` means that parallelisation occurs only in the Fourier directions: every processor will be given a full mesh of $256$ elements, but only $2$ planes (the minimum number); - `NZ = 16` means that parallelisation occurs only in the Fourier directions: every processor will be given a full mesh of $256$ elements, but only $2$ planes (the minimum number);
- `NZ = 2` means that parallelisation occurs in a mixture of directions; each process will get $32 / 2 = 16$ planes and $256 / 2 = 128$ elements. - `NZ = 2` means that parallelisation occurs in a mixture of directions; each process will get $32 / 2 = 16$ planes and $256 / 2 = 128$ elements.
Before running the solver in parallel we need to change the property `GlobalSysSoln` (defined in `SOLVERINFO` tag) from `DirectStaticCond` to `XxtMultiLevelStaticCond`, since direct solvers only work if the entire spectral/$hp$ element plane is contained on a single processor. Before running the solver in parallel we need to change the property `GlobalSysSoln` (defined in `SOLVERINFO` tag) from `DirectStaticCond` to `XxtMultiLevelStaticCond`, since direct solvers only work if the entire spectral/$hp$ element plane is contained on a single processor.
```xml ```xml
<I PROPERTY="GlobalSysSoln" VALUE="XxtMultiLevelStaticCond" /> <I PROPERTY="GlobalSysSoln" VALUE="XxtMultiLevelStaticCond" />
``` ```
`GlobalSysSoln` sets the approach we use to solve the global linear system $Ax=b$. `DirectStaticCond` uses direct static condensation technique, and can only be used for serial runs. `Xxt` solvers are still direct, but can be run in parallel (scaling to smaller core counts). `GlobalSysSoln` sets the approach we use to solve the global linear system $Ax=b$. `DirectStaticCond` uses direct static condensation technique, and can only be used for serial runs. `Xxt` solvers are still direct, but can be run in parallel (scaling to smaller core counts).
<div class="alert alert-warning"> <div class="alert alert-warning">
**Note:** Use of high number of processors might not work in the current online tutorial. Also, parallel run might be slower than the serial run due to resource constraints. **Note:** Use of high number of processors might not work in the current online tutorial. Also, parallel run might be slower than the serial run due to resource constraints.
</div> </div>
%% Cell type:code id:870b3f54 tags: %% Cell type:code id:31e67ace tags:
``` python ``` python
!mpirun -np 4 IncNavierStokesSolver --npz 2 cyl3DMesh.xml session3d.xml !mpirun -np 4 IncNavierStokesSolver --npz 2 cyl3DMesh.xml session3d.xml
``` ```
%% Cell type:markdown id:02038429 tags: %% Cell type:markdown id:561d10c7 tags:
## Extra tasks ## Extra tasks
<span style="color:Chocolate"><b> Exercise 4 </b> </span> <span style="color:Chocolate"><b> Exercise 4 </b> </span>
Observe the effect of solution stabilising (Spectral Vanishing Vicosity) parameters by changing SVVCutoffRatio and SVVDiffCoeff defined in `<PARAMETERS>`. Observe the effect of solution stabilising (Spectral Vanishing Vicosity) parameters by changing SVVCutoffRatio and SVVDiffCoeff defined in `<PARAMETERS>`.
<span style="color:Chocolate"><b> Exercise 5 </b> </span> <span style="color:Chocolate"><b> Exercise 5 </b> </span>
In this tutorial we have kept the flow domain artificially small (and the grid resolution very coarse) so that solutions can be obtained in few seconds without overwhelming the online system. The interested user can use appropriate flow domain and grid resolution to obtain the desired accuracy (see exercise 2 of this tutorial). In this tutorial we have kept the flow domain artificially small (and the grid resolution very coarse) so that solutions can be obtained in few seconds without overwhelming the online system. The interested user can use appropriate flow domain and grid resolution to obtain the desired accuracy (see exercise 2 of this tutorial).
%% Cell type:markdown id:928dcc78 tags: %% Cell type:markdown id:600d144e tags:
# Take-a-ways <a id='takeaways'> </a> # Take-a-ways <a id='takeaways'> </a>
From this tutorial, the following points are worth considering for your future simulations: From this tutorial, the following points are worth considering for your future simulations:
- Time step-size should be selected carefully to avoid issues associated with the CFL condition. - Time step-size should be selected carefully to avoid issues associated with the CFL condition.
- Reynolds number is not directly needed for the incompressible solver. It is the kinematic viscosity `Kinvis` that is read by the solver. - Reynolds number is not directly needed for the incompressible solver. It is the kinematic viscosity `Kinvis` that is read by the solver.
- Spectral-dealiasing should not be confused with Fourier-dealiasing, Quasi-3D (homogeneous) simulation uses Fourier-dealiasing together with regular spectral-dealiasing. - Spectral-dealiasing should not be confused with Fourier-dealiasing, Quasi-3D (homogeneous) simulation uses Fourier-dealiasing together with regular spectral-dealiasing.
- Number of Fourier modes (denoted by `HomModesZ`) used for z-direction must be an even number. - Number of Fourier modes (denoted by `HomModesZ`) used for z-direction must be an even number.
%% Cell type:markdown id:5dff8d55 tags: %% Cell type:markdown id:0f76e9c4 tags:
# References <a id='references'> </a> # References <a id='references'> </a>
1. Nektar++: Spectral/hp Element Framework, Version 5.0.1, [User guide](http://www.nektar.info/downloads), January, 2021. 1. Nektar++: Spectral/hp Element Framework, Version 5.0.1, [User guide](http://www.nektar.info/downloads), January, 2021.
2. Kirby R, Sherwin S. Stabilisation of spectral/hp methods through spectral vanishing viscosity: Application to fluid mechanics modelling. Comput. Methods Appl. Mech. Engrg. 2006; 195: 3128–3144 2. Kirby R, Sherwin S. Stabilisation of spectral/hp methods through spectral vanishing viscosity: Application to fluid mechanics modelling. Comput. Methods Appl. Mech. Engrg. 2006; 195: 3128–3144
3. G. Mengaldo, D. De Grazia, D. Moxey, P.E. Vincent, S. Sherwin, Dealiasing techniques for high-order spectral element methods on regular and irregular grids, J. Comput. Phys. 299 (2015) 56–81. 3. G. Mengaldo, D. De Grazia, D. Moxey, P.E. Vincent, S. Sherwin, Dealiasing techniques for high-order spectral element methods on regular and irregular grids, J. Comput. Phys. 299 (2015) 56–81.
4. A. Bolis, C.D. Cantwell, D. Moxey, D. Serson, S.J. Sherwin, An adaptable parallel algorithm for the direct numerical simulation of incompressible turbulent flows using a Fourier spectral/hp element method and MPI virtual topologies, Computer Physics Communications, 206 (2016) 17–25. 4. A. Bolis, C.D. Cantwell, D. Moxey, D. Serson, S.J. Sherwin, An adaptable parallel algorithm for the direct numerical simulation of incompressible turbulent flows using a Fourier spectral/hp element method and MPI virtual topologies, Computer Physics Communications, 206 (2016) 17–25.
%% Cell type:markdown id:c3b6f258 tags: %% Cell type:markdown id:3027a7c1 tags:
## End of the tutorial. ## End of the tutorial.
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment