Commit 700fa6b7 authored by Valentin Vallaeys's avatar Valentin Vallaeys Committed by Philippe Delandmeter
Browse files

interface3D: create output directory + update benchmarks

parent bfdcc864
import slim3d
from slimPre import make_directory, exit
output_directory = 'output'
make_directory(output_directory)
domain = slim3d.Domain('data/mesh3d.msh')
......
import slimPre
import numpy as np
global_dir = "/media/data/"
slimPre.set_global_data_directory(global_dir)
runFile = "data/"
fileName = "total"#"congoCoast-labels-clean"
dataDir = "slim_data/congo/"
meshFile = dataDir + "mesh/" + fileName + "_utm.msh"
bathDir = dataDir + "bathy/"
smoothBathy = "bathymetry_smooth"
latlon_pj = "+proj=latlong +ellps=WGS84"
utm_pj = "+proj=utm +ellps=WGS84 +zone=32"
tides_initial_time = '2010-12-01 00:00:00'
tides_end_time = '2012-06-01 00:00:00'
tides_time_step = 900
tides_physical_tags = ['ocean']
factor_show = 1000
nb_layers = 3
slimPre.make_directory(runFile)
slimPre.fetch_ftp(meshFile, runFile)
slimPre.fetch_ftp(bathDir+"bath_new.bin")
slimPre.fetch_ftp(bathDir+"dpp1Cor")
slimPre.fetch_ftp(bathDir+"dpp2Cor")
slimPre.fetch_ftp(bathDir+"dpp3Cor")
slimPre.fetch_ftp(bathDir+"dppc")
meshFile = runFile + fileName + "_utm.msh"
mesh = slimPre.Mesh(meshFile, mesh_proj=utm_pj)
region_global = slimPre.Region(mesh)
lonlat_global = slimPre.Coordinate_system(region_global, data_proj=latlon_pj)
lonlat_coordinates_deg = lonlat_global.coordinates * 180. / np.pi
print('Preprocessing Bathymetry')
bath1 = slimPre.interpolate_from_structured_file(lonlat_coordinates_deg, global_dir + bathDir + "dpp1Cor", transpose_matrix=True)
bath2 = slimPre.interpolate_from_structured_file(lonlat_coordinates_deg, global_dir + bathDir + "dpp2Cor", transpose_matrix=True)
bath3 = slimPre.interpolate_from_structured_file(lonlat_coordinates_deg, global_dir + bathDir + "dpp3Cor", transpose_matrix=True)
bathC = slimPre.interpolate_from_structured_file(lonlat_coordinates_deg, global_dir + bathDir + "dppc", transpose_matrix=True)
bathG = slimPre.interpolate_from_structured_file(lonlat_coordinates_deg, global_dir + bathDir + "bath_new.bin", transpose_matrix=True)
bath = np.maximum(bath1,np.maximum(bath2,np.maximum(bath3,np.maximum(bathC,np.maximum(bathG, 5.)))))
slimPre.write_file(runFile + "bath.nc", region=region_global, time=None, data=[('bath', bath)])
slimPre.netcdf_to_msh(meshFile, runFile + "bath.nc", "bath", runFile + "bathy_raw", time = 0)
slimPre.smooth_bathymetry(mesh, (runFile + "bath.nc", 'bath'), runFile+smoothBathy+'.nc', coefficient = 1./3., transform_p0=False)
slimPre.netcdf_to_msh(meshFile, runFile + smoothBathy+'.nc', "bath", runFile + "bathy_smooth", time = 0)
print('Extruding mesh')
slimPre.extrude(meshFile, (runFile+smoothBathy+'.nc','bath'), nb_layers=nb_layers, mesh_file_name_out=runFile+'mesh3d.msh', factor_show=factor_show)
print('Loading 3D mesh')
meshFile = runFile + "mesh3d.msh"
mesh = slimPre.Mesh(meshFile, mesh_proj=utm_pj)
region_global = slimPre.Region(mesh)
lonlat_global = slimPre.Coordinate_system(region_global, data_proj=latlon_pj)
print('Preprocessing Coriolis')
lonlat_coordinates = lonlat_global.coordinates
corio = 2 * 7.292e-5 * np.sin(lonlat_coordinates[:,1])
slimPre.write_file(runFile+'coriolis.nc', region=region_global, time=None, data=[('coriolis', corio)])
# check if correctly exported
slimPre.netcdf_to_msh(meshFile, runFile+'coriolis.nc', 'coriolis', runFile+'coriolis')
print('Preprocessing initial temperature and salinity')
slimPre.write_file(runFile+'initialCondition.nc', region=None, time=None, data=[('temperature', 10), ('salinity', 32)])
print('Preprocessing River Discharge')
slimPre.write_file(runFile+'constant_river_discharge.nc', region=None, time=None, data=[('constant_river_discharge', 40000)])
slimPre.write_file(runFile+'river_salinity.nc', region=None, time=None, data=[('river_salinity', 0)])
print('Preprocessing Open Sea Boundary')
timeTide = slimPre.Time(initial_time=tides_initial_time, final_time=tides_end_time, time_step=tides_time_step)
region = slimPre.Region(mesh, physical_tags=tides_physical_tags)
slimPre.tpxo_tide(region, timeTide, runFile+'tides.nc')
print('Preprocessing wind')
slimPre.write_file(runFile+'wind.nc', region=None, time=None, data=[('windx', 0), ('windy', 1)])
#print('Preprocessing wind')
#from netCDF4 import Dataset
#windFile = slimPre.fetch_ftp(dataDir+"wind/wind2011-2012.nc")
#f = Dataset(windFile,"r")
#
#t = f.variables['time']
#t_final = []
#t_final[:] = [x*3600 for x in t]
#print t[0]
#print t_final[0]
#u = np.array(f.variables['U_GRD_L103'][:])
#v = np.array(f.variables['V_GRD_L103'][:])
#x = f.variables['lon']
#y = f.variables['lat']
#
#nu_time = u.shape[0]
#nu_long = u.shape[2]
#nu_lat = u.shape[1]
#
#nx = x.shape[0]
#ny = y.shape[0]
#
#ox = x[0]
#oy = y[0]
#
#dx = x[1]-x[0]
#dy = y[1]-y[0]
#
#def interpolate_on_structured_grid(x, y, ox, oy, dx, dy, nx, ny, timej, data):
# if ((x < ox).any() or (x > ox + (nx-1)*dx).any()) and (dx > 0) :
# print('Fatal. (x,y) is outside the grid 1!!')
# exit(-1)
# if ((x > ox).any() or (x < ox + (nx-1)*dx).any()) and (dx < 0):
# print('Fatal. (x,y) is outside the grid 2!!')
# exit(-1)
# if ((y < oy).any() or (y > oy + (ny-1)*dy).any()) and (dy > 0):
# print('Fatal. (x,y) is outside the grid 3!!')
# exit(-1)
# if ((y > oy).any() or (y < oy + (ny-1)*dy).any()) and (dy < 0):
# print('Fatal. (x,y) is outside the grid 4!!')
# exit(-1)
# rx = x - ox
# ry = y - oy
# ix = (rx / dx).astype(int)
# iy = (ry / dy).astype(int)
# xsi = (rx - ix*dx) / dx
# eta = (ry - iy*dy) / dy
# val = (1-xsi)*(1-eta)*data[timej,iy,ix]+xsi*(1-eta)*data[timej,iy,(ix+1)%nx]+xsi*eta*data[timej,(iy+1)%ny,(ix+1)%nx]+(1-xsi)*eta*data[timej,(iy+1)%ny,ix]
# return val
#
#xyz = np.array(lonlat_coordinates) * 180./np.pi
#n = len(xyz)
#vit_u = np.empty((nu_time,n))
#vit_v = np.empty((nu_time,n))
#
#for j in range(nu_time):
# print j+1,"/",nu_time
# vit_u[j,:] = interpolate_on_structured_grid(xyz[:,0],xyz[:,1],ox,oy,dx,dy,nx,ny,j,u)
# vit_v[j,:] = interpolate_on_structured_grid(xyz[:,0],xyz[:,1],ox,oy,dx,dy,nx,ny,j,v)
#
#(vit_u, vit_v) = lonlat_global.rotate(vit_u, vit_v)
#wind_initial_time = '2011-01-01 00:00:00'
#timeWind = slimPre.Time(time_vector=t_final, initial_time=wind_initial_time)
#slimPre.write_file(runFile+'wind.nc', region=region_global, time=timeWind, data=[('vit_u', vit_u), ('vit_v', vit_v)])
print('preprocessing done')
slimPre.exit(0)
import slim3d
runFile = 'data/'
meshFile = runFile+'mesh3d.msh'
sim_Ti = "2011-01-01 00:00:00"
sim_Tf = "2011-01-02 00:00:00"
sim_export = 300
output_directory = "./output2"
domain = slim3d.Domain(meshFile)
equations = slim3d.Slim3d_equations(domain, salinity=True)
equations.set_coriolis((runFile+'coriolis.nc', 'coriolis'))
equations.set_horizontal_viscosity('smagorinsky')
equations.set_implicit_vertical_diffusion(True)
equations.set_vertical_diffusivity('gotm')
equations.set_vertical_viscosity('gotm')
equations.set_wind_stress("stress", wind_x=(runFile+'wind.nc', 'windx'), wind_y=(runFile+'wind.nc', 'windy'))
equations.set_initial_temperature('netcdf', temperature=(runFile+'initialCondition.nc', 'temperature'))
equations.set_initial_salinity('netcdf', salinity=(runFile+'initialCondition.nc', 'salinity'))
equations.set_boundary_coast('coast')
equations.set_boundary_open('source1', flux=(runFile+'constant_river_discharge.nc', 'constant_river_discharge'), salinity=('data/river_salinity.nc', 'river_salinity'))
equations.set_boundary_open('ocean', eta=(runFile+'tides.nc', 'h'), u=(runFile+'tides.nc', 'u'), v=(runFile+'tides.nc', 'v'), salinity=(runFile+'initialCondition.nc', 'salinity'), transport=True)
time_loop = slim3d.Loop(equations, 1, 50, sim_export, sim_Ti, sim_Tf, output_directory)
time_loop.export_elevation()
time_loop.export_salinity()
time_loop.export_uv()
time_loop.setup()
time_loop.loop()
......@@ -196,6 +196,10 @@ class Loop:
def setup(self):
slim_private.slim3d_setup(self)
slimSolver = self._slimSolver
if dgpy.Msg.GetCommRank() == 0:
if not slim_private.os.path.exists(self._odir):
slim_private.os.mkdir(self._odir)
dgpy.Msg.Barrier()
d = slimSolver.getDofs()
if self._export_uv: self._export.append(dgpy.dgIdxExporter(d.uvDof, self._odir + "/uv", True))
if self._export_w: self._export.append(dgpy.dgIdxExporter(d.wDof3d, self._odir + "/w"))
......
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