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

Add working land collision implementation

parent e39cb777
......@@ -5,6 +5,10 @@ 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
......@@ -22,3 +26,20 @@ class IBM:
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
import numpy as np
from ladim_ibm.chemicals import gridforce
import pytest
import pkg_resources
from ladim_ibm.chemicals import gridforce, IBM
class Test_nearest_unmasked:
......@@ -26,3 +28,52 @@ class Test_nearest_unmasked:
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]
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