=======================
The user options file
=======================

Simulation parameters are set in a user options file that is read in by *fluvial-particle* and evaluated as a Python script. See the examples section for a sample user options file.

Some parameters are required in the options file, while others are optional.

Required keyword arguments
============================

**file_name_2d**, str: Path to the 2D input mesh. Supported formats:

* ``.vts`` - VTK XML Structured Grid format (**recommended**). Modern binary format with compression, 5-10x smaller than legacy VTK.
* ``.vtk`` - VTK legacy format (ASCII or binary). Widely compatible with older tools.
* ``.npz`` - NumPy compressed archive. Python-specific format for custom workflows.

**file_name_3d**, str: Path to the 3D input mesh. Supports the same formats as file_name_2d. In the case of a 2D simulation, this entry will be ignored but it is still required.

**SimTime**, float: Maximum allowable simulation time, in seconds. The simulation will end before this time only if all of the tracked particles leave the domain before SimTime.

**dt**, float: The time step used in the simulation, in seconds.

**PrintAtTick**, float: The printing interval, in seconds. Data will be printed to stdout and written to the HDF5 files at intervals of less than or equal to PrintAtTick.

**Track3D**, bool (int 0 or 1): Indicates whether to run a 2D or 3D simulation. If set to 1, the simulation will be 3D and velocity vectors will be interpolated from the 3D mesh. Else, the simulation will be 2D with velocities from the 2D mesh.

**NumPart**, int: The number of particles to simulate, *per core*. This means that in parallel execution mode with *N* cores, the total number of simulated particles will be *N* * NumPart.

**StartLoc**, tuple or str: Indicates the starting location of the particles, can be given as a tuple of length 3 or as a string. Particles can be initialized from a single point, from an existing HDF5 file of particle locations, or from a time-delayed collection of points specified in a CSV file.

* tuple (x, y, z): all particles will be initiated from the same point
* str: path to a file from which these starting positions will be loaded; either an HDF5 or CSV file

 - ".h5" suffix: the starting locations will be loaded from the HDF5 file as a checkpoint file generated from a previous *fluvial-particle* simulation. The checkpoint file must have the same *total* number of particles as the current simulation, i.e. summed across all CPU cores. By default, data are loaded from the final entry in the HDF5 file, but a particular index to slice into along the time dimension (axis 0) can be provided with the *StartIdx* optional keyword argument.
 - ".csv" suffix: the starting locations will be loaded from a CSV file with 5 columns: start_time, x, y, z, numpart. For example, if a given row in the CSV file is "10.0, 6.14, 9.09, 10.3, 100", then 100 particles will be initiated from the point (6.14, 9.09, 10.3) starting at a simulation time of 10.0 seconds.

**ParticleType**: The type of particles to simulate, either the Particles class or a subclass (e.g. LarvalBotParticles). This argument should not be placed inside quotes or brackets of any kind.

**field_map_2d**, dict: A dictionary mapping standard internal field names to the model-specific names in your 2D mesh file. This allows *fluvial-particle* to work with output from different hydrodynamic models (e.g., Delft-FM, iRIC, HEC-RAS) that use different naming conventions.

Core required keys:

* ``bed_elevation``: bed/bottom elevation (e.g., "Elevation" in Delft-FM)
* ``velocity``: velocity vector (e.g., "Velocity")
* ``water_surface_elevation``: water surface elevation (e.g., "WaterSurfaceElevation")

Shear velocity source (at least one required):

* ``ustar``: direct shear velocity field [m/s] (preferred if available)
* ``shear_stress``: bed shear stress [Pa] (commonly available from hydrodynamic models)
* ``manning_n``: Manning's n field [-] (or use scalar option below)
* ``chezy_c``: Chézy C field [m^0.5/s] (or use scalar option)
* ``darcy_f``: Darcy-Weisbach f field [-] (or use scalar option)
* ``energy_slope``: energy slope field [-]
* ``tke``: turbulent kinetic energy field [m²/s²]

Optional key:

* ``wet_dry``: wet/dry indicator, 1=wet, 0=dry (e.g., "IBC" in Delft-FM). If ``wet_dry`` is omitted from ``field_map_2d``, this field will be computed automatically from depth using the ``min_depth`` threshold (defaults to 0.02m if not specified in the options file): cells with depth > min_depth are wet (1), otherwise dry (0).

Example with shear stress:

.. code-block:: python

    field_map_2d = {
        "bed_elevation": "Elevation",
        "shear_stress": "ShearStress (magnitude)",
        "velocity": "Velocity",
        "water_surface_elevation": "WaterSurfaceElevation",
    }

Example with direct ustar field:

.. code-block:: python

    field_map_2d = {
        "bed_elevation": "Elevation",
        "ustar": "FrictionVelocity",
        "velocity": "Velocity",
        "water_surface_elevation": "WaterSurfaceElevation",
    }

**field_map_3d**, dict: A dictionary mapping standard internal field names to the model-specific names in your 3D mesh file. Required keys:

* ``velocity``: velocity vector (e.g., "Velocity")

Example:

.. code-block:: python

    field_map_3d = {
        "velocity": "Velocity",
    }

.. note::
   For ``.npz`` files, the ``ibc`` field (wet/dry indicator) is optional. If not present in the npz file, ``wet_dry`` will be computed from depth using the ``min_depth`` threshold.


Optional keyword arguments
============================

These arguments can be specified in the options file. Otherwise, the default values will be used.

**beta**, float or tuple: Scales the 3D diffusion coefficients, can be specified as a scalar or a tuple of length 3. Default value: (0.067, 0.067, 0.067)

**lev**, float: Lateral eddy viscosity. Default value: 0.25

**min_depth**, float: The minimum depth that a particle may enter. If a depth update is less than min_depth, then the update is not permitted. Default value: 0.02

**StartIdx**, int: The index used to slice into the HDF5 particles checkpoint file along the 0th axis, i.e. the printing step axis. Only used if an HDF5 file is provided via the StartLoc keyword argument. Default value: -1

**output_vtp**, bool: If True, write particle output to VTK PolyData (.vtp) files in addition to HDF5. A PVD collection file (particles.pvd) is also created to reference all VTP files with timestamp information. VTP files are natively supported by ParaView without plugins and provide a simpler alternative to XDMF+HDF5 for visualization. Output is written to a ``vtp/`` subdirectory within the output directory. Default value: False

**startfrac**, float: If provided, will initialize particles to a vertical position as bed elevation + water depth multiplied with startfrac. startfrac should be between 0 and 1 -- values outside this range will initialize particles at the bed and water surface, respectively. A numpy array of length NumPart can also be used to vary the startfrac for every particle. Default value: None

**vertbound**, float: Bounds the particles in the fractional water column of [vertbound, 1-vertbound]. This prevents particles from moving out of the vertical domain, either by going below the channel bed or above the water surface. Default value: 0.01


Shear velocity (u*) configuration
====================================

*fluvial-particle* uses shear velocity (u*) to compute turbulent diffusion coefficients.
There are several methods to provide or compute u*, listed here in priority order.
If multiple sources are available, the highest priority method is automatically selected.

**Method 1: Direct u* field** (highest priority)

If your hydrodynamic model outputs shear velocity directly, map it in ``field_map_2d``:

.. code-block:: python

    field_map_2d = {
        "bed_elevation": "Elevation",
        "ustar": "FrictionVelocity",  # Direct u* [m/s]
        "velocity": "Velocity",
        "water_surface_elevation": "WaterSurfaceElevation",
    }

**Method 2: Bed shear stress**

Computes u* = √(τ_b / ρ), where τ_b is bed shear stress and ρ is water density.

.. code-block:: python

    field_map_2d = {
        "bed_elevation": "Elevation",
        "shear_stress": "ShearStress (magnitude)",  # Bed shear stress [Pa]
        "velocity": "Velocity",
        "water_surface_elevation": "WaterSurfaceElevation",
    }

    # Optional: customize water density (default 1000 kg/m³)
    water_density = 1025.0  # for seawater

**Method 3: Manning's n**

Computes u* = U · n · √g / h^(1/6). Can be provided as a field or scalar:

.. code-block:: python

    # As a scalar (uniform friction):
    manning_n = 0.03

    # Or as a spatially-varying field:
    field_map_2d = {
        "bed_elevation": "Elevation",
        "manning_n": "ManningN",
        "velocity": "Velocity",
        "water_surface_elevation": "WaterSurfaceElevation",
    }

**Method 4: Chézy C**

Computes u* = U · √g / C:

.. code-block:: python

    chezy_c = 50.0  # scalar value

**Method 5: Darcy-Weisbach f**

Computes u* = U · √(f / 8):

.. code-block:: python

    darcy_f = 0.02  # scalar value

**Method 6: Energy slope**

Computes u* = √(g · h · S), where S is the energy slope:

.. code-block:: python

    field_map_2d = {
        "bed_elevation": "Elevation",
        "energy_slope": "WaterSurfaceSlope",
        "velocity": "Velocity",
        "water_surface_elevation": "WaterSurfaceElevation",
    }

**Method 7: Turbulent kinetic energy (TKE)**

Computes u* = C_μ^(1/4) · √k, where C_μ = 0.09. Appropriate for RANS CFD outputs:

.. code-block:: python

    field_map_2d = {
        "bed_elevation": "Elevation",
        "tke": "TurbulentKineticEnergy",
        "velocity": "Velocity",
        "water_surface_elevation": "WaterSurfaceElevation",
    }

**Forcing a specific method**

If multiple u* sources are available, the highest-priority method is used by default.
To override this, use ``ustar_method``:

.. code-block:: python

    ustar_method = "manning"  # Force Manning's n method even if shear_stress is available


Time-varying grid settings
=============================

For simulations with unsteady (time-varying) flow conditions, *fluvial-particle* can load a sequence of pre-computed velocity fields and interpolate between them during the simulation. This is useful when:

* Flow conditions change significantly during the simulation period
* You have output from an unsteady CFD simulation
* You want to capture tidal, flood, or other transient flow effects

To enable time-varying grids, set ``time_dependent = True`` and provide the following settings:

**time_dependent**, bool: Enable time-varying grid mode. When True, velocity fields are loaded from a sequence of grid files and interpolated temporally. Default value: False

**file_pattern_2d**, str: Format string for 2D grid files with ``{}`` as placeholder for the file index. Example: ``"./data/flow_2d_{}.vts"`` will load ``flow_2d_0.vts``, ``flow_2d_1.vts``, etc.

**file_pattern_3d**, str: Format string for 3D grid files, same pattern as file_pattern_2d.

**grid_start_index**, int: Index of the first grid file in the sequence.

**grid_end_index**, int: Index of the last grid file in the sequence (inclusive).

**grid_dt**, float: Time interval between grid files in seconds.

**grid_start_time**, float: Simulation time corresponding to the first grid file. Default value: 0.0

**grid_interpolation**, str: Temporal interpolation method between grid timesteps. Options:

* ``"linear"`` (default): Linear interpolation between grid timesteps. Provides smooth velocity transitions.
* ``"nearest"``: Use the nearest grid timestep. Snaps to the closest available timestep.
* ``"hold"``: Use the most recent grid until the next timestep. Step-function behavior.

Example time-varying configuration:

.. code-block:: python

    # Enable time-varying mode
    time_dependent = True

    # File patterns (grid files: Result_2D_0.vts through Result_2D_10.vts)
    file_pattern_2d = "./data/unsteady/Result_2D_{}.vts"
    file_pattern_3d = "./data/unsteady/Result_3D_{}.vts"

    # Grid file indices
    grid_start_index = 0
    grid_end_index = 10

    # Timing: each grid file represents 60 seconds of flow
    grid_dt = 60.0
    grid_start_time = 0.0

    # Linear interpolation between timesteps
    grid_interpolation = "linear"

.. note::
   When using time-varying grids, ``file_name_2d`` and ``file_name_3d`` are ignored. The grid files are loaded based on the file patterns instead.

.. note::
   The ``TimeVaryingGrid`` class uses a sliding window approach, keeping only 2 grids in memory at a time. This allows simulations with many timesteps without excessive memory usage.


LarvalParticles optional keyword arguments
=============================================

**amp**, float or np.ndarray: The amplitude of larval sinusoidal swimming behavior. If a NumPy ndarray is provided, then it must be 1D and have length equal to NumPart. Only for the LarvalParticles subclasses. Default value: 0.2

**period**, float or np.ndarray: The temporal period of larval sinusoidal swimming behavior, in seconds. If a NumPy ndarray is provided, then it must be 1D and have length equal to NumPart. Only for the LarvalParticles subclasses. Default value: 60.0


FallingParticles optional keyword arguments
=============================================

**c1**, float or np.ndarray: The viscous drag coefficient, dimensionless. If a NumPy ndarray is provided, then it must be 1D and have length equal to NumPart. Only for the FallingParticles subclasses. Default value: 20.0

**c2**, float or np.ndarray: The turbulent wake drag coefficient, dimensionless. If a NumPy ndarray is provided, then it must be 1D and have length equal to NumPart. Only for the FallingParticles subclasses. Default value: 1.1

**radius**, float or np.ndarray: The particle radii, meters. If a NumPy ndarray is provided, then it must be 1D and have length equal to NumPart. Only for the FallingParticles subclasses. Default value: 0.0005

**rho**, float or np.ndarray: The particle density, kilograms per cubic meter. If a NumPy ndarray is provided, then it must be 1D and have length equal to NumPart. Only for the FallingParticles subclasses. Default value: 0.0005