slim3d: tanga2dxz adapt

benchmarks/swe3d/tanganyika/2dxz/prepro.py

 import slimPre
 import numpy as np
 
-data_dir = 'data/'
-slimPre.make_directory(data_dir)
+pre_data_dir_base = 'data_050_'+slimPre.partition_nb()+'/'
+pre_data_dir = pre_data_dir_base+slimPre.partition_id()+'/'
+slimPre.make_directory(pre_data_dir)
 
 
-mesh_file = slimPre.fetch_ftp('slim_data/tanganyika/tanga1dx.msh')
+mesh_file = slimPre.fetch_ftp('slim_data/tanganyika/tanga1dx_fine.msh', pre_data_dir_base)
+mesh_file = slimPre.partition_mesh(mesh_file, pre_data_dir_base+'mesh.msh')
 
 print('Extruding mesh')
 def shiftOperation(node, iPerBound) :
-  n = [node[0], node[1] - 1e4, node[2]]
+  n = [node[0], node[1] - 5e3, node[2]]
   return n
 cutTags = ["cut"]
 pasteTags = ["paste"]
-mapFilename = data_dir+"periodicMesh.txt"
+mapFilename = pre_data_dir_base+"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')
+bath_file_name = slimPre.fetch_ftp('slim_data/tanganyika/1d_bath.txt', pre_data_dir_base)
 f = open(bath_file_name, 'r');
 l = f.readline()
 data = []
@@ -46,14 +43,35 @@ for i in range(bath.shape[0]):
     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')
