# Tank filling 

This tutorial aims at running a simulation of the tank filling test case, see {numref}`fig:tankfilling`.
The test case deals with a 2D flow with Navier-Stokes and the equation for one constituent.

```{figure} /_static/FIGURES/tank2D.png
:class: custom-image-class
:name: fig:tankfilling
:alt: tank2D

Geometry of the 2D tank case
```

The following table summarizes the parameters of the simulation:

| **Fluid Properties** | **Value** |
|:---------------:|:---------------:|
| Dynamic viscosity ($\mu$) | $ 10^{-3}$ $kg \cdot m^{-1} \cdot s^{-1}$ |
| Density ($\rho$) | $1000$ $kg \cdot m^{-3}$ |
| Diffusion ($D$) | $10^{-9}$ $m^{2}\cdot s^{-1}$| 

| **Boundary Conditions** | **Value** |
|:---------------:|:---------------:|
| Inlet velocity ($V(t)$) | $\begin{cases} 1 -(y-0.025/0.005)^{2} & ,\; t \leq 0.5s \\ 0 & ,\; t>0.5s \end{cases}$ | 
| Concentration ($C$) | $ \begin{cases} 1 & ,\; t \leq 0.5s \\ 0 & ,\; t>0.5s \end{cases}$|
| Outlet pressure ($P_0$) | $0$ $Pa$ |
| Wall velocity | V= 0 |

| **Initial conditions** | **Value** |
|:---------------:|:---------------:|
| Velocity ($V(x,y,t=0)$)| $0 m \cdot s^{-1}$|
| Concentration ($C(x,y,t=0)$) | 0 |


You can copy the study named `diagonale`:
```bash
source $my_path_to_TRUST_installation/env_TRUST.sh
mkdir -p TRUST_tutorials
cd TRUST_tutorials
trust -copy diagonale
```

Now, edit the `diagonale.data` file.

Then, modify the fluid characteristics to the one given in the above table ($\mu, \rho, D$).

You then need to modify the geometry parameters, so that your geometry resembles {numref}`fig:tankblock`.

```{figure} /_static/FIGURES/tank2D_2.png
:width: 40%
:name: fig:tankblock
:alt: tank2D_mesh

Mesh blocks of the 2D tank filling case
```

To do so, you have to create three blocks, starting with $dx=dy=0.2cm$ which gives a total nodes number $Nx=51$ and $Ny=121$, see {numref}`fig:tankblock`.

- The first block,`Block1`, whose origin is (0, 0.03), $Nx=51$, $Ny=106$ (for $dx=dy=0.2cm$), $L=0.1 m$, $H=0.21 m$. Name the wall boundaries Left1, Outlet(=Top1) and Right1. (Don't forget the comma between block definitions.)

- The second block, `Block2`, whose origin is (0, 0.02), $Nx=51$, $Ny=6$ (for $dx=dy=0.2cm$), $L=0.1 m$, $H=0.01 m$. Name the wall boundaries Inlet(=Left2) and Right2.

- The third block, `Block3`, whose origin is (0, 0), $Nx=51$, $Ny=11$ (for $dx=dy=0.2cm$), $L=0.1 m$, $H=0.02 m$. Name the wall boundaries Left3, Bottom3 and Right3.

Then, define the boundary wall, using the keyword **{ref}`regroupebord`**.

You could also use **facteurs**, **symx** and **symy** keywords to define a refined mesh near the walls.

Once you are done with the geometry, change the values in the time scheme to stop the calculation at 1 second, and modify **dt\_min** and **dt\_max** values to let **TRUST** compute at least one time step.

Now, set the gravity value to $-9.81 m.s^{-2}$ along the y-axis. 

Note that the **beta\_co** keyword may be useful in order to have a Boussinesq coupling between momentum and concentration equations ($\beta C_0 g(C-C_0$) source term added to the Navier-Stokes equations).

## Boundary and initial conditions

You need to change the initial and boundary conditions for Navier-Stokes equations:

- for the Outlet boundary, impose $P=0$,

- for the Wall boundary, impose $V_x=V_y=0$ with **{ref}`paroi_fixe`** keyword,

- for the Inlet boundary, impose $(V_{x},V_{y})=(V(t),0)$ with:

   $V(t)= \begin{cases}
            1-(y-0.025/0.005)^{2} & ,\; t\leq0.5s\\
            0 & ,\; t>0.5s
            \end{cases}$.

   You will use the **{ref}`champ_front_fonc_txyz`** keyword for the velocity, to write something like: **Champ\_Front\_Fonc\_txyz $2$ $(1-((y-0.025)/0.005)^2)*(t<0.5)$ $0.$**
   ```{note}
   Use ($t[0.5)$ syntax if you prefer ($t<=0.5$)
   ```

For **Convection\_diffusion\_Concentration**, you need to use:

- For the Outlet, use the following keywords to ensure the external concentration is 0 :

   **Frontiere\_ouverte C\_ext Champ\_front\_uniforme 1 0.**

- For the Wall, the keyword for impermeable boundary condition for concentration is **{ref}`paroi`**.

- for the Inlet boundary condition of concentration, **{ref}`champ_front_fonc_txyz`** field

Then make sure to check that you have high-order schemes (i.e. **Quick** scheme) used in both equations to reduce numerical diffusion.

You can also neglect the diffusion term in concentration equation rather than using a small diffusion coefficient with:

   **Diffusion { negligeable }**

## Post-processing tools

To see the time evolution of velocity and concentration:

- Add a concentration probe near the inlet (e.g. at (0,0.025)) (period 0.01s).

- Add a velocity segment of probes (with 5 points between (0,0.021) and (0,0.029)) at the inlet boundary (period 0.01s).

Now, you are ready to run the study and follow the time evolution with the probes:
```bash
trust -evol diagonale
```

Press the `Start computation!` button and `Plot` or `Plot on same` for probes.

You can now check the flow rate at the inlet boundary in the `diagonale\_pb\_Debit.out file (plotted on the right of the `evol` window). You should find a value close to $6.8 \; 10^{-3} m^2.s^{-1}$.

Use VisIt to post-process the results at $t=0.2s$, $t=0.4s$ and $t=0.7s$. VisIt has some useful features for this study. It can, for example, display a concentration histogram to assess the numerical diffusion in the concentration equation. 

To do so, click on `Add` $\rightarrow$ `Histogram` $\rightarrow$ `CONCENTRATION\_ELEM\_dom`.

The volume of colored water (in $m^3$) is given by $Vol(t)= 6.66.10^{-3} t$ before $t=0.5s$ and $Vol(t)=3.33.10^{-3}$ after.

## VEF calculation
You will now create a variant of your test case. 

First, copy `diagonale.data` to `diagonale_VEF.data`.

In this new file, change the discretization from **VDF** to **VEFPreP1B**. Since the **VEF** discretization only works on simplices, you need to triangulate your mesh by adding the **trianguler ({ref}`triangulate`)** keyword in your `diagonale_VEF.data`.

Change the keyword **quick** (which is only available in VDF) to **muscl** in order to use a **MUSCL** scheme.

You can also switch **GCP** solver to **Cholesky** solver of the Petsc library: **GCP { precond ssor { omega 1.5 } seuil 1.e-6 } $\rightarrow$ Petsc Cholesky { }**. 
The Cholesky method is a direct method that works well on relatively small cases but that may need a large amount of RAM memory for larger problems.

Run the calculation. 

You should encounter an error, and **TRUST** will stop the calculation. You will find a `diagonale_VEF.decoupage_som` file in your working directory.

As **TRUST** indicates, to avoid this problem, you can:

- change the **trianguler** keyword to **{ref}`trianguler_h`**,

- or use the **{ref}`verifiercoin`** keyword. 

The first method is straightforward and works here due to the geometry of your domain. 

To use the second one, you will need the `diagonale_VEF.decoupage_som` file. Add the following: **VerifierCoin dom { read\_file diagonale\_VEF.decoupage\_som }** in your `diagonale\_VEF.data`, just after `trianguler dom`. This will subdivide any inconsistent 2D/3D cells.

Finally, run both of your `.data` files and compare the results between **VDF/quick** and **VEFPreP1B/muscl**.