Turbidity Current in Axisymmetric Configuration

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr  # noqa: F401

import xcompact3d_toolbox as x3d

Parameters

• Numerical precision

Use np.float64 if Xcompact3d was compiled with the flag -DDOUBLE_PREC, use np.float32 otherwise.

x3d.param["mytype"] = np.float32

• Xcompact3d’s parameters

For more information about them, checkout the API reference.

prm = x3d.Parameters(
    filename="input.i3d",
    # BasicParam
    itype=12,
    nx=501,
    ny=73,
    nz=501,
    xlx=20.0,
    yly=2.0,
    zlz=20.0,
    nclx1=1,
    nclxn=1,
    ncly1=2,
    nclyn=1,
    nclz1=1,
    nclzn=1,
    re=3500.0,
    init_noise=0.0125,
    dt=5e-4,
    ifirst=1,
    ilast=80000,
    numscalar=1,
    # ScalarParam
    nclxS1=1,
    nclxSn=1,
    nclyS1=1,
    nclySn=1,
    nclzS1=1,
    nclzSn=1,
    sc=[1.0],
    ri=[0.5],
    uset=[0.0],
    cp=[1.0],
)

Setup

Everything needed is in one Dataset (see API reference):

ds = x3d.init_dataset(prm)

Let’s see it, data and attributes are attached, try to interact with the icons:

ds

<xarray.Dataset> Size: 293MB
Dimensions:  (x: 501, y: 73, z: 501, n: 1)
Coordinates:
  * x        (x) float32 2kB 0.0 0.04 0.08 0.12 0.16 ... 19.88 19.92 19.96 20.0
  * y        (y) float32 292B 0.0 0.02778 0.05556 0.08333 ... 1.944 1.972 2.0
  * z        (z) float32 2kB 0.0 0.04 0.08 0.12 0.16 ... 19.88 19.92 19.96 20.0
  * n        (n) int64 8B 1
Data variables:
    ux       (x, y, z) float32 73MB 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    uy       (x, y, z) float32 73MB 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    uz       (x, y, z) float32 73MB 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    phi      (n, x, y, z) float32 73MB 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

xarray.Dataset

• Dimensions:
  • x: 501
  • y: 73
  • z: 501
  • n: 1
• Coordinates: (4)
  • x
    (x)
    float32
    0.0 0.04 0.08 ... 19.92 19.96 20.0

    name :

    Streamwise coordinate

    long_name :

    $x_1$

    array([ 0.  ,  0.04,  0.08, ..., 19.92, 19.96, 20.  ], dtype=float32)

  • y
    (y)
    float32
    0.0 0.02778 0.05556 ... 1.972 2.0

    name :

    Vertical coordinate

    long_name :

    $x_2$

    array([0.      , 0.027778, 0.055556, 0.083333, 0.111111, 0.138889, 0.166667,
           0.194444, 0.222222, 0.25    , 0.277778, 0.305556, 0.333333, 0.361111,
           0.388889, 0.416667, 0.444444, 0.472222, 0.5     , 0.527778, 0.555556,
           0.583333, 0.611111, 0.638889, 0.666667, 0.694444, 0.722222, 0.75    ,
           0.777778, 0.805556, 0.833333, 0.861111, 0.888889, 0.916667, 0.944444,
           0.972222, 1.      , 1.027778, 1.055556, 1.083333, 1.111111, 1.138889,
           1.166667, 1.194444, 1.222222, 1.25    , 1.277778, 1.305556, 1.333333,
           1.361111, 1.388889, 1.416667, 1.444444, 1.472222, 1.5     , 1.527778,
           1.555556, 1.583333, 1.611111, 1.638889, 1.666667, 1.694444, 1.722222,
           1.75    , 1.777778, 1.805556, 1.833333, 1.861111, 1.888889, 1.916667,
           1.944444, 1.972222, 2.      ], dtype=float32)

  • z
    (z)
    float32
    0.0 0.04 0.08 ... 19.92 19.96 20.0

    name :

    Spanwise coordinate

    long_name :

    $x_3$

    array([ 0.  ,  0.04,  0.08, ..., 19.92, 19.96, 20.  ], dtype=float32)

  • n
    (n)
    int64
    1

    name :

    Scalar fraction

    long_name :

    $\ell$

    array([1])