+slimPre.write_file(pre_data_dir+'bath_2d.nc', region=region_global, time=None, data=[('bath',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)
+
+def layersForZ(h) :
+    for i in range(len(zLayers)) :
+        if zLayers[i] > h : 
+            return i 
+    return len(zLayers)
+
+def layers_func(x, y, h):
+    #if h < 100:
+    #    if h < 30: print('bug')
+    #    zLay = zLayers[:3]
+    #    for i in range(1,11):
+    #        zLay.append(min(zLayers[i+3], zLayers[2] + (h-zLayers[2]) * i/10.) )
+    #else:
+    #    idLayMax = layersForZ(h)
+    #    zLay = [zLayers[i] for i in range(idLayMax)]
+    #return zLay
+    #h2 = max(h,100)
+    idLayMax = layersForZ(h)
+    zLay = [zLayers[i] for i in range(idLayMax)]
+    return zLay
+
+mesh3d_file = pre_data_dir_base+'tanga2dxz.msh'
+slimPre.extrude(mesh_file, (pre_data_dir+'bath_2d.nc','bath'), layers_function=layers_func, mesh_file_name_out=mesh3d_file, factor_show=200, periodicity=periodicity)
 
 print('Loading 3D mesh')
 mesh = slimPre.Mesh(mesh3d_file)
@@ -65,71 +83,116 @@ 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()
+slimPre.write_file(pre_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')
+from dgpy import *
+import slim3d
+domain = slim3d.Domain(mesh3d_file)
+slimSolver = domain._slimSolver
+
+layerDof = dgDofContainer(slimSolver.groups3d, 1)
+tmpDof = dgDofContainer(slimSolver.groups3d, 1)
+layerDof.setAll(0.)
+dof2d = dgDofContainer(slimSolver.groups2d, 1)
+dof2d.setAll(1)
+tmpDof.setAll(0.)
+nbLayers = 50
+for i in range(nbLayers):
+    slimSolver.copy2d3d.copyTo3D(dof2d, tmpDof, -2*i)
+    slimSolver.copy2d3d.copyTo3D(dof2d, tmpDof, -2*i+1)
+    layerDof.axpy(tmpDof, i)
+    tmpDof.setAll(0.)
+def TInitF(cm, lay) :
+    #return cm.get(lay)[:]
+    return 26 - 2.5 * (cm.get(lay)[:] > 7)
+TInitFunc = functorNumpy(lambda cm : TInitF(cm, layerDof.getFunction()))
+tmpDof.interpolate(TInitFunc)
+exporter = dgIdxExporter(tmpDof, pre_data_dir_base+'T_init')
+exporter.exportIdx()
+#print('Preprocessing Initial Temperature')
+#z = region_global.coordinates[:,2]
+#T0 = 26
+#Z0 = -49.
+#Z1 = -51.
+#T1 = 23.5
+#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)])
+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)
+ustress = np.zeros((n))
+wind_file_name = slimPre.fetch_ftp('slim_data/tanganyika/ty_Mpu_93_94.dat', pre_data_dir_base)
+f = open(wind_file_name, 'r') 
+for i in range(n):
+    ustress[i] = - float(f.readline()) * rho0 * 1 #putting ty in ustress is not a typo
+f.close()
+vstress = 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 wind')
+#initial_time = '1993-04-01 00:00:00'
+#time = slimPre.Time(initial_time=initial_time, time_vector=[900*i for i in range(1000000)])
+#uwind = -0.04-0.02 * np.sin(time._time * 2 * np.pi / (2*24*86400))
+#vwind = np.zeros(time._time.shape)
+#slimPre.write_file(data_dir+'wind_stress_sine.nc',region=None,time=time,data=[('u',uwind), ('v',vwind)])
+
 
 print( 'Preprocessing done ' )
 slimPre.exit(0)

benchmarks/swe3d/tanganyika/2dxz/run.py

 import slim3d
-from slimPre import make_directory
+import slimPre
 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')
+pythonFile = __file__
+#output_dir = '/scratch/ucl/mema/delandmeter/tanganyika/2dxz/output' + pythonFile[3:6] + '/'
+#output_dir = './output_' + pythonFile[3:6] + '/'
+output_dir = './output/'
+pre_data_dir_base = 'data_050_'+slimPre.partition_nb()+'/'
+pre_data_dir = pre_data_dir_base+slimPre.partition_id()+'/'
+slimPre.make_directory(output_dir)
+shutil.copyfile(pre_data_dir_base+'tanga2dxz.msh', output_dir+'mesh2dxz.msh')
+shutil.copyfile(pre_data_dir_base+'tanga2dxz_show.msh', output_dir+'mesh2dxz_show.msh')
 
 
-domain = slim3d.Domain(data_dir+'tanga2dxz.msh', data_dir+'periodicMesh.txt', reference_density=997.22441884)
+domain = slim3d.Domain(pre_data_dir_base+'tanga2dxz.msh', pre_data_dir_base+'periodicMesh.txt', reference_density=997.421366)
 
 equations = slim3d.Slim3d_equations(domain, temperature=True)
-equations.set_implicit_vertical_diffusion(True)
-#equations.set_implicit_vertical_advection(True)
-equations.set_vertical_viscosity('gotm')
-equations.set_vertical_diffusivity('gotm')
+equations.set_implicit_vertical(True)
+#equations.set_vertical_viscosity('gotm')
+#equations.set_vertical_diffusivity('gotm')
+equations.set_vertical_viscosity('constant', 1e-2)
+#equations.set_vertical_diffusivity('constant', 1e-2)
 equations.set_horizontal_viscosity('smagorinsky')
-equations.set_additional_artificial_horizontal_shear(data_dir+'additional_shear/additional_shear.idx')
+#equations.set_additional_artificial_horizontal_shear(data_dir+'additional_shear/additional_shear.idx')
 #equations.set_horizontal_diffusivity('okubo')
-equations.set_lax_friedrichs_factor(.5)
+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', (data_dir+'initial_temperature.nc','T_init'))
+#equations.set_initial_temperature('netcdf', (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=10*86400, belowBottom=100)
 equations.set_initial_salinity('vertical_gradient', None, 0, 0.)
-equations.set_wind_stress('stress', (data_dir+'windStress.nc','u'), (data_dir+'windStress.nc','v') )
-equations.set_coriolis((data_dir+'coriolis.nc', 'coriolis'))
+equations.set_wind_stress('stress', (pre_data_dir+'tanga_uniform_windStress.nc','ustress'), (pre_data_dir+'tanga_uniform_windStress.nc','vstress') )
+#equations.set_coriolis((data_dir+'coriolis.nc', 'coriolis'))
 equations.set_boundary_coast('coast')
-#equations.set_vertical_adaptation(tau=300, minimum_height=.5, maximum_height=1500, resize_factor=1.5, background_error=0.001)
+equations.set_vertical_adaptation(tau=100, uv_factor=0, rho_factor=1, minimum_height=.5, maximum_height=1500, resize_factor=1.5, background_error=0.00001, vertical_gradient_factor=.5)
 
 
 time_loop = slim3d.Loop(equations, 
         time_step=600,#1800,
-        export_time_step=3600*24*7,#86400, 
-        final_time=86400*365*3, 
+        export_time_step=3600*1, 
+        initial_time = '1993-04-01 00:00:00',
+        final_time = '1994-04-01 00:00:00',
         output_directory=output_dir)
 time_loop.export_elevation()
 time_loop.export_temperature()
@@ -42,6 +50,8 @@ time_loop.export_uv(False)
 time_loop.export_uv2d(False)
 time_loop.export_w()
 time_loop.export_rho()
+time_loop.export_vertical_viscosity()
+time_loop.export_vertical_diffusivity()
 
 #time_loop.loop()
 time_loop.setup()