In [1]:
import triflow as trf
import numpy as np
import pylab as pl
import scipy.signal
%matplotlib inline
In [2]:
# http://geoexamples.blogspot.fr/2014/03/shaded-relief-images-using-gdal-python.html
def hillshade(array, azimuth, angle_altitude):
x, y = np.gradient(array)
slope = np.pi/2. - np.arctan(np.sqrt(x*x + y*y))
aspect = np.arctan2(-x, y)
azimuthrad = azimuth*np.pi / 180.
altituderad = angle_altitude*np.pi / 180.
shaded = np.sin(altituderad) * np.sin(slope)\
+ np.cos(altituderad) * np.cos(slope)\
* np.cos(azimuthrad - aspect)
return 255*(shaded + 1)/2
Triflow is not able to handle 2D simution by itself. We have to help and provide our own routine. It will be slower than Triflow compiler, but enough for our computation. We have to provide:
In [3]:
class NonConservative_ShallowWater:
@staticmethod
def F(fields, pars):
Ffields = fields.copy()
x = fields["x"].values
h = fields["h"].values
u = fields["u"].values
v = fields["v"].values
H = fields["H"].values
delta_x = x.ptp() / (x.size - 1)
delta_y = y.ptp() / (y.size - 1)
def dx(U):
return (np.roll(U, -1, axis=0) - np.roll(U, 1, axis=0)) / (2 * delta_x)
def dy(U):
return (np.roll(U, -1, axis=1) - np.roll(U, 1, axis=1)) / (2 * delta_y)
def dxx(U):
return (np.roll(U, 1, axis=0) - 2 * U + np.roll(U, -1, axis=0)) / (delta_x**2)
def dyy(U):
return (np.roll(U, 1, axis=1) - 2 * U + np.roll(U, -1, axis=1)) / (delta_y**2)
eta = h + H
visc = lambda var: pars["nu"] * (dxx(var) + dyy(var))
dth = -(dx(u * eta) + dy(v * eta))
dtu = -(u * dx(u) + v * dy(u)) + pars["f"] * v - 9.81 * dx(h) + visc(u)
dtv = -(u * dx(v) + v * dy(v)) - pars["f"] * u - 9.81 * dy(h) + visc(v)
Ffields["h"][:] = dth
Ffields["u"][:] = dtu
Ffields["v"][:] = dtv
return Ffields.uflat
_indep_vars = ["x", "y"]
fields_template = trf.core.fields.BaseFields.factory(("x", "y"),
[("h", ("x", "y")),
("u", ("x", "y")),
("v", ("x", "y"))],
[("H", ("x", "y"))])
In [4]:
model = NonConservative_ShallowWater()
In [5]:
Nx = Ny = 350
x = np.linspace(0, 5, Nx)
y = np.linspace(0, 5, Ny)
u = np.zeros((Nx, Ny)) + 2
v = np.zeros((Nx, Ny)) + 2
h = np.zeros((Nx, Ny)) + (scipy.signal.gaussian(Nx, 5) * scipy.signal.gaussian(Ny, 5)[:, None]) * 1 # It's like a drop. Or close enough
H = np.ones((Nx, Ny)) * 10
In [6]:
pl.figure(figsize=(7, 7))
pl.contourf(hillshade(h + H, 45, 50), 200, cmap='Greys')
pl.grid(False)
pl.axis('off')
Out[6]:
In [7]:
init_fields = model.fields_template(x=x, y=y, h=h, u=u, v=v, H=H)
pars = {"f": 0, "nu": 1E-6} # no coriolis effect, water viscosity
In [8]:
simul = trf.Simulation(model, init_fields.copy(), pars, dt=.01, tmax=.4,
scheme=trf.schemes.scipy_ode, interpolator="dopri5",
id="wave_animation")
simul.attach_container(None)
Out[8]:
In [9]:
t, fields = simul.run()
In [19]:
pl.figure(figsize=(7, 7))
eta = fields["h"] + fields["H"]
norm_eta = ((eta - eta.min()) / (eta.max() - eta.min()))
pl.contourf(hillshade((eta - eta.min()) / (eta.max() - eta.min()), 45, 20), 200, cmap='Greys')
pl.grid(False)
pl.axis('off')
Out[19]:
In [ ]: