Skip to content
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
55 changes: 13 additions & 42 deletions simulator/flow_rock.py
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@
from pylops.utils.wavelets import ricker
from pylops.signalprocessing import Convolve1D
import sys
from PyGRDECL.GRDECL_Parser import GRDECL_Parser # https://github.com/BinWang0213/PyGRDECL/tree/master
#from PyGRDECL.GRDECL_Parser import GRDECL_Parser # https://github.com/BinWang0213/PyGRDECL/tree/master
from scipy.interpolate import interp1d
from scipy.interpolate import griddata
from pipt.misc_tools.analysis_tools import store_ensemble_sim_information
Expand Down Expand Up @@ -759,14 +759,13 @@ def run_fwd_sim(self, state, member_i, del_folder=True):

# The inherited simulator also has a run_fwd_sim. Call this.
self.ensemble_member = member_i
<<<<<<< Updated upstream
return super().run_fwd_sim(state, member_i, del_folder=del_folder)

=======
#return super().run_fwd_sim(state, member_i, del_folder=del_folder)


self.pred_data = super().run_fwd_sim(state, member_i, del_folder=del_folder)
return self.pred_data
>>>>>>> Stashed changes

def call_sim(self, folder=None, wait_for_proc=False, run_reservoir_model=None, save_folder=None):
# replace the sim2seis part (which is unusable) by avo based on Pylops

Expand Down Expand Up @@ -909,10 +908,7 @@ def get_avo_result(self, folder, save_folder):
def calc_velocities(self, folder, save_folder, grid, v, f_dim):
# The properties in pem are only given in the active cells
# indices of active cells:
<<<<<<< Updated upstream

=======
>>>>>>> Stashed changes
true_indices = np.where(grid['ACTNUM'])

# Alt 2
Expand Down Expand Up @@ -959,7 +955,7 @@ def calc_velocities(self, folder, save_folder, grid, v, f_dim):
save_dic = {'vp': vp, 'vs': vs, 'rho': rho}#, 'bulkmod': self.bulkmod, 'shearmod': self.shearmod,
#'Pov': self.poverburden, 'P': self.pressure, 'Peff': self.peff, 'por': porosity} # for debugging
else:
save_dic = {'vp': vp, 'vs': vs, 'rho': rho, 'por': porosity, 'sgas': sgas, 'Pd': pdyn}
save_dic = {'vp': vp, 'vs': vs, 'rho': rho}#, 'por': porosity, 'sgas': sgas, 'Pd': pdyn}

if save_folder is not None:
file_name = save_folder + os.sep + f"vp_vs_rho_vint{v}.npz" if save_folder[-1] != os.sep \
Expand Down Expand Up @@ -1033,12 +1029,12 @@ def _get_props(self, kw_file):
with np.load(file) as f:
exec(f'self.{kw} = f[ "{kw}" ]')
self.NX, self.NY, self.NZ = f['NX'], f['NY'], f['NZ']
else:
reader = GRDECL_Parser(filename=file)
reader.read_GRDECL()
exec(f"self.{kw} = reader.{kw}.reshape((reader.NX, reader.NY, reader.NZ), order='F')")
self.NX, self.NY, self.NZ = reader.NX, reader.NY, reader.NZ
eval(f'np.savez("./{kw}.npz", {kw}=self.{kw}, NX=self.NX, NY=self.NY, NZ=self.NZ)')
#else:
# reader = GRDECL_Parser(filename=file)
# reader.read_GRDECL()
# exec(f"self.{kw} = reader.{kw}.reshape((reader.NX, reader.NY, reader.NZ), order='F')")
# self.NX, self.NY, self.NZ = reader.NX, reader.NY, reader.NZ
# eval(f'np.savez("./{kw}.npz", {kw}=self.{kw}, NX=self.NX, NY=self.NY, NZ=self.NZ)')

def _calc_avo_props(self, dt=0.0005):
# dt is the fine resolution sampling rate
Expand Down Expand Up @@ -1249,6 +1245,7 @@ def _calc_avo_props_active_cells(self, grid, vp, vs, rho, dt=0.0005):
kind='nearest', fill_value='extrapolate')
trace_interp[ind, :] = f(t_interp)


if i == 0:
avo_data = trace_interp # 3D
elif i == 1:
Expand Down Expand Up @@ -2260,34 +2257,8 @@ def van_opstal(self, lambda_vals, z_res, z_base, poisson):

return value

def van_opstal_org(self, lambda_vals, z_res, z_base, poisson):
"""
Compute the Van Opstal transfer function.

Args:
lambda_vals -- Numpy array of lambda values.
z_res -- Depth to reservoir [m].
z_base -- Distance to the basement [m].
poisson -- Poisson's ratio.

Returns:
value -- Numpy array of computed values.
"""

term1 = np.exp(lambda_vals * z_res) * (2 * lambda_vals * z_base + 1)
term2 = np.exp(-lambda_vals * z_res) * (
4 * lambda_vals ** 2 * z_base ** 2 + 2 * lambda_vals * z_base + (3 - 4 * poisson) ** 2)

term3_numer = (3 - 4 * poisson) * (
np.exp(-lambda_vals * (2 * z_base + z_res)) - np.exp(-lambda_vals * (2 * z_base - z_res)))
term3_denom = 2 * ((1 - 2 * poisson) ** 2 + lambda_vals ** 2 * z_base ** 2 + (3 - 4 * poisson) * np.cosh(
lambda_vals * z_base) ** 2)

value = term1 - term2 - (term3_numer / term3_denom)

return value
def hankel_transform_order_0(self, f, r_max, num_points=1000):

def hankel_transform_order_0(f, r_max, num_points=1000):
"""
Computes the Hankel transform of order 0 of a function f(r).

Expand Down