Commit bec04ba2 authored by Philippe Delandmeter's avatar Philippe Delandmeter
Browse files

3D Interface. Should be working now!! with almost all open boundaries!!

(constant flow at open bnd not done in prepro. Bathymetry from 3d mesh
still unavailable)
This is a non checked commit.
3D interface is not veriried yet
parent fe97b80d
import slimPre
import numpy as np
pre_data_dir = 'data/'
mesh_file_name = 'data/mesh3d.msh'
mesh = slimPre.Mesh(mesh_file_name, mesh_proj='+proj=utm +ellps=WGS84 +zone=51')
region_global = slimPre.Region(mesh, '')
lonlat_global = slimPre.Coordinate_system(region_global, data_proj='+proj=latlong +ellps=WGS84')
initial_time = '2016-01-01 00:00:00'
final_time = '2016-01-05 00:00:00'
### Coriolis Process ###
print('Preprocessing Coriolis')
lonlat_coordinates = lonlat_global.coordinates
corio = 2 * 7.292e-5 * np.sin(lonlat_coordinates[:,1])
slimPre.write_file(pre_data_dir+'coriolis.nc', region=region_global, time=None, data=[('coriolis', corio)])
# check if correctly exported
slimPre.netcdf_to_msh(mesh_file_name, pre_data_dir+'coriolis.nc', 'coriolis', pre_data_dir+'coriolis')
### Initial conditions Process ###
print('Preprocessing Temperature and Salinity')
xyz_coordinates = region_global.coordinates
temp = 15 * np.ones((np.shape(lonlat_coordinates)[0]))
sal = 35 * np.ones((np.shape(lonlat_coordinates)[0]))
eta = 1- np.minimum(np.hypot(xyz_coordinates[:,0]-2.3e5, xyz_coordinates[:,1]-4.26e6)/5e4, 1)
slimPre.write_file(pre_data_dir+'initialCondition.nc', region=region_global, time=None, data=[('temperature', temp), ('salinity', sal), ('elevation', eta)])
slimPre.netcdf_to_msh(mesh_file_name, pre_data_dir+'initialCondition.nc', 'elevation', pre_data_dir+'elevationInit')
### River Discharge Process ###
print('Preprocessing River Discharge')
t = [3600*j for j in range(300)]
time = slimPre.Time(time_vector=t, initial_time=initial_time)
discharge = np.ones(len(t)) * 1000
discharge[0] = 0
slimPre.write_file(pre_data_dir+'river_discharge.nc', region=None, time=time, data=[('river_discharge', discharge)])
slimPre.write_file(pre_data_dir+'river_salinity.nc', region=None, time=None, data=[('river_salinity', 0)])
### Open Sea Process ###
print('Preprocessing Open Sea Boundary')
time = slimPre.Time(initial_time=initial_time, final_time=final_time, time_step=900.)
# Yellow Sea North
region = slimPre.Region(mesh, physical_tags=['yellow_sea_north'])
(ssh, ssu, ssv) = slimPre.tpxo_tide(region, time)
#flow = slimPre.residual_flow(region, 0.1*1e6, (pre_data_dir+'bathymetry.nc', 'bathymetry'))
#
#ssu += flow[None,:,0]
#ssv += flow[None,:,1]
slimPre.write_file(pre_data_dir+'yellow_sea_north.nc', region=region, time=time, data=[('h', ssh), ('u', ssu), ('v', ssv)])
# Yellow Sea Center
region = slimPre.Region(mesh, physical_tags=['yellow_sea_center'])
(ssh, ssu, ssv) = slimPre.tpxo_tide(region, time)
slimPre.write_file(pre_data_dir+'yellow_sea_center.nc', region=region, time=time, data=[('h', ssh), ('u', ssu), ('v', ssv)])
# Yellow Sea South
region = slimPre.Region(mesh, physical_tags=['yellow_sea_south'])
(ssh, ssu, ssv) = slimPre.tpxo_tide(region, time)
#flow = slimPre.residual_flow(region, -0.1*1e6, (pre_data_dir+'bathymetry.nc', 'bathymetry'))
#
#ssu += flow[None,:,0]
#ssv += flow[None,:,1]
slimPre.write_file(pre_data_dir+'yellow_sea_south.nc', region=region, time=time, data=[('h', ssh), ('u', ssu), ('v', ssv)])
print( 'Preprocessing done ' )
slimPre.exit(0)
import slimPre
pre_data_dir = 'data/'
slimPre.make_directory(pre_data_dir)
mesh_file_name = '../data/bohai_v1.msh'
bath_file_name = ('../data/bathymetry.nc', 'bathymetry')
nb_layers = 7
factor_show = 1000
z_layers = [0,0.5, 1,2]+[3*i for i in range(1,100)]
def layers_func(x, y, h):
h = max(h, 5.)
layers = [0, 1]
for z in range(nb_layers + 1 - 2) :
layers.append(2 + z * (h-2)/(nb_layers-2))
return layers
#slimPre.extrude(mesh_file_name, bath_file_name, nb_layers=nb_layers , mesh_file_name_out='data/mesh3d.msh' , factor_show=factor_show)
#slimPre.extrude(mesh_file_name, bath_file_name, z_layers=z_layers , mesh_file_name_out='data/mesh3dZ.msh' , factor_show=factor_show)
slimPre.extrude(mesh_file_name, bath_file_name, layers_function=layers_func, mesh_file_name_out='data/mesh3d.msh', factor_show=factor_show)
print('done')
slimPre.exit(0)
import slim3d
from slimPre import make_directory, exit
output_directory = 'output'
make_directory(output_directory)
domain = slim3d.Domain('data/mesh3d.msh')
equations = slim3d.Slim3d_equations(domain, salinity=True)
equations.set_coriolis(('data/coriolis.nc', 'coriolis'))
equations.set_initial_temperature('netcdf', temperature=('data/initialCondition.nc', 'temperature'))
equations.set_initial_salinity('netcdf', salinity=('data/initialCondition.nc', 'salinity'))
#equations.set_initial_elevation(('data/initialCondition.nc', 'elevation'))
equations.set_boundary_coast('coast')
#equations.set_boundary_coast('yellow_river')
#equations.set_boundary_coast('yellow_sea_north')
#equations.set_boundary_coast('yellow_sea_center')
#equations.set_boundary_coast('yellow_sea_south')
equations.set_boundary_open('yellow_river', flux=('data/river_discharge.nc', 'river_discharge'), temperature=('data/initialCondition.nc', 'temperature'), salinity=('data/river_salinity.nc', 'river_salinity'))
equations.set_boundary_open('yellow_sea_north', eta=('data/yellow_sea_north.nc', 'h'), u=('data/yellow_sea_north.nc', 'u'), v=('data/yellow_sea_north.nc', 'v'), temperature=('data/initialCondition.nc', 'temperature'), salinity=('data/initialCondition.nc', 'salinity'), transport=True)
equations.set_boundary_open('yellow_sea_center', eta=('data/yellow_sea_center.nc', 'h'), u=('data/yellow_sea_center.nc', 'u'), v=('data/yellow_sea_center.nc', 'v'), temperature=('data/initialCondition.nc', 'temperature'), salinity=('data/initialCondition.nc', 'salinity'), transport=True)
equations.set_boundary_open('yellow_sea_south', eta=('data/yellow_sea_south.nc', 'h'), u=('data/yellow_sea_south.nc', 'u'), v=('data/yellow_sea_south.nc', 'v'), temperature=('data/initialCondition.nc', 'temperature'), salinity=('data/initialCondition.nc', 'salinity'), transport=True)
initial_time = '2016-01-01 00:00:00'
final_time = '2016-01-05 00:00:00'
time_loop = slim3d.Loop(equations, 1, 50, 300, initial_time, final_time, output_directory)
time_loop.export_elevation()
time_loop.export_salinity()
time_loop.export_uv()
time_loop.export_uv2d()
time_loop.setup()
time_loop.loop()
......@@ -103,10 +103,10 @@ class Slim3d_equations :
def set_initial_elevation(self, elevation):
self._initial_elevation = elevation
def set_initial_salinity(self, mode, salinity, surface_salinity, vertical_gradient):
def set_initial_salinity(self, mode, salinity=None, surface_salinity=0, vertical_gradient=0):
self._initial_salinity = (mode, salinity, surface_salinity, vertical_gradient)
def set_initial_temperature(self, mode, temperature, surface_temperature, vertical_gradient):
def set_initial_temperature(self, mode, temperature=None, surface_temperature=0, vertical_gradient=0):
self._initial_temperature = (mode, temperature, surface_temperature, vertical_gradient)
def set_boundary_coast(self, physical_tags):
......@@ -116,13 +116,13 @@ class Slim3d_equations :
for tag in physical_tags:
self._boundary_coast.append(tag)
def set_boundary_open(self, physical_tags, uv=None, eta=None, salinity=None, temperature=None):
def set_boundary_open(self, physical_tags, eta=None, u=None, v=None, salinity=None, temperature=None, flux=None, transport=False):
if slim_private._is_string(physical_tags):
openBnd = slim_private.OpenBnd_3d(physical_tags, uv, eta, salinity, temperature)
openBnd = slim_private.OpenBnd_3d(physical_tags, eta, u, v, salinity, temperature, flux, transport)
self._boundary_open.append(openBnd)
else:
for tag in physical_tags:
openBnd = slim_private.OpenBnd_3d(tag, uv, eta, salinity, temperature)
openBnd = slim_private.OpenBnd_3d(tag, eta, u, v, salinity, temperature, flux, transport)
self._boundary_open.append(openBnd)
class Loop:
......
......@@ -21,12 +21,14 @@ def slim3d_setup(loop):
f = slimSolver.functions
f.initializeBath()
### TO DO
### TODO
#f.windFunc =
#f.windStressFunc =
if eq._initial_elevation:
f.etaInitFunc = slim_private._load_function(eq._initial_elevation, slimSolver.groups2d)
else:
f.etaInitFunc = dgpy.functionConstant(0.)
if eq._initial_salinity:
if eq._initial_salinity[0] == 'netcdf':
......@@ -38,7 +40,7 @@ def slim3d_setup(loop):
else:
dgpy.Msg.Fatal('this mode does not exist')
elif (not eq._linear_density) or (eq._linear_density[0]):
dgpy.Msg.Fatal('Initial temperature must be set if rhoFunc is not set')
dgpy.Msg.Fatal('Initial salinity must be set if rhoFunc is not set')
if eq._initial_temperature:
if eq._initial_temperature[0] == 'netcdf':
......@@ -115,6 +117,12 @@ def slim3d_setup(loop):
def etaOpenGivenUVNumpy(cmap, val, uv, bath, eta, normals) :
un = uv[:,0] * normals[:,0] + uv[:,1] * normals[:,1]
val[:] = numpy.sqrt( (bath[:]+eta[:]) / g) * un
def transport2UVNumpy(cmap, val, uv, bath, eta) :
val[:,0] = uv[:,0] / ( bath[:] + eta[:] )
val[:,1] = uv[:,1] / ( bath[:] + eta[:] )
def mergeUVNumpy(cmap, val, u, v) :
val[:,0] = u[:]
val[:,1] = v[:]
eta = f.eta2dFunc
uv = d.uvDof.getFunction()
......@@ -124,25 +132,67 @@ def slim3d_setup(loop):
eq._uv_open = {}
eq._uvInt2d_open = {}
eq._uvAv2d_open = {}
eq._uvTransport_open = {}
eq._uvAv2dTransport_open = {}
eq._u_open = {}
eq._uAv2d_open = {}
eq._v_open = {}
eq._vAv2d_open = {}
eq._flux_section = {}
eq._flux = {}
eq._flux2d = {}
eq._eta_open = {}
eq._S_open = {}
eq._T_open = {}
for openBnd in eq._boundary_open:
if openBnd.uv and openBnd.eta:
eq._uv_open[openBnd] = slim_private._load_function(openBnd.uv, slimSolver.groups3d)
eq._uvInt2d_open[openBnd] = slim_private._load_function(openBnd.uv, slimSolver.groups2d)
eq._uvAv2d_open[openBnd] = slim_private._load_function(openBnd.uv, slimSolver.groups2d)
if openBnd.u and openBnd.eta:
eq._u_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups3d)
eq._v_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups3d)
eq._uAv2d_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups2d)
eq._vAv2d_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups2d)
if openBnd.transport:
eq._uvTransport_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._u_open[openBnd], eq._v_open[openBnd]])
eq._uv_open[openBnd] = dgpy.functionNumpy(2, transport2UVNumpy, [eq._uvTransport_open[openBnd], f.bathFunc2d, eta])
eq._uvAv2dTransport_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]])
eq._uvAv2d_open[openBnd] = dgpy.functionNumpy(2, transport2UVNumpy, [eq._uvAv2dTransport_open[openBnd], f.bathFunc2d, eta])
eq._uvInt2d_open[openBnd] = eq._uvAv2d_open[openBnd]
else:
eq._uv_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._u_open[openBnd], eq._v_open[openBnd]])
eq._uvInt2d_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]])
eq._uvAv2d_open[openBnd] = eq._uvAv2d_open[openBnd]
eq._eta_open[openBnd] = slim_private._load_function(openBnd.eta, slimSolver.groups2d)
elif openBnd.uv:
eq._uv_open[openBnd] = slim_private._load_function(openBnd.uv, slimSolver.groups3d)
eq._uvInt2d_open[openBnd] = slim_private._load_function(openBnd.uv, slimSolver.groups2d)
eq._uvAv2d_open[openBnd] = slim_private._load_function(openBnd.uv, slimSolver.groups2d)
elif openBnd.u:
eq._u_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups3d)
eq._v_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups3d)
eq._uAv2d_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups2d)
eq._vAv2d_open[openBnd] = slim_private._load_function(openBnd.u, slimSolver.groups2d)
if openBnd.transport:
eq._uvTransport_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._u_open[openBnd], eq._v_open[openBnd]])
eq._uv_open[openBnd] = dgpy.functionNumpy(2, transport2UVNumpy, [eq._uvTransport_open[openBnd], f.bathFunc2d, eta])
eq._uvAv2dTransport_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]])
eq._uvAv2d_open[openBnd] = dgpy.functionNumpy(2, transport2UVNumpy, [eq._uvAv2dTransport_open[openBnd], f.bathFunc2d, eta])
eq._uvInt2d_open[openBnd] = eq._uvAv2d_open[openBnd]
else:
eq._uv_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._u_open[openBnd], eq._v_open[openBnd]])
eq._uvInt2d_open[openBnd] = dgpy.functionNumpy(2, mergeUVNumpy, [eq._uAv2d_open[openBnd], eq._vAv2d_open[openBnd]])
eq._uvAv2d_open[openBnd] = eq._uvAv2d_open[openBnd]
eq._eta_open[openBnd] = dgpy.functionNumpy(1, etaOpenGivenUVNumpy, [uvAv2d, f.bathFunc2d, eta, dgpy.function.getNormals()])
elif openBnd.eta:
eq._uv_open[openBnd] = dgpy.functionNumpy(2, uvOpenGivenEtaNumpy, [uv, dgpy.function.getNormals()])
eq._uvInt2d_open[openBnd] = dgpy.functionNumpy(2, uvOpenGivenEtaNumpy, [uvInt2d, dgpy.function.getNormals()])
eq._uvAv2d_open[openBnd] = dgpy.functionNumpy(2, uvOpenGivenEtaNumpy, [uvAv2d, dgpy.function.getNormals()])
eq._eta_open[openBnd] = slim_private._load_function(openBnd.eta, slimSolver.groups2d)
elif openBnd.flux:
integInterface = dgpy.dgFunctionIntegratorInterface(slimSolver.groups2d, f.bathFunc2d)
sectionFM = dgpy.fullMatrixDouble(1,1)
integInterface.compute(openBnd.tag, sectionFM)
eq._flux_section[openBnd] = sectionFM(0,0)
eq._flux[openBnd] = slim_private._load_function(openBnd.flux, slimSolver.groups3d)
eq._flux2d[openBnd] = slim_private._load_function(openBnd.flux, slimSolver.groups2d)
eq._uv_open[openBnd] = dgpy.slimFlowRateToVelocity(slimSolver.groups2d, eq._flux_section[openBnd], eq._flux[openBnd], f.bathFunc2d, eta)
eq._uvInt2d_open[openBnd] = dgpy.slimFlowRateToVelocity(slimSolver.groups2d, eq._flux_section[openBnd], eq._flux2d[openBnd], f.bathFunc2d, eta)
eq._uvAv2d_open[openBnd] = dgpy.slimFlowRateToVelocity(slimSolver.groups2d, eq._flux_section[openBnd], eq._flux2d[openBnd], f.bathFunc2d, eta)
eq._eta_open[openBnd] = eta
else:
eq._uv_open[openBnd] = dgpy.functionNumpy(2, uvOpenFreeNumpy, [f.bathFunc2d, eta, dgpy.function.getNormals()])
eq._uvInt2d_open[openBnd] = dgpy.functionNumpy(2, uvOpenFreeNumpy, [f.bathFunc2d, eta, dgpy.function.getNormals()])
......@@ -163,7 +213,7 @@ def slim3d_setup(loop):
horMomEq = e.horMomEq
horMomEq.setLaxFriedrichsFactor(0.0)
if slimSolver.getSolveUVImplicitVerticalDiffusion() or not f.windStressFunc :
#horMomEq.setBoundary0Flux(topTags) # CHECK WHICH TO DO !!!
#horMomEq.setBoundary0Flux(topTags) # CHECK WHICH TODO !!!
horMomEq.setBoundarySymmetry(topTags)
else :
horMomEq.addBoundaryCondition(topTags, horMomEq.newBoundarySurface(f.windStressFunc)) # zero for nonconst tracers!
......
......@@ -214,7 +214,6 @@ def _cpu_time(seconds):
h, m = divmod(m, 60)
return "%d:%02d:%02d" % (h, m, s)
<<<<<<< 158678f887f029c3ceab77a02a388f6a13e9187c
def _new_default_linear_system(domain, equation) :
try :
sys = dgpy.linearSystemPETScBlockDouble()
......@@ -230,7 +229,7 @@ def _new_default_linear_system(domain, equation) :
except :
dgpy.Msg.Fatal("No linear system available")
return sys
=======
def _findTopAndBottomTags(mesh_file_name):
bottomTags = []
topTags = []
......@@ -249,10 +248,14 @@ def _findTopAndBottomTags(mesh_file_name):
return (topTags, bottomTags)
class OpenBnd_3d:
def __init__(self, tag, uv, eta, salinity, temperature):
def __init__(self, tag, eta, u, v, salinity, temperature, flux, transport):
self.tag = tag
self.uv = uv
self.u = u
self.v = v
if (u and not v) or (v and not u):
dgpy.Msg.Fatal('Error. Either u and v or none of them must be set in open boundary conditions')
self.eta = eta
self.salinity = salinity
self.temperature = temperature
>>>>>>> slim3d Interface. Not all boundary conditions are ready yet
self.flux = flux
self.transport = transport
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