Commit 0bb0666b authored by Philippe Delandmeter's avatar Philippe Delandmeter Committed by Valentin Vallaeys
Browse files

slim3d: add 3d benchmarks

tanga 3d + 2dxz + upwelling + Tartinville + Ilicak
parent 34d4585a
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 = 1969;
General.GraphicsPositionY = 135;
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 = 2621;
General.OptionsPositionY = 523;
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";
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 = 70011.96365286827;
View[0].MaxX = 550000;
View[0].MaxY = 5000;
View[0].Min = -23049620.62291582;
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";
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 = 5721200.907856209;
View[1].MaxX = 550000;
View[1].MaxY = 5000;
View[1].Min = -2878507.451315449;
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";
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].Max = 315.2077599878262;
View[2].MaxVisible = 26.5;
View[2].MaxX = 550000;
View[2].MaxY = 5000;
View[2].Min = -28045.19971781372;
View[2].MinVisible = 24.3;
View[2].MinY = -5000;
View[2].MinZ = -280000;
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";
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 = 664454.69992007;
View[3].MaxX = 550000;
View[3].MaxY = 5000;
View[3].Min = -441810.6504328767;
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(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)
......@@ -21,7 +21,7 @@ equations.set_additional_artificial_horizontal_shear(data_dir+'additional_shear/
#equations.set_horizontal_diffusivity('okubo')
#equations.set_lax_friedrichs(.1)
equations.set_bottom_friction(True, z0B=0.005, z0S=0.02)
equations.set_limiter(False)
equations.set_limiter(True)
equations.set_initial_temperature('netcdf', (data_dir+'initial_temperature.nc','T_init'))
equations.set_initial_salinity('vertical_gradient', None, 0, 0.)
equations.set_wind_stress('stress', (data_dir+'windStress.nc','u'), (data_dir+'windStress.nc','v') )
......
cd output
#gmsh mesh2dxz_show.msh temperature/temperature.idx ../opt.opt
#gmsh mesh2dxz_show.msh uv/uv_COMP_0.idx ../opt.opt
#gmsh mesh2dxz_show.msh uv/uv_COMP_0.idx w/w.idx ../opt.opt
gmsh mesh2dxz_show.msh uv/uv_COMP_0.idx w/w.idx temperature/temperature.idx z/z.idx ../optAdapt.opt
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/tanganyika/tanganyika_v6.msh')
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')
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.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"
mesh = slimPre.Mesh(mesh_file_name, mesh_proj='+proj=utm +ellps=WGS84 +zone=35 +south')
region_global = slimPre.Region(mesh)
lonlat_global = slimPre.Coordinate_system(region_global, data_proj='+proj=latlong +ellps=WGS84')
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)])
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, pre_data_dir_base+'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(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)])
print( 'Preprocessing done ' )
slimPre.exit(0)
import slim3d
import slimPre
import shutil
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', reference_density=997.22441884)
equations = slim3d.Slim3d_equations(domain, temperature=True)
equations.set_implicit_vertical_diffusion(True)
equations.set_vertical_viscosity('gotm')
equations.set_vertical_diffusivity('gotm')
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_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_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_coriolis((pre_data_dir+'coriolis.nc', 'coriolis'))
equations.set_boundary_coast('coast')
time_loop = slim3d.Loop(equations,
time_step=600,#1800,
export_time_step=3600,#86400,
final_time=86400*50,
output_directory=output_dir)
time_loop.export_elevation()
time_loop.export_temperature()
time_loop.export_uv()
time_loop.export_uv2d()
time_loop.export_w()
time_loop.loop()
from extrude import *
import sys, os
nPartition = 4
#os.system("gmsh -2 squareSponge.geo -part %i" % (nPartition))
#
#mesh2d = "squareSponge.msh"
#mesh3d = "squareSponge_20layers.msh"
os.system("gmsh -2 square.geo -part %i" % (nPartition))
mesh2d = "square.msh"
mesh3d = "square_20layers.msh"
h = 20
nbLayers = 20
#zLayers = [0, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000]
def getLayers(e, v) :
return [z * h /nbLayers for z in range(nbLayers + 1)]
extrude(mesh(mesh2d), getLayers).write(mesh3d)
# .-------------------------------.
# | Description |
# '-------------------------------'
# Tartinville test-case
import numpy as np