diff --git a/PySDM/backends/impl_numba/methods/displacement_methods.py b/PySDM/backends/impl_numba/methods/displacement_methods.py index ac7b5fa44..7c8572462 100644 --- a/PySDM/backends/impl_numba/methods/displacement_methods.py +++ b/PySDM/backends/impl_numba/methods/displacement_methods.py @@ -14,12 +14,15 @@ @numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False}}) # pylint: disable=too-many-arguments def calculate_displacement_body_common( - dim, droplet, scheme, _l, _r, displacement, courant, position_in_cell, n_substeps + dim, droplet, scheme, _l, _r, displacement, courant, position_in_cell, cell_id, n_substeps, enable_monte_carlo, u01 ): displacement[dim, droplet] = scheme( position_in_cell[dim, droplet], + cell_id[droplet], courant[_l] / n_substeps, courant[_r] / n_substeps, + enable_monte_carlo, + u01 ) @@ -28,13 +31,14 @@ class DisplacementMethods(BackendMethods): @numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False, "cache": False}}) # pylint: disable=too-many-arguments def calculate_displacement_body_1d( - dim, scheme, displacement, courant, cell_origin, position_in_cell, n_substeps + dim, scheme, displacement, courant, cell_origin, position_in_cell, cell_id, n_substeps, enable_monte_carlo, rng ): length = displacement.shape[1] for droplet in numba.prange(length): # pylint: disable=not-an-iterable # Arakawa-C grid _l = cell_origin[0, droplet] _r = cell_origin[0, droplet] + 1 + u01 = rng[droplet] calculate_displacement_body_common( dim, droplet, @@ -44,14 +48,17 @@ def calculate_displacement_body_1d( displacement, courant, position_in_cell, + cell_id, n_substeps, + enable_monte_carlo, + u01, ) @staticmethod @numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False, "cache": False}}) # pylint: disable=too-many-arguments def calculate_displacement_body_2d( - dim, scheme, displacement, courant, cell_origin, position_in_cell, n_substeps + dim, scheme, displacement, courant, cell_origin, position_in_cell, cell_id, n_substeps, enable_monte_carlo, rng ): length = displacement.shape[1] for droplet in numba.prange(length): # pylint: disable=not-an-iterable @@ -64,6 +71,7 @@ def calculate_displacement_body_2d( cell_origin[0, droplet] + 1 * (dim == 0), cell_origin[1, droplet] + 1 * (dim == 1), ) + u01 = rng[droplet] calculate_displacement_body_common( dim, droplet, @@ -73,14 +81,17 @@ def calculate_displacement_body_2d( displacement, courant, position_in_cell, + cell_id, n_substeps, + enable_monte_carlo, + u01, ) @staticmethod @numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False, "cache": False}}) # pylint: disable=too-many-arguments def calculate_displacement_body_3d( - dim, scheme, displacement, courant, cell_origin, position_in_cell, n_substeps + dim, scheme, displacement, courant, cell_origin, position_in_cell, cell_id, n_substeps, enable_monte_carlo, rng ): n_sd = displacement.shape[1] for droplet in numba.prange(n_sd): # pylint: disable=not-an-iterable @@ -95,6 +106,7 @@ def calculate_displacement_body_3d( cell_origin[1, droplet] + 1 * (dim == 1), cell_origin[2, droplet] + 1 * (dim == 2), ) + u01 = rng[droplet] calculate_displacement_body_common( dim, droplet, @@ -104,11 +116,14 @@ def calculate_displacement_body_3d( displacement, courant, position_in_cell, + cell_id, n_substeps, + enable_monte_carlo, + u01, ) def calculate_displacement( - self, *, dim, displacement, courant, cell_origin, position_in_cell, n_substeps + self, *, dim, displacement, courant, cell_origin, position_in_cell, cell_id, n_substeps, enable_monte_carlo, rng ): n_dims = len(courant.shape) scheme = self.formulae.particle_advection.displacement @@ -120,7 +135,10 @@ def calculate_displacement( courant.data, cell_origin.data, position_in_cell.data, + cell_id.data, n_substeps, + enable_monte_carlo, + rng.uniform(0., 1., displacement.data.shape[1]), ) elif n_dims == 2: DisplacementMethods.calculate_displacement_body_2d( @@ -130,7 +148,10 @@ def calculate_displacement( courant.data, cell_origin.data, position_in_cell.data, + cell_id.data, n_substeps, + enable_monte_carlo, + rng.uniform(0., 1., displacement.data.shape[1]), ) elif n_dims == 3: DisplacementMethods.calculate_displacement_body_3d( @@ -140,7 +161,10 @@ def calculate_displacement( courant.data, cell_origin.data, position_in_cell.data, + cell_id.data, n_substeps, + enable_monte_carlo, + rng.uniform(0., 1., displacement.data.shape[1]), ) else: raise NotImplementedError() @@ -173,7 +197,8 @@ def body( # and crossed precip-counting level position_within_column < precipitation_counting_level_index ): - rainfall_mass += abs(water_mass[idx[i]]) * multiplicity[idx[i]] + rainfall_mass += abs(water_mass[idx[i]]) * \ + multiplicity[idx[i]] idx[i] = flag healthy[0] = 0 return rainfall_mass diff --git a/PySDM/backends/impl_thrust_rtc/methods/displacement_methods.py b/PySDM/backends/impl_thrust_rtc/methods/displacement_methods.py index 7c6758abd..0d340ec5a 100644 --- a/PySDM/backends/impl_thrust_rtc/methods/displacement_methods.py +++ b/PySDM/backends/impl_thrust_rtc/methods/displacement_methods.py @@ -24,8 +24,11 @@ def __calculate_displacement_body(self): "courant_shape_0", "courant_shape_1", "cell_origin", + "cell_id", "position_in_cell", "n_substeps", + "enable_monte_carlo", + "u01", ), name_iter="i", body=f""" @@ -52,8 +55,11 @@ def __calculate_displacement_body(self): displacement[i + n_sd * dim] = { self.formulae.particle_advection.displacement.c_inline( position_in_cell="position_in_cell[i + n_sd * dim]", + cell_id="cell_id[i + n_sd * dim]", c_l="courant[_l] / n_substeps", - c_r="courant[_r] / n_substeps" + c_r="courant[_r] / n_substeps", + enable_monte_carlo="false", + u01="u01" ) }; """.replace( @@ -96,7 +102,7 @@ def __flag_precipitated_body(self): @nice_thrust(**NICE_THRUST_FLAGS) def calculate_displacement( - self, *, dim, displacement, courant, cell_origin, position_in_cell, n_substeps + self, *, dim, displacement, courant, cell_origin, position_in_cell, cell_id, n_substeps, enable_monte_carlo, rng ): n_dim = len(courant.shape) n_sd = position_in_cell.shape[1] @@ -110,8 +116,11 @@ def calculate_displacement( trtc.DVInt64(courant.shape[0]), trtc.DVInt64(courant.shape[1] if n_dim > 2 else -1), cell_origin.data, + cell_id.data, position_in_cell.data, trtc.DVInt64(n_substeps), + trtc.DVBool(enable_monte_carlo), + trtc.DVDouble(rng.uniform(0., 1.)) ), ) diff --git a/PySDM/backends/impl_thrust_rtc/test_helpers/cpp2python.py b/PySDM/backends/impl_thrust_rtc/test_helpers/cpp2python.py index 8d3e7623b..b520f4fcb 100644 --- a/PySDM/backends/impl_thrust_rtc/test_helpers/cpp2python.py +++ b/PySDM/backends/impl_thrust_rtc/test_helpers/cpp2python.py @@ -199,7 +199,8 @@ def make(self): arcsinh as asinh, sinh, maximum, minimum, - where, # TODO #1295 + where, # TODO #1295 + absolute, sign ) import numba diff --git a/PySDM/dynamics/displacement.py b/PySDM/dynamics/displacement.py index f3e409f48..747410e00 100644 --- a/PySDM/dynamics/displacement.py +++ b/PySDM/dynamics/displacement.py @@ -21,12 +21,14 @@ class Displacement: # pylint: disable=too-many-instance-attributes def __init__( self, enable_sedimentation=False, + enable_monte_carlo=False, precipitation_counting_level_index: int = 0, adaptive=DEFAULTS.adaptive, rtol=DEFAULTS.rtol, ): # pylint: disable=too-many-arguments self.particulator = None self.enable_sedimentation = enable_sedimentation + self.enable_monte_carlo = enable_monte_carlo self.dimension = None self.grid = None self.courant = None @@ -101,11 +103,12 @@ def __call__(self): # TIP: not need all array only [idx[:sd_num]] cell_origin = self.particulator.attributes["cell origin"] position_in_cell = self.particulator.attributes["position in cell"] + cell_id = self.particulator.attributes["cell id"] self.precipitation_mass_in_last_step = 0.0 for _ in range(self._n_substeps): self.calculate_displacement( - self.displacement, self.courant, cell_origin, position_in_cell + self.displacement, self.courant, cell_origin, position_in_cell, cell_id ) self.update_position(position_in_cell, self.displacement) if self.enable_sedimentation: @@ -122,16 +125,23 @@ def __call__(self): self.particulator.attributes.mark_updated(key) def calculate_displacement( - self, displacement, courant, cell_origin, position_in_cell + self, displacement, courant, cell_origin, position_in_cell, cell_id ): + if self.enable_sedimentation and self.enable_monte_carlo: + dt = self.particulator.dt / self._n_substeps + dt_over_dz = dt / self.particulator.mesh.dz + courant += self.particulator.attributes["relative fall velocity"] / dt_over_dz + self.particulator.calculate_displacement( displacement=displacement, courant=courant, cell_origin=cell_origin, + cell_id=cell_id, position_in_cell=position_in_cell, n_substeps=self._n_substeps, + enable_monte_carlo=self.enable_monte_carlo, ) - if self.enable_sedimentation: + if self.enable_sedimentation and not self.enable_monte_carlo: displacement_z = displacement[self.dimension - 1, :] dt = self.particulator.dt / self._n_substeps dt_over_dz = dt / self.particulator.mesh.dz diff --git a/PySDM/particulator.py b/PySDM/particulator.py index 32ae0697d..b9fb209be 100644 --- a/PySDM/particulator.py +++ b/PySDM/particulator.py @@ -427,7 +427,7 @@ def flag_out_of_column(self): self.attributes.sanitize() def calculate_displacement( - self, *, displacement, courant, cell_origin, position_in_cell, n_substeps + self, *, displacement, courant, cell_origin, cell_id, position_in_cell, n_substeps, enable_monte_carlo ): for dim in range(len(self.environment.mesh.grid)): self.backend.calculate_displacement( @@ -435,8 +435,11 @@ def calculate_displacement( displacement=displacement, courant=courant[dim], cell_origin=cell_origin, + cell_id=cell_id, position_in_cell=position_in_cell, n_substeps=n_substeps, + enable_monte_carlo=enable_monte_carlo, + rng=self.Random(1,1).generator, ) def isotopic_fractionation(self, heavy_isotopes: tuple): diff --git a/PySDM/physics/particle_advection/explicit_in_space.py b/PySDM/physics/particle_advection/explicit_in_space.py index 120e64cb0..7bec7eeaa 100644 --- a/PySDM/physics/particle_advection/explicit_in_space.py +++ b/PySDM/physics/particle_advection/explicit_in_space.py @@ -1,12 +1,34 @@ """ basic explicit-in-space Euler scheme """ - +import numpy as np class ExplicitInSpace: # pylint: disable=too-few-public-methods def __init__(self, _): pass @staticmethod - def displacement(_, position_in_cell, c_l, c_r): - return c_l * (1 - position_in_cell) + c_r * position_in_cell + def displacement(_, position_in_cell, cell_id, c_l, c_r, enable_monte_carlo, u01): + return ( + position_in_cell + + ( + np.floor( + np.abs(max(c_l, c_r)) + ) + * np.sign( + np.abs(max(c_l, c_r)) + / max(c_l, c_r) + ) + ) + + ( + np.abs(max(c_l, c_r)) > u01 + ) + * np.sign( + np.abs(max(c_l, c_r)) + / max(c_l, c_r) + ) + ) if enable_monte_carlo else ( + c_l + * (1 - position_in_cell) + + c_r * position_in_cell + ) diff --git a/PySDM/physics/particle_advection/implicit_in_space.py b/PySDM/physics/particle_advection/implicit_in_space.py index 92ccec0fd..945ea0f12 100644 --- a/PySDM/physics/particle_advection/implicit_in_space.py +++ b/PySDM/physics/particle_advection/implicit_in_space.py @@ -1,12 +1,42 @@ """ eqs. 14-16 in [Arabas et al. 2015](https://doi.org/10.5194/gmd-8-1677-2015) """ - +import numpy as np class ImplicitInSpace: # pylint: disable=too-few-public-methods def __init__(self, _): pass @staticmethod - def displacement(_, position_in_cell, c_l, c_r): - return (c_l * (1 - position_in_cell) + c_r * position_in_cell) / (1 - c_r + c_l) + def displacement(_, position_in_cell, cell_id, c_l, c_r, enable_monte_carlo, u01): + return ( + position_in_cell + + ( + np.floor( + np.abs(max(c_l, c_r)) + ) + * np.sign( + np.abs(max(c_l, c_r)) + / max(c_l, c_r) + ) + ) + + ( + np.abs(max(c_l, c_r)) > u01 + ) + * np.sign( + np.abs(max(c_l, c_r)) + / max(c_l, c_r) + ) + ) if enable_monte_carlo else ( + ( + c_l + * (1 - position_in_cell) + + c_r + * position_in_cell + ) + / ( + 1 + - c_r + + c_l + ) + ) diff --git a/tests/unit_tests/dynamics/displacement/displacement_settings.py b/tests/unit_tests/dynamics/displacement/displacement_settings.py index e16d481ac..2dfee32e2 100644 --- a/tests/unit_tests/dynamics/displacement/displacement_settings.py +++ b/tests/unit_tests/dynamics/displacement/displacement_settings.py @@ -21,9 +21,10 @@ def __init__(self, n_sd=1, grid=None, positions=None, courant_field_data=None): self.sedimentation = False self.dt = None - def get_displacement(self, backend, scheme, adaptive=True): + def get_displacement(self, backend, scheme, adaptive=True, enable_monte_carlo=False): formulae = Formulae(particle_advection=scheme) - particulator = DummyParticulator(backend, n_sd=len(self.n), formulae=formulae) + particulator = DummyParticulator( + backend, n_sd=len(self.n), formulae=formulae) particulator.environment = DummyEnvironment( timestep=self.dt, grid=self.grid, courant_field_data=self.courant_field_data ) @@ -39,7 +40,8 @@ def get_displacement(self, backend, scheme, adaptive=True): "position in cell": position_in_cell, } particulator.build(attributes) - sut = Displacement(enable_sedimentation=self.sedimentation, adaptive=adaptive) + sut = Displacement(enable_sedimentation=self.sedimentation, + adaptive=adaptive, enable_monte_carlo=enable_monte_carlo) sut.register(particulator) sut.upload_courant_field(self.courant_field_data) diff --git a/tests/unit_tests/dynamics/displacement/test_advection.py b/tests/unit_tests/dynamics/displacement/test_advection.py index 14ff21caf..c62e85426 100644 --- a/tests/unit_tests/dynamics/displacement/test_advection.py +++ b/tests/unit_tests/dynamics/displacement/test_advection.py @@ -97,15 +97,16 @@ def test_calculate_displacement(backend_class): ) settings.positions = [[weight], [0]] sut, particulator = settings.get_displacement( - backend_class, scheme="ExplicitInSpace", adaptive=False + backend_class, scheme="ExplicitInSpace", adaptive=False, enable_monte_carlo=False ) # Act sut.calculate_displacement( - sut.displacement, - sut.courant, - particulator.attributes["cell origin"], - particulator.attributes["position in cell"], + displacement=sut.displacement, + courant=sut.courant, + cell_origin=particulator.attributes["cell origin"], + position_in_cell=particulator.attributes["position in cell"], + cell_id=particulator.attributes["cell id"], ) # Assert @@ -127,15 +128,16 @@ def test_calculate_displacement_dim1(backend_class): ) settings.positions = [[0], [weight]] sut, particulator = settings.get_displacement( - backend_class, scheme="ExplicitInSpace", adaptive=False + backend_class, scheme="ExplicitInSpace", adaptive=False, enable_monte_carlo=False ) # Act sut.calculate_displacement( - sut.displacement, - sut.courant, - particulator.attributes["cell origin"], - particulator.attributes["position in cell"], + displacement=sut.displacement, + courant=sut.courant, + cell_origin=particulator.attributes["cell origin"], + position_in_cell=particulator.attributes["position in cell"], + cell_id=particulator.attributes["cell id"], ) # Assert diff --git a/tests/unit_tests/dynamics/displacement/test_monte_carlo.py b/tests/unit_tests/dynamics/displacement/test_monte_carlo.py new file mode 100644 index 000000000..9b80e754d --- /dev/null +++ b/tests/unit_tests/dynamics/displacement/test_monte_carlo.py @@ -0,0 +1,326 @@ +# pylint: disable=missing-module-docstring,missing-class-docstring,missing-function-docstring +import numpy as np +import pytest + +from .displacement_settings import DisplacementSettings + + +class TestExplicitEulerWithInterpolation: + @staticmethod + @pytest.mark.parametrize( + "positions, courant_field, expected_positions, scheme", + ( + # ImplicitInSpace + # 1D + ( + [[0.5]], + (np.array([0.1, 0.2]),), + # (np.array([0.99, 0.99]),), + [[0.65827668]], + "ImplicitInSpace", + ), + # 2D + ( + [[0.5], [0.5]], + ( + np.array([0.1, 0.2]).reshape((2, 1)), + np.array([0.3, 0.4]).reshape((1, 2)), + ), + [[0.65827668], [0.86931224]], + "ImplicitInSpace", + ), + # 3D + ( + [[0.5], [0.5], [0.5]], + ( + np.array([0.1, 0.2]).reshape((2, 1, 1)), + np.array([0.3, 0.4]).reshape((1, 2, 1)), + np.array([0.2, 0.3]).reshape((1, 1, 2)), + ), + [ + [0.65827668], + [0.86931224], + [0.76379446], + ], + "ImplicitInSpace", + ), + # ExplicitInSpace + # 1D + ( + [[0.5]], + (np.array([0.1, 0.2]),), + # (np.array([0.99, 0.99]),), + [[0.657241]], + "ExplicitInSpace", + ), + # 2D + ( + [[0.5], [0.5]], + ( + np.array([0.1, 0.2]).reshape((2, 1)), + np.array([0.3, 0.4]).reshape((1, 2)), + ), + [[0.657241], [0.866895]], + "ExplicitInSpace", + ), + # 3D + ( + [[0.5], [0.5], [0.5]], + ( + np.array([0.1, 0.2]).reshape((2, 1, 1)), + np.array([0.3, 0.4]).reshape((1, 2, 1)), + np.array([0.2, 0.3]).reshape((1, 1, 2)), + ), + [ + [0.657241], + [0.866895], + [0.762068], + ], + "ExplicitInSpace", + ), + ), + ) + def test_single_cell( + backend_class, positions, expected_positions, scheme, courant_field: tuple + ): + # Arrange + settings = DisplacementSettings( + courant_field_data=courant_field, + positions=positions, + grid=tuple([1] * len(courant_field)), + n_sd=len(positions[0]), + ) + sut, particulator = settings.get_displacement( + backend_class, scheme=scheme + ) + + # Act + sut() + + # Assert + np.testing.assert_array_almost_equal( + np.asarray(expected_positions), + particulator.attributes["position in cell"].to_ndarray(), + ) + + @staticmethod + def test_advection(backend_class): + # Arrange + settings = DisplacementSettings() + settings.grid = (3, 3) + settings.courant_field_data = (np.ones((4, 3)), np.zeros((3, 4))) + settings.positions = [[1.5], [1.5]] + sut, particulator = settings.get_displacement( + backend_class, scheme="ImplicitInSpace" + ) + + # Act + sut() + + # Assert + dim_x = 0 + np.testing.assert_array_equal( + particulator.attributes["cell origin"][:, dim_x], np.array([2, 1]) + ) + + @staticmethod + def test_calculate_displacement(backend_class): + # Arrange + settings = DisplacementSettings() + value_a = 0.1 + value_b = 0.2 + weight = 0.125 + settings.courant_field_data = ( + np.array([[value_a, value_b]]).T, + np.array([[0, 0]]), + ) + settings.positions = [[weight], [0]] + sut, particulator = settings.get_displacement( + backend_class, scheme="ExplicitInSpace", adaptive=False, enable_monte_carlo=True + ) + # Act + sut.calculate_displacement( + displacement=sut.displacement, + courant=sut.courant, + cell_origin=particulator.attributes["cell origin"], + position_in_cell=particulator.attributes["position in cell"], + cell_id=particulator.attributes["cell id"], + ) + + # Assert + if backend_class.__name__ == "ThrustRTC": + np.testing.assert_equal( + sut.displacement[0, slice(0, 1)].to_ndarray(), + 0.1125 + ) + + if backend_class.__name__ == "Numba": + np.testing.assert_equal( + sut.displacement[0, slice(0, 1)].to_ndarray(), + 0.125 + ) + + @staticmethod + def test_calculate_displacement_2(backend_class): + # Arrange + settings = DisplacementSettings() + value_a = 1 + value_b = 1 + weight = 0.125 + settings.courant_field_data = ( + np.array([[value_a, value_b]]).T, + np.array([[0, 0]]), + ) + settings.positions = [[weight], [0]] + sut, particulator = settings.get_displacement( + backend_class, scheme="ExplicitInSpace", adaptive=False, enable_monte_carlo=True + ) + # Act + sut.calculate_displacement( + displacement=sut.displacement, + courant=sut.courant, + cell_origin=particulator.attributes["cell origin"], + position_in_cell=particulator.attributes["position in cell"], + cell_id=particulator.attributes["cell id"], + ) + + # Assert + if backend_class.__name__ == "ThrustRTC": + np.testing.assert_equal( + sut.displacement[0, slice(0, 1)].to_ndarray(), + 1. + ) + + if backend_class.__name__ == "Numba": + np.testing.assert_equal( + sut.displacement[0, slice(0, 1)].to_ndarray(), + 2.125 + ) + + @staticmethod + def test_calculate_displacement_dim1(backend_class): + # Arrange + settings = DisplacementSettings() + value_a = 0.1 + value_b = 0.2 + weight = 0.125 + settings.courant_field_data = ( + np.array([[0, 0]]).T, + np.array([[value_a, value_b]]), + ) + settings.positions = [[0], [weight]] + sut, particulator = settings.get_displacement( + backend_class, scheme="ExplicitInSpace", adaptive=False, enable_monte_carlo=True + ) + + # Act + sut.calculate_displacement( + displacement=sut.displacement, + courant=sut.courant, + cell_origin=particulator.attributes["cell origin"], + position_in_cell=particulator.attributes["position in cell"], + cell_id=particulator.attributes["cell id"], + ) + + # Assert + # np.testing.assert_equal( + # sut.displacement[1, slice(0, 1)].to_ndarray(), + # 0.1125 + # ) + + if backend_class.__name__ == "ThrustRTC": + np.testing.assert_equal( + sut.displacement[1, slice(0, 1)].to_ndarray(), + 0.1125 + ) + + if backend_class.__name__ == "Numba": + np.testing.assert_equal( + sut.displacement[1, slice(0, 1)].to_ndarray(), + 0.125 + ) + + @staticmethod + def test_update_position(backend_class): + # Arrange + settings = DisplacementSettings() + p_x = 0.1 + p_y = 0.2 + settings.positions = [[p_x], [p_y]] + sut, particulator = settings.get_displacement( + backend_class, scheme="ImplicitInSpace" + ) + + droplet_id = slice(0, 1) + sut.displacement[:] = particulator.backend.Storage.from_ndarray( + np.asarray( + [ + [ + 0.1, + ], + [0.2], + ] + ) + ) + + # Act + sut.update_position( + particulator.attributes["position in cell"], sut.displacement + ) + + # Assert + for d in range(2): + np.testing.assert_array_almost_equal( + particulator.attributes["position in cell"][d, droplet_id].to_ndarray(), + settings.positions[d][droplet_id] + + sut.displacement[d, droplet_id].to_ndarray(), + ) + + @staticmethod + def test_update_cell_origin(backend_class): + # Arrange + settings = DisplacementSettings() + sut, particulator = settings.get_displacement( + backend_class, scheme="ImplicitInSpace" + ) + + droplet_id = 0 + state = particulator.attributes + state["position in cell"][:] = particulator.backend.Storage.from_ndarray( + np.asarray([[1.1], [1.2]]) + ) + + # Act + sut.update_cell_origin(state["cell origin"], state["position in cell"]) + + # Assert + for d in range(2): + assert ( + state["cell origin"][d, droplet_id] + == settings.positions[d][droplet_id] + 1 + ) + assert state["position in cell"][d, droplet_id] == ( + state["position in cell"][d, droplet_id] + - np.floor(state["position in cell"][d, droplet_id]) + ) + + @staticmethod + def test_boundary_condition(backend_class): + # Arrange + settings = DisplacementSettings() + sut, particulator = settings.get_displacement( + backend_class, scheme="ImplicitInSpace" + ) + + droplet_id = 0 + state = particulator.attributes + state["cell origin"][:] = particulator.backend.Storage.from_ndarray( + np.asarray([[1], [1]]) + ) + + # Act + sut.boundary_condition(state["cell origin"]) + + # Assert + assert state["cell origin"][0, droplet_id] == 0 + assert state["cell origin"][1, droplet_id] == 0