Post-processing with Matplotlib

In this example we show how to use matplotlib as a pre and post-treatment tool to easily plot meshes and solutions. River example is used as a basis for this tutorial. For advanced visualizations it is strongly advised to use a post-treatment software like paraview or mayavi.

import matplotlib.pyplot as plt
import matplotlib.tri as mtri
import numpy as np
import freshkiss3d as fk

Quick case set up:

simutime = fk.SimuTime(final_time=1500., time_iteration_max=100000)

triangular_mesh = fk.TriangularMesh.from_msh_file('../simulations/inputs/river.msh')

NL = 8
layer = fk.Layer(NL, triangular_mesh, topography=0.)

primitives = fk.Primitive(triangular_mesh, layer, height=1.0)

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)]

slides = [fk.Slide(ref=1)]
fluvial_heights = [fk.FluvialHeight(ref=3, height=1.)]

Plot boundary conditions on mesh:

To check boundary labels:

colors = {0:'blue', 1:'red', 2:'green', 3:'purple'}

fig = plt.figure()
ax = fig.add_subplot(111)
ax.triplot(triangular_mesh.triangulation.x,
           triangular_mesh.triangulation.y,
           triangular_mesh.triangulation.trivtx)

# Plot edges.
for B in range(triangular_mesh.boundary_edges.size):
    i0 = triangular_mesh.boundary_edges.vertices[B, 0]
    i1 = triangular_mesh.boundary_edges.vertices[B, 1]
    x0,y0 = triangular_mesh.triangulation.x[i0], triangular_mesh.triangulation.y[i0]
    x1,y1 = triangular_mesh.triangulation.x[i1], triangular_mesh.triangulation.y[i1]
    ax.plot([x0, x1], [y0, y1], color=colors[triangular_mesh.boundary_edges.label[B]], linewidth=2.)
    offset = 0.04
    x_edge = 0.5*( x0 + x1) + offset
    y_edge = 0.5*( y0 + y1) + offset
    ax.text(x_edge, y_edge, triangular_mesh.boundary_edges.label[B],
            color= colors[triangular_mesh.boundary_edges.label[B]],fontsize=10)

plt.grid()
plt.axis('scaled')
plt.show()

Problem solving:

river = fk.Problem(simutime, triangular_mesh, layer, primitives,
                   slides=slides,
                   fluvial_flowrates=fluvial_flowrates,
                   fluvial_heights=fluvial_heights)

river.solve()

