Commit 92e9f57e authored by Sævik Pål Næverlid's avatar Sævik Pål Næverlid
Browse files

Merge branch 'feature/collision_treatment'

parents e05f91bf 40246c09
# Chemicals
The module represents passive particles that have a vertical migration rate due
to upwelling/downwelling and vertical turbulence.
## Usage
Modify `ladim.yaml` and `particles.rls` to suit your needs.
Common changes applied to `ladim.yaml`:
- Start date of simulation (`time_control.start_time`)
- Stop date of simulation (`time_control.stop_time`)
- Forcing input file (`gridforce.input_file`)
- Horizontal diffusion parameter (`numerics.diffusion`)
- Vertical diffusion parameter (`ibm.vertical_mixing`)
- Time step length (`numerics.dt`)
- Output frequency (`output_variables.outper`)
The file `particles.rls` is a tab-delimited text file containing particle
release time and location, as well as particle attributes at the release time.
The order of the columns is given by the entry `particle_release.variables`
within `ladim.yaml`.
Finally, copy `ladim.yaml` and `particles.rls` to a separate directory and
run `ladim` here.
## Output
The simulation result is stored in a file specified by the `files.output_file`
entry in `ladim.yaml`. The output variables are specified by the
`output_variables` entries.
from .gridforce import Grid, Forcing
from .ibm import IBM
This diff is collapsed.
import numpy as np
class IBM:
def __init__(self, config):
self.D = config["ibm"].get('vertical_mixing', 0) # Vertical mixing [m*2/s]
self.dt = config['dt']
self.x = np.array([])
self.y = np.array([])
self.pid = np.array([])
self.land_collision = config["ibm"].get('land_collision', 'reposition')
def update_ibm(self, grid, state, forcing):
# Vertical advection velocity
W = forcing.forcing.wvel(state.X, state.Y, state.Z)
# Vertical diffusion velocity
rand = np.random.normal(size=len(state.X))
W += rand * (2 * self.D / self.dt) ** 0.5
# Update vertical position, using reflexive boundary condition at top
state.Z += W * self.dt
state.Z[state.Z < 0] *= -1
# Reflexive boundary condition at bottom
H = grid.sample_depth(state.X, state.Y) # Water depth
below_seabed = state.Z > H
state.Z[below_seabed] = 2*H[below_seabed] - state.Z[below_seabed]
if self.land_collision == "reposition":
# If particles have not moved: Assume they ended up on land.
# If that is the case, reposition them within the cell.
pid, pidx_old, pidx_new = np.intersect1d(self.pid, state.pid, return_indices=True)
onland = ((self.x[pidx_old] == state.X[pidx_new]) &
(self.y[pidx_old] == state.Y[pidx_new]))
num_onland = np.count_nonzero(onland)
pidx_new_onland = pidx_new[onland]
x_new = np.round(state.X[pidx_new_onland]) - 0.5 + np.random.rand(num_onland)
y_new = np.round(state.Y[pidx_new_onland]) - 0.5 + np.random.rand(num_onland)
state.X[pidx_new_onland] = x_new
state.Y[pidx_new_onland] = y_new
self.x = state.X
self.y = state.Y
self.pid = state.pid
# Configuration file for sedimentation model
time_control:
# Start and stop of simulation
start_time : 2015-09-07 01:00:00
stop_time : 2015-09-07 01:05:00
files:
particle_release_file : particles.rls
output_file : out.nc
gridforce:
# Module name is mandatory.
# The other items are module dependent
module: ladim_ibm.chemicals
input_file: forcing.nc # Use wildcards (*) to select a range of files
ibm:
module: ladim_ibm.chemicals
vertical_mixing: 0.0001 # [m*2/s]
particle_release:
variables: [release_time, X, Y, Z]
release_time: time # np.datetime64[s]
# Mark variables as time independent
particle_variables: [release_time]
output_variables:
# Frequency of output
outper: [60, s]
# Variables included in output
particle: [release_time]
instance: [pid, X, Y, Z]
# --- Output format for standard variables ---
# Output format for the particle release time
release_time:
ncformat: f8
long_name: particle release time
units: seconds since reference_time
# Output format for the particle identifier
pid:
ncformat: i4
long_name: particle identifier
# Output format for the X coordinate
X:
ncformat: f4
long_name: particle X-coordinate
# Output format for the Y coordinate
Y:
ncformat: f4
long_name: particle Y-coordinate
# Output format for the particle depth
Z:
ncformat: f4
long_name: particle depth
standard_name: depth_below_surface
units: m
positive: down
numerics:
# Model time step
dt: [60, s] # Format = [value, unit]
# Advection method, EF, RK2, or RK4
advection: RK4
diffusion: 1 # [m/s**2]
2015-09-07T01:00:00 3 3 31.793
2015-09-07T01:01:00 2 2 0
2015-09-07T01:01:00 7 5.5001 5
import numpy as np
import pytest
import pkg_resources
from ladim_ibm.chemicals import gridforce, IBM
class Test_nearest_unmasked:
def test_correct_when_all_unmasked(self):
mask = np.zeros((4, 3))
i = np.array([1, 1, 2])
j = np.array([2, 3, 3])
ii, jj = gridforce.nearest_unmasked(mask, i, j)
assert ii.tolist() == i.tolist()
assert jj.tolist() == j.tolist()
def test_correct_when_south_edge(self):
mask_south = np.array([[0, 0, 0, 0], [1, 1, 1, 1], [1, 1, 1, 1]])
i = np.array([0, 1, 2.51])
j = np.array([0, 1, 1.49])
ii, jj = gridforce.nearest_unmasked(mask_south, i, j)
assert ii.tolist() == [0, 1, 3]
assert jj.tolist() == [0, 0, 0]
def test_correct_when_corner(self):
mask = np.array([[0, 0, 0, 0], [0, 1, 1, 1], [0, 1, 1, 0]])
i = np.array([0.51, 0.51, 0.99, 1.49, 1.51, 2.00, 3.00])
j = np.array([0.52, 0.98, 0.52, 1.01, 1.01, 1.01, 1.01])
ii, jj = gridforce.nearest_unmasked(mask, i, j)
assert ii.tolist() == [0, 0, 1, 1, 2, 2, 3]
assert jj.tolist() == [1, 1, 0, 0, 0, 0, 2]
class Test_ibm_land_collision:
@pytest.fixture()
def config(self):
from ladim.configuration import configure
pkg = 'ladim_ibm.chemicals'
with pkg_resources.resource_stream(pkg, 'ladim.yaml') as config_file:
return configure(config_file)
@pytest.fixture()
def grid(self, config):
from ladim.gridforce import Grid
return Grid(config)
@pytest.fixture()
def forcing(self, config, grid):
from ladim.gridforce import Forcing
return Forcing(config, grid)
@pytest.fixture()
def state(self, config, grid):
from ladim.state import State
return State(config, grid)
@pytest.fixture()
def ibm_chemicals(self, config):
return IBM(config)
def test_land_collision(self, ibm_chemicals, grid, state, forcing):
np.random.seed(0)
state.X = np.float32([1, 1, 1])
state.Y = np.float32([1, 1, 1])
state.Z = np.float32([1, 1, 1])
state.pid = np.int32([0, 1, 2])
ibm_chemicals.update_ibm(grid, state, forcing)
assert state.X.tolist() == [1, 1, 1]
assert state.Y.tolist() == [1, 1, 1]
state.X = np.float32([1, 2, 3, 4])
state.Y = np.float32([1, 1, 1, 1])
state.Z = np.float32([1, 1, 1, 1])
state.pid = np.int32([1, 2, 3, 4])
ibm_chemicals.update_ibm(grid, state, forcing)
assert state.X.tolist() == [1.4636627435684204, 2, 3, 4]
assert state.Y.tolist() == [0.8834415078163147, 1, 1, 1]
......@@ -10,5 +10,5 @@ setup(
author='Pål Næverlid Sævik',
author_email='a5606@hi.no',
description='IBMs for LADiM', install_requires=['numpy', 'pytest', 'xarray',
'PyYAML'],
'PyYAML', 'netCDF4'],
)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment