Commit 6c0a7085 authored by Philippe Delandmeter's avatar Philippe Delandmeter

slim3d: Tanganyika 3D

parent 50e3ed1a
Pipeline #2410 failed with stage
in 1 minute and 40 seconds
import slimPre
import numpy as np
pre_data_dir_base = 'data_'+slimPre.partition_nb()+'/'
pre_data_dir_base = 'data_051_'+slimPre.partition_nb()+'/'
pre_data_dir = pre_data_dir_base+slimPre.partition_id()+'/'
slimPre.make_directory(pre_data_dir)
mesh_file_name = slimPre.fetch_ftp('slim_data/tanganyika/tanganyika_v6.msh')
mesh_file_name = slimPre.fetch_ftp('slim_data/tanganyika/tanganyika_v6.msh', pre_data_dir_base)
initial_time = '2007-01-01 00:00:00'
final_time = '2008-02-01 00:00:00'
mesh_file_name = slimPre.partition_mesh(mesh_file_name, pre_data_dir_base+'mesh.msh')
mesh = slimPre.Mesh(mesh_file_name, mesh_proj='+proj=utm +ellps=WGS84 +zone=35 +south')
......@@ -16,11 +14,15 @@ region_global = slimPre.Region(mesh, '')
lonlat_global = slimPre.Coordinate_system(region_global, data_proj='+proj=latlong +ellps=WGS84')
### Extrusion ###
bath_file_name = slimPre.fetch_ftp('slim_data/tanganyika/bath_v6smooth_COMP_0.msh')
print('Extruding mesh')
zLayers = [0, 2, 5]
zLayers += [10 * h for h in range(1, 11)]
zLayers += [150, 200, 300, 500, 700, 950, 1200, 1400, 1500]
slimPre.make_directory(pre_data_dir_base+'bath_v6smooth/')
slimPre.fetch_ftp('slim_data/tanganyika/bath_v6smooth/bath_v6smooth_COMP_0.msh', pre_data_dir_base+'bath_v6smooth/')
bath_file_name = slimPre.fetch_ftp('slim_data/tanganyika/bath_v6smooth/bath_v6smooth.idx', pre_data_dir_base+'bath_v6smooth/')
slimPre.extrude(mesh_file_name, bath_file_name, z_layers=zLayers, mesh_file_name_out=pre_data_dir_base+'mesh3d.msh', factor_show=200)
mesh_file_name = pre_data_dir_base + "mesh3d.msh"
......@@ -72,22 +74,268 @@ for i in range(len(shear)-1):
idxExp = dgpy.dgIdxExporter(dof, pre_data_dir_base+'additional_shear')
idxExp.exportIdx()
print('Preprocessing Wind')
from netCDF4 import Dataset
import numpy as np
from scipy import ma
iFile = slimPre.fetch_ftp('slim_data/tanganyika/FLake_2000-2005.nc', pre_data_dir_base)
iF = Dataset(iFile, 'r')
lonVarIn = iF.variables['lon'][:]
latVarIn = iF.variables['lat'][:]
timeVarIn = iF.variables['time']
uVarIn = iF.variables['U_10M'][:]
vVarIn = iF.variables['V_10M'][:]
minLon = np.amin(lonVarIn)
maxLon = np.amax(lonVarIn)
minLat = np.amin(latVarIn)
maxLat = np.amax(latVarIn)
t = np.size(uVarIn,0)
print('min/max lon %g/%g' % (minLon, maxLon) )
print('min/max lat %g/%g' % (minLat, maxLat) )
def interpolate(lonloc, latloc):
i = (latloc < latVarIn[:,0]).argmax() -1
j = (lonloc < lonVarIn[i,:]).argmax() -1
i = (latloc < latVarIn[:,j]).argmax() -1
j = (lonloc < lonVarIn[i,:]).argmax() -1
x0 = 0.5 * ( lonVarIn[i,j] + lonVarIn[i+1,j] )
x1 = 0.5 * ( lonVarIn[i,j+1] + lonVarIn[i+1,j+1] )
y0 = 0.5 * ( latVarIn[i,j] + latVarIn[i,j+1] )
y1 = 0.5 * ( latVarIn[i+1,j] + latVarIn[i+1,j+1] )
xsi = (lonloc - x0) / (x1 - x0)
eta = (latloc - y0) / (y1 - y0)
if i < 0: xsi = 0
if j < 0: eta = 0
#if xsi < 0 or xsi > 1:
# print('warning xsi !! (i,j) = (%d,%d) xsi = %g' % (i,j,xsi))
#if eta < 0 or eta > 1:
# print('warning eta !! (i,j) = (%d,%d) eta = %g' % (i,j,eta))
xsi = min(1,max(0,xsi))
eta = min(1,max(0,eta))
uloc = ( (1-xsi) * (1-eta) * uVarIn[:,0,i ,j ] +
( xsi) * (1-eta) * uVarIn[:,0,i ,j+1] +
( xsi) * ( eta) * uVarIn[:,0,i+1,j+1] +
(1-xsi) * ( eta) * uVarIn[:,0,i+1,j ] )
vloc = ( (1-xsi) * (1-eta) * vVarIn[:,0,i ,j ] +
( xsi) * (1-eta) * vVarIn[:,0,i ,j+1] +
( xsi) * ( eta) * vVarIn[:,0,i+1,j+1] +
(1-xsi) * ( eta) * vVarIn[:,0,i+1,j ] )
mask = ma.is_masked( uloc[0])
if mask:
uloc = -999 * np.ones(np.size(uloc))
vloc = -999 * np.ones(np.size(uloc))
return (mask, uloc, vloc)
mesh = slimPre.Mesh(mesh_file_name, mesh_proj='+proj=utm +ellps=WGS84 +zone=35 +south')
region_top = slimPre.Region(mesh, 'top_lake')
lonlat_top = slimPre.Coordinate_system(region_top, data_proj='+proj=latlong +ellps=WGS84')
lonlat_coordinates = lonlat_top.coordinates * 180 / np.pi
n = np.shape(lonlat_coordinates)[0]
uData = np.empty((t, n))
vData = np.empty((t, n))
dltab = [.02, .02, .05, .05, .07, .07, .1, .1, .2, .2]
itMax = 10
for i in range(n):
x = lonlat_coordinates[i,0]
y = lonlat_coordinates[i,1]
(mask, uloc, vloc) = interpolate(x, y)
it = 0
origMask = mask
while mask and it < itMax :
dl = dltab[it]
if (it % 2 == 0) :
(mask1, uloc1, vloc1) = interpolate(x+dl, y)
(mask2, uloc2, vloc2) = interpolate(x-dl, y)
(mask3, uloc3, vloc3) = interpolate(x, y+dl)
(mask4, uloc4, vloc4) = interpolate(x, y-dl)
else :
(mask1, uloc1, vloc1) = interpolate(x+dl, y+dl)
(mask2, uloc2, vloc2) = interpolate(x-dl, y-dl)
(mask3, uloc3, vloc3) = interpolate(x-dl, y+dl)
(mask4, uloc4, vloc4) = interpolate(x+dl, y-dl)
mask = (mask1+mask2+mask3+mask4 == 4)
it = it + 1
if origMask and not mask:
uloc = ( (uloc1 * (not mask1) + uloc2 * (not mask2) +
uloc3 * (not mask3) + uloc4 * (not mask4) ) / (4-(mask1+mask2+mask3+mask4)) )
vloc = ( (vloc1 * (not mask1) + vloc2 * (not mask2) +
vloc3 * (not mask3) + vloc4 * (not mask4) ) / (4-(mask1+mask2+mask3+mask4)) )
uData[:,i] = uloc
vData[:,i] = vloc
(uData, vData) = lonlat_top.rotate(uData, vData)
### TIME REINTERPOLATION ###
nbIter = 5832-1
nbIter = int( 8*365*5.9)
timeData = [timeVarIn[0] + 3600*3*i for i in range(nbIter)]
uData2 = np.empty((nbIter, n))
vData2 = np.empty((nbIter, n))
j = 0
for tt in range(nbIter):
tloc = timeData[tt]
if timeVarIn[j+1] < tloc:
j = j+1
t0 = timeVarIn[j]
t1 = timeVarIn[j+1]
xsi = (tloc - t0) / (t1 - t0)
if xsi < 0 or xsi > 1:
print('WROOOOOOOONG')
exit(-1)
uData2[tt,:] = (1-xsi) * uData[j,:] + xsi * uData[j+1,:]
vData2[tt,:] = (1-xsi) * vData[j,:] + xsi * vData[j+1,:]
time = slimPre.Time(time_vector=timeData, initial_time = '1996-01-01 00:00:00')
slimPre.write_file(pre_data_dir+'wind2000-2005.nc', region=region_top, time=time, data=[('u', uData2), ('v', vData2)])
#print('Preprocessing wind')
#rho0 = 1000
#n = int(35040/2)
#t = [1800*i for i in range(int(3*n))]
#t = np.array(t)
#initial_time = '1993-04-01 00:00:00'
#time = slimPre.Time(time_vector=t, initial_time=initial_time)
#vstress = np.zeros((n))
#f = open('./windDataJaya/ty_Mpu_93_94.dat', 'r')
#for i in range(n):
# vstress[i] = float(f.readline()) * rho0 * 1
#f.close()
#ustress = np.zeros((n))
#
#ustress3 = np.concatenate((ustress, ustress, ustress))
#vstress3 = np.concatenate((vstress, vstress, vstress))
#slimPre.write_file(pre_data_dir+'tanga_uniform_windStress.nc', region=None, time=time, data=[('ustress', ustress3), ('vstress', vstress3)])
print('Preprocessing Call back temperature')
from netCDF4 import Dataset
import numpy as np
from scipy import ma
iF = Dataset(iFile, 'r')
lonVarIn = iF.variables['lon'][:]
latVarIn = iF.variables['lat'][:]
timeVarIn = iF.variables['time']
TVarIn = iF.variables['T_S_LAKE'][:]
minLon = np.amin(lonVarIn)
maxLon = np.amax(lonVarIn)
minLat = np.amin(latVarIn)
maxLat = np.amax(latVarIn)
t = np.size(uVarIn,0)
print('min/max lon %g/%g' % (minLon, maxLon) )
print('min/max lat %g/%g' % (minLat, maxLat) )
def interpolateT(lonloc, latloc):
i = (latloc < latVarIn[:,0]).argmax() -1
j = (lonloc < lonVarIn[i,:]).argmax() -1
i = (latloc < latVarIn[:,j]).argmax() -1
j = (lonloc < lonVarIn[i,:]).argmax() -1
x0 = 0.5 * ( lonVarIn[i,j] + lonVarIn[i+1,j] )
x1 = 0.5 * ( lonVarIn[i,j+1] + lonVarIn[i+1,j+1] )
y0 = 0.5 * ( latVarIn[i,j] + latVarIn[i,j+1] )
y1 = 0.5 * ( latVarIn[i+1,j] + latVarIn[i+1,j+1] )
xsi = (lonloc - x0) / (x1 - x0)
eta = (latloc - y0) / (y1 - y0)
if i < 0: xsi = 0
if j < 0: eta = 0
#if xsi < 0 or xsi > 1:
# print('warning xsi !! (i,j) = (%d,%d) xsi = %g' % (i,j,xsi))
#if eta < 0 or eta > 1:
# print('warning eta !! (i,j) = (%d,%d) eta = %g' % (i,j,eta))
xsi = min(1,max(0,xsi))
eta = min(1,max(0,eta))
Tloc = ( (1-xsi) * (1-eta) * TVarIn[:,i ,j ] +
( xsi) * (1-eta) * TVarIn[:,i ,j+1] +
( xsi) * ( eta) * TVarIn[:,i+1,j+1] +
(1-xsi) * ( eta) * TVarIn[:,i+1,j ] )
mask = ma.is_masked( Tloc[0])
if mask:
Tloc = -999 * np.ones(np.size(Tloc))
return (mask, Tloc)
mesh = slimPre.Mesh(mesh_file_name, mesh_proj='+proj=utm +ellps=WGS84 +zone=35 +south')
region = slimPre.Region(mesh)
lonlat_top = slimPre.Coordinate_system(region, data_proj='+proj=latlong +ellps=WGS84')
lonlat_coordinates = lonlat_top.coordinates * 180 / np.pi
n = np.shape(lonlat_coordinates)[0]
TData = np.empty((t, n))
dltab = [.02, .02, .05, .05, .07, .07, .1, .1, .2, .2]
itMax = 10
for i in range(n):
x = lonlat_coordinates[i,0]
y = lonlat_coordinates[i,1]
(mask, Tloc) = interpolateT(x, y)
it = 0
origMask = mask
while mask and it < itMax :
dl = dltab[it]
if (it % 2 == 0) :
(mask1, Tloc1) = interpolateT(x+dl, y)
(mask2, Tloc2) = interpolateT(x-dl, y)
(mask3, Tloc3) = interpolateT(x, y+dl)
(mask4, Tloc4) = interpolateT(x, y-dl)
else :
(mask1, Tloc1) = interpolateT(x+dl, y+dl)
(mask2, Tloc2) = interpolateT(x-dl, y-dl)
(mask3, Tloc3) = interpolateT(x-dl, y+dl)
(mask4, Tloc4) = interpolateT(x+dl, y-dl)
mask = (mask1+mask2+mask3+mask4 == 4)
it = it + 1
if origMask and not mask:
Tloc = ( (Tloc1 * (not mask1) + Tloc2 * (not mask2) +
Tloc3 * (not mask3) + Tloc4 * (not mask4) ) / (4-(mask1+mask2+mask3+mask4)) )
TData[:,i] = Tloc -273.15
### TIME REINTERPOLATION ###
nbIter = 5832-1
nbIter = int( 8*365*5.9)
timeData = [timeVarIn[0] + 3600.*3*i for i in range(nbIter)]
TData2 = np.empty((nbIter, n))
j = 0
for tt in range(nbIter):
tloc = timeData[tt]
if timeVarIn[j+1] < tloc:
j = j+1
t0 = timeVarIn[j]
t1 = timeVarIn[j+1]
xsi = (tloc - t0) / (t1 - t0)
if xsi < 0 or xsi > 1:
print('WROOOOOOOONG')
exit(-1)
TData2[tt,:] = (1-xsi) * TData[j,:] + xsi * TData[j+1,:]
time = slimPre.Time(time_vector=timeData, initial_time = '1996-01-01 00:00:00')
slimPre.write_file(pre_data_dir+'surface_temperature.nc', region=region, time=time, data=[('temp',TData2)])
#slimPre.write_file(pre_data_dir+'surface_init_temperature.nc', region=region, data=[('temp',TData2[0,:])])
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(pre_data_dir+'initial_temperature.nc', region=region_global, time=None, data=[('T_init', temp)])
print('Preprocessing wind')
time = slimPre.Time(time_vector=[900*i for i in range(10000)])
uwind = np.zeros(time._time.shape)
vwind = 0.2 * np.sin(time._time * 2 * np.pi / (4*86400))
slimPre.write_file(pre_data_dir+'windStress.nc',region=None,time=time,data=[('u',uwind), ('v',vwind)])
from dgpy import *
import slim3d
domain = slim3d.Domain(mesh_file_name)
slimSolver = domain._slimSolver
layerDof = dgDofContainer(slimSolver.groups3d, 1)
tmpDof = dgDofContainer(slimSolver.groups3d, 1)
slimPre.make_directory(pre_data_dir_base+'Tinit_v6/')
for iProc in range(8):
slimPre.fetch_ftp(('slim_data/tanganyika/Tinit_v6/temperature-000699_COMP_0.msh_%g' % iProc), pre_data_dir_base+'Tinit_v6/')
Tinit_v6_file_name = slimPre.fetch_ftp('slim_data/tanganyika/Tinit_v6/temperature-000699.idx', pre_data_dir_base+'Tinit_v6/')
print( 'Preprocessing done ' )
slimPre.exit(0)
......@@ -2,40 +2,71 @@ import slim3d
import slimPre
import shutil
pythonFile = __file__
print('running file %s\n' % pythonFile)
output_dir = 'output/'
pre_data_dir_base = 'data_'+slimPre.partition_nb()+'/'
pre_data_dir_base = 'data_051_'+slimPre.partition_nb()+'/'
pre_data_dir = pre_data_dir_base+slimPre.partition_id()+'/'
slimPre.make_directory(output_dir)
shutil.copyfile(pre_data_dir_base+'mesh3d.msh', output_dir+'mesh3d.msh')
shutil.copyfile(pre_data_dir_base+'mesh3d_show.msh', output_dir+'mesh3d_show.msh')
shutil.copyfile(pythonFile, output_dir+'run.py')
domain = slim3d.Domain(pre_data_dir_base+'mesh3d.msh', reference_density=997.22441884)
domain = slim3d.Domain(pre_data_dir_base+'mesh3d.msh', reference_density=997.421366)
equations = slim3d.Slim3d_equations(domain, temperature=True)
equations.set_implicit_vertical_diffusion(True)
equations.set_implicit_vertical(True)
equations.set_vertical_viscosity('gotm')
#equations.set_vertical_viscosity('pacanowski_philander')
equations.set_vertical_diffusivity('gotm')
#equations.set_vertical_diffusivity('pacanowski_philander')
equations.set_horizontal_viscosity('smagorinsky')
equations.set_additional_artificial_horizontal_shear(pre_data_dir_base+'additional_shear/additional_shear.idx')
#equations.set_horizontal_diffusivity('okubo')
equations.set_lax_friedrichs_factor(.1)
equations.set_bottom_friction(True, z0B=0.005, z0S=0.02)
equations.set_limiter(True)
equations.set_initial_temperature('netcdf', (pre_data_dir+'initial_temperature.nc', 'T_init'))
#equations.set_initial_temperature('netcdf', (pre_data_dir+'initial_temperature.nc', 'T_init'))
equations.set_initial_temperature('netcdf', pre_data_dir_base+'T_init/T_init.idx')
#equations.set_temperature_boundary_relaxation('netcdf', pre_data_dir_base+'T_init/T_init.idx', tau=4*86400, belowBottom=100)
equations.set_temperature_boundary_relaxation('netcdf', (pre_data_dir+'surface_temperature.nc', 'temp'), tau=10*86400, belowBottom=100)
equations.set_initial_salinity('vertical_gradient', None, 0, 0.)
equations.set_wind_stress('stress', (pre_data_dir+'windStress.nc','u'), (pre_data_dir+'windStress.nc','v') )
equations.set_wind_stress('speed', (pre_data_dir+'wind2000-2005.nc','u'), (pre_data_dir+'wind2000-2005.nc','v') )
equations.set_coriolis((pre_data_dir+'coriolis.nc', 'coriolis'))
equations.set_boundary_coast('coast')
equations.set_vertical_adaptation(tau=3600, uv_factor=0, rho_factor=1, minimum_height=1.5, maximum_height=1500, resize_factor=1.5, background_error=0.00001)
time_loop = slim3d.Loop(equations,
time_step=600,#1800,
export_time_step=3600,#86400,
final_time=86400*50,
output_directory=output_dir)
export_time_step=3600,
ratio_full_export=-1,
initial_time = '2000-12-01 00:00:00',
final_time = '2005-12-01 00:00:00',
output_directory=output_dir,
full_export_directory=output_dir[:-1]+'_full/')
time_loop.export_elevation()
time_loop.export_temperature()
time_loop.export_uv()
time_loop.export_uv2d()
time_loop.export_z_coordinate()
time_loop.export_uv(True)
time_loop.export_uv2d(True)
time_loop.export_w()
time_loop.export_vertical_viscosity()
time_loop.export_vertical_diffusivity()
time_loop.loop()
#time_loop.loop()
time_loop.setup()
while time_loop.get_time() < time_loop.final_time:
time_loop.advance_one_time_step()
time_loop.check_sanity()
#print('******************')
#print(time_loop.equations._slimSolver.getUseConservativeALE())
if time_loop.check_export():
time_loop.print_iter_info()
[relVolErr, relSErr, relTErr] = time_loop.check_mass_conservation()
[absSDev, absTDev] = time_loop.check_tracer_consistency()
time_loop.export_fields()
if time_loop.check_full_export():
time_loop.export_full()
time_loop.terminate()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment