Skip to content

Commit

Permalink
adding impedance option to simpeg_2d
Browse files Browse the repository at this point in the history
  • Loading branch information
@kujaku11
kujaku11 committed Oct 23, 2024
1 parent 62635aa commit 6040c3e
Show file tree
Hide file tree
Showing 3 changed files with 136 additions and 53 deletions.
179 changes: 130 additions & 49 deletions mtpy/modeling/simpeg/data_2d.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,10 +32,38 @@ def __init__(self, dataframe, **kwargs):
self.include_elevation = True
self.invert_te = True
self.invert_tm = True
self.invert_zxy = False
self.invert_zyx = False

self.invert_impedance = False

for key, value in kwargs.items():
setattr(self, key, value)

@property
def invert_impedance(self):
return self._invert_impedance

@invert_impedance.setter
def invert_impedance(self, value):
if not isinstance(value, bool):
raise TypeError(
f"invert_impedance must be a boolean, not type{type(value)}"
)

if value:
self.invert_zxy = True
self.invert_zyx = True
self.invert_te = False
self.invert_tm = False
self._invert_impedance = True
if not value:
self.invert_zxy = False
self.invert_zyx = False
self.invert_te = True
self.invert_tm = True
self._invert_impedance = False

@property
def station_locations(self):
"""
Expand Down Expand Up @@ -91,21 +119,41 @@ def _get_mode_sources(self, simpeg_mode):
"""
rx_locs = self.station_locations.copy()
rx_list = [
nsem.receivers.PointNaturalSource(
rx_locs,
orientation=simpeg_mode,
component="apparent_resistivity",
),
nsem.receivers.PointNaturalSource(
rx_locs, orientation=simpeg_mode, component="phase"
),
]

src_list = [
nsem.sources.Planewave(rx_list, frequency=f)
for f in self.frequencies
]

if not self.invert_impedance:
rx_list = [
nsem.receivers.PointNaturalSource(
rx_locs,
orientation=simpeg_mode,
component="apparent_resistivity",
),
nsem.receivers.PointNaturalSource(
rx_locs, orientation=simpeg_mode, component="phase"
),
]

src_list = [
nsem.sources.Planewave(rx_list, frequency=f)
for f in self.frequencies
]
else:
rx_list = [
nsem.receivers.PointNaturalSource(
rx_locs,
orientation=simpeg_mode,
component="real",
),
nsem.receivers.PointNaturalSource(
rx_locs,
orientation=simpeg_mode,
component="imag",
),
]

src_list = [
nsem.sources.Planewave(rx_list, frequency=f)
for f in self.frequencies
]
return nsem.Survey(src_list)

@property
Expand All @@ -132,7 +180,7 @@ def tm_survey(self):

return self._get_mode_sources(self.component_map["tm"]["simpeg"])

