Commit 7dd562c2 authored by Philippe Delandmeter's avatar Philippe Delandmeter Committed by Valentin Vallaeys
Browse files

slim3d: adaptive coord. laxFriedrichs proper option

parent 04b71712
......@@ -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',
......
......@@ -21,8 +21,8 @@ General.FieldPositionY = 303;
General.FieldHeight = 601;
General.FieldWidth = 903;
General.GraphicsHeight = 1002;
General.GraphicsPositionX = 1969;
General.GraphicsPositionY = 135;
General.GraphicsPositionX = 49;
General.GraphicsPositionY = 24;
General.GraphicsWidth = 1871;
General.MaxX = 550000;
General.MaxY = 5000;
......@@ -32,8 +32,8 @@ General.MenuPositionY = 624;
General.MessageHeight = 140;
General.MinY = -5000;
General.MinZ = -280000;
General.OptionsPositionX = 2621;
General.OptionsPositionY = 523;
General.OptionsPositionX = 1438;
General.OptionsPositionY = 550;
General.PluginPositionX = 635;
General.PluginPositionY = 314;
General.PluginHeight = 641;
......@@ -60,7 +60,7 @@ 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].FileName = "uv/uv-000013/uv-000013_COMP_0.msh_1";
View[0].GeneralizedRaiseX = "0";
View[0].GeneralizedRaiseY = "0";
View[0].GeneralizedRaiseZ = "v0";
......@@ -71,10 +71,10 @@ View[0].CustomMin = 24.3;
View[0].GeneralizedRaiseFactor = 200;
View[0].GeneralizedRaiseView = 3;
View[0].IntervalsType = 3;
View[0].Max = 70011.96365286827;
View[0].Max = 6.227853439825584;
View[0].MaxX = 550000;
View[0].MaxY = 5000;
View[0].Min = -23049620.62291582;
View[0].Min = -3.622001857057784;
View[0].MinY = -5000;
View[0].MinZ = -280000;
View[0].NbIso = 21;
......@@ -84,7 +84,7 @@ 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].FileName = "w/w-000013/w-000013_COMP_0.msh_1";
View[1].GeneralizedRaiseX = "0";
View[1].GeneralizedRaiseY = "0";
View[1].GeneralizedRaiseZ = "v0";
......@@ -92,10 +92,12 @@ View[1].Name = "w";
View[1].ColormapNumber = 10;
View[1].GeneralizedRaiseFactor = 200;
View[1].GeneralizedRaiseView = 3;
View[1].Max = 5721200.907856209;
View[1].Max = 1.742882729259567e+97;
View[1].MaxVisible = 9.731221743813752e-312;
View[1].MaxX = 550000;
View[1].MaxY = 5000;
View[1].Min = -2878507.451315449;
View[1].Min = -8.714413646297006e+96;
View[1].MinVisible = -4.865610871899465e-312;
View[1].MinY = -5000;
View[1].MinZ = -280000;
View[1].NbTimeStep = 14;
......@@ -104,7 +106,7 @@ 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].FileName = "temperature/temperature-000013/temperature-000013_COMP_0.msh_1";
View[2].GeneralizedRaiseX = "0";
View[2].GeneralizedRaiseY = "0";
View[2].GeneralizedRaiseZ = "v0";
......@@ -112,21 +114,23 @@ View[2].Name = "temperature";
View[2].ColormapNumber = 10;
View[2].GeneralizedRaiseFactor = 200;
View[2].GeneralizedRaiseView = 3;
View[2].Max = 315.2077599878262;
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 = -28045.19971781372;
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";
View[3].FileName = "z/z-000013/z-000013_COMP_0.msh_1";
View[3].GeneralizedRaiseX = "0";
View[3].GeneralizedRaiseY = "0";
View[3].GeneralizedRaiseZ = "v0";
......@@ -134,10 +138,10 @@ View[3].Name = "z";
View[3].ColormapNumber = 10;
View[3].GeneralizedRaiseFactor = 200;
View[3].GeneralizedRaiseView = 3;
View[3].Max = 664454.69992007;
View[3].Max = 0.07048042688673509;
View[3].MaxX = 550000;
View[3].MaxY = 5000;
View[3].Min = -441810.6504328767;
View[3].Min = -1400;
View[3].MinVisible = -1400;
View[3].MinY = -5000;
View[3].MinZ = -280000;
......
......@@ -125,8 +125,9 @@ 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))
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)])
......
......@@ -19,7 +19,7 @@ equations.set_vertical_diffusivity('gotm')
equations.set_horizontal_viscosity('smagorinsky')
equations.set_additional_artificial_horizontal_shear(data_dir+'additional_shear/additional_shear.idx')
#equations.set_horizontal_diffusivity('okubo')
#equations.set_lax_friedrichs(.1)
equations.set_lax_friedrichs_factor(.5)
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'))
......@@ -27,13 +27,13 @@ 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_boundary_coast('coast')
equations.set_conservative_ale(False)#True)
#equations.set_vertical_adaptation(tau=300, minimum_height=.5, maximum_height=1500, resize_factor=1.5, background_error=0.001)
time_loop = slim3d.Loop(equations,
time_step=600,#1800,
export_time_step=3600,#86400,
final_time=86400*50,
export_time_step=3600*24*7,#86400,
final_time=86400*365*3,
output_directory=output_dir)
time_loop.export_elevation()
time_loop.export_temperature()
......
cd output
#cd output_adapt
cd output_adaptLong
#cd output_fixLong
#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
......
......@@ -31,7 +31,8 @@ xyz = region_global.coordinates
bath = 40 - 35 * xyz[:,0] * xyz[:,0] / 1e8
slimPre.write_file(data_dir+'bath_2d.nc', region=region_global, time=None, data=[('bath',bath)])
slimPre.extrude(mesh_file, (data_dir+'bath_2d.nc','bath'), nb_layers=50, mesh_file_name_out=data_dir+'mesh3d.msh', factor_show=200, periodicity=periodicity)
#slimPre.extrude(mesh_file, (data_dir+'bath_2d.nc','bath'), nb_layers=int(50), mesh_file_name_out=data_dir+'mesh3d.msh', factor_show=200, periodicity=periodicity)
slimPre.extrude(mesh_file, (data_dir+'bath_2d.nc','bath'), nb_layers=int(10), mesh_file_name_out=data_dir+'mesh3d.msh', factor_show=200, periodicity=periodicity)
##zLayers = [0, .25, .5, .75, 1, 1.5, 2, 3, 4, 6, 8, 10]
##zLayers += [12 + 2*i for i in range(15)]
#zLayers = [.5 * i for i in range(1000)]
......
......@@ -29,13 +29,26 @@ mesh2d = slimPre.Mesh(mesh_file)
region_global = slimPre.Region(mesh2d)
xyz = region_global.coordinates
bath = 40 - 35 * xyz[:,0] * xyz[:,0] / 1e8
slimPre.write_file(data_dir+'bath_2d.nc', region=region_global, time=None, data=[('bath',bath)])
#slimPre.extrude(mesh_file, (data_dir+'bath_2d.nc','bath'), nb_layers=50, mesh_file_name_out=data_dir+'mesh3d.msh', factor_show=200, periodicity=periodicity)
#zLayers = [0, .25, .5, .75, 1, 1.5, 2, 3, 4, 6, 8, 10]
#zLayers += [12 + 2*i for i in range(15)]
#zLayers = [.5 * i for i in range(1000)]
zLayers = [1 * i for i in range(1000)]
zLayersFine = [1 * i for i in range(1000)]
zLayersCoarse = [1 * i for i in range(5)]
zLayersCoarse += [5+4 * i for i in range(1000)]
#def layersForZ(h) :
# for i in range(len(zLayersCoarse)) :
# if zLayersCoarse[i] > h :
# return zLayersCoarse[i-1]
# return zLayersCoarse[-1]
#
#for i in range(len(bath)):
# bath[i] = layersForZ(bath[i])
slimPre.write_file(data_dir+'bath_2d.nc', region=region_global, time=None, data=[('bath',bath)])
zLayers = zLayersCoarse
slimPre.extrude(mesh_file, (data_dir+'bath_2d.nc','bath'), z_layers=zLayers, mesh_file_name_out=data_dir+'mesh3d.msh', factor_show=200, periodicity=periodicity)
print('Loading 3D mesh')
......
......@@ -17,23 +17,28 @@ equations.set_vertical_viscosity('gotm')
equations.set_vertical_diffusivity('gotm')
equations.set_boundary_coast(['coastL', 'coastR'])
equations.set_bottom_friction(True, z0B=0.005, z0S=0.02)
equations.set_lax_friedrichs_factor(0.05)
equations.set_limiter(False)#True)
equations.set_initial_temperature('vertical_gradient', None, 20, 0.25)
equations.set_initial_salinity('vertical_gradient', None, 35, 0.)
equations.set_wind_stress('stress', (data_dir+'wind.nc','u'), (data_dir+'wind.nc','v') )
equations.set_coriolis((data_dir+'coriolis.nc', 'coriolis'))
equations.set_vertical_adaptation(tau=.5, uv_factor=1, rho_factor=0, minimum_height=.1, maximum_height=15, resize_factor=3, background_error=1e-4)
#equations.set_vertical_adaptation(tau=5, minimum_height=.1, maximum_height=15, resize_factor=3, background_error=1e-3)
time_loop = slim3d.Loop(equations,
time_step=20,
time_step=100,
export_time_step=1800,
final_time=86400*2,
final_time=86400,#*2,
output_directory=output_dir)
time_loop.export_elevation()
time_loop.export_temperature()
time_loop.export_uv(False)
time_loop.export_w()
time_loop.export_z_coordinate()
time_loop.export_uv2d(False)
#time_loop.loop()
......
......@@ -17,23 +17,26 @@ equations.set_vertical_viscosity('gotm')
equations.set_vertical_diffusivity('gotm')
equations.set_boundary_coast(['coastL', 'coastR'])
equations.set_bottom_friction(True, z0B=0.005, z0S=0.02)
equations.set_limiter(False)
equations.set_initial_temperature('vertical_gradient', None, 20, 0.25)
equations.set_lax_friedrichs_factor(0.1)
equations.set_limiter(True)#False)
equations.set_initial_temperature('vertical_gradient', None, 10, 0.25)
equations.set_initial_salinity('vertical_gradient', None, 35, 0.)
equations.set_wind_stress('stress', (data_dir+'wind.nc','u'), (data_dir+'wind.nc','v') )
equations.set_coriolis((data_dir+'coriolis.nc', 'coriolis'))
equations.set_vertical_adaptation(tau=100, minimum_height=.1, maximum_height=15, resize_factor=3, minimum_error=0.00000000002)#1e-5)
time_loop = slim3d.Loop(equations,
time_step=20,
export_time_step=3600,
final_time=86400*2,
time_step=100,
export_time_step=1800,
final_time=86400,#*2,
output_directory=output_dir)
time_loop.export_elevation()
time_loop.export_temperature()
time_loop.export_uv(False)
time_loop.export_w()
time_loop.export_z_coordinate()
time_loop.export_uv2d(False)
#time_loop.loop()
......
General.ExecutableFileName = "/home/delandmeter/sources/gmsh/projects/dg/build/gmsh/gmsh";
General.FileName = "data/mesh3d_show.msh";
General.RecentFile0 = "data/mesh3d_show.msh";
General.RecentFile1 = "output/mesh3d_show.msh";
General.RecentFile2 = "mesh3d_show.msh";
General.RecentFile3 = "../006/meshShow.msh";
General.RecentFile4 = "data_2/0/mesh.msh";
General.RecentFile5 = "../../data/geo/tanganyika_v5.msh";
General.RecentFile6 = "t4.geo";
General.RecentFile7 = "data/mesh3d.msh";
General.RecentFile8 = "box.msh";
General.RecentFile9 = "../boxShow.msh";
General.BoundingBoxSize = 20397.61015413325;
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 = 811;
General.GraphicsPositionX = 1969;
General.GraphicsPositionY = 332;
General.GraphicsWidth = 1371;
General.MaxX = 10000;
General.MaxY = 125;
General.MenuWidth = 179;
General.MenuPositionX = 665;
General.MenuPositionY = 624;
General.MessageHeight = 1;
General.MinX = -10000;
General.MinY = -125;
General.MinZ = -4000;
General.OptionsPositionX = 3347;
General.OptionsPositionY = 487;
General.PluginPositionX = 635;
General.PluginPositionY = 314;
General.PluginHeight = 641;
General.PluginWidth = 1284;
General.RotationX = 270;
General.StatisticsPositionX = 347;
General.StatisticsPositionY = 561;
General.TrackballQuaternion0 = 0.7071067811865475;
General.TrackballQuaternion3 = 0.7071067811865476;
General.TranslationX = 681.2080536912745;
General.TranslationY = 3523.489932885905;
General.VisibilityPositionX = 1340;
General.VisibilityPositionY = 540;
Mesh.ChacoHypercubeDim = 1;
Mesh.ChacoMeshDim1 = 1;
Mesh.NbHexahedra = 240;
Mesh.NbNodes = 648;
Mesh.NbQuadrangles = 646;
PostProcessing.NbViews = 4;
View[0].FileName = "output/eta/eta-000011/eta-000011_COMP_0.msh";
View[0].Name = "eta";
View[0].Max = 0.09859503264417313;
View[0].MaxX = 10000;
View[0].MaxY = 125;
View[0].Min = -0.09859443862851336;
View[0].MinX = -10000;
View[0].MinY = -125;
View[0].MinZ = -4000;
View[0].NbTimeStep = 12;
View[0].OffsetZ = 3000;
View[0].RaiseZ = 10000;
View[0].Time = 0;
View[0].TransformZY = 1;
View[1].FileName = "output/uvAv2d/uvAv2d-000011/uvAv2d-000011_COMP_0.msh";
View[1].Name = "uvAv2d";
View[1].Max = 0.06905589734740324;
View[1].MaxX = 10000;
View[1].MaxY = 125;
View[1].Min = -0.06898112634799605;
View[1].MinX = -10000;
View[1].MinY = -125;
View[1].MinZ = -4000;
View[1].NbTimeStep = 12;
View[1].OffsetZ = 1500;
View[1].Time = 0;
View[1].TransformZY = 1;
View[2].FileName = "output/uv/uv-000011/uv-000011_COMP_0.msh";
View[2].Name = "uv";
View[2].Max = 0.06949757297468448;
View[2].MaxX = 10000;
View[2].MaxY = 125;
View[2].Min = -0.06922376129169332;
View[2].MinX = -10000;
View[2].MinY = -125;
View[2].MinZ = -4000;
View[2].NbTimeStep = 12;
View[2].Time = 0;
View[3].FileName = "output/w/w-000011/w-000011_COMP_0.msh";
View[3].Name = "w";
View[3].ColormapNumber = 10;
View[3].CustomMax = 0.001;
View[3].CustomMin = -0.001;
View[3].IntervalsType = 3;
View[3].Max = 0.01456612224776507;
View[3].MaxVisible = 0.007096412783709145;
View[3].MaxX = 10000;
View[3].MaxY = 125;
View[3].Min = -0.01480159060422249;
View[3].MinVisible = -0.003548206391854618;
View[3].MinX = -10000;
View[3].MinY = -125;
View[3].MinZ = -4000;
View[3].NbIso = 21;
View[3].NbTimeStep = 12;
View[3].OffsetZ = -5000;
View[3].RangeType = 2;
View[3].Time = 0;
......@@ -34,6 +34,7 @@ time = slimPre.Time(time_vector=tvec, periodic = True)
eta = np.ones_like(tvec)
eta[:] = 0.1*np.sin(2*np.pi*tvec[:]/tau)
slimPre.write_file(run_dir+'eta.nc', region=None, time=time, data=[('eta', eta)])
slimPre.write_file(run_dir+'uv.nc', region=None, time=time, data=[('u', eta), ('v', 0*eta)])
print('preprocessing done')
slimPre.exit(0)
......@@ -11,16 +11,20 @@ output_directory = "./output"
domain = slim3d.Domain(meshFile, periodic_map_file=runFile+"periodicMesh.txt")
equations = slim3d.Slim3d_equations(domain, temperature=True)
equations.set_horizontal_viscosity('smagorinsky')
#equations.set_bottom_friction(False)
#equations.set_horizontal_viscosity('smagorinsky')
equations.set_bottom_friction(False)
equations.set_initial_salinity('vertical_gradient', salinity=None, surface_salinity=0, vertical_gradient=0)
equations.set_initial_temperature('vertical_gradient', temperature=None, surface_temperature=30., vertical_gradient=0.1)
equations.set_initial_temperature('vertical_gradient', temperature=None, surface_temperature=0., vertical_gradient=0.)
equations.set_linear_density('temperature', 0, 0)
equations.set_boundary_open("coastL", eta=(runFile+"eta.nc","eta"))
#equations.set_boundary_open("coastL", u=(runFile+"uv.nc","u"), v=(runFile+"uv.nc", "v"))
equations.set_boundary_open("coastR")
equations.set_limiter(False)
time_loop = slim3d.Loop(equations, time_step=20, 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_uv(False)
time_loop.export_uv2d(False)
time_loop.export_w()
time_loop.loop()
gmsh data/mesh3d_show.msh output/eta/eta.idx output/uvAv2d/uvAv2d_COMP_0.idx output/uv/uv_COMP_0.idx output/w/w.idx opt.opt
......@@ -397,13 +397,17 @@ void dgFaceTermAssembler::mapFromInterface (const dgGroupOfFaces &faces, int nFi
}
else {
for (size_t i=0; i < faces.size(); i++) {
int nbConnectionsInFlux = v.size2() / (nFields * faces.size());
for (int iConnection = 0; iConnection < nbConnections; iConnection++) {
const std::vector<int> &closure = faces.closure(i, iConnection);
int nFieldsInProxy = proxies[iConnection]->size2() / faces.elementGroup(iConnection).getNbElements();
const int nFieldsInProxy = proxies[iConnection]->size2() / faces.elementGroup(iConnection).getNbElements();
const int eid = faces.elementId(i, iConnection);
for (int iField=0; iField<nFieldsInProxy; iField++){
for(size_t j =0; j < closure.size(); j++)
(*proxies[iConnection])(closure[j], faces.elementId(i, iConnection)*nFieldsInProxy + iField) +=
v(j, (i*nbConnections+iConnection)*nFields + iField);
double *data = &(*proxies[iConnection])(0,eid*nFields + iField);
for(size_t j =0; j < closure.size(); j++){
double flux = v(j, (i*nbConnectionsInFlux+iConnection)*nFields + iField);
data[closure[j]] += flux;
}
}
}
}
......
......@@ -202,7 +202,8 @@ import re
def read_msh(filename, vals={}) :
try :
with open(filename, "r") as f:
for l in f :
while True:
l = f.readline()
if l.strip() == "$ElementNodeData" :
break
if l.strip() != "$ElementNodeData" :
......
......@@ -47,6 +47,8 @@ class Slim3d_equations :
self._initial_elevation = None
self._initial_salinity = None
self._initial_temperature = None
self._salinity_boundary_relaxation = None
self._temperature_boundary_relaxation = None
self._wind_stress = None
self._boundary_coast = []
self._boundary_open = []
......@@ -54,7 +56,7 @@ class Slim3d_equations :
self._z0 = [0.005, 0.02]
self._linear2d = False
self._linear_solver_2D = True
self._LFF = 0.
self._LFF = 1.
self._additional_shear = None
self._areaFunc = None
self._nuh = None
......@@ -72,6 +74,11 @@ class Slim3d_equations :
if slimSolver.getSolveT() and self._vertical_diffusivity:
slimSolver.setSolveTImplicitVerticalDiffusion(True)
def set_implicit_vertical_advection(self, flag):
slimSolver = self._slimSolver
self._verticalAdvImplicit = flag
slimSolver.setSolveSImplicitVerticalAdvection(flag)
def set_linear_2d_equations(self, flag):
self._linear2d = flag
......@@ -86,15 +93,18 @@ class Slim3d_equations :
if slimSolver.getSolveT():
slimSolver.setFlagTLimiter(flag)
def set_vertical_adaptation(self, tau, minimum_height=1, maximum_height=1e10, minimum_error=0, resize_factor=1e10) :
def set_vertical_adaptation(self, tau, uv_factor, rho_factor, minimum_height=1, maximum_height=1e10, background_error=0, resize_factor=1e10, smoothing_factor=1) :
self._slimSolver.setUseAdaptiveVerticalGrid(True)
self._verticalAdaptation['tau'] = tau
self._verticalAdaptation['uv_factor'] = uv_factor
self._verticalAdaptation['rho_factor'] = rho_factor
self._verticalAdaptation['minimum_height'] = minimum_height
self._verticalAdaptation['maximum_height'] = maximum_height
self._verticalAdaptation['minimum_error'] = minimum_error
self._verticalAdaptation['background_error'] = background_error
self._verticalAdaptation['resize_factor'] = resize_factor
self._verticalAdaptation['smoothing_factor'] = smoothing_factor
def set_lax_friedrichs(self, factor=1.):
def set_lax_friedrichs_factor(self, factor=1.):
self._LFF = factor
def set_conservative_ale(self, flag):
......@@ -180,9 +190,15 @@ class Slim3d_equations :
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_salinity_boundary_relaxation(self, mode, salinity=None, surface_salinity=0, vertical_gradient=0):
self._salinity_boundary_relaxation = (mode, salinity, surface_salinity, 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_temperature_boundary_relaxation(self, mode, temperature=None, surface_temperature=0, vertical_gradient=0, tau=0, belowBottom=0):
self._temperature_boundary_relaxation = (mode, temperature, surface_temperature, vertical_gradient, tau, belowBottom)
# TODO check value of air_density
def set_wind_stress(self, mode, wind_x, wind_y, density_air=1.25):
"""Add a wind stress term in the shallow water equation.
......@@ -225,8 +241,10 @@ class Slim3d_equations :
groups3 = self._slimSolver.groups3d
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')
def functorMerge(allf) :
return dgpy.functorNumpy(lambda cm : slim3d_private.numpy.hstack([cm.get(f) for f in allf]))
uvF = functorMerge([slim_private._load_function(u, groups3), slim_private._load_function(v, groups3)]) if u else None
uv2F = functorMerge([slim_private._load_function(u, groups2), slim_private._load_function(u, groups2)]) if u else None
uv2F = functorMerge([slim_private._load_function(u, groups2), slim_private._load_function(v, groups2)]) if u else None
etaF = slim_private._load_function(eta, groups3) if eta else None
eta2F = slim_private._load_function(eta,groups2) if eta else None
flowF = slim_private._load_function(flux, groups3) if flux else None
......@@ -296,6 +314,7 @@ class Loop:
def export_nuh(self): self._export_nuh = True
def export_uv_nc(self): self._export_uv_nc = True
def export_salinity_nc(self): self._export_S_nc = True
def export_temperature_nc(self): self._export_T_nc = True
def export_elevation_nc(self): self._export_eta_nc = True
def export_z_nc(self): self._export_z_nc = True
......@@ -370,14 +389,19 @@ class Loop:
self._ncWriter2d = dgpy.slimNetCDFIO(slimSolver.groups2d)
self._ncWriter2d.setInitTime(self._initial_date.year,self._initial_date.month,self._initial_date.day,self._initial_date.hour,self._initial_date.minute,self._initial_date.second,0,0)
self._ncWriter3d = dgpy.slim3dNetCDFIO(slimSolver.columnInfo,self._ncWriter2d)
if self._export_eta_nc: self._ncExport.append([d.etaDof2d, self._odir + '/eta','eta'])
if self._export_eta_nc: self._ncExport.append([d.etaDof2dExport, self._odir + '/eta','eta'])
if self._export_uv_nc: self._ncExport.append([d.uvDof, self._odir + '/uv','uv'])
if self._export_z_nc: self._ncExport.append([d.zDof, self._odir + '/z','z'])
if self._export_z_nc: self._ncExport.append([d.zDofExport, self._odir + '/z','z'])
if self._export_S_nc:
if slimSolver.getSolveS():
self._ncExport.append([d.SDof, self._odir + '/S','S'])
else:
dgpy.Msg.Warning("Salinity is not solved. It won't be exported")
dgpy.Msg.Warning("Salinity is not solved. It won't be exported")
if self._export_T_nc:
if slimSolver.getSolveT():
self._ncExport.append([d.TDof, self._odir + '/T','T'])
else:
dgpy.Msg.Warning("Temperature is not solved. It won't be exported")
else:
self._ncWriter2d = None
self._ncWriter3d = None
......
......@@ -27,7 +27,7 @@ class OpenBoundary : # etaF can be None
groups2 = self.solver.groups2d
self.hF = self.solver.functions.bathFunc2d
if self.flowF :
self.section = dgpy.dgFunctionIntegratorInterface(groups2, hF).compute(self.tag)
self.section = dgpy.dgFunctionIntegratorInterface(groups2, self.hF).compute(self.tag)
def compute_uv_ext(self, cm, uvin, etain) :
dim = cm.elementGroup().getDimUVW()
if self.uvF and not self.transport :
......@@ -42,32 +42,35 @@ class OpenBoundary : # etaF can be None
return UV/(h+eta)
elif self.flowF :
flowF = self.flowF if dim == 3 else self.flow2F
return -cm.get(flow)/self.section * h/(h+eta) * cm.get(dgpy.function.getNormals())
return -cm.get(flowF)/self.section * h/(h+eta) * cm.get(dgpy.function.getNormals())
else :
n = cm.get(dgpy.function.getNormals())
un = uvin[:,[0]] * n[:,[0]] + uvin[:,[1]] * n[:,[1]]
#todo : this is not correct
return un * n
def compute_eta_ext(self, cm, uvin, etain) : #only 2d
return un * n[:,0:2]
def compute_eta_ext(self, cm, uvin, etain) :
dim = cm.elementGroup().getDimUVW()
if self.etaF :