# Source code for opendrift.models.eulerdrift.simulation

import logging
logger = logging.getLogger(__name__)
from .grid import RegularGrid
from . import srs
from .interp import vec_nearest
from abc import abstractmethod
import numpy as np
from datetime import timedelta, datetime

[docs]
class Simulation:
r"""
The simulation.

It contains the problem to be simulated, means to read necessary input
variables, and the physics for modeling the convection of the initial
conditions.

Convection:

Convection consists of advection and diffusion.

Diffusion is given by:

.. math::

\frac{\partial U}{\partial t} = D \left( \frac{\partial^2 U}{\partial x^2} + \frac{\partial^2 U}{\partial y^2} \right)

The convection equation is (wiki
<https://en.wikipedia.org/wiki/Convection%E2%80%93diffusion_equation>_):

.. math::

\frac{\partial c}{\partial t} = ...

with the assumptions that:

* the diffusion constant D is constant for the field,
* and that the flow u is incompressible (i.e. has *no divergence*).

the equation simplifies to:

.. math::

\frac{\partial c}{\partial t} = D \nabla^2 c - \mathbf{v} \cdot \nabla T

where :math:\nabla^2 = \triangle is the Laplacian.
"""

"""The model grid."""
grid = None

"""Reference time datetime.Datetime"""
t0 = None

"""Current time offset after t0"""
t = 0.

"""
Diffusivity (:math:m^2/s). E.g. between 0.01 and 0.1 for oil on the
surface of the ocean (Matsuzakia et. al., 2017
<https://www.sciencedirect.com/science/article/pii/S0025326X16308426>_).

Decreasing diffusivity places stricter stability criteria on time step.
"""
D = 0.1

"""Porosity, rate of liquid volume to total volume (fraction of flux)"""
rho = 1.

def __init__(self, grid):
self.grid = grid
self.t0 = datetime.utcnow()

[docs]
@classmethod
def new(cls, lon0 = 10., lat0 = 65., res = 10., shape=(100, 100)):
"""
New simulation on a :class:.grid.RegularGrid.

Args:

lon0: lower-left corner longitude

lat0: lower-left corner latitude

res: resolution (dx and dy)

shape: shape (size) of grid
"""
s = cls(grid=RegularGrid.new(lon0, lat0, res, shape))
logger.info("New %s with Regular Grid" % (cls.__name__))

return s

[docs]
def source(self, lon, lat, X):
"""
Source X onto grid with lower-left corner lon, lat.
"""
x0, y0 = self.grid.srs(lon, lat)
assert self.grid.contains(x0, y0)

ix0 = vec_nearest(self.grid.x, x0)
iy0 = vec_nearest(self.grid.y, y0)
print(ix0)
ix1 = ix0 + X.shape
iy1 = iy0 + X.shape

assert ix1 < self.grid.grid.shape
assert iy1 < self.grid.grid.shape

self.grid.grid[ix0:ix1, iy0:iy1] = X

[docs]
def source_gaussian_blob(self, lon, lat, A=1., N=10, sigma=10.):
r"""
Source a Gaussian blob (2D normal distribution) at lon and lat with
sigma (standard deviation, meters) radius.

Args:

lon, lat: Center coordinates, or :math:\bar \mu.

A: Amplitude.

N: Kernel size.

sigma: standard deviation (:math:\sigma) in meters.
"""
def gaussian2d(n, std):
"""
Gaussian 2D window with size n and standard deviation std.
"""

from scipy import signal
gk = signal.windows.gaussian(n, std=std)
return np.outer(gk, gk)

# center of kernel
x, y = self.grid.srs(lon, lat)
assert self.grid.contains(x, y)

ix0 = vec_nearest(self.grid.x, x) - N // 2
iy0 = vec_nearest(self.grid.y, y) - N // 2
ix1 = ix0 + N
iy1 = iy0 + N

assert ix0 > 0 and iy0 > 0

S = gaussian2d(N, sigma / self.grid.res)

self.grid.grid[ix0:ix1, iy0:iy1] = S

[docs]
def U(self, t):
"""
Get U (ocean current) for t0 + t
"""
Ux, Uy = None, None
t = self.t0 + timedelta(seconds = t)

vv = r.variables()
if 'x_sea_water_velocity' in vv and 'y_sea_water_velocity' in vv:
Ux, Uy = r.read_grid(self.grid, ['x_sea_water_velocity', 'y_sea_water_velocity'], t)

if Ux is not None and Uy is not None:
break

return Ux, Uy

[docs]
@abstractmethod
def step(self, dt=None):
"""
Step the simulation.

Stepping the simulation involves applying diffusion and advection to the field.

Args:

dt: time delta (or use automatic).
"""
pass

[docs]
@abstractmethod
def integrate(self, dt=None, max_t=None, max_steps=None, observer = None):
"""
Run simulation until termination condition is met.

Args:

dt: override time step

observer: function to call after each step taking simulation object as first argument. the function may return False to stop the integration.

Termination conditions:

max_t: max time
max_steps: max iterations
"""

if max_t is not None and max_steps is not None:
logger.warning(
"no termination condition supplied, using max_steps = 1000")
max_steps = 1000

logger.info("integrating.. (dt = %s)", dt)

n = 0

while True:
logger.debug("step: n = %d, t = %s (dt = %s)" % (n, self.t, dt))

if max_t is not None and self.t >= max_t:
logger.info("max_t reached.")
break

if max_steps is not None and n >= max_steps:
logger.info("max_steps reached.")
break

self.step(dt=dt)

n += 1

[docs]
class ExplSimulation(Simulation):
r"""
A simple explicit scheme for integrating the convection-equation.

* Forward difference in time
* ndimage.laplace and np.gradient for spatial differences.

https://en.wikipedia.org/wiki/Numerical_solution_of_the_convection%E2%80%93diffusion_equation#Solving_the_convection%E2%80%93diffusion_equation_using_the_finite_difference_method

.. seealso::

:class:Simulation.

"""

[docs]
def stability(self, dx, D, umax):
"""
https://en.wikipedia.org/wiki/Numerical_solution_of_the_convection%E2%80%93diffusion_equation#Solving_the_convection%E2%80%93diffusion_equation_using_the_finite_difference_method
"""
h = 2 * D / (self.rho * umax)
dt = dx**2 / (2 * D)

return h, dt

[docs]
def step(self, dt=None):
from scipy import ndimage

dx, dy = self.grid.res, self.grid.res
D = self.D
# Ux = np.sqrt(50 * 1.) * np.ones(self.grid.grid.shape)
# Uy = Ux
Ux, Uy = self.U(self.t)
maxu = np.max(np.sqrt(Ux**2 + Uy**2).ravel())
logger.debug("maxu = %s" % maxu)

h, ddt = self.stability(dx, D, maxu)

if h < dx:
logger.warning("dx too big, dx = %s > h = %s" % (dx, h))

if dt is None:
dt = ddt
elif ddt < dt:
logger.warning("dt too big, dt = %s > ddt = %s" % (dt, ddt))

# diffusion
diff = D * 1. / dx**2 * ndimage.laplace(self.grid.grid)