Note
Click here to download the full example code
River¶
In this example a flow in a bended river is considered. A simple configuration is chosen with water viscosity, Navier friction and flat topography.
This example purpose is to show the basic layout of frashkiss3d scripts.
import os, sys
import argparse
import matplotlib.pyplot as plt
from matplotlib.pyplot import cm
import matplotlib.tri as mtri
import numpy as np
import freshkiss3d as fk
import freshkiss3d.extra.plots as fk_plt
os.system('rm -r outputs')
parser = argparse.ArgumentParser()
parser.add_argument('--nographics', action='store_true')
args = parser.parse_args()
#sphinx_gallery_thumbnail_number = 2
Parameters:¶
SECOND_ORDER=False
NUM_PARAMS={'ipres':True,
'space_second_order':SECOND_ORDER,
'flux_type':1}
PHY_PARAMS={'friction_law':'Navier',
'friction_coeff':0.01,
'horizontal_viscosity':0.0013,
'vertical_viscosity':0.0013}
Time loop:¶
simutime = fk.SimuTime(final_time=500., time_iteration_max=100000,
second_order=SECOND_ORDER)
Mesh:¶
Mesh is imported from mesh file located in /inputs directory.
dir_path = os.path.abspath(os.path.dirname(sys.argv[0]))
triangular_mesh = fk.TriangularMesh.from_msh_file(dir_path + '/inputs/river.msh')
if not args.nographics:
fk_plt.plot_mesh(triangular_mesh, plot_labels=False)
Layers:¶
In this examlpe topography is flat all along the river.
NL = 8
layer = fk.Layer(NL, triangular_mesh, topography=0.)
Primitives:¶
primitives = fk.Primitive(triangular_mesh, layer, height=10.0)
Boundary conditions:¶
FluvialFlowRate is set for inlet boundary condition and FluvialHeight for outlet.
On side walls a Slide boundary condition is used with wall_friction_coeff.
flowrate = fk.TimeDependentFlowRate(times=[0.0, 1000.0, 10000],
flowrates=[0.0, 10.0, 10.0])
fluvial_flowrates = [fk.FluvialFlowRate(ref=2,
time_dependent_flowrate=flowrate,
x_flux_direction=1.0,
y_flux_direction=0.0)]
fluvial_heights = [fk.FluvialHeight(ref=3, height=10.)]
slides = [fk.Slide(ref=1, horizontal_viscosity=PHY_PARAMS['horizontal_viscosity'],
wall_friction_coeff=PHY_PARAMS['friction_coeff'])]
Writter:¶
vtk_writer = fk.VTKWriter(triangular_mesh, scheduler=fk.schedules(count=10), scale_h=10.)
Problem definition:¶
Problem is called with basic parameters (simutime, mesh, layer, primitives and
boundary conditions). See API for other useful fk.Problem parameters.
river = fk.Problem(simutime, triangular_mesh, layer, primitives,
slides=slides,
fluvial_flowrates=fluvial_flowrates,
fluvial_heights=fluvial_heights,
numerical_parameters=NUM_PARAMS,
physical_parameters=PHY_PARAMS,
vtk_writer=vtk_writer)
===================================================================
| INITIALIZATION |
===================================================================
Problem size: Nodes=655, Layers=8, Triangles=1154,
Iter = 0 ; Dt = 0.0000s ; Time = 0.00s ; ETA = 0.00s
Problem solving:¶
river.solve()
===================================================================
| TIME LOOP |
===================================================================
Iter = 36 ; Dt = 0.2908s ; Time = 10.47s ; ETA = 7.88s
Iter = 71 ; Dt = 0.2908s ; Time = 20.64s ; ETA = 7.47s
Iter = 106 ; Dt = 0.2908s ; Time = 30.82s ; ETA = 7.63s
Iter = 141 ; Dt = 0.2908s ; Time = 41.00s ; ETA = 7.45s
Iter = 176 ; Dt = 0.2908s ; Time = 51.17s ; ETA = 7.10s
Iter = 211 ; Dt = 0.2908s ; Time = 61.35s ; ETA = 7.05s
Iter = 246 ; Dt = 0.2908s ; Time = 71.53s ; ETA = 6.95s
Iter = 281 ; Dt = 0.2908s ; Time = 81.70s ; ETA = 6.78s
Iter = 316 ; Dt = 0.2908s ; Time = 91.88s ; ETA = 6.41s
Iter = 351 ; Dt = 0.2908s ; Time = 102.06s ; ETA = 6.26s
Iter = 387 ; Dt = 0.2908s ; Time = 112.52s ; ETA = 6.36s
Iter = 422 ; Dt = 0.2908s ; Time = 122.70s ; ETA = 6.07s
Iter = 457 ; Dt = 0.2908s ; Time = 132.88s ; ETA = 5.90s
Iter = 492 ; Dt = 0.2908s ; Time = 143.05s ; ETA = 5.71s
Iter = 527 ; Dt = 0.2908s ; Time = 153.23s ; ETA = 5.44s
Iter = 562 ; Dt = 0.2908s ; Time = 163.40s ; ETA = 5.32s
Iter = 597 ; Dt = 0.2908s ; Time = 173.58s ; ETA = 5.13s
Iter = 632 ; Dt = 0.2908s ; Time = 183.76s ; ETA = 4.85s
Iter = 667 ; Dt = 0.2908s ; Time = 193.93s ; ETA = 4.63s
Iter = 702 ; Dt = 0.2908s ; Time = 204.11s ; ETA = 4.43s
Iter = 737 ; Dt = 0.2908s ; Time = 214.29s ; ETA = 4.32s
Iter = 773 ; Dt = 0.2908s ; Time = 224.75s ; ETA = 4.32s
Iter = 808 ; Dt = 0.2908s ; Time = 234.93s ; ETA = 4.13s
Iter = 843 ; Dt = 0.2908s ; Time = 245.11s ; ETA = 3.91s
Iter = 878 ; Dt = 0.2908s ; Time = 255.28s ; ETA = 3.83s
Iter = 913 ; Dt = 0.2908s ; Time = 265.46s ; ETA = 3.55s
Iter = 948 ; Dt = 0.2908s ; Time = 275.64s ; ETA = 3.50s
Iter = 983 ; Dt = 0.2908s ; Time = 285.81s ; ETA = 3.31s
Iter = 1018 ; Dt = 0.2908s ; Time = 295.99s ; ETA = 3.10s
Iter = 1053 ; Dt = 0.2908s ; Time = 306.17s ; ETA = 3.00s
Iter = 1088 ; Dt = 0.2908s ; Time = 316.34s ; ETA = 2.79s
Iter = 1124 ; Dt = 0.2908s ; Time = 326.81s ; ETA = 2.62s
Iter = 1159 ; Dt = 0.2908s ; Time = 336.99s ; ETA = 2.48s
Iter = 1194 ; Dt = 0.2908s ; Time = 347.16s ; ETA = 2.36s
Iter = 1229 ; Dt = 0.2908s ; Time = 357.34s ; ETA = 2.18s
Iter = 1264 ; Dt = 0.2908s ; Time = 367.51s ; ETA = 2.01s
Iter = 1299 ; Dt = 0.2908s ; Time = 377.69s ; ETA = 1.87s
Iter = 1334 ; Dt = 0.2908s ; Time = 387.87s ; ETA = 1.71s
Iter = 1369 ; Dt = 0.2908s ; Time = 398.04s ; ETA = 1.54s
Iter = 1404 ; Dt = 0.2908s ; Time = 408.22s ; ETA = 1.41s
Iter = 1439 ; Dt = 0.2908s ; Time = 418.40s ; ETA = 1.26s
Iter = 1474 ; Dt = 0.2908s ; Time = 428.57s ; ETA = 1.11s
Iter = 1510 ; Dt = 0.2908s ; Time = 439.04s ; ETA = 0.95s
Iter = 1545 ; Dt = 0.2908s ; Time = 449.22s ; ETA = 0.78s
Iter = 1580 ; Dt = 0.2908s ; Time = 459.39s ; ETA = 0.63s
Iter = 1615 ; Dt = 0.2908s ; Time = 469.57s ; ETA = 0.47s
Iter = 1650 ; Dt = 0.2908s ; Time = 479.75s ; ETA = 0.32s
Iter = 1685 ; Dt = 0.2908s ; Time = 489.92s ; ETA = 0.15s
Iter = 1720 ; Dt = 0.1919s ; Time = 500.00s ; ETA = 0.00s
===================================================================
| END |
===================================================================
Problem.solve() completed in 9.125890731811523s (wall time)
Basic plots:¶
if not args.nographics:
x = np.asarray(triangular_mesh.triangulation.x)
y = np.asarray(triangular_mesh.triangulation.y)
trivtx = np.asarray(triangular_mesh.triangulation.trivtx)
triang = mtri.Triangulation(x, y, trivtx)
LAYER_IDX = int(NL/2)
plt.rcParams["figure.figsize"] = [16, 6]
fig = plt.figure()
# Subplot 1: x velocity
ax1 = fig.add_subplot(121)
ax1.triplot(x, y, trivtx, color='k', lw=0.5)
im1 = ax1.tricontourf(x, y, trivtx, river.primitives.U[:, LAYER_IDX], 30, cmap=plt.cm.jet)
fig.colorbar(im1, ax=ax1)
ax1.set_title("Velocity x component on layer {}".format(LAYER_IDX))
# Subplot 2: y velocity
ax2 = fig.add_subplot(122)
ax2.triplot(x, y, trivtx, color='k', lw=0.5)
im2 = ax2.tricontourf(x, y, trivtx, river.primitives.V[:, LAYER_IDX], 30, cmap=plt.cm.jet)
fig.colorbar(im2, ax=ax2)
ax2.set_title("Velocity y component on layer {}".format(LAYER_IDX))
plt.show()
Streamlines plot:¶
if not args.nographics:
plt.rcParams["figure.figsize"] = [8, 6]
fig = plt.figure()
ax1 = fig.add_subplot(111)
xi, yi = np.meshgrid(np.linspace(0, 2000, 500),
np.linspace(0, 2000, 500))
U_interp = mtri.LinearTriInterpolator(triang, river.primitives.U[:, LAYER_IDX])
V_interp = mtri.LinearTriInterpolator(triang, river.primitives.V[:, LAYER_IDX])
U_i = U_interp(xi, yi)
V_i = V_interp(xi, yi)
speed = np.sqrt(U_i**2 + V_i**2)
ax1.set_xlim(0, 2000)
ax1.set_ylim(0, 2000)
cmap = cm.get_cmap(name='jet', lut=None)
strm = ax1.streamplot(xi, yi, U_i, V_i, density=(2, 2), color=speed, linewidth=2*speed/speed.max(), cmap=cmap)
fig.colorbar(strm.lines)
ax1.set_title("Velocity & Streamlines on layer {}".format(LAYER_IDX))
plt.show()
Total running time of the script: ( 0 minutes 11.286 seconds)