Commit 4cc571eb authored by Jonathan Lambrechts's avatar Jonathan Lambrechts Committed by Valentin Vallaeys
Browse files

add thacker test case

parent 87c517a3
R = 430620;
Point(1) = {0,0,0};
Point(2) = {0,-R,0};
Point(3) = {R,0,0};
Point(4) = {0,R,0};
Point(5) = {-R,0,0};
Circle(1)= {2,1,3};
Circle(2)= {3,1,4};
Circle(3)= {4,1,5};
Circle(4)= {5,1,2};
Line Loop(1) = {1,2,3,4};
Plane Surface(1) = {1};
Physical Line("Wall") = {1,2,3,4};
Physical Surface("Basin") = {1};
Field[1] = MathEval;
Field[1].F = Sprintf("Sqrt(x*x+y*y+z*z)/%g",R);
Field[2] = MathEval;
Field[2].F = "10e3 + (1-F1)*90e3";
//Mesh.Algorithm = 6;
Background Field = 2;
Mesh.CharacteristicLengthExtendFromBoundary = 0;
Mesh.CharacteristicLengthFromPoints = 0;
import slimPre
import slim
import numpy
#########################
### global parameters ###
#########################
meshfile = "thacker.msh"
tend = 4*12*3600
dt = 50
h0 = 50
eta0 = 2
R = 430.62e3
g = 9.81
def ssebnd(x,y) :
r = numpy.sqrt(x**2 + y**2)
A = ((h0+eta0)**2 - h0**2)/((h0+eta0)**2 + h0**2)
o = numpy.sqrt(8*g*h0/(R**2))
t = numpy.pi/o
eta = h0*(numpy.sqrt(1-A**2)/(1-A*numpy.cos(o*t)) - 1 - (r**2/R**2)*((1-A**2)/(1-A*numpy.cos(o*t))**2 - 1) )
return eta
def bathf(x,y) :
r2 = x**2 + y**2
return h0 * (R**2 - r2)/(R**2)
hlimt = 1
hlimd = 0.1
linDrag = 0.01
###################
### pre process ###
###################
### mesh ###
mesh = slimPre.Mesh(meshfile)
region = slimPre.Region(mesh)
xyz = region.coordinates
### bathymetry ###
bath = numpy.zeros((xyz.shape[0]))
for i in range(xyz.shape[0]) :
bath[i] = bathf(xyz[i,0], xyz[i,1])
slimPre.write_file('bath.nc', region=region, time=None, data=[('bath', bath)])
slimPre.netcdf_to_msh(meshfile, "bath.nc", "bath", "output/bath")
### init elev ###
etaInit = numpy.zeros((xyz.shape[0]))
for i in range(xyz.shape[0]) :
etaInit[i] = numpy.maximum(ssebnd(xyz[i,0],xyz[i,1]), -bath[i]+.05)
zeros = numpy.zeros_like(etaInit)
slimPre.write_file('init.nc', region=region, time=None, data=[('eta', etaInit),('u', zeros),('v', zeros)])
slimPre.write_file('g.nc', region=None, time=None, data=[('g', g)])
###########
### run ###
###########
domain = slim.Domain(meshfile,('bath.nc','bath'), g=("g.nc","g"))
eq = slim.ShallowWater2dWD(domain, "implicit", final_time=tend, hlimt=hlimt, hlimd=hlimd, linDrag=linDrag)
eq.set_initial_condition(('init.nc','eta'),('init.nc','u'),('init.nc','v'))
eq.set_boundary_coast('Wall')
loop = slim.Loop(maximum_time_step=dt, export_time=600)
loop.add_equation(eq)
loop.run()
Merge "thacker.msh";
Merge "output/sw2d/eta/eta.idx";
Merge "output/bath/bath.idx";
Merge "view.opt";
General.ExecutableFileName = "/home/jonathan/code/gmsh/build/gmsh";
General.FileName = "thacker.msh";
General.RecentFile0 = "thacker.msh";
General.RecentFile1 = "thacker.geo";
General.RecentFile2 = "thaker.msh";
General.RecentFile3 = "circle.geo";
General.RecentFile4 = "florida2.msh";
General.RecentFile5 = "test.pos";
General.RecentFile6 = "mesh.msh";
General.RecentFile7 = "mesh.geo";
General.RecentFile8 = "mod_neuron.msh";
General.RecentFile9 = "mod_neuron.geo";
General.BoundingBoxSize = 1217977.288458204;
General.ColorScheme = 2;
General.FieldPositionX = 667;
General.FieldPositionY = 647;
General.FieldHeight = 680;
General.FieldWidth = 923;
General.GraphicsHeight = 1398;
General.GraphicsPositionX = 430;
General.GraphicsPositionY = 460;
General.GraphicsWidth = 2880;
General.MaxX = 430620;
General.MaxY = 430620;
General.MenuWidth = 303;
General.MenuPositionX = 665;
General.MenuPositionY = 59;
General.MessageHeight = 585;
General.MinX = -430620;
General.MinY = -430620;
General.OptionsPositionX = 349;
General.OptionsPositionY = 613;
General.PluginHeight = 680;
General.PluginWidth = 1004;
General.RotationX = 270;
General.RotationZ = 89.99999999999999;
General.TrackballQuaternion0 = 0.5;
General.TrackballQuaternion1 = -0.4999999999999999;
General.TrackballQuaternion2 = -0.5;
General.TrackballQuaternion3 = 0.5000000000000001;
General.TranslationX = -56060.68669527897;
General.TranslationY = 395874.6952789699;
Mesh.ChacoHypercubeDim = 1;
Mesh.ChacoMeshDim1 = 1;
Mesh.NbNodes = 1495;
Mesh.NbTriangles = 2716;
Mesh.SurfaceEdges = 0;
Mesh.VolumeEdges = 0;
Solver.Executable1 = "test.py";
Solver.Executable2 = "mesh.py";
Solver.Name1 = "test";
Solver.Name2 = "mesh";
PostProcessing.NbViews = 2;
View[0].FileName = "output/sw2d/eta/eta-000103/eta-000103_COMP_0.msh";
View[0].Name = "eta";
View[0].Clip = 1;
View[0].CustomMax = 2;
View[0].CustomMin = -2;
View[0].IntervalsType = 3;
View[0].Light = 0;
View[0].Max = 2.319530919607473;
View[0].MaxVisible = 1.849112426035505;
View[0].MaxX = 430620;
View[0].MaxY = 430620;
View[0].Min = -2.58368102279384;
View[0].MinVisible = -1.92307692307692;
View[0].MinX = -430620;
View[0].MinY = -430620;
View[0].NbIso = 21;
View[0].NbTimeStep = 104;
View[0].RaiseZ = 100000;
View[0].RangeType = 2;
View[0].Time = 0;
View[1].FileName = "output/bath/bath_COMP_0.msh";
View[1].Name = "comp_0";
View[1].Clip = 1;
View[1].IntervalsType = 3;
View[1].Light = 0;
View[1].Max = 50;
View[1].MaxVisible = 50;
View[1].MaxX = 430620;
View[1].MaxY = 430620;
View[1].Min = -2.468612540474626e-14;
View[1].MinVisible = -2.468612540474626e-14;
View[1].MinX = -430620;
View[1].MinY = -430620;
View[1].NbIso = 100;
View[1].RaiseZ = -100000;
View[1].Time = 0;
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