Running the models#
We understand plume evolution and melt generation as it ascends using the accompanying parameter file, plume_with_melt.prm:
In the cells below, it is assumed that the ASPECT directory is installed system-wide or that you are using the HubZero ASPECT installation. If not, modify the ASPECT executable to the location where it is installed.
ASPECT can be configured in either “release” mode or “debug” mode. Release mode is more optimized than debug mode, and will usually run much faster but lacks debugging symbols and more descriptive error messages. If you are making changes to ASPECT or altering parameter files, it might be a good idea to try running in debug mode in order to resolve any problems. The code below is calling the release version using
aspect-release. If you want to try running debug mode, callaspect-debuginstead.
! aspect-release ./prm_files/plume_with_melt.prm
import os
os.environ['VTK_DEFAULT_OPENGL_WINDOW'] = 'vtkEGLRenderWindow' # Suppresses a warning with Pyvista
# Load the relevant libraries
import pyvista as pv
import glob
from IPython.display import HTML
pv.global_theme.allow_empty_mesh = True
# pv.start_xvfb()
def plot_deformation(output_dir, field, gif_file, cmap='RdBu', ncontours=5):
'''
Process a series of VTU files in the given output directory, returning the 'field'
data for plotting.
Inputs:
-----------
output_dir : str
Path to the folder containing solution VTU files (expected in solution/solution-*.vtu)
field : str
Name of the field we want to look at in the solution
label : str
Label for the field being plotted
'''
solution_file_names = sorted(glob.glob(output_dir + 'solution/solution-*.pvtu'))
mesh = pv.read(solution_file_names[0])
plotter = pv.Plotter()
# Create the first frame
plotter.add_mesh(mesh, scalars=field, show_scalar_bar=False)
plotter.add_axes()
plotter.hide_axes()
# We know that the model is XY plane so the following
# lines just sets up the camera
plotter.view_xy()
plotter.camera.roll = 0
plotter.camera.parallel_projection = True
sargs_1 = dict(height=0.5, vertical=True, position_x=0.05, position_y=0.05)
sargs_2 = dict(height=0.5, vertical=True, position_x=0.85, position_y=0.05)
plotter.open_gif(gif_file, fps=4)
# Update function for each frame
def update(frame):
frame_int = int(frame)
new_mesh = pv.read(solution_file_names[frame_int])
# Generate contour lines
contours = new_mesh.contour(isosurfaces=ncontours, scalars='porosity')
plotter.add_mesh(new_mesh, clim=[0.0, 1.0], scalars=field, cmap=cmap, \
scalar_bar_args=sargs_1)
if contours.n_points > 0:
plotter.add_mesh(contours, clim=[0.0, 0.5], line_width=2, cmap='viridis_r', scalar_bar_args=sargs_2)
plotter.add_text(f"Time step: {frame}", name="tstep")
plotter.add_text(field, position='upper_right', font_size=14)
n_frames = len(solution_file_names)
for i in range(n_frames):
update(i)
plotter.write_frame()
plotter.close()
plot_deformation('plume-with-melt/', 'peridotite', 'images/peridotite.gif')
Now, you would have the gifs saved in the current directory this cell visualizes the model outputs side by side for comparison
Exercises#
How does the plume evolve if you increase/decrease the shear viscosity? How is the melt generated affected?
The melt volume generated in the above example is ~ 20 \(km^3\), much smaller than what is expected for continental flood basalts. What parameters could you change to increase the melt volume? Try changing them and see how the plume evolves.
Appendix#
The example file provided takes about 20 minutes on 6 cores. When the model is complete, you should expect the following output on your terminal: