ppv.py

import itasca as it
it.command('python-reset-state false')
import numpy as np
np.set_printoptions(threshold=20)
from itasca import gridpointarray as gpa

ppv = None

def store_ppv(*args):
    global ppv
    if ppv is None:
        ppv = np.zeros(it.gridpoint.count())
    ppv = np.maximum(np.linalg.norm(gpa.vel(), axis=1), ppv)

def ppv_save(*args):
    global ppv
    print("saving PPV")
    if ppv.any():
        gpa.set_extra(1,ppv)

def ppv_restore(*args):
    global ppv
    print("restoring PPV")
    try:
        print(it.zone.count())
        ppv = gpa.extra(1)
        print("max ppv during restore", ppv.max())
    except: pass

def ppv_new(*args):
    global ppv
    print("resetting PPV")
    ppv = None

def transit_time():
    z = it.zone.find(1)
    K, G, rho = z.prop("bulk"), z.prop("shear"), z.density()
    vp = np.sqrt((K + 4.0/3.0*G)/rho)
    d = np.sqrt(3*100**2)
    return d/vp

it.set_callback("store_ppv", 0.1)
it.set_callback("ppv_save", "save")
it.set_callback("ppv_restore", "restore")
it.set_callback("ppv_new", "new")