import itasca as it
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")