Note
Click here to download the full example code
Ravaflow24 Mixture Model
This example shows how to simulate mass flow on a topography using
Ravaflow24Mixture.
First, import the class Ravaflow24Mixture by
from psimpy.simulator import Ravaflow24Mixture
To create an instance of this class, we must specify the parameter dir_sim.
It represents the directory in which output files generated by r.avaflow will be
saved.
import os
# Here we create a folder called `temp_Ravaflow24Mixture_example` to save output
# files generated by r.avaflow
cwd = os.getcwd()
if not os.path.exists('temp_Ravaflow24Mixture_example'):
os.mkdir('temp_Ravaflow24Mixture_example')
dir_sim = os.path.join(cwd, 'temp_Ravaflow24Mixture_example')
Now, we can create an instance of Ravaflow24Mixture by given
dir_sim and leaving other parameters to their default values
(Other parameters include time_step, time_end, curvature_control,
entrainment_control, stopping_control, etc.)
voellmy_model = Ravaflow24Mixture(dir_sim)
- To run a simulation using above voellmy_model, one needs to specify
elevation– Name of elevation raster file (including its path).hrelease– Name of release height raster file (including its path).prefix– Prefix required by r.avaflow to name output files.If
elevationis not a georeferenced file which can be used to create a GRASS Location,EPSGmust be provided.
Other optional parameters include internal_friction, basal_friction,
turbulent_friction, entrainment_coef etc.
The synthetic topography synthetic_topo.tif and release mass synthetic_rel.tif
are used here for illustration. They are located at the /tests/data/ folder.
dir_data = os.path.abspath('../../../tests/data/')
elevation = os.path.join(dir_data, 'synthetic_topo.tif')
hrelease = os.path.join(dir_data, 'synthetic_rel.tif')
prefix = 'synthetic'
Note
You may need to modify dir_data according to where you save them
on your local machine.
Given above inputs, we can create a GRASS Location and a shell file by calling
Ravaflow24Mixture.preprocess() method, and run the simulation using
Ravaflow24Mixture.run() method
grass_location, sh_file = voellmy_model.preprocess(
prefix=prefix, elevation=elevation, hrelease=hrelease, EPSG='2326')
voellmy_model.run(grass_location, sh_file)
Once the simulation is finished, the results are located in the folder /dir_sim/your_prefix_results. In this case, they are located in /temp_Ravaflow24Mixture_example/synthetic_results. We can extract desired outputs using respective method and visualize them. For example, we can extract the overall impact area, maximum flow velocity, or maximum flow velocity at specific locations:
import numpy as np
# overall impact area
impact_area = voellmy_model.extract_impact_area(prefix)
print(f"Impact_area is {impact_area} m^2")
# overall maximum flow velocity
v_mmax = voellmy_model.extract_qoi_max(prefix, 'v', aggregate=True)
print(f"Maximum flow velocity is {v_mmax} m/s^2")
# maximum flow velocity at specific locations
loc = np.array([[500, 2000], [1500, 2000]])
v_max_loc = voellmy_model.extract_qoi_max_loc(prefix, loc, 'v')
print(f"Maximum flow velocity at location {loc[0]} is {v_max_loc[0]} m/s^2")
print(f"Maximum flow velocity at location {loc[1]} is {v_max_loc[1]} m/s^2")
Impact_area is 1394800.0 m^2
Maximum flow velocity is 41.39 m/s^2
Maximum flow velocity at location [ 500 2000] is 34.264 m/s^2
Maximum flow velocity at location [1500 2000] is 36.689 m/s^2
We can visulize point-wise maximum flow height and velocity using heatmaps:
import matplotlib.pyplot as plt
import linecache
# extract head information of result ascii files and compute x and y coordinates
hmax_asc = os.path.join(dir_sim, f'{prefix}_results', f'{prefix}_ascii', f'{prefix}_hflow_max.asc')
header = [linecache.getline(hmax_asc, i) for i in range(1,6)]
header_values = [float(h.split()[-1].strip()) for h in header]
ncols, nrows, xll, yll, cellsize = header_values
ncols = int(ncols)
nrows = int(nrows)
x = np.arange(xll, xll+(cellsize*ncols), cellsize)
y = np.arange(yll, yll+(cellsize*nrows), cellsize)
# point-wise maximum flow height
h_max = voellmy_model.extract_qoi_max(prefix, 'h', aggregate=False)
# point-wise maximum flow velocity
v_max = voellmy_model.extract_qoi_max(prefix, 'v', aggregate=False)
fig, ax = plt.subplots(2, 1, figsize=(10,6))
fig0 = ax[0].contourf(x, y, h_max, levels=20)
ax[0].set_xlabel('x')
ax[0].set_ylabel('y')
ax[0].set_title('maximum flow height' + r' $(m/s)$')
cbar0 = plt.colorbar(fig0, ax=ax[0], format='%.1f', orientation='vertical', fraction=0.1)
fig1 = ax[1].contourf(x, y, v_max, levels=20)
ax[1].set_xlabel('x')
ax[1].set_ylabel('y')
ax[1].set_title('maximum flow velocity' + r' ($m/s^2$)')
cbar1 = plt.colorbar(fig1, ax=ax[1], format='%.1f', orientation='vertical', fraction=0.1)
plt.tight_layout()

Here we delete the folder temp_Ravaflow24Mixture_example and all files therein.
import shutil
shutil.rmtree(dir_sim)
Total running time of the script: ( 5 minutes 56.157 seconds)