===================================================================
|                          INITIALIZATION                         |
===================================================================
Problem size: Nodes=655, Layers=8, Triangles=1154,
Iter =        0 ; Dt = 0.0000s ; Time =     0.00s ; ETA =     0.00s
===================================================================
|                            TIME LOOP                            |
===================================================================
Iter =       34 ; Dt = 0.9194s ; Time =    31.26s ; ETA =     4.38s
Iter =       67 ; Dt = 0.9194s ; Time =    61.60s ; ETA =     4.30s
Iter =      100 ; Dt = 0.9193s ; Time =    91.94s ; ETA =     4.30s
Iter =      134 ; Dt = 0.9193s ; Time =   123.19s ; ETA =     4.18s
Iter =      167 ; Dt = 0.9192s ; Time =   153.53s ; ETA =     4.02s
Iter =      200 ; Dt = 0.9192s ; Time =   183.86s ; ETA =     3.96s
Iter =      234 ; Dt = 0.9191s ; Time =   215.11s ; ETA =     3.82s
Iter =      267 ; Dt = 0.9191s ; Time =   245.45s ; ETA =     3.81s
Iter =      300 ; Dt = 0.9191s ; Time =   275.78s ; ETA =     3.74s
Iter =      334 ; Dt = 0.9190s ; Time =   307.02s ; ETA =     3.66s
Iter =      367 ; Dt = 0.9190s ; Time =   337.35s ; ETA =     3.49s
Iter =      400 ; Dt = 0.9189s ; Time =   367.67s ; ETA =     3.48s
Iter =      433 ; Dt = 0.9189s ; Time =   398.00s ; ETA =     3.35s
Iter =      467 ; Dt = 0.9188s ; Time =   429.24s ; ETA =     3.24s
Iter =      500 ; Dt = 0.9188s ; Time =   459.56s ; ETA =     3.11s
Iter =      533 ; Dt = 0.9188s ; Time =   489.88s ; ETA =     3.10s
Iter =      567 ; Dt = 0.9187s ; Time =   521.12s ; ETA =     2.97s
Iter =      600 ; Dt = 0.9187s ; Time =   551.43s ; ETA =     2.88s
Iter =      633 ; Dt = 0.9186s ; Time =   581.75s ; ETA =     2.81s
Iter =      667 ; Dt = 0.9186s ; Time =   612.98s ; ETA =     2.70s
Iter =      700 ; Dt = 0.9185s ; Time =   643.29s ; ETA =     2.67s
Iter =      733 ; Dt = 0.9185s ; Time =   673.60s ; ETA =     2.58s
Iter =      767 ; Dt = 0.9184s ; Time =   704.83s ; ETA =     2.43s
Iter =      800 ; Dt = 0.9184s ; Time =   735.14s ; ETA =     2.39s
Iter =      833 ; Dt = 0.9184s ; Time =   765.45s ; ETA =     2.30s
Iter =      867 ; Dt = 0.9183s ; Time =   796.67s ; ETA =     2.11s
Iter =      900 ; Dt = 0.9183s ; Time =   826.97s ; ETA =     2.09s
Iter =      933 ; Dt = 0.9182s ; Time =   857.27s ; ETA =     1.95s
Iter =      967 ; Dt = 0.9182s ; Time =   888.49s ; ETA =     1.83s
Iter =     1000 ; Dt = 0.9181s ; Time =   918.79s ; ETA =     1.78s
Iter =     1033 ; Dt = 0.9181s ; Time =   949.09s ; ETA =     1.67s
Iter =     1067 ; Dt = 0.9180s ; Time =   980.30s ; ETA =     1.58s
Iter =     1100 ; Dt = 0.9180s ; Time =  1010.60s ; ETA =     1.49s
Iter =     1133 ; Dt = 0.9180s ; Time =  1040.89s ; ETA =     1.41s
Iter =     1167 ; Dt = 0.9180s ; Time =  1072.11s ; ETA =     1.33s
Iter =     1200 ; Dt = 0.9180s ; Time =  1102.40s ; ETA =     1.21s
Iter =     1233 ; Dt = 0.9180s ; Time =  1132.69s ; ETA =     1.12s
Iter =     1267 ; Dt = 0.9180s ; Time =  1163.91s ; ETA =     1.04s
Iter =     1300 ; Dt = 0.9180s ; Time =  1194.20s ; ETA =     0.97s
Iter =     1333 ; Dt = 0.9180s ; Time =  1224.50s ; ETA =     0.84s
Iter =     1367 ; Dt = 0.9180s ; Time =  1255.71s ; ETA =     0.76s
Iter =     1400 ; Dt = 0.9180s ; Time =  1286.01s ; ETA =     0.66s
Iter =     1434 ; Dt = 0.9180s ; Time =  1317.22s ; ETA =     0.57s
Iter =     1467 ; Dt = 0.9181s ; Time =  1347.51s ; ETA =     0.47s
Iter =     1500 ; Dt = 0.9181s ; Time =  1377.81s ; ETA =     0.38s
Iter =     1534 ; Dt = 0.9181s ; Time =  1409.02s ; ETA =     0.28s
Iter =     1567 ; Dt = 0.9181s ; Time =  1439.32s ; ETA =     0.19s
Iter =     1600 ; Dt = 0.9177s ; Time =  1469.61s ; ETA =     0.09s
Iter =     1634 ; Dt = 0.1105s ; Time =  1500.00s ; ETA =     0.00s
===================================================================
|                               END                               |
===================================================================
Problem.solve() completed in 4.604763746261597s (wall time)

Basic plots of velocity components:

First velocity components are plotted on the middle layer (see LAYER_IDX).

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)
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)
fig.colorbar(im2, ax=ax2)
ax2.set_title("Velocity y component on layer {}".format(LAYER_IDX))

plt.show()

Quivers plot:

plt.rcParams["figure.figsize"] = [6, 6]
fig = plt.figure()
ax1 = fig.add_subplot(111)
ax1.quiver(triang.x, triang.y, river.primitives.U[:, LAYER_IDX],
                               river.primitives.V[:, LAYER_IDX])
ax1.set_title("Quivers on layer {}".format(LAYER_IDX))

plt.show()

Streamlines plot:

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)
strm = ax1.streamplot(xi, yi, U_i, V_i, density=(2, 2), color=speed, linewidth=2*speed/speed.max())
fig.colorbar(strm.lines)
ax1.set_title("Velocity & Streamlines on layer {}".format(LAYER_IDX))

plt.show()