wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-08-07 11:28:42 +08:00
Code Issues Packages Projects Releases Wiki Activity

breaking up tests, also pass in u to df_dm

Browse Source
This commit is contained in:
Lindsey Heagy
2015-12-10 08:51:06 -08:00
parent 9e69455d6e
commit 22a6310a59
9 changed files with 446 additions and 296 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/EM/FDEM/FDEM.py
+2 -2
View File
@@ -78,10 +78,10 @@ class BaseFDEMProblem(BaseEMProblem):
for rx in src.rxList:
df_duFun = getattr(f, '_%sDeriv_u'%rx.projField, None)
df_dudu_dm = df_duFun(src, du_dm, adjoint=False)
df_dudu_dm = df_duFun(src, du_dm, u_src, adjoint=False)
df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None)
df_dm = df_dmFun(src, v, adjoint=False)
df_dm = df_dmFun(src, v, u_src, adjoint=False)
Df_Dm = np.array(df_dudu_dm + df_dm,dtype=complex)
SimPEG/EM/FDEM/FieldsFDEM.py
+36 -38
View File
@@ -32,10 +32,9 @@ class Fields_e(Fields):
self._edgeCurl = self.survey.prob.mesh.edgeCurl
self._aveE2CCV = self.survey.prob.mesh.aveE2CCV
self._aveF2CCV = self.survey.prob.mesh.aveF2CCV
self._sigma = self.survey.prob.curModel.sigma
self._sigmaDeriv = self.survey.prob.curModel.sigmaDeriv
self._mui = self.survey.prob.curModel.mui
self._nC = self.survey.prob.mesh.nC
self._MeSigma = self.survey.prob.MeSigma
self._MeSigmaDeriv = self.survey.prob.MeSigmaDeriv
def _GLoc(self,fieldType):
if fieldType == 'e':
@@ -60,10 +59,10 @@ class Fields_e(Fields):
def _e(self, eSolution, srcList):
return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList)
def _eDeriv_u(self, src, v, adjoint = False):
def _eDeriv_u(self, src, v, eSolution, adjoint = False):
return Identity()*v
def _eDeriv_m(self, src, v, adjoint = False):
def _eDeriv_m(self, src, v, eSolution, adjoint = False):
# assuming primary does not depend on the model
return Zero()
@@ -83,13 +82,13 @@ class Fields_e(Fields):
b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * S_m
return b
def _bSecondaryDeriv_u(self, src, v, adjoint = False):
def _bSecondaryDeriv_u(self, src, v, eSolution, adjoint = False):
C = self._edgeCurl
if adjoint:
return - 1./(1j*omega(src.freq)) * (C.T * v)
return - 1./(1j*omega(src.freq)) * (C * v)
def _bSecondaryDeriv_m(self, src, v, adjoint = False):
def _bSecondaryDeriv_m(self, src, v, eSolution, adjoint = False):
S_mDeriv, _ = src.evalDeriv(self.prob, adjoint)
S_mDeriv = S_mDeriv(v)
return 1./(1j * omega(src.freq)) * S_mDeriv
@@ -97,27 +96,50 @@ class Fields_e(Fields):
def _b(self, eSolution, srcList):
return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList)
def _bDeriv_u(self, src, v, adjoint=False):
def _bDeriv_u(self, src, v, eSolution, adjoint = False):
# Primary does not depend on u
return self._bSecondaryDeriv_u(src, v, adjoint)
def _bDeriv_m(self, src, v, adjoint=False):
def _bDeriv_m(self, src, v, eSolution, adjoint = False):
# Assuming the primary does not depend on the model
return self._bSecondaryDeriv_m(src, v, adjoint)
def _j(self, eSolution, srcList):
sigma = self._sigma
aveE2CCV = self._aveE2CCV
n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy
# Sigma = sdiag(np.kron(np.ones(n), sigma))
Sigma = self.prob.MeSigma
Sigma = self._MeSigma
VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol))
e = self._e(eSolution, srcList)
return VI * (aveE2CCV * (Sigma *e) )
def _h(self, eolution, srcList):
def _jDeriv_u(self, src, eSolution, v, adjoint = False):
aveE2CCV = self._aveE2CCV
n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy
Sigma = self._MeSigma
VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol))
if not adjoint:
return VI * (aveE2CCV * (Sigma * (self._eDeriv_u(src, v, adjoint) ) ) )
return self._eDeriv_u(src, Sigma.T * (aveE2CCV.T * (VI.T * v) ), adjoint)
def _jDeriv_m(self, src, v, eSolution, adjoint = False):
aveE2CCV = self._aveE2CCV
Sigma = self._MeSigma
SigmaDeriv = self._MeSigmaDeriv
e = self._e(eSolution, [src])
n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy
VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol))
if not adjoint:
return VI * (aveE2CCV * ( SigmaDeriv(e) * v + self._eDeriv_m(src, v, adjoint) ))
return SigmaDeriv(aveE2CCV.T * (VI.T * e), adjoint) * v + self._eDeriv_m(src, aveE2CCV.T * (VI.T * v), adjoint)
def _h(self, eSolution, srcList):
b = self._b(eSolution, srcList)
Mui = self.survey.prob.MfMui
aveF2CCV = self._aveF2CCV
@@ -129,30 +151,6 @@ class Fields_e(Fields):
return VI * (aveF2CCV * (Mui * b))
def _jDeriv_u(self, src, v, adjoint=False):
raise NotImplementedError
sigma = self._sigma
aveE2CCV = self._aveE2CCV
n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy
Sigma = sdiag(sp.kron(np.ones(n), sigma))
if not adjoint:
return Sigma * (aveE2CCV * (v + self._eDeriv_u(src, v, adjoint)))
return aveE2CCV.T * Sigma.T * v
def _jDeriv_m(self, src, v, adjoint=False):
raise NotImplementedError
sigma = self._sigma
aveE2CCV = self._aveE2CCV
n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy
Sigma = sdiag(sp.kron(np.ones(n), sigma))
if not adjoint:
dsigma_dm = self._sigmaDeriv(v)
dSigma_dm = sdiag(sp.kron(np.ones(n), dsigma_dm))
class Fields_b(Fields):
knownFields = {'bSolution':'F'}
SimPEG/EM/Utils/testingUtils.py
+47 -1
View File
@@ -4,6 +4,13 @@ from SimPEG import EM
import sys
from scipy.constants import mu_0
FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
CONDUCTIVITY = 1e1
MU = mu_0
freq = 5e-1
addrandoms = False
def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False):
cs = 10.
ncx, ncy, ncz = 0, 0, 0
@@ -72,4 +79,43 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False):
except ImportError, e:
prb.Solver = SolverLU
return prb
return prb
def crossCheckTest(SrcList, fdemType1, fdemType2, comp, TOL=1e-5, verbose=False):
l2norm = lambda r: np.sqrt(r.dot(r))
prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, verbose)
mesh = prb1.mesh
print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp)
m = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
mu = np.log(np.ones(mesh.nC)*MU)
if addrandoms is True:
m = m + np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1
mu = mu + np.random.randn(mesh.nC)*MU*1e-1
# prb1.PropMap.PropModel.mu = mu
# prb1.PropMap.PropModel.mui = 1./mu
survey1 = prb1.survey
d1 = survey1.dpred(m)
if verbose:
print ' Problem 1 solved'
prb2 = getFDEMProblem(fdemType2, comp, SrcList, freq, verbose)
# prb2.mu = mu
survey2 = prb2.survey
d2 = survey2.dpred(m)
if verbose:
print ' Problem 2 solved'
r = d2-d1
l2r = l2norm(r)
tol = np.max([TOL*(10**int(np.log10(l2norm(d1)))),FLR])
print l2norm(d1), l2norm(d2), l2r , tol, l2r < tol
return l2r < tol