Source code for fluidsim.util.scripts.turb_trandom_anisotropic

"""Simulations with the solver ns3d.strat and the forcing tcrandom_anisotropic.

.. autofunction:: create_parser

.. autofunction:: main

"""

from math import pi, asin, sin
import argparse
import sys

import matplotlib.pyplot as plt

from fluiddyn.util import mpi
from fluidsim.util.scripts import parse_args

doc = """Launcher for simulations with the solver ns3d.strat and
the forcing tcrandom_anisotropic.

Examples
--------

```
./run_simul.py --only-print-params
./run_simul.py --only-print-params-as-code

./run_simul.py -F 0.3 --delta-F 0.1 --ratio-kfmin-kf 0.8 --ratio-kfmax-kf 1.5 -opf
./run_simul.py -F 1.0 --delta-F 0.1 --ratio-kfmin-kf 0.8 --ratio-kfmax-kf 1.5 -opf

mpirun -np 2 ./run_simul.py
```

Notes
-----

This script is designed to study stratified turbulence forced with an
anisotropic forcing in toroidal or poloidal modes.

The regime depends on the value of the horizontal Froude number Fh and buoyancy
Reynolds numbers R and R4:

- Fh = epsK / (Uh^2 N)
- R = epsK / (N^2 nu)
- R4 = epsK Uh^2 / (nu4 N^4)

Fh has to be very small to be in a strongly stratified regime. R and R4 has to
be "large" to be in a turbulent regime.

For this forcing, we fix the injection rate P (very close to epsK). We will
work at P = 1.0, such that N, nu and nu4 determine the non dimensional numbers.

Note that Uh is not directly fixed by the forcing but should be a function of
the other input parameters. Dimensionally, we can write Uh = (P Lfh)**(1/3).

For simplicity, we'd like to have Lfh=1.0. We want to force at "large
horizontal scale" (compared to the size of the numerical domain). This length
(params.oper.Lx = params.oper.Ly) is computed with this condition.

"""

keys_versus_kind = {
    "toro": "vt_fft",
    "polo": "vp_fft",
    "buoyancy": "b_fft",
    "vert": "vz_fft",
}


