Commit 56868c17 authored by Philippe Delandmeter's avatar Philippe Delandmeter
Browse files

dgpy interface 3D is now running

parent 4eec5991
$MeshFormat
2.2 0 8
$EndMeshFormat
$PhysicalNames
5
1 2 "open"
1 3 "river"
1 4 "cut"
1 5 "paste"
2 1 "sea"
$EndPhysicalNames
$Nodes
52
1 -7500 -250 0
2 7500 -250 0
3 7500 250 0
4 -7500 250 0
5 -6899.999999999091 -250 0
6 -6299.999999998182 -250 0
7 -5699.999999997273 -250 0
8 -5099.999999996365 -250 0
9 -4499.999999995456 -250 0
10 -3899.999999994981 -250 0
11 -3299.999999996367 -250 0
12 -2699.999999999097 -250 0
13 -2100.000000001827 -250 0
14 -1500.000000004557 -250 0
15 -900.0000000072878 -250 0
16 -300.0000000100181 -250 0
17 299.9999999890724 -250 0
18 899.9999999899846 -250 0
19 1499.999999990894 -250 0
20 2099.999999991804 -250 0
21 2699.999999992715 -250 0
22 3299.999999993624 -250 0
23 3899.999999994536 -250 0
24 4499.999999995447 -250 0
25 5099.999999996357 -250 0
26 5699.999999997268 -250 0
27 6299.999999998179 -250 0
28 6899.999999999089 -250 0
29 6899.999999999089 250 0
30 6299.999999998179 250 0
31 5699.999999997268 250 0
32 5099.999999996357 250 0
33 4499.999999995447 250 0
34 3899.999999994536 250 0
35 3299.999999993624 250 0
36 2699.999999992715 250 0
37 2099.999999991804 250 0
38 1499.999999990894 250 0
39 899.9999999899846 250 0
40 299.9999999890724 250 0
41 -300.0000000100181 250 0
42 -900.0000000072878 250 0
43 -1500.000000004557 250 0
44 -2100.000000001827 250 0
45 -2699.999999999097 250 0
46 -3299.999999996367 250 0
47 -3899.999999994981 250 0
48 -4499.999999995456 250 0
49 -5099.999999996365 250 0
50 -5699.999999997273 250 0
51 -6299.999999998182 250 0
52 -6899.999999999091 250 0
$EndNodes
$Elements
77
1 1 2 5 1 1 5
2 1 2 5 1 5 6
3 1 2 5 1 6 7
4 1 2 5 1 7 8
5 1 2 5 1 8 9
6 1 2 5 1 9 10
7 1 2 5 1 10 11
8 1 2 5 1 11 12
9 1 2 5 1 12 13
10 1 2 5 1 13 14
11 1 2 5 1 14 15
12 1 2 5 1 15 16
13 1 2 5 1 16 17
14 1 2 5 1 17 18
15 1 2 5 1 18 19
16 1 2 5 1 19 20
17 1 2 5 1 20 21
18 1 2 5 1 21 22
19 1 2 5 1 22 23
20 1 2 5 1 23 24
21 1 2 5 1 24 25
22 1 2 5 1 25 26
23 1 2 5 1 26 27
24 1 2 5 1 27 28
25 1 2 5 1 28 2
26 1 2 2 2 2 3
27 1 2 4 3 3 29
28 1 2 4 3 29 30
29 1 2 4 3 30 31
30 1 2 4 3 31 32
31 1 2 4 3 32 33
32 1 2 4 3 33 34
33 1 2 4 3 34 35
34 1 2 4 3 35 36
35 1 2 4 3 36 37
36 1 2 4 3 37 38
37 1 2 4 3 38 39
38 1 2 4 3 39 40
39 1 2 4 3 40 41
40 1 2 4 3 41 42
41 1 2 4 3 42 43
42 1 2 4 3 43 44
43 1 2 4 3 44 45
44 1 2 4 3 45 46
45 1 2 4 3 46 47
46 1 2 4 3 47 48
47 1 2 4 3 48 49
48 1 2 4 3 49 50
49 1 2 4 3 50 51
50 1 2 4 3 51 52
51 1 2 4 3 52 4
52 1 2 3 4 4 1
53 3 2 1 1 2 3 29 28
54 3 2 1 1 28 29 30 27
55 3 2 1 1 52 4 1 5
56 3 2 1 1 51 52 5 6
57 3 2 1 1 50 51 6 7
58 3 2 1 1 49 50 7 8
59 3 2 1 1 48 49 8 9
60 3 2 1 1 24 25 32 33
61 3 2 1 1 25 26 31 32
62 3 2 1 1 26 27 30 31
63 3 2 1 1 47 48 9 10
64 3 2 1 1 23 24 33 34
65 3 2 1 1 16 17 40 41
66 3 2 1 1 46 47 10 11
67 3 2 1 1 22 23 34 35
68 3 2 1 1 43 44 13 14
69 3 2 1 1 19 20 37 38
70 3 2 1 1 45 46 11 12
71 3 2 1 1 44 45 12 13
72 3 2 1 1 41 42 15 16
73 3 2 1 1 17 18 39 40
74 3 2 1 1 20 21 36 37
75 3 2 1 1 21 22 35 36
76 3 2 1 1 42 43 14 15
77 3 2 1 1 18 19 38 39
$EndElements
Plugin(MathEval).Expression0 = "v0";
Plugin(MathEval).View = 1;
Plugin(MathEval).Run;
Plugin(MathEval).Expression0 = "v0";
Plugin(MathEval).View = 2;
Plugin(MathEval).Run;
Plugin(MathEval).Expression0 = "v0+w0";
Plugin(MathEval).View = 0;
Plugin(MathEval).OtherView = 4;
//Plugin(MathEval).OtherTimeStep = 0;
Plugin(MathEval).Run;
Plugin(MathEval).Expression0 = "v0*w0";
Plugin(MathEval).View = 7;
Plugin(MathEval).OtherView = 5;
Plugin(MathEval).Run;
General.ExecutableFileName = "/home/delandmeter/sources/gmsh/projects/dg/build/gmsh/gmsh";
General.FileName = "boxShow.msh";
General.RecentFile0 = "sigmaShow.msh";
General.RecentFile1 = "sigmaShow.msh";
General.RecentFile2 = "sigmaShow.msh";
General.RecentFile3 = "sigmaShow.msh";
General.RecentFile4 = "sigmaShow.msh";
General.RecentFile5 = "box.msh";
General.RecentFile6 = "boxSigma.msh";
General.RecentFile7 = "sigmaShow.msh";
General.RecentFile8 = "sigmaShow.msh";
General.RecentFile9 = "sigmaShow.msh";
General.BoundingBoxSize = 15819.29201955637;
General.ClipPositionX = 2239;
General.ClipPositionY = 285;
General.ColorScheme = 0;
General.ContextPositionX = 1485;
General.FieldPositionX = 614;
General.FieldPositionY = 300;
General.FieldHeight = 601;
General.FieldWidth = 903;
General.GraphicsHeight = 1002;
General.GraphicsPositionX = 1920;
General.GraphicsPositionY = 108;
General.GraphicsWidth = 1871;
General.MaxX = 7500;
General.MaxY = 250;
General.MenuWidth = 167;
General.MessageHeight = 836;
General.MinX = -7500;
General.MinY = -250;
General.MinZ = -5000;
General.OptionsPositionX = 86;
General.OptionsPositionY = 537;
General.PluginPositionX = 339;
General.PluginPositionY = 183;
General.PluginHeight = 641;
General.PluginWidth = 1284;
General.RotationX = 270;
General.ScaleX = 0.738123063193709;
General.ScaleY = 0.738123063193709;
General.ScaleZ = 0.738123063193709;
General.StatisticsPositionX = 347;
General.StatisticsPositionY = 561;
General.TrackballQuaternion0 = 0.7071067811865475;
General.TrackballQuaternion3 = 0.7071067811865476;
General.TranslationX = 56.58661247778832;
General.TranslationY = 3566.951958200472;
General.VisibilityPositionX = 175;
General.VisibilityPositionY = 533;
Mesh.ChacoHypercubeDim = 1;
Mesh.ChacoMeshDim1 = 1;
Mesh.NbHexahedra = 300;
Mesh.NbNodes = 676;
Mesh.NbQuadrangles = 674;
PostProcessing.NbViews = 9;
View[0].FileName = "output/eta/eta-000117/eta-000117_COMP_0.msh";
View[0].Name = "eta";
View[0].CenterGlyphs = 1;
View[0].Max = 0.6094007853020991;
View[0].MaxX = 7500;
View[0].MaxY = 250;
View[0].Min = -0.6092751716751482;
View[0].MinX = -7500;
View[0].MinY = -250;
View[0].MinZ = -5000;
View[0].NbTimeStep = 118;
View[0].OffsetZ = 2000;
View[0].RaiseZ = 300;
View[0].RangeType = 3;
View[0].Time = 0;
View[0].TransformYY = 0;
View[0].TransformZY = 0.3;
View[0].Color.Background2D = {255,255,255};
View[1].FileName = "output/uvAv2d/uvAv2d-000117/uvAv2d-000117.msh";
View[1].Name = "uvAv2d";
View[1].CenterGlyphs = 1;
View[1].Max = 0.1905701320914563;
View[1].MaxX = 7500;
View[1].MaxY = 250;
View[1].MinX = -7500;
View[1].MinY = -250;
View[1].MinZ = -5000;
View[1].NbTimeStep = 118;
View[1].Time = 0;
View[1].Visible = 0;
View[1].Color.Background2D = {255,255,255};
View[2].FileName = "output/uv/uv-000117/uv-000117.msh";
View[2].Name = "uv";
View[2].CenterGlyphs = 1;
View[2].Max = 0.1957123578257393;
View[2].MaxX = 7500;
View[2].MaxY = 250;
View[2].MinX = -7500;
View[2].MinY = -250;
View[2].MinZ = -5000;
View[2].NbTimeStep = 118;
View[2].Time = 0;
View[2].Visible = 0;
View[2].Color.Background2D = {255,255,255};
View[3].FileName = "output/w/w-000117/w-000117_COMP_0.msh";
View[3].Name = "w";
View[3].CenterGlyphs = 1;
View[3].IntervalsType = 3;
View[3].Max = 0.007010609259835515;
View[3].MaxX = 7500;
View[3].MaxY = 250;
View[3].Min = -0.006080809138238502;
View[3].MinX = -7500;
View[3].MinY = -250;
View[3].MinZ = -5000;
View[3].NbIso = 20;
View[3].NbTimeStep = 118;
View[3].OffsetZ = -6000;
View[3].RangeType = 3;
View[3].Time = 0;
View[3].Color.Background2D = {255,255,255};
View[4].FileName = "bath2d_COMP_0.msh";
View[4].Name = "h";
View[4].CenterGlyphs = 1;
View[4].Max = 100;
View[4].MaxVisible = 100;
View[4].MaxX = 7500;
View[4].MaxY = 250;
View[4].Min = 100;
View[4].MinVisible = 100;
View[4].MinX = -7500;
View[4].MinY = -250;
View[4].MinZ = -5000;
View[4].OffsetZ = 1000;
View[4].RangeType = 3;
View[4].Time = 0;
View[4].TransformYY = 0;
View[4].TransformZY = 0.3;
View[4].Visible = 0;
View[4].Color.Background2D = {255,255,255};
View[5].FileName = "uvAv2d_MathEval.pos";
View[5].Name = "u2d";
View[5].CenterGlyphs = 1;
View[5].Max = 0.190354967371036;
View[5].MaxX = 7500;
View[5].MaxY = 250;
View[5].Min = -0.1905701320914563;
View[5].MinX = -7500;
View[5].MinY = -250;
View[5].NbTimeStep = 118;
View[5].OffsetZ = 1000;
View[5].RaiseZ = 200;
View[5].RangeType = 3;
View[5].Time = 0;
View[5].TransformZY = 0.3;
View[5].Color.Background2D = {255,255,255};
View[6].FileName = "uv_MathEval.pos";
View[6].Name = "u";
View[6].MaxRecursionLevel = 1;
View[6].Max = 0.1955898729076838;
View[6].MaxX = 7500;
View[6].MaxY = 250;
View[6].Min = -0.1957123578257393;
View[6].MinX = -7500;
View[6].MinY = -250;
View[6].MinZ = -5000;
View[6].NbTimeStep = 118;
View[6].RangeType = 3;
View[6].TargetError = -0.0001;
View[6].Time = 0;
View[7].FileName = "eta_MathEval.pos";
View[7].Name = "H";
View[7].Max = 100.6094007853021;
View[7].MaxVisible = 100;
View[7].MaxX = 7500;
View[7].MaxY = 250;
View[7].Min = 99.39072482832485;
View[7].MinVisible = 100;
View[7].MinX = -7500;
View[7].MinY = -250;
View[7].NbTimeStep = 118;
View[7].OffsetZ = 9000;
View[7].RangeType = 3;
View[7].Time = 0;
View[7].TransformZY = 0.3;
View[7].Visible = 0;
View[8].FileName = "eta_MathEval_MathEval.pos";
View[8].Name = "Hu2d";
View[8].Max = 19.15123429581403;
View[8].MaxX = 7500;
View[8].MaxY = 250;
View[8].Min = -18.94113245423248;
View[8].MinX = -7500;
View[8].MinY = -250;
View[8].NbTimeStep = 118;
View[8].OffsetZ = 500;
View[8].RaiseZ = 20;
View[8].RangeType = 3;
View[8].Time = 0;
View[8].TransformZY = 0.3;
View[8].TransformZZ = 0;
View[8].Visible = 0;
import slimPre
import numpy as np
mesh_file_name = 'box.msh'
### Bathymetry Process ###
bath = 50
slimPre.write_file('bathymetry.nc', data=[('bathymetry', bath)])
### Extrusion Process ###
bath_file_name = ('bathymetry.nc', 'bathymetry')
nb_layers = 7
factor_show = 100
def shiftOperation(node, iPerBound) :
n = [node[0], node[1] - 2e3/4., node[2]]
return n
cutTags = ["cut"]
pasteTags = ["paste"]
mapFilename = "periodicMesh.txt"
periodicity = (shiftOperation, cutTags, pasteTags, mapFilename)
slimPre.extrude(mesh_file_name, bath_file_name, nb_layers=nb_layers, factor_show=factor_show, periodicity=periodicity)
### Initial Elevation Process ###
mesh_3d = mesh_file_name[0:-4] + '_3d.msh'
mesh = slimPre.Mesh(mesh_3d)
region_global = slimPre.Region(mesh, '')
x = region_global.coordinates[:,0]
eta = 1 - np.minimum((np.absolute(x)/2000), 1)
slimPre.write_file('initialCondition.nc', region=region_global, time=None, data=[('elevation', eta)])
### Initial Elevation Process ###
t = [10*j for j in range(1000)]
t = np.array(t)
time = slimPre.Time(time_vector=t)
eta = 0.6 * np.sin(t/100)
slimPre.write_file('wave.nc', region=None, time=time, data=[('eta', eta)])
print('done')
slimPre.exit(0)
import slim3d
from slimPre import make_directory, exit
output_directory = 'outputNew'
make_directory(output_directory)
domain = slim3d.Domain('box_3d.msh', 'periodicMesh.txt')
equations = slim3d.Slim3d_equations(domain, temperature=True)
equations.set_linear_density(temperature=True, constant_coefficient=10, linear_coefficient=-1000*1.94552591e-4)
equations.set_initial_temperature('vertical_gradient', None, 30, 0.1)
#equations.set_initial_elevation(('initialCondition.nc', 'elevation'))
#equations.set_boundary_coast(['river'])
equations.set_boundary_open('open', temperature='vertical_gradient')
equations.set_boundary_open('river', eta=('wave.nc', 'eta'), temperature='vertical_gradient')
equations.set_bottom_friction(False)
equations.set_limiter(False)
final_time = 86400*2
time_loop = slim3d.Loop(equations, 1, 1, 20, final_time=final_time, output_directory=output_directory)
time_loop.export_elevation()
time_loop.export_temperature()
time_loop.export_uv()
time_loop.export_uv2d()
time_loop.setup()
time_loop.loop()
#gmsh box_3d_show.msh output/eta/eta.idx output/uvAv2d/uvAv2d.idx output/uv/uv.idx output/T/T.idx bath2d_COMP_0.msh checkPluginTransport optTransport.opt
gmsh box_3d_show.msh outputNew/eta/eta.idx outputNew/uvAv2d/uvAv2d.idx outputNew/uv/uv.idx outputNew/temperature/temperature.idx bath2d_COMP_0.msh checkPluginTransport optTransport.opt
import dgpy
import slim_private
class Domain :
"""Create the slim3d Domain"""
def __init__(self, mesh_file_name, periodic_map_file='') :
"""keyword arguments:
* mesh_file_name
path to the mesh file (.msh format). The partitioned mesh will automatically be loaded in multiprocessing
"""
(topTags, bottomTags) = slim_private._findTopAndBottomTags(mesh_file_name)
self._topTags = topTags
self._bottomTags = bottomTags
self._slimSolver = dgpy.slim3dSolver(mesh_file_name, bottomTags, topTags, periodic_map_file)
class Slim3d_equations :
"""Create the slim3d Equation"""
def __init__(self, Domain, salinity=False, temperature=False) :
self._domain = Domain
self._slimSolver = Domain._slimSolver
slimSolver = self._slimSolver
slimSolver.setSolveS(salinity)
slimSolver.setSolveT(temperature)
slimSolver.setSolveUVImplicitVerticalDiffusion(False)
slimSolver.setSolveSImplicitVerticalDiffusion(False)
slimSolver.setSolveTImplicitVerticalDiffusion(False)
slimSolver.setFlagUVLimiter(True)
if salinity:
slimSolver.setFlagSLimiter(True)
if temperature:
slimSolver.setFlagTLimiter(True)
slimSolver.setComputeBottomFriction(True)
self._horizontal_viscosity = None
self._horizontal_diffusivity = None
self._vertical_viscosity = None
self._vertical_diffusivity = None
self._coriolis = None
self._linear_density = None
self._reference_density = 1027
self._initial_elevation = None
self._initial_salinity = None
self._initial_temperature = None
self._boundary_coast = []
self._boundary_open = []
def set_implicit_vertical_diffusion(self, flag):
slimSolver = self._slimSolver
slimSolver.setSolveUVImplicitVerticalDiffusion(flag)
if slimSolver.getSolveS():
slimSolver.setSolveSImplicitVerticalDiffusion(flag)
if slimSolver.getSolveT():
slimSolver.setSolveTImplicitVerticalDiffusion(flag)
def set_limiter(self, flag):
slimSolver = self._slimSolver
slimSolver.setFlagUVLimiter(flag)
if slimSolver.getSolveS():
slimSolver.setFlagSLimiter(flag)
if slimSolver.getSolveT():
slimSolver.setFlagTLimiter(flag)
def set_horizontal_viscosity(self, mode, constant=1, factor=0.02, maximum=1e9):
self._horizontal_viscosity = mode
self._hor_visc_const = constant
self._hor_visc_fact = factor
self._hor_visc_max = maximum
def set_horizontal_diffusivity(self, mode, constant=1, factor=0.0, maximum=1e9):
self._horizontal_diffusivity = mode
self._hor_diff_const = constant
self._hor_diff_fact = factor
self._hor_diff_max = maximum
def set_vertical_viscosity(self, viscosity):
self._vertical_viscosity = viscosity
def set_vertical_diffusivity(self, diffusivity):
self._vertical_diffusivity = diffusivity
def set_coriolis(self, coriolis):
self._coriolis = coriolis
def set_bottom_friction(self, flag):
self._slimSolver.setComputeBottomFriction(flag)
def set_linear_density(self, salinity=False, temperature=False, constant_coefficient=0, linear_coefficient=0):
if (temperature and salinity) or (not temperature and not salinity):
dgpy.Msg.Fatal('density is a linear function of either salinity or temperature, not both, not none')
self._linear_density = (salinity, constant_coefficient, linear_coefficient)
def set_reference_density(self, density):
self._reference_density = density
def set_initial_elevation(self, elevation):
self._initial_elevation = elevation
def set_initial_salinity(self, mode, salinity, surface_salinity, vertical_gradient):
self._initial_salinity = (mode, salinity, surface_salinity, vertical_gradient)
def set_initial_temperature(self, mode, temperature, surface_temperature, vertical_gradient):
self._initial_temperature = (mode, temperature, surface_temperature, vertical_gradient)
def set_boundary_coast(self, physical_tags):
if slim_private._is_string(physical_tags):
self._boundary_coast.append(physical_tags)
else:
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):
if slim_private._is_string(physical_tags):
openBnd = slim_private.OpenBnd_3d(physical_tags, uv, eta, salinity, temperature)
self._boundary_open.append(openBnd)
else:
for tag in physical_tags:
openBnd = slim_private.OpenBnd_3d(tag, uv, eta, salinity, temperature)
self._boundary_open.append(openBnd)
class Loop:
"""Temporal solver"""
def __init__(self, equations, maximum_2d_time_step=3600, ratio_3dvs2d_time_step=1, export_time_step=3600, initial_time='1970-01-01 00:00:00', final_time=0, output_directory="./output"):
self._slimSolver = equations._slimSolver
self._equations = equations
self._max_2d_dt = maximum_2d_time_step
self._dt2d3d_ratio = ratio_3dvs2d_time_step
self._export_dt = export_time_step
fmt = '%Y-%m-%d %H:%M:%S'
date = slim_private.datetime.datetime.strptime(initial_time, fmt)
self._initial_time = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
if slim_private._is_string(final_time):
date = slim_private.datetime.datetime.strptime(final_time, fmt)
self._final_time = slim_private.calendar.timegm([date.year, date.month, date.day, date.hour, date.minute, date.second])
else:
self._final_time = self._initial_time + final_time
self._odir = output_directory
self._export_uv = False
self._export_w = False
self._export_S = False
self._export_T = False
self._export_uvAv2d = False
self._export_eta = False
self._export_kappav = False
self._export_nuv = False
self._timeIntegrator = None
self._export = []
def export_uv(self): self._export_uv = True
def export_w(self): self._export_w = True
def export_salinity(self): self._export_S = True
def export_temperature(self): self._export_T = True
def export_uv2d(self): self._export_uvAv2d = True
def export_elevation(self): self._export_eta = True
def export_vertical_diffusivity(self): self._export_kappav = True
def export_vertical_viscosity(self): self._export_nuv = True
def setup(self):
slim_private.slim3d_setup(self)
slimSolver = self._slimSolver
d = slimSolver.getDofs()
if self._export_uv: self._export.append(dgpy.dgIdxExporter(d.uvDof, self._odir + "/uv", True))
if self._export_w: self._export.append(dgpy.dgIdxExporter(d.wDof3d, self._odir + "/w"))
if self._export_uvAv2d: self._export.append(dgpy.dgIdxExporter(d.uvAvDof2d, self._odir + "/uvAv2d", True))
if self._export_eta: self._export.append(dgpy.dgIdxExporter(d.etaDof2d, self._odir + "/eta"))
if self._export_S:
if slimSolver.getSolveS():
self._export.append(dgpy.dgIdxExporter(d.SDof, self._odir + "/salinity"))
else:
dgpy.Msg.Warning("Salinity is not solved. It won't be exported")
if self._export_T:
if slimSolver.getSolveT():
self._export.append(dgpy.dgIdxExporter(d.TDof, self._odir + "/temperature"))
else:
dgpy.Msg.Warning("Temperature is not solved. It won't be exported")
if self._export_kappav: