rough.py 2.04 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
import slimPre
import slim
import numpy
import scipy.interpolate

#########################
### global parameters ###
#########################

meshfile = "balzano.msh"
tend = 10*12*3600

dt = 900

def ssebnd(t) :
    h = 21-t/(3600)
    return max(h,3)

N=10
pts = numpy.ndarray((N+4,2))
vals = numpy.zeros((N+4,))
L = 13800
l = 3600
pts[:4,:] = [[0,-l],[0,l],[L,-l],[L,l]]
vals[0:2] = -20
vals[2:4]=0
numpy.random.seed(0)
pts[4:,0] = numpy.random.rand(N)*0.79*L
pts[4:,1] = numpy.random.rand(N)*2*l-l
vals[4:] = numpy.random.rand(N)*-20+1
interp = scipy.interpolate.LinearNDInterpolator(pts,vals,0)
def bathf(x,y) :
    return interp(x,y)

hlimt = 1
hlimd = 0.5
linDrag = 0.01

###################
### pre process ###
###################

### mesh ###
mesh = slimPre.Mesh(meshfile)
region = slimPre.Region(mesh)
xyz = region.coordinates

### bathymetry ###
bath = numpy.zeros((xyz.shape[0]))
for i in range(xyz.shape[0]) :
    bath[i] = bathf(xyz[i,0], xyz[i,1])
slimPre.write_file('bath.nc', region=region, time=None, data=[('bath', bath)])
slimPre.netcdf_to_msh(meshfile, "bath.nc", "bath", "output/bath")

### open boundary ###
time_vector = numpy.arange(0,tend+2*dt,dt,numpy.float64)
time = slimPre.Time(time_vector)
eta = numpy.zeros_like(time_vector)
for i in range(eta.shape[0]) :
    eta[i] = ssebnd(time_vector[i])
slimPre.write_file('eta.nc', region=None, time=time, data=[('eta', eta)])

### initial condition ###
zeros = numpy.zeros_like(xyz)
etaInit = numpy.ones((xyz.shape[0],1))*ssebnd(0)
etaInit = numpy.maximum(-bath,etaInit)
slimPre.write_file('init.nc', region=region, time=None, data=[('eta', etaInit),('u', zeros),('v', zeros)])

###########
### run ###
###########

domain = slim.Domain(meshfile,('bath.nc','bath'))
eq = slim.ShallowWater2dWD(domain, "implicit", final_time=tend, hlimt=hlimt, hlimd=hlimd)
eq.set_initial_condition(('init.nc','eta'),('init.nc','u'),('init.nc','v'))    
eq.set_boundary_open('OpenSea', sse=('eta.nc','eta'))
eq.set_boundary_coast('Wall')
loop = slim.Loop(maximum_time_step=dt, export_time=3600)
loop.add_equation(eq)
loop.run()