[docs] def create_parser(): """Create the argument parser with default arguments""" parser = argparse.ArgumentParser( description=doc, formatter_class=argparse.RawDescriptionHelpFormatter ) # parameters to study the input parameters without running the simulation parser.add_argument( "-oppac", "--only-print-params-as-code", action="store_true", help="Only run initialization phase and print params as code", ) parser.add_argument( "-opp", "--only-print-params", action="store_true", help="Only run initialization phase and print params", ) parser.add_argument( "-opf", "--only-plot-forcing", action="store_true", help="Only run initialization phase and plot forcing", ) # grid parameter parser.add_argument( "-nz", "--nz", type=int, default=32, help="Number of numerical points over one vertical axis", ) # physical parameters parser.add_argument( "--ratio-nh-nz", type=int, default=4, help=( "Ratio nh / nz (has to be a positive integer, " "fixes the aspect ratio of the numerical domain)" ), ) parser.add_argument("--t_end", type=float, default=10.0, help="End time") parser.add_argument( "-N", type=float, default=10.0, help="Brunt-Väisälä frequency" ) parser.add_argument("-nu", type=float, default=None, help="Viscosity") parser.add_argument( "--nu4", type=float, default=None, help="Order-4 hyper viscosity", ) parser.add_argument( "-Rb", type=float, default=None, help="Input buoyancy Reynolds number (injection_rate / (nu * N^2))", ) parser.add_argument( "--Rb4", type=float, default=None, help="Input order-4 buoyancy Reynolds number (injection_rate / (nu_4 * N^4))", ) parser.add_argument( "--coef-nu4", type=float, default=1.0, help="Coefficient used to compute the order-4 hyper viscosity", ) parser.add_argument( "--forced-field", type=str, default="polo", help='Forced field (can be "polo", "toro", ...)', ) parser.add_argument( "--init-velo-max", type=float, default=0.01, help="params.init_fields.noise.max", ) # shape of the forcing region in spectral space parser.add_argument( "--nh-forcing", type=int, default=3, help="Dimensionless horizontal wavenumber forced", ) parser.add_argument( "--ratio-kfmin-kf", type=float, default=0.5, help="", ) parser.add_argument( "--ratio-kfmax-kf", type=float, default=2.0, help="", ) parser.add_argument( "-F", type=float, default=0.3, help="Ratio omega_f / N, fixing the mean angle between the vertical and the forced wavenumber", ) parser.add_argument( "--delta-F", type=float, default=0.1, help="delta F, fixing angle_max - angle_min", ) # Output parameters parser.add_argument( "--spatiotemporal-spectra", action="store_true", help="Activate the output spatiotemporal_spectra", ) # Other parameters parser.add_argument( "--projection", type=str, default=None, help="Type of projection (changes the equations solved)", ) parser.add_argument( "--disable-no-vz-kz0", action="store_true", help="Disable params.no_vz_kz0", ) parser.add_argument( "--disable-NO-SHEAR-MODES", action="store_true", help="Disable params.oper.NO_SHEAR_MODES", ) parser.add_argument( "--max-elapsed", type=str, default="23:50:00", help="Max elapsed time", ) parser.add_argument( "--coef-dealiasing", type=float, default=0.8, help="params.oper.coef_dealiasing", ) parser.add_argument( "--periods_save-phys_fields", type=float, default=1.0, help="params.output.periods_save.phys_fields", ) parser.add_argument( "--sub-directory", type=str, default="aniso", help="Sub directory where the simulation will be saved", ) parser.add_argument( "--modify-params", type=str, default=None, help="Code modifying the `params` object.", ) return parser
def create_params(args): """Create the params object from the script arguments""" from fluidsim.solvers.ns3d.strat.solver import Simul params = Simul.create_default_params() params.output.sub_directory = args.sub_directory params.short_name_type_run = args.forced_field if args.projection is not None: params.projection = args.projection params.short_name_type_run += "_proj" params.no_vz_kz0 = not args.disable_no_vz_kz0 params.oper.NO_SHEAR_MODES = not args.disable_NO_SHEAR_MODES params.oper.coef_dealiasing = args.coef_dealiasing params.oper.truncation_shape = "no_multiple_aliases" params.oper.nz = nz = args.nz nh = nz * args.ratio_nh_nz if args.ratio_nh_nz < 1: raise ValueError("args.ratio_nh_nz < 1") params.oper.nx = params.oper.ny = nh Lfh = 1.0 injection_rate = 1.0 Uh = (injection_rate * Lfh) ** (1 / 3) Lh = args.nh_forcing * Lfh params.oper.Lx = params.oper.Ly = Lh params.oper.Lz = Lh / args.ratio_nh_nz delta_kz = 2 * pi / params.oper.Lz # Brunt Vaisala frequency params.N = args.N Fh = Uh / (args.N * Lfh) params.short_name_type_run += f"_Fh{Fh:.3e}" mpi.printby0(f"Input horizontal Froude number: {Fh:.3g}") nu = None Rb = None if args.nu is not None and args.Rb is not None: raise ValueError("args.nu is not None and args.Rb is not None") if args.nu is not None: nu = args.nu if nu != 0.0: Rb = injection_rate / (nu * args.N**2) params.short_name_type_run += f"_nu{nu:.3e}" mpi.printby0(f"Input viscosity: {nu:.3e}") elif args.Rb is not None: Rb = args.Rb nu = injection_rate / (Rb * args.N**2) params.short_name_type_run += f"_Rb{Rb:.3g}" mpi.printby0(f"Input buoyancy Reynolds number: {Rb:.3g}") else: Rb = 5.0 nu = injection_rate / (Rb * args.N**2) params.nu_2 = nu if args.nu4 is not None and args.Rb4 is not None: raise ValueError("args.nu4 is not None and args.Rb4 is not None") if args.nu4 is not None: params.nu_4 = args.nu4 params.short_name_type_run += f"_nuh{params.nu_4:.3e}" mpi.printby0(f"Input order-4 hyper viscosity: {params.nu_4:.3e}") elif args.Rb4 is not None: Rb_4 = args.Rb4 params.nu_4 = injection_rate / (Rb_4 * args.N**4) params.short_name_type_run += f"_Rbh{Rb_4:.3g}" mpi.printby0(f"Input order-4 buoyancy Reynolds number: {Rb_4:.3g}") else: # compute nu_4 from injection_rate and dx # Kolmogorov length scale eta = (nu**3 / injection_rate) ** 0.25 k_max = params.oper.coef_dealiasing * delta_kz * nz / 2 mpi.printby0(f"{eta * k_max = :.3e}") if eta * k_max > 1: mpi.printby0("Well resolved simulation, no need for nu_4") params.nu_4 = 0.0 else: # TODO: more clever injection_rate_4 (tends to 0 when eta*k_max = 1) injection_rate_4 = injection_rate # only valid if R4 >> 1 (isotropic turbulence at small scales) params.nu_4 = ( args.coef_nu4 * injection_rate_4 ** (1 / 3) / k_max ** (10 / 3) ) Rb_4 = injection_rate / (params.nu_4 * args.N**4) mpi.printby0( f"Resolution too coarse, we add order-4 hyper viscosity nu_4={params.nu_4:.3e}." ) params.init_fields.type = "noise" params.init_fields.noise.length = params.oper.Lz / 2 params.init_fields.noise.velo_max = args.init_velo_max params.forcing.enable = True params.forcing.type = "tcrandom_anisotropic" params.forcing.forcing_rate = injection_rate params.forcing.key_forced = keys_versus_kind[args.forced_field] """ Since args.ratio_nh_nz > 1, - kf_min = delta_kz * nkmin_forcing - kf_max = delta_kz * nkmax_forcing """ def round3(number): return round(number, 3) angle = asin(args.F) mpi.printby0(f"angle = {angle / pi * 180:.2f}°") delta_angle = asin(args.delta_F) kfh = 2 * pi / Lfh kf = kfh / sin(angle) kf_min = kf * args.ratio_kfmin_kf kf_max = kf * args.ratio_kfmax_kf params.forcing.nkmin_forcing = max(0, round3(kf_min / delta_kz)) params.forcing.nkmax_forcing = min(nz // 2, round3(kf_max / delta_kz)) mpi.printby0( f"{params.forcing.nkmin_forcing = }\n{params.forcing.nkmax_forcing = }" ) period_N = 2 * pi / args.N omega_l = args.N * args.F # time_stepping fixed to follow waves params.time_stepping.USE_T_END = True params.time_stepping.t_end = args.t_end params.time_stepping.max_elapsed = args.max_elapsed params.time_stepping.deltat_max = min(0.1, period_N / 16) # time_correlation is fixed to forced wave period params.forcing.tcrandom.time_correlation = 2 * pi / omega_l params.forcing.tcrandom_anisotropic.angle = round3(angle) params.forcing.tcrandom_anisotropic.delta_angle = round3(delta_angle) params.forcing.tcrandom_anisotropic.kz_negative_enable = True params.output.periods_print.print_stdout = 1e-1 params.output.periods_save.phys_fields = args.periods_save_phys_fields params.output.periods_save.spatial_means = 0.02 params.output.periods_save.spectra = 0.05 params.output.periods_save.spect_energy_budg = 0.1 params.output.spectra.kzkh_periodicity = 1 if args.spatiotemporal_spectra: params.output.periods_save.spatiotemporal_spectra = period_N / 8 params.output.spatiotemporal_spectra.file_max_size = 80.0 # (Mo) # probes_region in nondimensional units (mode indices). ikzmax = 16 ikhmax = ikzmax * args.ratio_nh_nz params.output.spatiotemporal_spectra.probes_region = (ikhmax, ikhmax, ikzmax) if args.modify_params is not None: exec(args.modify_params) return params
[docs] def main(args=None, **defaults): """Main function for the scripts based on turb_trandom_anisotropic""" parser = create_parser() if defaults: parser.set_defaults(**defaults) args = parse_args(parser, args) params = create_params(args) if ( args.only_plot_forcing or args.only_print_params_as_code or args.only_print_params ): params.output.HAS_TO_SAVE = False sim = None if args.only_print_params_as_code: params._print_as_code() return params, sim if args.only_print_params: print(params) return params, sim from fluidsim.solvers.ns3d.strat.solver import Simul sim = Simul(params) if args.only_plot_forcing: sim.forcing.forcing_maker.plot_forcing_region() plt.show() return params, sim sim.time_stepping.start() mpi.printby0( f""" # To visualize the output with Paraview, create a file states_phys.xmf with: fluidsim-create-xml-description {sim.output.path_run} # To visualize with fluidsim: fluidsim-ipy-load {sim.output.path_run} # in IPython: sim.output.phys_fields.set_equation_crosssection('x={params.oper.Lx/2}') sim.output.phys_fields.animate('b') """ ) return params, sim
if "sphinx" in sys.modules: from textwrap import indent from unittest.mock import patch with patch.object(sys, "argv", ["run_simul.py"]): parser = create_parser() __doc__ += """ Example of help message ----------------------- .. code-block:: """ + indent( parser.format_help(), " " ) if __name__ == "__main__": params, sim = main()