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
elevation
is not a georeferenced file which can be used to create a GRASS Location,EPSG
must 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)