• Data variables: (4)
  • ux
    (x, y, z)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    file_name :

    ux

    name :

    Initial Condition for Streamwise Velocity

    long_name :

    $u_1 (x_1,x_2,x_3,t=0)$

    BC :

    {'x': {'ncl1': 1, 'ncln': 1, 'npaire': 0}, 'y': {'ncl1': 2, 'ncln': 1, 'npaire': 1, 'istret': <Istret.NO_REFINEMENT: 0>, 'beta': 1.0}, 'z': {'ncl1': 1, 'ncln': 1, 'npaire': 1}}

    array([[[0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            ...,
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.]],
    
           [[0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            ...,
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.]],
    
           [[0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            ...,
    ...
            ...,
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.]],
    
           [[0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            ...,
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.]],
    
           [[0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            ...,
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.]]], dtype=float32)

  • uy
    (x, y, z)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    file_name :

    uy

    name :

    Initial Condition for Vertical Velocity

    long_name :

    $u_2 (x_1,x_2,x_3,t=0)$

    BC :

    {'x': {'ncl1': 1, 'ncln': 1, 'npaire': 1}, 'y': {'ncl1': 2, 'ncln': 1, 'npaire': 0, 'istret': <Istret.NO_REFINEMENT: 0>, 'beta': 1.0}, 'z': {'ncl1': 1, 'ncln': 1, 'npaire': 1}}

    array([[[0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            ...,
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.]],
    
           [[0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            ...,
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.]],
    
           [[0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            ...,
    ...
            ...,
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.]],
    
           [[0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            ...,
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.]],
    
           [[0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            ...,
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.]]], dtype=float32)

  • uz
    (x, y, z)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    file_name :

    uz

    name :

    Initial Condition for Spanwise Velocity

    long_name :

    $u_3 (x_1,x_2,x_3,t=0)$

    BC :

    {'x': {'ncl1': 1, 'ncln': 1, 'npaire': 1}, 'y': {'ncl1': 2, 'ncln': 1, 'npaire': 1, 'istret': <Istret.NO_REFINEMENT: 0>, 'beta': 1.0}, 'z': {'ncl1': 1, 'ncln': 1, 'npaire': 0}}

    array([[[0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            ...,
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.]],
    
           [[0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            ...,
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.]],
    
           [[0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            ...,
    ...
            ...,
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.]],
    
           [[0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            ...,
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.]],
    
           [[0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            ...,
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.],
            [0., 0., 0., ..., 0., 0., 0.]]], dtype=float32)

  • phi
    (n, x, y, z)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    file_name :

    phi

    name :

    Initial Condition for Scalar field(s)

    long_name :

    $\varphi (x_1,x_2,x_3,n,t=0)$

    BC :

    {'x': {'ncl1': 1, 'ncln': 1, 'npaire': 1}, 'y': {'ncl1': 1, 'ncln': 1, 'npaire': 1, 'istret': <Istret.NO_REFINEMENT: 0>, 'beta': 1.0}, 'z': {'ncl1': 1, 'ncln': 1, 'npaire': 1}}

    array([[[[0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.],
             ...,
             [0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.]],
    
            [[0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.],
             ...,
             [0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.]],
    
            [[0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.],
             ...,
    ...
             ...,
             [0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.]],
    
            [[0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.],
             ...,
             [0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.]],
    
            [[0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.],
             ...,
             [0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.],
             [0., 0., 0., ..., 0., 0., 0.]]]], dtype=float32)

• Indexes: (4)
  • x
    PandasIndex

    PandasIndex(Index([                0.0, 0.03999999910593033, 0.07999999821186066,
           0.11999999731779099,  0.1599999964237213, 0.20000000298023224,
           0.23999999463558197,  0.2800000011920929,  0.3199999928474426,
           0.36000001430511475,
           ...
            19.639999389648438,   19.68000030517578,  19.719999313354492,
            19.760000228881836,  19.799999237060547,   19.84000015258789,
              19.8799991607666,  19.920000076293945,  19.959999084472656,
                          20.0],
          dtype='float32', name='x', length=501))

  • y
    PandasIndex

    PandasIndex(Index([                0.0, 0.02777777798473835,  0.0555555559694767,
            0.0833333358168602,  0.1111111119389534,  0.1388888955116272,
            0.1666666716337204,  0.1944444477558136,  0.2222222238779068,
                          0.25,  0.2777777910232544,  0.3055555522441864,
            0.3333333432674408,  0.3611111044883728,  0.3888888955116272,
            0.4166666567325592,  0.4444444477558136,  0.4722222089767456,
                           0.5,  0.5277777910232544,  0.5555555820465088,
            0.5833333134651184,  0.6111111044883728,  0.6388888955116272,
            0.6666666865348816,  0.6944444179534912,  0.7222222089767456,
                          0.75,  0.7777777910232544,  0.8055555820465088,
            0.8333333134651184,  0.8611111044883728,  0.8888888955116272,
            0.9166666865348816,  0.9444444179534912,  0.9722222089767456,
                           1.0,  1.0277777910232544,  1.0555555820465088,
            1.0833333730697632,  1.1111111640930176,  1.1388888359069824,
            1.1666666269302368,  1.1944444179534912,  1.2222222089767456,
                          1.25,  1.2777777910232544,  1.3055555820465088,
            1.3333333730697632,  1.3611111640930176,  1.3888888359069824,
            1.4166666269302368,  1.4444444179534912,  1.4722222089767456,
                           1.5,  1.5277777910232544,  1.5555555820465088,
            1.5833333730697632,  1.6111111640930176,  1.6388888359069824,
            1.6666666269302368,  1.6944444179534912,  1.7222222089767456,
                          1.75,  1.7777777910232544,  1.8055555820465088,
            1.8333333730697632,  1.8611111640930176,  1.8888888359069824,
            1.9166666269302368,  1.9444444179534912,  1.9722222089767456,
                           2.0],
          dtype='float32', name='y'))

  • z
    PandasIndex

    PandasIndex(Index([                0.0, 0.03999999910593033, 0.07999999821186066,
           0.11999999731779099,  0.1599999964237213, 0.20000000298023224,
           0.23999999463558197,  0.2800000011920929,  0.3199999928474426,
           0.36000001430511475,
           ...
            19.639999389648438,   19.68000030517578,  19.719999313354492,
            19.760000228881836,  19.799999237060547,   19.84000015258789,
              19.8799991607666,  19.920000076293945,  19.959999084472656,
                          20.0],
          dtype='float32', name='z', length=501))

  • n
    PandasIndex

    PandasIndex(Index([1], dtype='int64', name='n'))

• Attributes: (0)

Initial Condition

A random noise will be applied to the initial velocity field, we start creating a modulation function mod, to apply it just near the fresh/turbidity interface:

# Position of the initial interface in the polar coordinate
r0 = 4.0

# Random noise with fixed seed,
# important for reproducibility, development and debugging
if prm.iin == 2:
    np.random.seed(seed=67)

radius = np.sqrt(ds.x**2 + ds.z**2.0)

mod = np.exp(-25.0 * (radius - r0) ** 2.0)

# This attribute will be shown at the colorbar
mod.attrs["long_name"] = "Noise modulation"

mod.plot()

<matplotlib.collections.QuadMesh at 0x7f947c174fe0>

Now we reset velocity fields ds[key] *= 0.0, just to guarantee consistency in the case of multiple executions of this cell. Notice that ds[key] = 0.0 may overwrite all the metadata contained in the array, so it should be avoided.

We then add a random number array with the right shape and multiply by the noise amplitude at the initial condition init_noise and multiply again by our modulation function mod, defined previously.

Plotting a xarray.DataArray is as simple as da.plot() (see its user guide), I’m adding extra options just to exemplify how easily we can slice the vertical coordinate and produce multiple plots:

for key in "ux uy uz".split():
    #
    print(ds[key].attrs["name"])
    #
    ds[key] *= 0.0
    #
    ds[key] += prm.init_noise * (np.random.random(ds[key].shape) - 0.5)
    ds[key] *= mod
    #
    ds[key].sel(y=slice(None, None, ds.y.size // 3)).plot(x="x", y="z", col="y", col_wrap=2)
    plt.show()

plt.close("all")

Initial Condition for Streamwise Velocity

Initial Condition for Vertical Velocity

Initial Condition for Spanwise Velocity

A smooth transition at the interface fresh/turbidity fluid is used for the initial concentration(s) field(s), it is defined in the polar coordinates as:

\[ \varphi_n = c_{0,n} \dfrac{1}{2} \left( 1 - \tanh \left( (r - r_0) \sqrt{Sc_n Re} \right) \right). \]

The code block includes the same procedures, we reset the scalar field ds['phi'] *= 0.0, just to guarantee consistency, we compute the equation above, we add it to the array and make a plot:

# Concentration

print(ds["phi"].attrs["name"])

ds["phi"] *= 0.0

for n in range(prm.numscalar):
    #
    fun = 0.5 * prm.cp[n] * (1.0 - np.tanh((radius - r0) * (prm.sc[n] * prm.re) ** 0.5))
    #
    ds["phi"][{"n": n}] += fun
    #
    ds.phi.isel(n=n).sel(y=prm.yly / 2.0).T.plot()
    plt.show()

plt.close("all")

Initial Condition for Scalar field(s)

Writing to disc

is as simple as:

prm.dataset.write(ds)

prm.write()

Running the Simulation

It was just to show the capabilities of xcompact3d_toolbox.sandbox, keep in mind the aspects of numerical stability of our Navier-Stokes solver. It is up to the user to find the right set of numerical and physical parameters.

Make sure that the compiling flags and options at Makefile are what you expect. Then, compile the main code at the root folder with make.

And finally, we are good to go:

mpirun -n [number of cores] ./xcompact3d |tee log.out