Deterministic simulation

Let's simulate a volcanic eruption at Etna using GFS forecasts. The following input files are required for a single FALL3D run:

  • etna.inp: FALL3D configuration file
  • etna.gfs.nc: GFS weather data in netCDF format
  • GFS.tbl: Dictionary file with GFS variable names (included in the FALL3D distribution under the folder Other/Meteo/Tables)
  • JOB.cmd: A job script with a series of directives to inform the batch system about the characteristics of the job
  • Fall3d.r8.x: FALL3D executable program

Copy the required files to your projects folder:

cd /gpfs/projects/nct01/$USER
cp -r /gpfs/projects/nct00/nct00014/FALL3D_material/hands-on-2 .

Enter to the new folder and create a symlink to the FALL3D executable file. In other words, if the FALL3D installation path is $FALL3D_PATH, run the commands:

cd hands-on-2
ln -s $FALL3D_PATH/bin/Fall3d.r8.x .

Now you can run FALL3D!

Submitting jobs

The job script JOB.cmd contains a series of directives to inform the batch system about the characteristics of the job:

#!/bin/bash
#SBATCH --job-name=FALL3D
#SBATCH --output=%x_%j.out
#SBATCH --error=%x_%j.err
#SBATCH --nodes=1
#SBATCH --ntasks=20
#SBATCH --time=00:10:00
#SBATCH --qos=training
#SBATCH --reservation=Computational24

module purge
module load intel/2017.4
module load impi/2017.4
module load netcdf

INPFILE="etna.inp"
FALLTASK="all"
NX=5
NY=2
NZ=2

srun ./Fall3d.r8.x ${FALLTASK} ${INPFILE} ${NX} ${NY} ${NZ}

You can submit the job with sbatch:

sbatch JOB.cmd

In this case, we are requesting 20 tasks for the job FALL3D.

Visualizing and analyzing model outputs

An output file etna.res.nc will be generated if the run has been completed successfully. The ncdump program generates an ASCII representation of a netCDF file and can be used to explore the content of the FALL3D output file. First, load the netcdf module in the cluster:

module load netcdf

and run the following command to show the list of variables and some metadata information:

ncdump -h etna.res.nc

Quick visualization using ncview

Next, we use the ncview tool to generate quick view of the model output. Load the required modules:

module load netcdf
module load udunits
module load gsl
module load nco
module load ncview

You can execute ncview from the cluster now:

ncview etna.res.nc

Notes:

  • To run your X apps remotely, log in to the remote server over SSH with the -X option, which will enable X forwarding on the client end.
ssh -X {username}@server

Visualization using Python

The model results can be plotted using the python package Cartopy. First, activate the following Anaconda environment:

module load anaconda
source activate volcanology

and run the python script under the folder POSTPROCESSING:

python plot_map.py
Content of file plot_map.py. Click to expand!
import numpy as np
import xarray as xr
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.colors import BoundaryNorm
import cartopy.crs as crs
import cartopy.feature as cfeature

###
### Parameters
###
minval          = 0.1
key             = "tephra_col_mass"
fname           = "../etna.res.nc"
levels          = np.arange(0.0,4,0.25)
vlon, vlat      = 15.0, 37.75
cmap            = plt.cm.RdYlBu_r

###
### Set mininmum level
###
if minval>0: levels[0] = minval

###
### Read file
###
ds = xr.open_dataset(fname)

###
### Generate map
###
proj = crs.PlateCarree()
fig, ax = plt.subplots( subplot_kw={'projection': proj} )

###
### Add map features
###
BORDERS = cfeature.NaturalEarthFeature(
        scale     = '10m',
        category  = 'cultural',
        name      = 'admin_0_countries',
        edgecolor = 'gray',
        facecolor = 'none'
        )
LAND = cfeature.NaturalEarthFeature(
        'physical', 'land', '10m',
        edgecolor = 'none',
        facecolor = 'lightgrey',
        alpha     = 0.8
        )

ax.add_feature(LAND,zorder=0)
ax.add_feature(BORDERS, linewidth=0.4)

###
### Add grid lines
###
gl = ax.gridlines(
	crs         = crs.PlateCarree(),
    draw_labels = True,
	linewidth   = 0.5,
	color       = 'gray',
	alpha       = 0.5,
	linestyle   = '--')
gl.top_labels    = False
gl.right_labels  = False
gl.ylabel_style  = {'rotation': 90}

###
### Add vent location
###
ax.plot(vlon,vlat,color='red',marker='^')

###
### Plot contours
###
cbar = None
for it in range(ds.time.size):
	time_fmt = ds.isel(time=it)['time'].dt.strftime("%d/%m/%Y %H:%M").item()
	ax.set_title(time_fmt, loc='right')
	fc = ax.contourf(
		ds.lon,ds.lat,ds.isel(time=it)[key],
		levels    = levels,
		norm      = BoundaryNorm(levels,cmap.N),
		cmap      = cmap,
		extend    = 'max',
		transform = crs.PlateCarree()
		)

	###
	### Generate colorbar
	###
	if not cbar:
		cbar=fig.colorbar(fc,
			orientation = 'horizontal',
			label       = r'Tephra column mass [$g~m^{-2}$]',
			)

	###
	### Output plot
	###
	fname_plt = f"map_{it:03d}.png"
	plt.savefig(fname_plt,dpi=300,bbox_inches='tight')

	###
	### Clear contours
	###
	for item in fc.collections: item.remove()

You can create an animated gif using the following command:

convert -delay 10 -loop 0 *.png animation.gif

Notes:

  • The first time you run the python script, the Cartopy package will try to download some data, including coastlines information. Since we have no internet connection in the cluster, you can copy this data tou your home folder:
cp -r /gpfs/projects/nct00/nct00014/share ~/.local