prepro.py 4.82 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
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
import slimPre
import numpy as np

data_dir = 'data/'
slimPre.make_directory(data_dir)


mesh_file = slimPre.fetch_ftp('slim_data/tanganyika/tanga1dx.msh')

print('Extruding mesh')
def shiftOperation(node, iPerBound) :
  n = [node[0], node[1] - 1e4, node[2]]
  return n
cutTags = ["cut"]
pasteTags = ["paste"]
mapFilename = data_dir+"periodicMesh.txt"
periodicity = (shiftOperation, cutTags, pasteTags, mapFilename)

nPart = 2
if nPart > 1:
    slimPre.dgpy.dgMeshPartition(mesh_file, nPart)
    mesh_file = mesh_file[:-4] + '_' + str(nPart) + '.msh'

mesh2d = slimPre.Mesh(mesh_file)
region_global = slimPre.Region(mesh2d)
xyz = region_global.coordinates

bath_file_name = slimPre.fetch_ftp('slim_data/tanganyika/1d_bath.txt')
f = open(bath_file_name, 'r');
l = f.readline()
data = []
while l != '':
    (x1d, x, y, h) = [float(d) for d in l.split(',')]
    data.append([x1d, x, y, h])
    l = f.readline()
f.close()
data = np.array(data)

bath = np.zeros((xyz.shape[0]))
for i in range(bath.shape[0]):
    x = xyz[i,0]
    ind = np.where(x < data[:,0])[0][0]
    if ind == 0:
        print('bug')
        exit(-1)
    xsi = (x-data[ind-1,0]) / (data[ind,0]-data[ind-1,0])
    bath[i] = (1-xsi) * data[ind-1,3] + xsi * data[ind,3]

slimPre.write_file(data_dir+'bath_2d.nc', region=region_global, time=None, data=[('bath',bath)])
slimPre.netcdf_to_msh(mesh_file, data_dir+'bath_2d.nc','bath',data_dir+'bath')

zLayers  = [0, 2, 5]
zLayers += [10 * h for h in range(1, 11)]
zLayers += [150, 200, 300, 500, 700, 950, 1200, 1400, 1500]
mesh3d_file = data_dir+'tanga2dxz.msh'
slimPre.extrude(mesh_file, (data_dir+'bath_2d.nc','bath'), z_layers=zLayers, mesh_file_name_out=mesh3d_file, factor_show=200, periodicity=periodicity)

print('Loading 3D mesh')
mesh = slimPre.Mesh(mesh3d_file)
region_global = slimPre.Region(mesh)

print('Preprocessing Coriolis')
Tday = 0.99726968*24*60*60; # sidereal time of Earth revolution
OmegaEarth = 2*np.pi/Tday
latDeg = -6.2
phi = (np.pi/180)*latDeg # latitude in radians
corio = 2*OmegaEarth*np.sin(phi) # // [rad/s] Coriolis parameter
slimPre.write_file(data_dir+'coriolis.nc',None,None,[('coriolis',corio)])

print('Fixed viscosity')
import dgpy
shear = np.array([5, 5, 5, 1])/1e3
groups = mesh._groups
data = dgpy.fullMatrixDouble()
dof = dgpy.dgDofContainer(groups, 1)
dof_copy = dgpy.dgDofContainer(groups, 1)
dof.setAll(0)
for iFaceGroup in range(groups.getNbFaceGroups()):
    faceGroup = groups.getFaceGroup(iFaceGroup)
    if faceGroup.physicalTag() == 'vertical_bottom' or faceGroup.physicalTag() == 'bottom_sea':
        iGroup = groups.getElementGroupId(faceGroup.elementGroup(0))
        for iFace in range(faceGroup.size()):
            iElem = faceGroup.elementId(iFace, 0)
            dof.getElementProxy(iGroup, iElem, data)
            data.setAll(shear[0])
data0 = dgpy.fullMatrixDouble()
data1 = dgpy.fullMatrixDouble()
for i in range(len(shear)-1):
    dof_copy.copy(dof)
    for iFaceGroup in range(groups.getNbFaceGroups()):
        faceGroup = groups.getFaceGroup(iFaceGroup)
        if faceGroup.nConnection() == 2:
            iGroup0 = groups.getElementGroupId(faceGroup.elementGroup(0))
            iGroup1 = groups.getElementGroupId(faceGroup.elementGroup(1))
            for iFace in range(faceGroup.size()):
                iElem0 = faceGroup.elementId(iFace, 0)
                iElem1 = faceGroup.elementId(iFace, 1)
                dof_copy.getElementProxy(iGroup0, iElem0, data0)
                dof.getElementProxy(iGroup1, iElem1, data1)
                if data0.get(0,0) > 1e-10 and data1.get(0,0) < 1e-10:
                    data1.setAll(shear[i+1])
                dof_copy.getElementProxy(iGroup1, iElem1, data0)
                dof.getElementProxy(iGroup0, iElem0, data1)
                if data0.get(0,0) > 1e-10 and data1.get(0,0) < 1e-10:
                    data1.setAll(shear[i+1])
idxExp = dgpy.dgIdxExporter(dof, data_dir+'additional_shear')
idxExp.exportIdx()

print('Preprocessing Initial Temperature')
z = region_global.coordinates[:,2]
T0 = 26.5
Z0 = -40.
T1 = 24.3
Z1 = -100
temp = np.where(z > Z0, T0,
       np.where(z > Z1, T0 + (z-Z0) * ( (T1-T0)/ (Z1-Z0) ),
                        T1) ) 
slimPre.write_file(data_dir+'initial_temperature.nc', region=region_global, time=None, data=[('T_init', temp)])

print('Preprocessing Initial Elevation')
x = region_global.coordinates[:,0]
eta = np.maximum(0,(1-(x-2e5)*(x-2e5)/1e10))
slimPre.write_file(data_dir+'initial_eta.nc', region=region_global, time=None, data=[('eta_init', eta)])
slimPre.netcdf_to_msh(mesh3d_file, data_dir+'initial_eta.nc', 'eta_init', 'eta')


print('Preprocessing wind')
time = slimPre.Time(time_vector=[900*i for i in range(10000)])
uwind = 0.2 * np.sin(time._time * 2 * np.pi / (4*86400))
vwind = np.zeros(time._time.shape)
slimPre.write_file(data_dir+'windStress.nc',region=None,time=time,data=[('u',uwind), ('v',vwind)])

print( 'Preprocessing done ' )
slimPre.exit(0)