Lab simulation#
This jupyter notebook is a complete example of how to use Watlab. It will guide you through the different steps of a classic Watlab usage: the mesh creation with GMSH, the model description with one of the modules of Watlab, and the visualisation with the visualisation tools of Watlab.
Now make sure you have installed Watlab correctly, and dig in with us!
In this example, we will perform a dam break on the Uclouvain’s lab geometry.
All you need is GMSH!#
First, lets create the mesh :
[1]:
import gmsh
import sys
gmsh.initialize()
gmsh.model.add("canalLab")
#characteristic size of cells
lc = 0.08
# create points
gmsh.model.geo.addPoint(-7.8, 1.8, 0.155, lc, 1)
gmsh.model.geo.addPoint(-7.8, 1.46, 0, lc, 2)
gmsh.model.geo.addPoint(-7.8, -1.46, 0, lc, 3)
gmsh.model.geo.addPoint(-7.8, -1.8, 0.155, lc, 4)
gmsh.model.geo.addPoint(-0.8, 1.8, 0.155, lc, 5)
gmsh.model.geo.addPoint(-0.8, 1.46, 0, lc, 6)
gmsh.model.geo.addPoint(-0.8, 0.5, 0, lc, 7)
gmsh.model.geo.addPoint(-0.8, -0.5, 0, lc, 8)
gmsh.model.geo.addPoint(-0.8, -1.46, 0, lc, 9)
gmsh.model.geo.addPoint(-0.8, -1.8, 0.155, lc, 10)
gmsh.model.geo.addPoint(0, 1.8, 0.155, lc, 11)
gmsh.model.geo.addPoint(0, 1.46, 0, lc, 12)
gmsh.model.geo.addPoint(0, 0.5, 0, lc, 13)
gmsh.model.geo.addPoint(0, -0.5, 0, lc, 14)
gmsh.model.geo.addPoint(0, -1.46, 0, lc, 15)
gmsh.model.geo.addPoint(0, -1.8, 0.155, lc, 16)
gmsh.model.geo.addPoint(28.10, 1.8, 0.155, lc, 17)
gmsh.model.geo.addPoint(28.10, 1.46, 0 , lc, 18)
gmsh.model.geo.addPoint(28.10, -1.46, 0 , lc, 19)
gmsh.model.geo.addPoint(28.10, -1.8, 0.155, lc, 20)
#create lines between those points
gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(2, 3, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(5, 6, 4)
gmsh.model.geo.addLine(6, 7, 5)
gmsh.model.geo.addLine(7, 8, 6)
gmsh.model.geo.addLine(8, 9, 7)
gmsh.model.geo.addLine(9, 10, 8)
gmsh.model.geo.addLine(11, 12, 9)
gmsh.model.geo.addLine(12, 13, 10)
gmsh.model.geo.addLine(13, 14, 11)
gmsh.model.geo.addLine(14, 15, 12)
gmsh.model.geo.addLine(15, 16, 13)
gmsh.model.geo.addLine(17, 18, 14)
gmsh.model.geo.addLine(18, 19, 15)
gmsh.model.geo.addLine(19, 20, 16)
gmsh.model.geo.addLine(1, 5, 17)
gmsh.model.geo.addLine(5, 11, 18)
gmsh.model.geo.addLine(11, 17, 19)
gmsh.model.geo.addLine(2, 6, 20)
gmsh.model.geo.addLine(6, 12, 21)
gmsh.model.geo.addLine(12, 18, 22)
gmsh.model.geo.addLine(7, 13, 23)
gmsh.model.geo.addLine(8, 14, 24)
gmsh.model.geo.addLine(3, 9, 25)
gmsh.model.geo.addLine(9, 15, 26)
gmsh.model.geo.addLine(15, 19, 27)
gmsh.model.geo.addLine(4, 10, 28)
gmsh.model.geo.addLine(10, 16, 29)
gmsh.model.geo.addLine(16, 20, 30)
#add curve loop : must be created at first before creating a surface
# it is defined as an ordered list of connected curves
gmsh.model.geo.addCurveLoop([1, 20, -4, -17], 1)
gmsh.model.geo.addCurveLoop([2, 25, -7, -6, -5, -20], 2)
gmsh.model.geo.addCurveLoop([3, 28, -8, -25], 3)
gmsh.model.geo.addCurveLoop([4, 21, -9, -18], 4)
gmsh.model.geo.addCurveLoop([5, 23, -10, -21], 5)
gmsh.model.geo.addCurveLoop([6, 24, -11, -23], 6)
gmsh.model.geo.addCurveLoop([7, 26, -12, -24], 7)
gmsh.model.geo.addCurveLoop([8, 29, -13, -26], 8)
gmsh.model.geo.addCurveLoop([9, 22, -14, -19], 9)
gmsh.model.geo.addCurveLoop([10, 11, 12, 27, -15, -22], 10)
gmsh.model.geo.addCurveLoop([13, 30, -16, -27], 11)
# create the surface : it is a list of curve loops
gmsh.model.geo.addPlaneSurface([1], 1)
gmsh.model.geo.addPlaneSurface([2], 2)
gmsh.model.geo.addPlaneSurface([3], 3)
gmsh.model.geo.addPlaneSurface([6], 6)
gmsh.model.geo.addPlaneSurface([9], 9)
gmsh.model.geo.addPlaneSurface([10], 10)
gmsh.model.geo.addPlaneSurface([11], 11)
# it is necessary to synchronise those informations with the Gmsh model, creating the relevant data structure
gmsh.model.geo.synchronize()
#entities as domain,boundary,material properties are called physical groups
gmsh.model.addPhysicalGroup(1, [14,15,16],name="Qout")
gmsh.model.addPhysicalGroup(1,[1,3,4,5,7,8,9,10,12,13,17,18,19,23,24,28,29,30],name="Walls")
gmsh.model.addPhysicalGroup(1, [2],name="Qin")
gmsh.model.addPhysicalGroup(1, [11],name="The gate")
gmsh.model.addPhysicalGroup(2, [1,2, 3, 6],name="Before the gate") # 2 = dimension, [1] = tagsof the plane surface
gmsh.model.addPhysicalGroup(2,[9,11],name="Canal slopes")
gmsh.model.addPhysicalGroup(2, [10],name="After the gate")
#Mesh on the surface
gmsh.model.mesh.generate(2)
gmsh.model.mesh.optimize('Laplace2D')
gmsh.write("canalLabMesh.msh")
if '-nopopup' not in sys.argv:
gmsh.fltk.run()
gmsh.finalize()
Watlab helps you master the model#
Now that the geometry of your domain is done, you can focus on solving your hydraulic problem. For this, you need a solver (hopefully, Watlab provides you one!). You also need to tell Watlab what mesh represents the domain, and what conditions are set to the cells and the boundaries.
The example here is only hydrodynamics, so the solver is hydroflow.
In this example, we use a the HydroFlow’s API to model another geometrical case. It consists of a dambreak happening in the UCLouvain’s lab
First, import Watlab, import the mesh, and create a Watlab model!
[2]:
import watlab
mesh = watlab.Mesh("canalLabMesh.msh")
model = watlab.HydroflowModel(mesh)
Several properties can be given to the model (see the documentation for a complete set of properties).
[3]:
model.name = "Canal_Simulation"
model.ending_time = 20
model.Cfl_number = 0.95
Initial conditions and physical properties, such as friction coefficient or water level can be given as follows:
[4]:
model.set_friction_coefficient("Before the gate",0.01)
model.set_friction_coefficient("After the gate",0.01)
model.set_friction_coefficient("Canal slopes",0.04)
model.set_initial_water_level("Before the gate", 0.4)
model.set_initial_water_level("After the gate", 0.165)
model.set_initial_water_level("Canal slopes",0.165)
We add walls, a transmissive boundary condition at the end of the canal and a constant input discharge:
[5]:
model.set_wall_boundaries("Walls")
model.set_transmissive_boundaries("Qout")
model.set_boundary_water_discharge('Qin', 1)
Provide the number of required results files:
[6]:
model.set_picture_times(n_pic = 21)
Finally, we can specify output parameters, generate the data files and run this specific model:
[7]:
model.export_data()
print("Input files generated!")
model.solve()
Input files generated!
hydroflow.exe exists
Launching the executable ...
Total execution time: 3994 sec. for 20 sec. of simulation.
Visualize the results with the Plotter class combined with matplotlib#
[4]:
import matplotlib.pyplot as plt
#create the mesh and plotter object
plotter = watlab.Plotter(mesh)
#time of the plot
time0 = "0_00"
time1 = "1_00"
time2 = "4_00"
time3 = "6_00"
time4 = "10_00"
time5 = "19_00"
#path of the pic at the corresponding times
myPic0 = "output\\pic_"+ time0 +".txt"
myPic1 = "output\\pic_"+ time1 +".txt"
myPic2 = "output\\pic_"+ time2 +".txt"
myPic3 = "output\\pic_"+ time3 +".txt"
myPic4 = "output\\pic_"+ time4 +".txt"
myPic5 = "output\\pic_"+ time5 +".txt"
plotter.plot(myPic0, "h", colorbar_values=[0,0.5])
plt.title("Time: " + time0 +" s")
plt.show()
plotter.plot(myPic1, "h", colorbar_values=[0,0.5])
plt.title("Time: " + time1 +" s")
plt.show()
plotter.plot(myPic2, "h", colorbar_values=[0,0.6])
plt.title("Time: " + time2 +" s")
plt.show()
plotter.plot(myPic3, "h", colorbar_values=[0,0.55])
plotter.show_velocities(myPic3, scale=150, velocity_ds=0.2)
plt.title("Time: " + time3 +" s")
plotter.plot(myPic4, "h", colorbar_values=[0,0.65])
plt.title("Time: " + time4 +" s")
plt.show()
plotter.plot(myPic5, "h", colorbar_values=[0,0.7])
plt.title("Time: " + time5 +" s")
plt.show()
plotter.plot(myPic3, "h", colorbar_values=[0,0.5])
plotter.show_velocities(myPic3, scale=150, velocity_ds=0.2)
plt.title("Zoom at time: " + time3+" s")
plt.xlim([0, 5])
plt.ylim([-1,1])
# #cross section plot
plotter.plot_profile_along_line(myPic0, "h", x_coordinate=[-10, 10], y_coordinate=[0.2,0.5], new_fig = True, label = time0)
plotter.plot_profile_along_line(myPic1, "h", x_coordinate=[-10, 10], y_coordinate=[0.2,0.5], label = time1)
plotter.plot_profile_along_line(myPic2, "h", x_coordinate=[-10, 10], y_coordinate=[0.2,0.5], label = time2)
plotter.plot_profile_along_line(myPic3, "h", x_coordinate=[-10, 10], y_coordinate=[0.2,0.5], label = time3)
plotter.plot_profile_along_line(myPic4, "h", x_coordinate=[-10, 10], y_coordinate=[0.2,0.5], label = time4)
plotter.plot_profile_along_line(myPic4, "h", x_coordinate=[-10, 10], y_coordinate=[0.2,0.5], label = time5)
[ ]: