Note
Click here to download the full example code
Thacker¶
Pseudo-1D thacker simulation.
import os, sys
import argparse
import matplotlib.pyplot as plt
import numpy as np
import freshkiss3d as fk
import freshkiss3d.extra.plots as fk_plt
#sphinx_gallery_thumbnail_number = 2
parser = argparse.ArgumentParser()
parser.add_argument('--nographics', action='store_true')
args = parser.parse_args()
Parameters:¶
A = 1.
H = 0.1
#Numerical scheme:
SECOND_ORDER = True
Time loop:¶
simutime = fk.SimuTime(final_time=4.6, time_iteration_max=2000,
second_order=SECOND_ORDER)
Mesh:¶
dir_path = os.path.abspath(os.path.dirname(sys.argv[0]))
triangular_mesh = fk.TriangularMesh.from_msh_file(dir_path + '/inputs/thacker.msh')
if not args.nographics:
fk_plt.plot_mesh(triangular_mesh)
Topography:¶
Layers:¶
NL = 1
layer = fk.Layer(NL, triangular_mesh, topography_funct=Topography)
Primitives:¶
def eta_0(x, y):
h = -0.5*x + 4.0
return h
def T_0(x, y):
if x < -1.5:
t = 35.0
else:
t = 0.0
return t
#Plot initial free surface in (x,z) slice plane:
if not args.nographics:
fk_plt.plot_init_1d_slice(triangular_mesh, surf_xy=eta_0, topo_xy=Topography, primitive_xy=T_0)
primitives = fk.Primitive(triangular_mesh, layer, free_surface_funct=eta_0)
tracers = [fk.Tracer(triangular_mesh, layer, primitives, Tinit_funct=T_0)]
Boundary conditions:¶
slides = [fk.Slide(ref=1), fk.Slide(ref=3)]
torrential_outs = [fk.TorrentialOut(ref=2)]
tube = fk.RectangularTube(xlim=(0.0, 4.0), ylim=(0.0, 6.0), zlim=(0.0, 6.0))
flowrate = fk.TimeDependentFlowRate(times=[0.0, 0.0], flowrates=[0.0, 0.0])
fluvial_flowrates = [fk.FluvialFlowRate(ref=4,
time_dependent_flowrate=flowrate,
rectangular_tube=tube,
x_flux_direction=1.0,
y_flux_direction=0.0)]
Problem definition:¶
thacker = fk.Problem(simutime, triangular_mesh, layer, primitives,
tracers=tracers,
slides=slides,
numerical_parameters={'space_second_order':SECOND_ORDER},
torrential_outs=torrential_outs,
fluvial_flowrates=fluvial_flowrates)
===================================================================
| INITIALIZATION |
===================================================================
Problem size: Nodes=441, Layers=1, Triangles=800,
Iter = 0 ; Dt = 0.0000s ; Time = 0.00s ; ETA = 0.00s
Problem solving:¶
if not args.nographics:
# Plot triangular mesh.
TG = thacker.triangular_mesh.triangulation
x = np.asarray(TG.x)
y = np.asarray(TG.y)
trivtx = np.asarray(TG.trivtx)
# Plot y=5 line
y5 = np.argwhere(np.isclose(thacker.triangular_mesh.vertices.y, 5))
y5 = y5.flatten()
# Plot topography
plt.rcParams["figure.figsize"] = [10, 8]
topo = np.asarray(thacker.layer.topography)
plot_H = fk.TimeScheduler(np.linspace(0, 4.5, 10))
H_figure, H_axes = plt.subplots(5, 2, sharex='col', sharey='row')
H_axes = H_axes.T.flatten()
# Time steping:
while not thacker.simutime.is_finished:
thacker.forward()
if plot_H.now(thacker.simutime):
H = np.asarray(thacker.primitives.H)
ax = H_axes[plot_H.iteration]
ax.plot(x[y5], topo[y5] + H[y5])
ax.set_title("time = {:.1f}".format(thacker.simutime.time))
ax.plot(x[y5], topo[y5], color='grey')
ax.fill_between(x[y5], -0.5, topo[y5], facecolor='lightgrey', interpolate=True)
ax.set_xlim(-10, 10)
ax.set_ylim(-0.5, 10)
H_figure.subplots_adjust(hspace=0.5)
plt.show()
plt.rcParams["figure.figsize"] = [8, 6]
plt.figure()
plt.triplot(x, y, trivtx)
plt.plot(x[y5], y[y5], 'ro')
Iter = 11 ; Dt = 0.0086s ; Time = 0.10s ; ETA = 1.09s
Iter = 22 ; Dt = 0.0084s ; Time = 0.19s ; ETA = 1.09s
Iter = 34 ; Dt = 0.0082s ; Time = 0.29s ; ETA = 1.09s
Iter = 45 ; Dt = 0.0080s ; Time = 0.38s ; ETA = 1.10s
Iter = 57 ; Dt = 0.0078s ; Time = 0.47s ; ETA = 1.11s
Iter = 69 ; Dt = 0.0077s ; Time = 0.56s ; ETA = 1.15s
Iter = 82 ; Dt = 0.0076s ; Time = 0.66s ; ETA = 1.16s
Iter = 94 ; Dt = 0.0075s ; Time = 0.75s ; ETA = 1.12s
Iter = 107 ; Dt = 0.0074s ; Time = 0.85s ; ETA = 1.14s
Iter = 119 ; Dt = 0.0074s ; Time = 0.94s ; ETA = 1.14s
Iter = 132 ; Dt = 0.0074s ; Time = 1.04s ; ETA = 1.11s
Iter = 145 ; Dt = 0.0074s ; Time = 1.13s ; ETA = 1.10s
Iter = 157 ; Dt = 0.0074s ; Time = 1.22s ; ETA = 1.10s
Iter = 170 ; Dt = 0.0074s ; Time = 1.32s ; ETA = 1.04s
Iter = 183 ; Dt = 0.0075s ; Time = 1.41s ; ETA = 1.01s
Iter = 195 ; Dt = 0.0075s ; Time = 1.50s ; ETA = 0.95s
Iter = 208 ; Dt = 0.0076s ; Time = 1.60s ; ETA = 0.96s
Iter = 220 ; Dt = 0.0077s ; Time = 1.69s ; ETA = 0.91s
Iter = 232 ; Dt = 0.0076s ; Time = 1.79s ; ETA = 0.92s
Iter = 245 ; Dt = 0.0074s ; Time = 1.88s ; ETA = 0.87s
Iter = 257 ; Dt = 0.0072s ; Time = 1.97s ; ETA = 0.89s
Iter = 272 ; Dt = 0.0059s ; Time = 2.07s ; ETA = 1.03s
Iter = 288 ; Dt = 0.0059s ; Time = 2.16s ; ETA = 1.00s
Iter = 304 ; Dt = 0.0063s ; Time = 2.26s ; ETA = 0.91s
Iter = 317 ; Dt = 0.0084s ; Time = 2.35s ; ETA = 0.64s
Iter = 328 ; Dt = 0.0084s ; Time = 2.44s ; ETA = 0.60s
Iter = 340 ; Dt = 0.0082s ; Time = 2.54s ; ETA = 0.60s
Iter = 351 ; Dt = 0.0081s ; Time = 2.63s ; ETA = 0.58s
Iter = 364 ; Dt = 0.0070s ; Time = 2.73s ; ETA = 0.62s
Iter = 378 ; Dt = 0.0060s ; Time = 2.82s ; ETA = 0.68s
Iter = 396 ; Dt = 0.0047s ; Time = 2.91s ; ETA = 0.82s
Iter = 419 ; Dt = 0.0034s ; Time = 3.01s ; ETA = 3.09s
Iter = 454 ; Dt = 0.0021s ; Time = 3.10s ; ETA = 1.68s
Iter = 471 ; Dt = 0.0076s ; Time = 3.20s ; ETA = 0.43s
Iter = 483 ; Dt = 0.0075s ; Time = 3.29s ; ETA = 0.40s
Iter = 496 ; Dt = 0.0075s ; Time = 3.39s ; ETA = 0.37s
Iter = 508 ; Dt = 0.0075s ; Time = 3.48s ; ETA = 0.35s
Iter = 521 ; Dt = 0.0075s ; Time = 3.57s ; ETA = 0.32s
Iter = 533 ; Dt = 0.0076s ; Time = 3.66s ; ETA = 0.29s
Iter = 546 ; Dt = 0.0076s ; Time = 3.76s ; ETA = 0.26s
Iter = 558 ; Dt = 0.0077s ; Time = 3.85s ; ETA = 0.23s
Iter = 570 ; Dt = 0.0078s ; Time = 3.95s ; ETA = 0.20s
Iter = 582 ; Dt = 0.0079s ; Time = 4.04s ; ETA = 0.17s
Iter = 594 ; Dt = 0.0080s ; Time = 4.14s ; ETA = 0.14s
Iter = 605 ; Dt = 0.0082s ; Time = 4.23s ; ETA = 0.11s
Iter = 617 ; Dt = 0.0082s ; Time = 4.33s ; ETA = 0.08s
Iter = 629 ; Dt = 0.0066s ; Time = 4.41s ; ETA = 0.07s
Iter = 644 ; Dt = 0.0065s ; Time = 4.51s ; ETA = 0.09s
Iter = 658 ; Dt = 0.0041s ; Time = 4.60s ; ETA = 0.00s
Total running time of the script: ( 0 minutes 2.681 seconds)