Skip to content

Instantly share code, notes, and snippets.

@jkfurtney
Created May 3, 2022 22:49
Show Gist options
  • Save jkfurtney/c2bd32a974ffc90e9abf8437688b72ac to your computer and use it in GitHub Desktop.
Save jkfurtney/c2bd32a974ffc90e9abf8437688b72ac to your computer and use it in GitHub Desktop.
Map porepressure from a csv file with only positive values to FLAC3D gridpoints.
import itasca as it
from itasca import gpa
import numpy as np
from scipy.constants import foot
from collections import defaultdict
import os
from scipy.spatial import cKDTree
import scipy
if os.path.exists("data.npy"):
data = np.load("data.npy")
else:
data = np.loadtxt("My_pore_pressure_data.csv",
delimiter=",", usecols=(1,2,3,4,5), skiprows=1)
np.save("data.npy", data)
# write for Paraview so we can visualize
data *= [foot,foot,foot,1,1000*9.81*foot] # convert to m,m,m,s,Pa
X,Y,Z,time,head = data.T
f3_gridpoints = gpa.pos()
def draw_fragment(gname, points, faces):
f = open(gname, "w")
print("ITASCA GEOMETRY3D", file=f)
print("NODES", file=f)
for j, p in enumerate(points):
x,y,z = p
print(j+1, x, y, z, "EXTRA 1 0", file=f)
print("POLYS", file=f)
for j, t in enumerate(faces):
n0, n1, n2 = t
print(j+1, "NODE", n0+1, n1+1, n2+1, "EXTRA 1 0", file=f)
f.close()
# this type of data has repeating x,y coordinates
# we bin them first
times = list(set(time))
for t in times:
mask = time==t
print("running t=", t)
gridded = defaultdict(list)
for row in data[mask]:
x,y,z,_,head = row
gridded[(x,y)].append((z,head))
new_points, new_head = [], []
free_surface = []
for x_y, z_head in gridded.items():
z_head = list(set(z_head))
z_head.sort(key=lambda _:_[0])
z0, h0 = z_head[-2]
z1, h1 = z_head[-1]
assert not z1==z0
assert not h1==h0
# look at the vertical pore pressure gradient and project upwards to find the pressure = zero height
slope = (h1-h0)/(z1-z0)
# we can calculate the slope or use the (known) hydrostatic gradient
slope = -1_000*9.81
assert slope < 0
delta_h = -h1/slope
new_points.append((*x_y, z1 + delta_h*2))
# because we are using nearest neighbor, we need to double delta h so we define the zero surface correctly.
free_surface.append((*x_y, z1 + delta_h))
new_head.append(0)
all_points = np.vstack((data[mask][:,:3], np.array(new_points)))
all_head = np.hstack((data[mask][:,-1], np.array(new_head)))
#all_head
distances, indices = cKDTree(all_points).query(f3_gridpoints)
print(np.histogram(distances))
gp_heads = all_head[indices]
print(f"saving pp_t_{int(t)}.npy")
np.save(f"gp_pp_t_{int(t)}.npy", gp_heads)
new_points = np.array(new_points)
mesh = scipy.spatial.Delaunay(np.array(free_surface)[:,:2])
draw_fragment(f"free_surface_t_{int(t)}.geom", np.array(free_surface), mesh.simplices)
@jkfurtney
Copy link
Author

This is for a special case of porepressure data that is kind of common. All the points are aligned on a repeating x,y grid but they can be at any height. There are no zero values so the nearest neighbor interpolation does not directly work. The procedure is to identify all the columns of points with the same x and y coordinates, sort them by height, look at the highest one, extrapolate upwards to find the phreatic surface, and add those as zero pressure points to the data. Then we do the nearest neighbor interpolation.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment