Commit 6dae3b12 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

Merge branch simplified consistant-conservative temporal scheme

parents 753bc118 c0b7ceca
Pipeline #1758 passed with stage
in 34 minutes and 34 seconds
......@@ -18,6 +18,7 @@ option(ENABLE_PETSC4PY "Enable PETSc4Py" ON)
option(ENABLE_FUNCTIONC "ENABLE C CALLBACK" ON)
project(dg CXX)
set( CMAKE_EXPORT_COMPILE_COMMANDS 1 )
set(CMAKE_CXX_STANDARD 11)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF)
......
......@@ -38,8 +38,7 @@ equations.set_boundary_open('burdekinMouth',
salinity=(pre_data_dir+'river_salinity.nc', 'river_salinity'))
time_loop = slim3d.Loop(equations,
maximum_2d_time_step=1,
ratio_3dvs2d_time_step=20,
time_step=20,
export_time_step=1800,
initial_time='2007-01-01 00:00:00',
final_time='2007-02-01 00:00:00',
......
import slimPre
import numpy as np
pre_data_dir_base = 'data_'+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/columbia/meshes/columbia_v1.msh')
mesh_file_name = slimPre.partition_mesh(mesh_file_name, pre_data_dir_base+"mesh2d.msh")
tides_initial_time = '2006-05-1 07:45:00'
tides_final_time = '2006-05-14 00:00:00'
tides_time_step = 900
tides_physical_tags = ["sea_west", "sea_north", "sea_south"]
time_zone = -8.*3600
nb_layers = 10
print('Preprocessing Bathymetry')
bath_file_name = slimPre.fetch_ftp('slim_data/columbia/bath/bath_smooth_v1/bath_smooth_v1.idx', pre_data_dir_base)
bath_file = slimPre.fetch_ftp('slim_data/columbia/bath/bath_smooth_v1/bath_smooth_v1_COMP_0.msh', pre_data_dir_base)
zLayers = [0,1,3,5,10,20,50,100,300,800]
def layers_func(x, y, h):
n = nb_layers-1
f = nb_layers-1
det = (1-pow(f, (n+1.)/n))/(1-pow(f, 1./n))
layers = []
for z in range(n + 1):
a = 0
for i in range(z):
a+= pow(f, float(i)/n) * h /det
if (z < len(zLayers)):
a = min(a, zLayers[z])
layers.append(a)
layers.append(h)
return layers
print('Extruding mesh')
slimPre.extrude(mesh_file_name, bath_file_name, layers_function=layers_func, mesh_file_name_out=pre_data_dir_base+'mesh3d.msh', factor_show=500)
print('Loading 3D mesh')
mesh_file_name = pre_data_dir_base + 'mesh3d.msh'
mesh = slimPre.Mesh(mesh_file_name, mesh_proj="+proj=lcc +datum=NAD27 +lon_0=-120d30 +lat_1=46 +lat_2=44d20 +lat_0=43d40 +x_0=609601.2192024384 +y_0=0")
region_global = slimPre.Region(mesh)
lonlat_global = slimPre.Coordinate_system(region_global, data_proj="+proj=latlong +ellps=WGS84 +nadgrids=WO")
lonlat_coordinates = lonlat_global.coordinates
lonlat_coordinates_deg = lonlat_global.coordinates * 180. / np.pi
print('Preprocessing Coriolis')
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)])
#slimPre.netcdf_to_msh(mesh_file_name, pre_data_dir+'coriolis.nc', 'coriolis', 'coriolis')
print('Preprocessing Open Sea Boundary')
timeTide = slimPre.Time(initial_time=tides_initial_time, final_time=tides_final_time, time_step=tides_time_step)
region = slimPre.Region(mesh, physical_tags=tides_physical_tags)
(ssh,ssu,ssv) = slimPre.tpxo_tide(region, timeTide)
for i in range(len(timeTide._time)):
# from GMT to PST
timeTide._time[i] += time_zone
delta = (timeTide._time[i] - timeTide._time[0]) / (12 * 60. * 60.)
if delta > 1:
continue
if len(ssh) > 1:
ssh[i,:] *= max(0,min(delta,1))
ssu[i,:] *= max(0,min(delta,1))
ssv[i,:] *= max(0,min(delta,1))
slimPre.write_file(pre_data_dir+'tides.nc', region=region, time=timeTide, data=[('h', ssh),('u',ssu),('v',ssv)])
print('Preprocessing river discharge')
bath_file_name = slimPre.fetch_ftp('slim_data/columbia/forcings/beaver_ramp2.nc', pre_data_dir_base)
print('Preprocessing river salt')
slimPre.write_file(pre_data_dir+'river_salinity.nc', region=None, time=None, data=[('river_salinity', 0)])
print('Preprocessing open bnd salt')
slimPre.write_file(pre_data_dir+'open_salinity.nc', region=None, time=None, data=[('salt', 34)])
print('Preprocessing temp')
slimPre.write_file(pre_data_dir+'temp.nc', region=None, time=None, data=[('temp', 10)])
print('preprocessing done')
slimPre.exit(0)
import slim3d
import slimPre
import shutil
sim_Ti = '2006-05-1 00:00:00'
sim_Tf = '2006-05-14 00:30:00'
dt = 30.
sim_export = 900.
sim_exportFullRatio = -1#4*24 #each day
output_dir = 'output/'
pre_data_dir_base = 'data_'+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')
domain = slim3d.Domain(pre_data_dir_base+'mesh3d.msh')
equations = slim3d.Slim3d_equations(domain, salinity=True)
equations.set_coriolis((pre_data_dir+'coriolis.nc', 'coriolis'))
equations.set_horizontal_viscosity('smagorinsky', factor=0.5)
equations.set_implicit_vertical_diffusion(True)
equations.set_vertical_diffusivity('gotm')
equations.set_vertical_viscosity('gotm')
equations.set_limiter(True)
equations.set_bottom_friction(True, z0B=0.005, z0S=0.02)
equations.set_initial_temperature('netcdf', temperature=(pre_data_dir+'temp.nc', 'temp'))
equations.set_initial_salinity('netcdf', salinity=(pre_data_dir+'open_salinity.nc', 'salt'))
coast_tags = ["Land", "land", "island_0", "island_1", "island_2", "island_3",
"island_4", "island_5", "island_6", "island_7", "island_8", "island_9", "island_10"]
for coast in coast_tags:
equations.set_boundary_coast(coast)
equations.set_boundary_open('river',
flux=(pre_data_dir_base+'beaver_ramp2.nc', 'flux'),
salinity=(pre_data_dir+'river_salinity.nc', 'river_salinity'),
temperature=(pre_data_dir+'temp.nc', 'temp'))
open_tags = ["sea_west", "sea_north", "sea_south"]
for tag in open_tags:
equations.set_boundary_open(tag,
eta=(pre_data_dir+'tides.nc', 'h'),
u=(pre_data_dir+'tides.nc', 'u'),
v=(pre_data_dir+'tides.nc', 'v'),
temperature=(pre_data_dir+'temp.nc', 'temp'),
salinity=(pre_data_dir+'open_salinity.nc', 'salt'),
transport=True)
time_loop = slim3d.Loop(equations, time_step=dt,
export_time_step=sim_export,
ratio_full_export=sim_exportFullRatio,
initial_time=sim_Ti, final_time=sim_Tf,
output_directory=output_dir)
time_loop.export_elevation()
time_loop.export_salinity()
time_loop.export_uv()
time_loop.export_uv2d()
time_loop.loop()
\ No newline at end of file
import slimPre
import numpy as np
run_dir = 'data/'
slimPre.make_directory(run_dir)
slimPre.write_file(run_dir+'bath.nc',None,None,[('bath',1000)])
mesh_file = 'square.msh'
print('Extruding mesh')
def shiftOperation(node, iPerBound) :
n = [node[0] - 1.6e5, node[1], node[2]]
return n
cutTags = ["cut"]
pasteTags = ["paste"]
mapFilename = "periodicMesh.txt"
periodicity = (shiftOperation, cutTags, pasteTags, mapFilename)
nPart = 2
slimPre.dgpy.dgMeshPartition(mesh_file, nPart)
mesh_file = mesh_file[:-4] + '_' + str(nPart)+'.msh'
slimPre.extrude(mesh_file, (run_dir+'bath.nc','bath'), nb_layers=20, mesh_file_name_out=run_dir+'mesh3d.msh', periodicity=periodicity)
print('Loading 3D mesh')
mesh_file = run_dir + "mesh3d.msh"
mesh = slimPre.Mesh(mesh_file)
region_global = slimPre.Region(mesh)
print('Preprocessing coriolis')
xyz = region_global.coordinates
coriolis = -1.2e-4
slimPre.write_file(run_dir+'coriolis.nc', region=None, time=None, data=[('coriolis', coriolis)])
slimPre.netcdf_to_msh(mesh_file, run_dir+'coriolis.nc', 'coriolis', 'coriolis')
print('Preprocessing initial temperature')
xyz = region_global.coordinates
theta0 = 10.1 + 3*(-975.-xyz[:,2])/(-975.)
yw = 2.5e5-4e4*np.sin(2*np.pi*3*xyz[:,0]/1.6e5)
fact = 1. - (xyz[:,1]-yw[:])/4.e4
indx = fact < 0.
fact[indx] = 0.
indx = fact > 1.
fact[indx] = 1.
theta = theta0[:]-1.2*fact[:]
yw_t = 2.5e5-2.0e4*np.sin(np.pi*(xyz[:,0]-1.1e5)/2.0e4)
theta_t = 0.3*(1.0 - (xyz[:,1] - yw_t)/2.0e4)
indx = (xyz[:,0] > 1.1e5) * (xyz[:,0] < 1.3e5) * (xyz[:,1] > yw_t[:] - 2.0e4) * (xyz[:,1] < yw_t[:] + 2.0e4)
theta[indx] += theta_t[indx]
slimPre.write_file(run_dir+'initialTemp.nc', region=region_global, time=None, data=[('temp', theta)])
slimPre.netcdf_to_msh(mesh_file, run_dir+'initialTemp.nc', 'temp', 'temp')
print('preprocessing done')
slimPre.exit(0)
import slim3d
runFile = 'data/'
meshFile = runFile+'mesh3d.msh'
sim_Ti = "2012-03-24 00:00:00"
sim_Tf = "2012-10-07 00:00:00"
sim_export = 86400/2.
sim_exportFullRatio = -1
output_directory = "./output"
domain = slim3d.Domain(meshFile, "periodicMesh.txt", reference_density=1000)
equations = slim3d.Slim3d_equations(domain,temperature=True)
equations.set_coriolis((runFile+'coriolis.nc', 'coriolis'))
equations.set_horizontal_viscosity('smagorinsky',factor=0.5)
equations.set_vertical_viscosity("constant", 1.0e-4)
equations.set_implicit_vertical_diffusion(True)
equations.set_lax_friedrichs(0.1)
equations.set_linear_density(mode="temperature", constant_coefficient=5., linear_coefficient=-0.2)
equations.set_initial_temperature('netcdf', temperature=(runFile+'initialTemp.nc','temp'))
equations.set_boundary_coast(["Wall"])
time_loop = slim3d.Loop(equations, time_step=500, export_time_step=sim_export, ratio_full_export=sim_exportFullRatio, initial_time=sim_Ti, final_time=sim_Tf, output_directory=output_directory)
time_loop.export_elevation()
time_loop.export_temperature()
time_loop.export_uv()
time_loop.export_uv2d()
time_loop.export_rho()
time_loop._slimSolver.functions.bottomFrictionDragCoeff2d = slim3d.dgpy.functionConstant(0.01)
time_loop.loop()
L = 1.6e5;
l = 5e5;
Point(1) = {0, 0, 0};
Point(2) = {L, 0, 0};
Point(3) = {L, l, 0};
Point(4) = {0, l, 0};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop(1) = {3, 4, 1, 2};
Plane Surface(1) = {1};
Transfinite Line{1} = 41;
Transfinite Line{2} = 126;
Transfinite Line{3} = 41;
Transfinite Line{4} = 126;
Transfinite Surface{1} = {1,2,3,4};
Recombine Surface{1};
Physical Surface("Surface")={1};
Physical Line("Wall") = {1,3};
Physical Line("cut") = {2};
Physical Line("paste") = {4};
General.ExecutableFileName = "/home/delandmeter/sources/gmsh/projects/dg/build/gmsh/gmsh";
General.FileName = "mesh2dxz_show.msh";
General.GraphicsFontEngine = "StringTexture";
General.RecentFile0 = "mesh3d_show.msh";
General.RecentFile1 = "data/mesh3d_show.msh";
General.RecentFile2 = "data/mesh3d.msh";
General.RecentFile3 = "box.msh";
General.RecentFile4 = "../boxShow.msh";
General.RecentFile5 = "mesh.msh";
General.RecentFile6 = "square.msh";
General.RecentFile7 = "meshShow.msh";
General.RecentFile8 = "temperature.pos";
General.RecentFile9 = "bath.pos";
General.BoundingBoxSize = 617251.9744804386;
General.ClipPositionX = 144;
General.ClipPositionY = 645;
General.ColorScheme = 0;
General.ContextPositionX = 1485;
General.FieldPositionX = 2006;
General.FieldPositionY = 303;
General.FieldHeight = 601;
General.FieldWidth = 903;
General.GraphicsHeight = 1002;
General.GraphicsPositionX = 49;
General.GraphicsPositionY = 24;
General.GraphicsWidth = 1871;
General.MaxX = 550000;
General.MaxY = 5000;
General.MenuWidth = 179;
General.MenuPositionX = 665;
General.MenuPositionY = 624;
General.MessageHeight = 140;
General.MinY = -5000;
General.MinZ = -280000;
General.OptionsPositionX = 1438;
General.OptionsPositionY = 550;
General.PluginPositionX = 635;
General.PluginPositionY = 314;
General.PluginHeight = 641;
General.PluginWidth = 1284;
General.RotationX = 270;
General.ScaleX = 1.324867563680734;
General.ScaleY = 1.324867563680734;
General.ScaleZ = 1.324867563680734;
General.SmallAxes = 0;
General.StatisticsPositionX = 347;
General.StatisticsPositionY = 561;
General.TrackballQuaternion0 = 0.7071067811865475;
General.TrackballQuaternion3 = 0.7071067811865476;
General.TranslationX = -71747.50086237336;
General.TranslationY = 7673.625977357318;
General.VisibilityPositionX = 1340;
General.VisibilityPositionY = 540;
Mesh.ChacoHypercubeDim = 1;
Mesh.ChacoMeshDim1 = 1;
Mesh.NbHexahedra = 872;
Mesh.NbNodes = 1906;
Mesh.NbQuadrangles = 1904;
Mesh.RemeshParametrization = 0;
Mesh.SurfaceEdges = 0;
Mesh.VolumeEdges = 0;
PostProcessing.NbViews = 4;
View[0].FileName = "uv/uv-000013/uv-000013_COMP_0.msh_1";
View[0].GeneralizedRaiseX = "0";
View[0].GeneralizedRaiseY = "0";
View[0].GeneralizedRaiseZ = "v0";
View[0].Name = "u";
View[0].ColormapNumber = 10;
View[0].CustomMax = 26.5;
View[0].CustomMin = 24.3;
View[0].GeneralizedRaiseFactor = 200;
View[0].GeneralizedRaiseView = 3;
View[0].IntervalsType = 3;
View[0].Max = 6.227853439825584;
View[0].MaxX = 550000;
View[0].MaxY = 5000;
View[0].Min = -3.622001857057784;
View[0].MinY = -5000;
View[0].MinZ = -280000;
View[0].NbIso = 21;
View[0].NbTimeStep = 14;
View[0].RangeType = 3;
View[0].ShowElement = 1;
View[0].Time = 0;
View[0].TransformZZ = 0;
View[0].UseGeneralizedRaise = 1;
View[1].FileName = "w/w-000013/w-000013_COMP_0.msh_1";
View[1].GeneralizedRaiseX = "0";
View[1].GeneralizedRaiseY = "0";
View[1].GeneralizedRaiseZ = "v0";
View[1].Name = "w";
View[1].ColormapNumber = 10;
View[1].GeneralizedRaiseFactor = 200;
View[1].GeneralizedRaiseView = 3;
View[1].Max = 1.742882729259567e+97;
View[1].MaxVisible = 9.731221743813752e-312;
View[1].MaxX = 550000;
View[1].MaxY = 5000;
View[1].Min = -8.714413646297006e+96;
View[1].MinVisible = -4.865610871899465e-312;
View[1].MinY = -5000;
View[1].MinZ = -280000;
View[1].NbTimeStep = 14;
View[1].RangeType = 3;
View[1].Time = 0;
View[1].TransformZZ = 0;
View[1].UseGeneralizedRaise = 1;
View[1].Visible = 0;
View[2].FileName = "temperature/temperature-000013/temperature-000013_COMP_0.msh_1";
View[2].GeneralizedRaiseX = "0";
View[2].GeneralizedRaiseY = "0";
View[2].GeneralizedRaiseZ = "v0";
View[2].Name = "temperature";
View[2].ColormapNumber = 10;
View[2].GeneralizedRaiseFactor = 200;
View[2].GeneralizedRaiseView = 3;
View[2].IntervalsType = 3;
View[2].Max = 26.74644988945463;
View[2].MaxVisible = 26.5;
View[2].MaxX = 550000;
View[2].MaxY = 5000;
View[2].Min = 24.11740178250944;
View[2].MinVisible = 24.3;
View[2].MinY = -5000;
View[2].MinZ = -280000;
View[2].NbIso = 21;
View[2].NbTimeStep = 14;
View[2].RangeType = 3;
View[2].Time = 0;
View[2].TransformZZ = 0;
View[2].UseGeneralizedRaise = 1;
View[2].Visible = 0;
View[3].FileName = "z/z-000013/z-000013_COMP_0.msh_1";
View[3].GeneralizedRaiseX = "0";
View[3].GeneralizedRaiseY = "0";
View[3].GeneralizedRaiseZ = "v0";
View[3].Name = "z";
View[3].ColormapNumber = 10;
View[3].GeneralizedRaiseFactor = 200;
View[3].GeneralizedRaiseView = 3;
View[3].Max = 0.07048042688673509;
View[3].MaxX = 550000;
View[3].MaxY = 5000;
View[3].Min = -1400;
View[3].MinVisible = -1400;
View[3].MinY = -5000;
View[3].MinZ = -280000;
View[3].NbTimeStep = 14;
View[3].RangeType = 3;
View[3].Time = 0;
View[3].TransformZZ = 0;
View[3].UseGeneralizedRaise = 1;
View[3].Visible = 0;
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(200000)])
uwind = 0.2 * np.sin(time._time * 2 * np.pi / (8*86400))
uwind = 0.4 * np.sin(time._time * 2 * np.pi / (16*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)
import slim3d
from slimPre import make_directory
import shutil
output_dir = 'output/'
data_dir = './data/'
make_directory(output_dir)
shutil.copyfile(data_dir+'tanga2dxz.msh', output_dir+'mesh2dxz.msh')
shutil.copyfile(data_dir+'tanga2dxz_show.msh', output_dir+'mesh2dxz_show.msh')
domain = slim3d.Domain(data_dir+'tanga2dxz.msh', data_dir+'periodicMesh.txt', reference_density=997.22441884)
equations = slim3d.Slim3d_equations(domain, temperature=True)
equations.set_implicit_vertical_diffusion(True)
#equations.set_implicit_vertical_advection(True)
equations.set_vertical_viscosity('gotm')