Commit f53d5199 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts

tmp

parent 6c72d461
Pipeline #2415 failed with stage
in 32 minutes and 20 seconds
//*********** box2d.geo *************//
Cx = 10000;
Cy = 2500;
lc = 250;
Point(1) = {-Cx, -Cy, 0, lc};
Point(2) = {Cx , -Cy, 0, lc};
Point(3) = {Cx , Cy , 0, lc};
Point(4) = {-Cx, Cy , 0, lc};
ex = 1.7*lc;
Cx = Cx-ex;
Cy = Cy-ex;
Point(11) = {-Cx, -Cy, 0, lc};
Point(12) = {Cx , -Cy, 0, lc};
Point(13) = {Cx , Cy , 0, lc};
Point(14) = {-Cx, Cy , 0, lc};
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,1};
Line(11) = {11,12};
Line(12) = {12,13};
Line(13) = {13,14};
Line(14) = {14,11};
Physical Line("Coast") = {1,2,3,4};
Line Loop(5) = {1,2,3,4};
Line Loop(15) = {11,12,13,14};
Plane Surface(6) = {5,15};
Plane Surface(7) = {15};
Physical Surface("In") = {7};
Physical Surface("Bnd") = {6};
Field[1]=Attractor;
Field[1].EdgesList={1,2,3,4,11,12,13,14};
Field[1].NNodesByEdge = 100;
Field[2]=Threshold;
Field[2].IField = 1;
Field[2].LcMin = lc;
Field[2].LcMax = 2*lc;
Field[2].DistMin = lc;
Field[2].DistMax = 4*lc;
Background Field = 2;
Mesh.CharacteristicLengthExtendFromBoundary = 0;
Mesh.CharacteristicLengthFromPoints = 0;
import os
import dgpy
import numpy
import time
groups2d = dgpy.dgGroupCollection("box.msh",1,None,2,["top_Sea"])
uvdof = dgpy.dgDofContainer(groups2d, 2)
vorticitydof = dgpy.dgDofContainer(groups2d,1)
idxw = dgpy.dgIdxExporter(vorticitydof,"vorticity")
F = 1
def vorticityNoJumpF(cm) :
ug = cm.get(uvdof.getFunctionGradient())
f = F if cm.elementGroup().getPhysicalTag() == "top_Bnd" else 1
return -ug[:,[1]] + ug[:,[4]]
vorticityNoJump = dgpy.functorNumpy(vorticityNoJumpF)
integralNoJump = dgpy.dgFunctionIntegrator(groups2d, vorticityNoJump)
oldct = 0;
i = 0
while(True):
f = "./output/sw2d/uv/uv-%06d.idx"%i
if i != 0 :
while not os.path.isfile(f) or oldct > os.path.getctime(f):
time.sleep(1)
oldct = os.path.getctime(f)
uvdof.importIdx(f)
vorticitydof.interpolate(vorticityNoJump)
idxw.exportIdx(i,i)
F = 0
i0 = integralNoJump.compute()
F = 1
i1 = integralNoJump.compute()
print("integral ",i,i0,i1)
i += 1
import slim
import slimPre
Ek = 0.05
nu = 0.01
f = 1e-4
H = (nu/(f*Ek))**0.5
D = 1e-5
C = 1e-5
### Prepro ###
mesh2d = slimPre.Mesh("box.msh",dimension=2,tags=["top_Sea"])
slimPre.write_file('bath_2d.nc', time=None, data=[('bath',H)])
slimPre.write_file('coriolis.nc', time=None, data=[('f',f)])
region_global = slimPre.Region(mesh2d)
xyz = region_global.coordinates
taux = (D/2*xyz[:,0]-C/2*xyz[:,1])
tauy = (C/2*xyz[:,0]+D/2*xyz[:,1])
slimPre.write_file('windstress.nc', region=region_global, time=None, data=[('taux',taux),('tauy',tauy)])
### Run ###
domain = slim.Domain('box.msh',("bath_2d.nc","bath"),dimension=2,tags=["top_Sea"])
equations = slim.ShallowWater2d(domain,"DIRK2", final_time=100*24*3600)
equations.set_boundary_coast(['Coast'],slip=True)
equations.set_wind_stress('stress',('windstress.nc','taux'),('windstress.nc','tauy'))
equations.set_coriolis(('coriolis.nc','f'))
equations.set_dissipation("quadratic")
equations.set_viscosity("smagorinsky")
time_loop = slim.Loop(
maximum_time_step=4*3600,
export_time=12*3600,
path='output')
time_loop.add_equation(equations)
time_loop.run()
import slim3d
import slimPre
import numpy
Ek = 0.05
nu = 0.01
f = 1e-4
H = (nu/(f*Ek))**0.5
D = 1e-5
C = 1e-5
nb_layers = 10
print("H %.16g" % H)
print("Ek %.16g" % Ek)
print("nu %.16g" % nu)
print("f %.16g" % f)
print("D %.16g" % D)
print("C %.16g" % C)
print("nb_layers %d" % nb_layers)
### Prepro ###
mesh_file = 'box2d.msh'
mesh2d = slimPre.Mesh(mesh_file)
slimPre.write_file('bath_2d.nc', time=None, data=[('bath',H)])
slimPre.write_file('coriolis.nc', time=None, data=[('f',f)])
slimPre.write_file('eta.nc', time=None, data=[('eta',0)])
slimPre.extrude(mesh_file, ('bath_2d.nc','bath'), nb_layers = nb_layers, mesh_file_name_out='box.msh',factor_show=500)
mesh_file = 'box.msh'
mesh = slimPre.Mesh(mesh_file)
region_global = slimPre.Region(mesh)
xyz = region_global.coordinates
taux = (D/2*xyz[:,0]-C/2*xyz[:,1])
tauy = (C/2*xyz[:,0]+D/2*xyz[:,1])
slimPre.write_file('windstress.nc', region=region_global, time=None, data=[('taux',taux),('tauy',tauy)])
slimPre.netcdf_to_msh('box.msh','windstress.nc','taux', 'taux')
slimPre.netcdf_to_msh('box.msh','windstress.nc','tauy', 'tauy')
### Run ###
domain = slim3d.Domain('box.msh', reference_density=1027)
equations = slim3d.Slim3d_equations(domain, temperature=False, salinity=False)
equations.set_implicit_vertical(True)
equations.set_lax_friedrichs_factor(1.)
equations.set_bottom_friction(True)
equations.set_boundary_coast(['Coast'])
#equations.set_boundary_open(['Coast'])#,eta=("eta.nc","eta"))
equations.set_limiter(False)
#equations.set_initial_elevation("./run2/output/sw2d/eta/eta-000025.idx")
#equations.set_initial_uv_2d("./run2/output/sw2d/uv/uv-000025.idx")
equations.set_wind_stress('stress',('windstress.nc','taux'),('windstress.nc','tauy'))
equations.set_linear_density('z_coordinate', 0, 0)
equations.set_coriolis(('coriolis.nc','f'))
equations.set_vertical_viscosity(mode="constant", constant_value=nu)
equations.set_horizontal_viscosity(mode="smagorinsky")
#equations.set_horizontal_viscosity(mode="constant",constant=0)
time_loop = slim3d.Loop(equations,
time_step=450,
export_time_step=0.5*3600,
final_time=200*24*3600,
output_directory='output')
time_loop.export_elevation()
time_loop.export_uv(True)
time_loop.export_w()
time_loop.export_uv2d(True)
################
# vorticity 2d #
################
import dgpy
groups2d = dgpy.dgGroupCollection("box.msh", 1, None, 2, ["top_In","top_Bnd"])
uvdof2d = dgpy.dgDofContainer(groups2d, 2)
etadof = dgpy.dgDofContainer(groups2d,1)
sys = dgpy.linearSystemPETScDouble()
dof = dgpy.dgDofManager.newCG(groups2d,2,sys)
l2p = dgpy.L2ProjectionContinuous(dof)
ucg = dgpy.dgDofContainer(groups2d,2)
etagradf = dgpy.dgFEGradient(etadof.getFunction())
etagradf.setBoundarySymmetry("Coast")
etagrad = dgpy.dgDofContainer(groups2d,3)
def ugeof(cm) :
eg = cm.get(etagrad.getFunction())
return 9.81/f*numpy.hstack([-eg[:,[1]], eg[:,[0]]])
ugeo = dgpy.functorNumpy(ugeof)
ugeodof = dgpy.dgDofContainer(groups2d,2)
ugeow = dgpy.dgIdxExporter(ugeodof,"output/ugeo",True)
def vorticityNoJumpF(cm,field,F) :
ug = cm.get(field.getFunctionGradient())
FL = F if cm.elementGroup().getPhysicalTag() == "top_Bnd" else 1
return (-ug[:,[1]] + ug[:,[4]])*FL
vfin = dgpy.functorNumpy(lambda cm : vorticityNoJumpF(cm,ucg,0))
vf = dgpy.functorNumpy(lambda cm : vorticityNoJumpF(cm,ucg,1))
viin = dgpy.dgFunctionIntegrator(groups2d, vfin)
vitot = dgpy.dgFunctionIntegrator(groups2d, vf)
def computeVorticity(u,v) :
global F
l2p.apply(ucg,u.getFunction())
v.interpolate(vf)
return viin.compute(), vitot.compute()
u2vdof = dgpy.dgDofContainer(groups2d,1)
ugvdof = dgpy.dgDofContainer(groups2d,1)
u2vw = dgpy.dgIdxExporter(u2vdof,"output/vort_u2")
ugvw = dgpy.dgIdxExporter(ugvdof,"output/vort_ug")
def compute_vorticity(i) :
global F
uvdof2d.importIdx("./output/uvAv2d/uvAv2d-%06d.idx"%i)
etadof.importIdx("./output/eta/eta-%06d.idx"%i)
etagradf.computeAllTerms(0,etagrad,etagrad)
etagrad.multiplyByInvMassMatrix()
ugeodof.interpolate(ugeo)
ugeow.exportIdx(i,i)
i2, j2 = computeVorticity(uvdof2d,u2vdof)
ig, jg = computeVorticity(ugeodof,ugvdof)
u2vw.exportIdx(i,i)
ugvw.exportIdx(i,i)
print("2d ",i,i2,j2)
print("geo ",i,ig,jg)
return i2
################
time_loop.setup()
iout = 1
iold = 1.e9
while time_loop.get_time() < time_loop.final_time:
time_loop.advance_one_time_step()
time_loop.check_sanity()
if time_loop.check_export():
time_loop.print_iter_info()
time_loop.export_fields()
inew = compute_vorticity(iout)
if abs(iold-inew) < 1 :
break
iold = inew
iout += 1
time_loop.terminate()
This diff is collapsed.
This diff is collapsed.
General.ExecutableFileName = "/home/jonathan/code/gmsh/build/gmsh";
General.FileName = "box_show.msh";
General.RecentFile0 = "box_show.msh";
General.RecentFile1 = "box.msh";
General.RecentFile2 = "box2d.msh";
General.RecentFile3 = "square_regular.msh";
General.RecentFile4 = "test.pos";
General.RecentFile5 = "square_regular.geo";
General.RecentFile6 = "../torres8_meadow_utm.msh";
General.RecentFile7 = "torres8_meadow_utm_1.msh";
General.RecentFile8 = "outputGoutte4/fluid_00010.msh";
General.RecentFile9 = "outputGoutte4/fluid_00000.msh";
General.BoundingBoxSize = 30413.81265149109;
General.ClipPositionX = 0;
General.ClipPositionY = 632;
General.ColorScheme = 3;
General.FieldHeight = 680;
General.FieldWidth = 923;
General.GraphicsHeight = 1925;
General.GraphicsPositionX = 1;
General.GraphicsPositionY = 129;
General.GraphicsWidth = 1920;
General.MaxX = 10000;
General.MaxY = 2500;
General.MenuWidth = 303;
General.MinX = -10000;
General.MinY = -2500;
General.MinZ = -22360.67977499789;
General.OptionsPositionX = 1716;
General.OptionsPositionY = 830;
General.PluginPositionX = 2008;
General.PluginPositionY = 482;
General.PluginHeight = 680;
General.PluginWidth = 1004;
General.RotationX = 270;
General.RotationZ = 180;
General.TrackballQuaternion0 = 4.329780281177466e-17;
General.TrackballQuaternion1 = -0.7071067811865475;
General.TrackballQuaternion2 = -0.7071067811865476;
General.TrackballQuaternion3 = 4.329780281177467e-17;
Mesh.ChacoHypercubeDim = 0;
Mesh.ChacoMeshDim1 = 1;
Mesh.NbNodes = 6204;
Mesh.NbPrisms = 10260;
Mesh.NbQuadrangles = 1000;
Mesh.NbTriangles = 2052;
Mesh.SurfaceEdges = 0;
Mesh.VolumeEdges = 0;
PostProcessing.NbViews = 2;
View[0].FileName = "output/uv/uv-000002/uv-000002.msh";
View[0].Name = "uv";
View[0].Max = 0.005825227745814053;
View[0].MaxX = 10000;
View[0].MaxY = 2500;
View[0].MinX = -10000;
View[0].MinY = -2500;
View[0].MinZ = -22360.67977499789;
View[0].NbTimeStep = 3;
View[0].Time = 0;
View[0].Visible = 0;
View[1].FileName = "output/uv/uv-000002/uv-000002.msh_Levelset.pos";
View[1].Name = "uv_Levelset";
View[1].Max = 0.005715740961671321;
View[1].MaxVisible = 0.004092803493164712;
View[1].MaxX = 10000;
View[1].MaxY = -2000;
View[1].MinX = -10000;
View[1].MinY = -2000;
View[1].MinZ = -22360.67977499789;
View[1].NbTimeStep = 3;
View[1].TimeStep = 1;
View[1].Time = 3600;
General.ExecutableFileName = "/home/jonathan/code/gmsh/build/gmsh";
General.FileName = "box_show.msh";
General.RecentFile0 = "box_show.msh";
General.RecentFile1 = "box.msh";
General.RecentFile2 = "box2d.msh";
General.RecentFile3 = "box.mshoutput/uvAv2d/uvAv2d-000003.idx";
General.RecentFile4 = "output/uvAv2d/uvAv2d-000003.idx";
General.RecentFile5 = "vorticity/vorticity-000003.idx";
General.RecentFile6 = "test.msh";
General.RecentFile7 = "adapt/mesh.msh";
General.RecentFile8 = "mesh.msh";
General.RecentFile9 = "adapt/fluid_00000.msh";
General.BoundingBoxSize = 30413.81265149109;
General.ClipPositionX = 158;
General.ClipPositionY = 536;
General.ColorScheme = 0;
General.FieldHeight = 680;
General.FieldWidth = 923;
General.GraphicsHeight = 1925;
General.GraphicsPositionX = 0;
General.GraphicsPositionY = 129;
General.GraphicsWidth = 3840;
General.MaxX = 10000;
General.MaxY = 2500;
General.MenuWidth = 303;
General.MinX = -10000;
General.MinY = -2500;
General.MinZ = -22360.67977499789;
General.OptionsPositionX = 238;
General.OptionsPositionY = 1089;
General.PluginPositionX = 2008;
General.PluginPositionY = 482;
General.PluginHeight = 680;
General.PluginWidth = 1004;
General.TranslationX = 15.83262651964833;
General.TranslationY = 3087.362171331638;
General.VisibilityPositionX = 158;
General.VisibilityPositionY = 1027;
Mesh.ChacoHypercubeDim = 0;
Mesh.ChacoMeshDim1 = 1;
Mesh.NbNodes = 1668;
Mesh.NbPrisms = 2196;
Mesh.NbQuadrangles = 300;
Mesh.NbTriangles = 1464;
Mesh.SurfaceEdges = 0;
Mesh.VolumeEdges = 0;
PostProcessing.NbViews = 2;
View[0].FileName = "output/vort_u2/vort_u2-000056/vort_u2-000056_COMP_0.msh";
View[0].Name = "vort_u2";
View[0].CustomMax = 5e-05;
View[0].CustomMin = -5e-05;
View[0].Max = 0.0001446561212635729;
View[0].MaxVisible = 1.76085858125056e-06;
View[0].MaxX = 10000;
View[0].MaxY = 2500;
View[0].Min = -0.0001319136519110398;
View[0].MinVisible = -1.366886649708295e-06;
View[0].MinX = -10000;
View[0].MinY = -2500;
View[0].MinZ = -22360.67977499789;
View[0].NbTimeStep = 57;
View[0].RangeType = 2;
View[0].SaturateValues = 1;
View[0].TimeStep = 1;
View[0].Time = 1;
View[1].FileName = "output/vort_ug/vort_ug-000056/vort_ug-000056_COMP_0.msh";
View[1].Name = "vort_ug";
View[1].CustomMax = 5e-05;
View[1].CustomMin = -5e-05;
View[1].Max = 0.000103490089339881;
View[1].MaxVisible = 7.567041085906988e-06;
View[1].MaxX = 10000;
View[1].MaxY = 2500;
View[1].Min = -9.213736117725958e-05;
View[1].MinVisible = -3.994424432704968e-06;
View[1].MinX = -10000;
View[1].MinY = -2500;
View[1].MinZ = -22360.67977499789;
View[1].NbTimeStep = 57;
View[1].OffsetY = -6000;
View[1].RangeType = 2;
View[1].SaturateValues = 1;
View[1].TimeStep = 1;
View[1].Time = 1;
General.ExecutableFileName = "/home/jonathan/code/gmsh/build/gmsh";
General.FileName = "box.msh";
General.GraphicsFontEngine = "Cairo";
General.RecentFile0 = "box.msh";
General.RecentFile1 = "box2d.geo";
General.RecentFile2 = "box2d.msh";
General.RecentFile3 = "vorticity.msh";
General.RecentFile4 = "box_show.msh";
General.RecentFile5 = "square_regular.msh";
General.RecentFile6 = "test.pos";
General.RecentFile7 = "square_regular.geo";
General.RecentFile8 = "../torres8_meadow_utm.msh";
General.RecentFile9 = "torres8_meadow_utm_1.msh";
General.BackgroundGradient = 0;
General.BoundingBoxSize = 20615.57663515624;
General.ClipPositionX = 1;
General.ClipPositionY = 632;
General.ColorScheme = 0;
General.FieldHeight = 680;
General.FieldWidth = 923;
General.GraphicsFontSize = 30;
General.GraphicsFontSizeTitle = 40;
General.GraphicsHeight = 786;
General.GraphicsPositionX = 630;
General.GraphicsPositionY = 260;
General.GraphicsWidth = 2883;
General.MaxX = 10000;
General.MaxY = 2500;
General.MenuWidth = 468;
General.MinX = -10000;
General.MinY = -2500;
General.MinZ = -44.72135954999579;
General.OptionsPositionX = 1832;
General.OptionsPositionY = 1459;
General.PluginPositionX = 2008;
General.PluginPositionY = 482;
General.PluginHeight = 680;
General.PluginWidth = 1004;
General.ScaleX = 1.393574381816881;
General.ScaleY = 1.393574381816881;
General.ScaleZ = 1.393574381816881;
General.SmallAxes = 0;
General.TranslationX = -1.295017368273875;
General.TranslationY = 567.977869479571;
General.VisibilityPositionX = 158;
General.VisibilityPositionY = 1027;
General.Color.Background = {255,255,255};
Mesh.ChacoHypercubeDim = 1;
Mesh.ChacoMeshDim1 = 1;
Mesh.ColorCarousel = 0;
Mesh.Lines = 1;
Mesh.LineWidth = 3;
Mesh.NbNodes = 17988;
Mesh.NbPrisms = 27470;
Mesh.NbQuadrangles = 2500;
Mesh.NbTriangles = 10988;
Mesh.RemeshParametrization = 0;
Mesh.SurfaceEdges = 0;
Mesh.VolumeEdges = 0;
PostProcessing.NbViews = 1;
View[0].FileName = "vorticity/vorticity-000024/vorticity-000024_COMP_0.msh";
View[0].Name = "vorticity";
View[0].CustomMax = 1.189555032086522e-06;
View[0].Light = 0;
View[0].Max = 7.407563789331132e-05;
View[0].MaxVisible = 1.189555032086522e-06;
View[0].MaxX = 10000;
View[0].MaxY = 2500;
View[0].MinVisible = 7.138565323781235e-10;
View[0].MinX = -10000;
View[0].MinY = -2500;
View[0].MinZ = -44.72135954999579;
View[0].NbTimeStep = 25;
View[0].RangeType = 2;
View[0].TimeStep = 24;
View[0].Time = 24;
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment