Add equation example

[tomsail]

Hi @lucduron,

I see you added lots of tools to calculate parameters (in variable).
Could you provide a simple example in order to calculate CHEZY BOTTOM FRICTION and add it in a geo file?

I think I could add it in the README then.

Thanks

[lucduron]

Hi @tomsail,

In 2D, W is the friction variable (see 2D variables), its unit depends on the friction law (Chézy, Strickler...).
It can not be computed (from which variables would you want to do it?) in PyTelTools, but can be used to derive other variables.
For example, it can be used to compute US and/or TAU for example (see equations here).

The handling of variables and equations is optimized and quite complex in PyTelTools, it requires some more files (I am not the original developper of this part).
If you include also variables.py, you should be able to use this kind of script:

import xarray as xr

from xarray_selafin.variable.variables_2d import get_US_equation, NIKURADSE_ID
from xarray_selafin.variables import do_calculation, get_necessary_equations
from xarray_selafin.xarray_backend import SelafinBackendEntrypoint


slf_in = "tests/data/r2d_tidal_flats.slf"
ds = xr.open_dataset(slf_in, engine=SelafinBackendEntrypoint)
# https://gitlab.pam-retd.fr/otm/telemac-mascaret/-/blob/main/examples/telemac3d/tidal_flats/t3d_tidal_flats.cas
# LAW OF BOTTOM FRICTION              = 5
# FRICTION COEFFICIENT FOR THE BOTTOM = 0.01
ds = ds.assign(W=0.01)  # Add a constant Nikuradse of 0.01 (Is it the best way to assign it?)


us_equation = get_US_equation(NIKURADSE_ID)
is_2d = "plan" not in ds.dims
necessary_equations = get_necessary_equations(
    ds.attrs["variables"].keys(),
    ["TAU"],  # Some new variables to add
    is_2d=is_2d,
    us_equation=us_equation
)
for equation in necessary_equations:
    input_var_IDs = list(map(lambda x: x.ID(), equation.input))

    if is_2d:
        # handle the special case for US (user-specified equation)
        if equation.output.ID() == 'US':
            ds['US'] = xr.Variable(
                dims=("time", "node"),
                data=do_calculation(us_equation, [ds['W'], ds['H'], ds['M']])
            )
            # Clean US values in case of negative or null water depth
            xr.where(ds["H"] > 0, ds["US"], 0.0)

    # handle the normal case (if not already done)
    if equation.output.ID() not in ds.data_vars:
        ds[equation.output.ID()] = xr.Variable(
            dims=("time", "node"),
            data=do_calculation(equation, [ds[var_ID] for var_ID in input_var_IDs])
        )

print(ds)
ds.selafin.write("r2d_with_new_variables.slf")

This script will rewrite original variables but add some computed variables:

• M (from U and V)
• US (from W, H and M)
• TAU (from US)

Sorry I do not have much time to integrate it properly, but I am available to help if required.

Hope it helps,
Luc

[tomsail]

Ok thanks for taking the time.

I was interested in varying bottom friction.
I had implemented something in my pipeline already:

manning = 0.027 # 
bf = (abs(ds['B']) ** (1 / 6)) / manning
ds['W'] = xr.Variable(dims=["node"], data=bf)

(for the complete code)
but since you had a lot implemented already, I wanted to know if I could "plug-in" some of your equations.

No rush for this, it can be done down the line.

[lucduron]

I am not used to Chézy but is it not H (or better: the hydraulic radius) instead of B?
Source: https://en.wikipedia.org/wiki/Ch%C3%A9zy_formula#Ch%C3%A9zy_formula_vs_Manning_formula

[tomsail]

Indeed you use the hydraulic radius in manning formula.

However when it comes to ocean modeling, you are in the configuration of a stream whose width is much greater than the depth.
Hence why I use the depth in my models (I exclusively do ocean modeling).

[lucduron]

Ok, I understand.