simpegEM1D
simpegEM1D copied to clipboard
Half space option causes TDsurvey.dpred() fail
I may not be doing this correctly...
Ive attached a notebook cell, and the error messages below.
The only difference is half_switch = True
in the instantiation of TDsurvey. Looks like the ExpMap is confused about the 19 layer model, and the request for only the results due to one "layer".
Perhaps a separate issue, but related. I have not been able to create a single layer mesh to represent a half space instead of a multi-layer mesh with a half_switch logical.
It would be a good idea to be able to pass a single layer model to the forward and sensitivity methods. Inside simply use if (half_switch or n_layer == 1):
instead of forcing the user to specify the half_switch argument?
Here's a notebook cell
from SimPEG import Mesh, Maps, Utils
%matplotlib notebook
%load_ext autoreload
%autoreload 2
import matplotlib.pyplot as plt
from simpegEM1D import (
EM1D, EM1DSurveyTD, EM1DAnalytics, piecewise_pulse, set_mesh_1d,
skytem_HM_2015, skytem_LM_2015, get_vertical_discretization_time
)
import numpy as np
from scipy import io
from scipy.interpolate import interp1d
from geobipy import StatArray
wave_HM = skytem_HM_2015()
wave_LM = skytem_LM_2015()
time_HM = wave_HM.time_gate_center[0::2]
time_LM = wave_LM.time_gate_center[0::2]
hz = get_vertical_discretization_time(
np.unique(np.r_[time_HM, time_LM]), facter_tmax=0.5, factor_tmin=10.,
n_layer=19
)
mesh1D = set_mesh_1d(hz)
depth = -mesh1D.gridN[:-1]
LocSigZ = -mesh1D.gridCC
time_input_currents_HM = wave_HM.current_times[-7:]
input_currents_HM = wave_HM.currents[-7:]
time_input_currents_LM = wave_LM.current_times[-13:]
input_currents_LM = wave_LM.currents[-13:]
TDsurvey = EM1DSurveyTD(
rx_location=np.array([0., 0., 0.]),
src_location=np.array([0., 0., 0.]),
topo=np.r_[0., 0., 0.],
depth=depth,
rx_type='dBzdt',
wave_type='general',
src_type='CircularLoop',
a=13.,
I=1.,
time=time_HM,
time_input_currents=time_input_currents_HM,
input_currents=input_currents_HM,
n_pulse=2,
base_frequency=25.,
use_lowpass_filter=True,
high_cut_frequency=7e4,
moment_type='dual',
time_dual_moment=time_LM,
time_input_currents_dual_moment=time_input_currents_LM,
input_currents_dual_moment=input_currents_LM,
base_frequency_dual_moment=210,
half_switch = True
)
sig_half=1e-2
chi_half=0.
expmap = Maps.ExpMap(mesh1D)
m_1D = np.log(np.arange(TDsurvey.n_layer)*sig_half)
chi = np.zeros(TDsurvey.n_layer)
prob = EM1D(
mesh1D, sigmaMap=expmap, chi=chi
)
prob.unpair()
prob.pair(TDsurvey)
dBzdtTD = TDsurvey.dpred(m_1D)
This causes
<ipython-input-7-6d55b4039c6a> in <module>()
1 prob.unpair()
2 prob.pair(TDsurvey)
----> 3 dBzdtTD = TDsurvey.dpred(m_1D)
4 err = np.linalg.norm(p_saved-dBzdtTD)
5 print(err, err < 1e-12)
/Users/nfoks/anaconda/lib/python3.5/site-packages/SimPEG/Utils/codeutils.py in requiresVarWrapper(self, *args, **kwargs)
214 if getattr(self, var, None) is None:
215 raise Exception(extra)
--> 216 return f(self, *args, **kwargs)
217
218 doc = requiresVarWrapper.__doc__
/Users/nfoks/anaconda/lib/python3.5/site-packages/simpegEM1D/Survey.py in dpred(self, m, f)
535 """
536 if f is None:
--> 537 f = self.prob.fields(m)
538
539 return self._pred
/Users/nfoks/anaconda/lib/python3.5/site-packages/simpegEM1D/EM1D.py in fields(self, m)
422 # @profile
423 def fields(self, m):
--> 424 f = self.forward(m, output_type='response')
425 self.survey._pred = Utils.mkvc(self.survey.projectFields(f))
426 return f
/Users/nfoks/anaconda/lib/python3.5/site-packages/simpegEM1D/EM1D.py in forward(self, m, output_type)
301
302 # TODO: potentially store
--> 303 sig = self.sigma_cole()
304
305 if output_type == 'response':
/Users/nfoks/anaconda/lib/python3.5/site-packages/simpegEM1D/EM1D.py in sigma_cole(self)
228 )
229
--> 230 sigma = np.tile(self.sigma.reshape([-1, 1]), (1, n_frequency))
231 if np.isscalar(self.eta):
232 eta = self.eta
/Users/nfoks/anaconda/lib/python3.5/site-packages/SimPEG/Props.py in fget(self)
211 )
212 )
--> 213 return mapping * self.model
214
215 def fset(self, value):
/Users/nfoks/anaconda/lib/python3.5/site-packages/SimPEG/Maps.py in __mul__(self, val)
192 raise ValueError(
193 'Dimension mismatch in {0!s} and np.ndarray{1!s}.'.format(
--> 194 str(self), str(val.shape)
195 )
196 )
ValueError: Dimension mismatch in ExpMap(19,19) and np.ndarray(1,).