Process Complex Rift Inversion Model

This notebook contains example code for processing a production-level ASPECT model of rift inversion using GDTchron. Because the source files are ~100 Gb, this code is not fully reproducible without a local copy of these files. Nevertheless, it demonstrates the workflow.

The ASPECT parameter files needed to reproduce the model results locally are available in the ri_prms directory in the GDTchron GitHub repository.

# Imports
import os

import pyvista as pv
from joblib import Parallel, delayed
from tqdm import tqdm

import gdtchron as gdt

The cell below gathers the paths to the .pvtu files output by ASPECT from the local source for processing

# Get full output files from local source
part_path = (
    '/home/vaseylab/glaucophane/fastscape_testing/060524_rif/output_ri_rift/particles'
)
files = [
    os.path.join(part_path, file) for file in os.listdir(part_path)
    if file.endswith('.pvtu')
]
files.sort()

The cell below trims these large files by reducing their bounds to the upper and middle part of the model and by only including the ID and temperature fields needed for GDTchron. This makes subsequent processing of the files faster.

# Isolate temperature and particle id and trim files for size
clip_dir = 'ri_clipped'
bounds = [200e3, 800e3, 550e3, 620e3, 0, 0]  # m (X,Y,Z)
os.makedirs(clip_dir, exist_ok=True)


def trim_file(file, bounds, clip_dir):
    """Trims large VTK files to reduce size."""
    basename = os.path.basename(file)[:-5] + '.vtu'
    filename = os.path.join(clip_dir, basename)

    if os.path.exists(filename):
        return

    mesh = pv.read(file)
    temps = mesh['T']
    ids = mesh['id']
    mesh.point_data.clear()
    mesh['T'] = temps
    mesh['id'] = ids
    clipped = mesh.clip_box(bounds=bounds, invert=False)
    clipped.save(filename)
    return


output = Parallel(n_jobs=os.cpu_count() / 2)(
    delayed(trim_file)(file, bounds, clip_dir) for file in tqdm(files, position=0)
    )

The cell below gets the paths to the newly trimmed files for GDTchron processing

# Get path to new trimmed files
files = [
    os.path.join(clip_dir, file) for file in os.listdir(clip_dir)
    if file.endswith('.vtu')
]
files.sort()

The next few cells below run GDTchron to predict AHe ages and create a new set of VTK files with these ages, then add AFT and ZHe ages to that same set of VTK files.

# Add AHe ages and create new meshes.
gdt.run_vtk(files=files, system='AHe', time_interval=0.1, processes=os.cpu_count() / 2,
            batch_size=1000)

# Set path to reuse files from AHe
new_path = 'meshes_tchron'
files = [
    os.path.join(new_path, file) for file in os.listdir(new_path)
    if file.endswith('.vtu')
]
files.sort()

# Add AFT
gdt.run_vtk(files=files, system='AFT', time_interval=0.1, processes=os.cpu_count() / 2,
            batch_size=1000)

# Add ZHe
gdt.run_vtk(files=files, system='ZHe', time_interval=0.1, processes=os.cpu_count() / 2,
            batch_size=1000)

This final cell re-runs AHe and creates a new set of VTK files without interpolation of thermal histories for new particles, so particles that were added during the ASPECT run will have no AHe ages. This is to demonstrate the effects of the interpolation scheme on GDTchron results.

# Run AHe without interpolation in a new set of meshes
gdt.run_vtk(files=files, system='AHe', time_interval=0.1, processes=os.cpu_count() / 2,
            batch_size=1000, file_prefix='meshes_ahe_nointerp', interpolate_vals=False)