def _get_data_observations(self, mode):
def _get_data_observations(self, mode, impedance=False):
"""
get data
Expand All @@ -154,12 +202,12 @@ def _get_data_observations(self, mode):
phase = []
for ff in self.frequencies:
f_df = self.dataframe[self.dataframe.period == 1.0 / ff]
res.append(f_df[f"res_{comp}"])
# flip into the appropriate coordinate system
if comp in ["xy"]:
phase.append(1 * f_df[f"phase_{comp}"])
elif comp in ["yx"]:
phase.append(1 * f_df[f"phase_{comp}"])
if not self.invert_impedance:
res.append(f_df[f"res_{comp}"])
phase.append(f_df[f"phase_{comp}"])
else:
res.append(f_df[f"z_{comp}"].values.real)
phase.append(f_df[f"z_{comp}"].values.imag)

return np.hstack((res, phase)).flatten()

Expand Down Expand Up @@ -193,8 +241,12 @@ def _get_data_errors(self, mode):
# there is probably a more efficient method here using pandas
for ff in np.sort(self.frequencies):
f_df = self.dataframe[self.dataframe.period == 1.0 / ff]
res.append(f_df[f"res_{comp}_model_error"])
phase.append(f_df[f"phase_{comp}_model_error"])
if not self.invert_impedance:
res.append(f_df[f"res_{comp}_model_error"])
phase.append(f_df[f"phase_{comp}_model_error"])
else:
res.append(f_df[f"z_{comp}_model_error"])
phase.append(f_df[f"z_{comp}_model_error"])

return np.hstack((res, phase)).flatten()

Expand Down Expand Up @@ -260,32 +312,61 @@ def plot_response(self, **kwargs):
tm_data = self.tm_data.dobs.reshape(
(self.n_frequencies, 2, self.n_stations)
)
for ii in range(self.n_stations):
ax_xy_res.loglog(
1.0 / self.frequencies,
te_data[:, 0, ii],
color=(0.5, 0.5, ii / self.n_stations),
)
ax_xy_phase.semilogx(
1.0 / self.frequencies,
te_data[:, 1, ii],
color=(0.25, 0.25, ii / self.n_stations),
)
ax_yx_res.loglog(
1.0 / self.frequencies,
tm_data[:, 0, ii],
color=(0.5, ii / self.n_stations, 0.75),
)
ax_yx_phase.semilogx(
1.0 / self.frequencies,
tm_data[:, 1, ii],
color=(0.25, ii / self.n_stations, 0.75),
)

ax_xy_phase.set_xlabel("Period (s)")
ax_yx_phase.set_xlabel("Period (s)")
ax_xy_res.set_ylabel("Apparent Resistivity")
ax_xy_phase.set_ylabel("Phase")
if not self.invert_impedance:
for ii in range(self.n_stations):
ax_xy_res.loglog(
1.0 / self.frequencies,
te_data[:, 0, ii],
color=(0.5, 0.5, ii / self.n_stations),
)
ax_xy_phase.semilogx(
1.0 / self.frequencies,
te_data[:, 1, ii],
color=(0.25, 0.25, ii / self.n_stations),
)
ax_yx_res.loglog(
1.0 / self.frequencies,
tm_data[:, 0, ii],
color=(0.5, ii / self.n_stations, 0.75),
)
ax_yx_phase.semilogx(
1.0 / self.frequencies,
tm_data[:, 1, ii],
color=(0.25, ii / self.n_stations, 0.75),
)

ax_xy_phase.set_xlabel("Period (s)")
ax_yx_phase.set_xlabel("Period (s)")
ax_xy_res.set_ylabel("Apparent Resistivity")
ax_xy_phase.set_ylabel("Phase")
else:
for ii in range(self.n_stations):
ax_xy_res.loglog(
1.0 / self.frequencies,
te_data[:, 0, ii],
color=(0.5, 0.5, ii / self.n_stations),
)
ax_xy_phase.loglog(
1.0 / self.frequencies,
te_data[:, 1, ii],
color=(0.25, 0.25, ii / self.n_stations),
)
ax_yx_res.loglog(
1.0 / self.frequencies,
tm_data[:, 0, ii],
color=(0.5, ii / self.n_stations, 0.75),
)
ax_yx_phase.loglog(
1.0 / self.frequencies,
tm_data[:, 1, ii],
color=(0.25, ii / self.n_stations, 0.75),
)

ax_xy_phase.set_xlabel("Period (s)")
ax_yx_phase.set_xlabel("Period (s)")
ax_xy_res.set_ylabel("Real Impedance [Ohms]")
ax_xy_phase.set_ylabel("Imag Impedance [Ohms]")

for ax in [ax_xy_res, ax_xy_phase, ax_yx_res, ax_yx_phase]:
ax.grid(
Expand Down
8 changes: 5 additions & 3 deletions mtpy/modeling/simpeg/recipes/inversion_2d.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -322,7 +322,7 @@ def regularization(self):
reference_model=self.reference_model,
alpha_s=self.alpha_s,
alpha_x=self.alpha_y,
alpha_z=self.alpha_z,
alpha_y=self.alpha_z,
mapping=maps.IdentityMap(nP=self.mesh.number_of_active_cells),
)

Expand Down Expand Up @@ -425,7 +425,7 @@ def directives(self):
self.starting_beta,
self.beta_schedule,
self.saved_model_outputs,
self.target_misfit,
# self.target_misfit,
]

def run_inversion(self):
Expand Down Expand Up @@ -540,7 +540,9 @@ def plot_tikhonov_curve(self):
ax.plot(
xlim,
np.ones(2)
* (self.data.te_observations.size + self.data.tm_observations.size),
* (
self.data.te_observations.size + self.data.tm_observations.size
),
"--",
)
ax.set_xlim(xlim)
Expand Down
2 changes: 1 addition & 1 deletion tests/modeling/simpeg/test_simpeg_2d_inversion_recipe.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -145,7 +145,7 @@ def target_misfit(self):
)

def test_directives(self):
self.assertEqual(4, len(self.simpeg_inversion.directives))
self.assertEqual(3, len(self.simpeg_inversion.directives))


class TestSimpeg2DRecipeRun(unittest.TestCase):
Expand Down

0 comments on commit 6040c3e

Please sign in to comment.