mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-30 08:27:46 +08:00
FDEM problems: e,b,h,j up and running with Jvec and Jtvec within the re-factored framework. Whenever a new element is created, we create its derive wrt u (the computed field) and wrt m (the model). Jvec and Jtvec then stitch these pieces together using chain rule
This commit is contained in:
Lindsey Heagy
+2
-1
@@ -130,7 +130,8 @@ class BaseEMProblem(Problem.BaseProblem):
|
||||
|
||||
# TODO: This should take a vector
|
||||
def MfRhoDeriv(self,u):
|
||||
return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * self.curModel.rhoDeriv
|
||||
return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * -Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv
|
||||
# self.curModel.rhoDeriv
|
||||
|
||||
@property
|
||||
def MfRhoI(self):
|
||||
|
||||
+8
-27
@@ -36,7 +36,7 @@ class BaseFDEMProblem(BaseEMProblem):
|
||||
|
||||
def Jvec(self, m, v, f=None):
|
||||
if f is None:
|
||||
f = self.fields(m) # rename to f?
|
||||
f = self.fields(m)
|
||||
|
||||
self.curModel = m
|
||||
|
||||
@@ -238,8 +238,6 @@ class ProblemFDEM_e(BaseFDEMProblem):
|
||||
C = self.mesh.edgeCurl
|
||||
MfMui = self.MfMui
|
||||
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
|
||||
# # evalDeriv(MfMui.T* (C * v), adjoint)
|
||||
# raise Exception('Not implemented')
|
||||
|
||||
if adjoint:
|
||||
dRHS = MfMui * (C * v)
|
||||
@@ -303,19 +301,14 @@ class ProblemFDEM_b(BaseFDEMProblem):
|
||||
|
||||
MeSigmaIDeriv = MeSigmaIDeriv(vec)
|
||||
|
||||
# dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(vec)
|
||||
|
||||
if adjoint:
|
||||
if self._makeASymmetric is True:
|
||||
v = MfMui * v
|
||||
return MeSigmaIDeriv.T * (C.T * v)
|
||||
# dsig_dm.T * ( dMe_dsig.T * ( dMeSigmaI_dI.T * ( C.T * v ) ) )
|
||||
|
||||
if self._makeASymmetric is True:
|
||||
# return MfMui.T * ( C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) ) )
|
||||
return MfMui.T * ( C * ( MeSigmaIDeriv * v ) )
|
||||
return C * ( MeSigmaIDeriv * v )
|
||||
# C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) )
|
||||
|
||||
|
||||
def getRHS(self, freq):
|
||||
@@ -456,22 +449,13 @@ class ProblemFDEM_j(BaseFDEMProblem):
|
||||
C = self.mesh.edgeCurl
|
||||
MfRhoDeriv_m = self.MfRhoDeriv(u)
|
||||
|
||||
# sigi = self.curModel.rho
|
||||
# dsig_dm = self.curModel.rhoDeriv
|
||||
# dsigi_dsig = -Utils.sdiag(sigi)**2
|
||||
# dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(u)
|
||||
# MfRhoDeriv_m = (dMf_dsigi * ( dsigi_dsig * ( dsig_dm)))
|
||||
|
||||
if adjoint:
|
||||
if self._makeASymmetric is True:
|
||||
v = MfRho * v
|
||||
# return dsig_dm.T * ( dsigi_dsig.T *( dMf_dsigi.T * ( C * ( MeMuI.T * ( C.T * v ) ) ) ) )
|
||||
return MfRhoDeriv_m.T * (C * (MeMuI.T * v))
|
||||
return MfRhoDeriv_m.T * (C * (MeMuI.T * (C.T * v)))
|
||||
|
||||
if self._makeASymmetric is True:
|
||||
# return MfRho.T * ( C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) )
|
||||
return MfRho.T * (C * ( MeMuI * (C.T * (MfRhoDeriv_m * v) )))
|
||||
# return C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) )
|
||||
return C * (MeMuI * (C.T * (MfRhoDeriv_m * v)))
|
||||
|
||||
|
||||
@@ -501,6 +485,7 @@ class ProblemFDEM_j(BaseFDEMProblem):
|
||||
|
||||
if adjoint:
|
||||
if self._makeASymmetric:
|
||||
MfRho = self.MfRho
|
||||
v = MfRho*v
|
||||
S_mDerivv = S_mDeriv(MeMuI.T * (C.T * v))
|
||||
S_eDerivv = S_eDeriv(v)
|
||||
@@ -525,6 +510,7 @@ class ProblemFDEM_j(BaseFDEMProblem):
|
||||
return None
|
||||
|
||||
if self._makeASymmetric:
|
||||
MfRho = self.MfRho
|
||||
return MfRho.T * RHSDeriv
|
||||
return RHSDeriv
|
||||
|
||||
@@ -580,17 +566,10 @@ class ProblemFDEM_h(BaseFDEMProblem):
|
||||
|
||||
MeMu = self.MeMu
|
||||
C = self.mesh.edgeCurl
|
||||
# sigi = self.sigmai
|
||||
# dsig_dm = self.curModel.transformDeriv
|
||||
# dsigi_dsig = -Utils.sdiag(sigi)**2
|
||||
|
||||
MfRhoDeriv_m = self.MfRhoDeriv(C*u)
|
||||
# dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(C*u)
|
||||
|
||||
if adjoint:
|
||||
# return (dsig_dm.T * (dsigi_dsig.T * (dMf_dsigi.T * (C * v))))
|
||||
return MfRhoDeriv_m.T * (C * v)
|
||||
# return (C.T * (dMf_dsigi * (dsigi_dsig * (dsig_dm * v))))
|
||||
return C.T * (MfRhoDeriv_m * v)
|
||||
|
||||
def getRHS(self, freq):
|
||||
@@ -623,8 +602,10 @@ class ProblemFDEM_h(BaseFDEMProblem):
|
||||
|
||||
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
|
||||
|
||||
RHSDeriv = C.T * (MfRhoDeriv * v)
|
||||
# = S_m + C.T * ( MfRho * S_e )
|
||||
if not adjoint:
|
||||
RHSDeriv = C.T * (MfRhoDeriv * v)
|
||||
elif adjoint:
|
||||
RHSDeriv = MfRhoDeriv.T * (C * v)
|
||||
|
||||
S_mDeriv = S_mDeriv(v)
|
||||
S_eDeriv = S_eDeriv(v)
|
||||
|
||||
+24
-81
@@ -206,11 +206,10 @@ class FieldsFDEM_j(FieldsFDEM):
|
||||
self._edgeCurl = self.survey.prob.mesh.edgeCurl
|
||||
self._MeMuI = self.survey.prob.MeMuI
|
||||
self._MfRho = self.survey.prob.MfRho
|
||||
self._curModel = self.survey.prob.curModel
|
||||
self._MfRhoDeriv = self.survey.prob.MfRhoDeriv
|
||||
|
||||
def _jPrimary(self, jSolution, srcList):
|
||||
jPrimary = np.zeros_like(jSolution)
|
||||
jPrimary = np.zeros_like(jSolution,dtype = complex)
|
||||
for i, src in enumerate(srcList):
|
||||
jp = src.jPrimary(self.survey.prob)
|
||||
if jp is not None:
|
||||
@@ -257,7 +256,7 @@ class FieldsFDEM_j(FieldsFDEM):
|
||||
if not adjoint:
|
||||
return -1./(1j*omega(src.freq)) * MeMuI * (C.T * (MfRho * v) )
|
||||
elif adjoint:
|
||||
return -1./(1j*omega(src.freq)) * MfRho * (C * ( MeMuI .T * v))
|
||||
return -1./(1j*omega(src.freq)) * MfRho.T * (C * ( MeMuI.T * v))
|
||||
|
||||
def _hSecondaryDeriv_m(self, src, v, adjoint=False):
|
||||
jSolution = self[[src],'jSolution']
|
||||
@@ -281,36 +280,18 @@ class FieldsFDEM_j(FieldsFDEM):
|
||||
S_mDeriv = S_mDeriv(MeMuI.T * v)
|
||||
if S_mDeriv is not None:
|
||||
hDeriv_m += 1./(1j*omega(src.freq)) * S_mDeriv
|
||||
return h
|
||||
return hDeriv_m
|
||||
|
||||
|
||||
def _h(self, jSolution, srcList):
|
||||
return self._hPrimary(jSolution, srcList) + self._hSecondary(jSolution, srcList)
|
||||
|
||||
# raise NotImplementedError('Fields Derivs Not Implemented Yet')
|
||||
# sig = self._curModel.transform
|
||||
# sigi = 1/sig
|
||||
# dsig_dm = self._curModel.transformDeriv
|
||||
# dsigi_dsig = -Utils.sdiag(sigi)**2
|
||||
# dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j)
|
||||
# sigi = self._MfRho
|
||||
|
||||
# S_mDeriv,_ = src.getSourceDeriv(self.survey.prob, v, adjoint)
|
||||
|
||||
# if not adjoint:
|
||||
# h_Deriv= -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) )
|
||||
# else:
|
||||
# h_Deriv= -(1./(1j*omega(freq))) * dsig_dm.T * ( dsigi_dsig.T * ( dMf_dsigi.T * ( C * ( MeMuI.T * v ) ) ) )
|
||||
|
||||
# if S_mDeriv is not None:
|
||||
# return 1./(1j*omega(src.freq)) * S_mDeriv + h_Deriv
|
||||
|
||||
def _hDeriv_u(self, src, v, adjoint=False):
|
||||
return _hSecondaryDeriv_u(self, src, v, adjoint)
|
||||
return self._hSecondaryDeriv_u(src, v, adjoint)
|
||||
|
||||
def _hDeriv_m(self, src, v, adjoint=False):
|
||||
# assuming the primary doesn't depend on the model
|
||||
return _hSecondaryDeriv_u(self, src, v, adjoint)
|
||||
return self._hSecondaryDeriv_m(src, v, adjoint)
|
||||
|
||||
|
||||
class FieldsFDEM_h(FieldsFDEM):
|
||||
@@ -333,7 +314,7 @@ class FieldsFDEM_h(FieldsFDEM):
|
||||
self._MfRho = self.survey.prob.MfRho
|
||||
|
||||
def _hPrimary(self, hSolution, srcList):
|
||||
hPrimary = np.zeros_like(hSolution)
|
||||
hPrimary = np.zeros_like(hSolution,dtype = complex)
|
||||
for i, src in enumerate(srcList):
|
||||
hp = src.hPrimary(self.survey.prob)
|
||||
if hp is not None:
|
||||
@@ -369,63 +350,25 @@ class FieldsFDEM_h(FieldsFDEM):
|
||||
j[:,i] += -S_e
|
||||
return j
|
||||
|
||||
def _jSecondaryDeriv_u(self, src, v, adjoint=False):
|
||||
if not adjoint:
|
||||
return self._edgeCurl*v
|
||||
elif adjoint:
|
||||
return self._edgeCurl.T*v
|
||||
|
||||
def _jSecondaryDeriv_m(self, src, v, adjoint=False):
|
||||
_,S_eDeriv = src.evalDeriv(self.survey.prob, adjoint)
|
||||
S_eDeriv = S_eDeriv(v)
|
||||
if S_eDeriv is not None:
|
||||
return -S_eDeriv
|
||||
return None
|
||||
|
||||
def _j(self, hSolution, srcList):
|
||||
return self._jPrimary(hSolution, srcList) + self._jSecondary(hSolution, srcList)
|
||||
|
||||
def _jDeriv(self, hSolution, srcList, v, adjoint=False):
|
||||
raise NotImplementedError('Fields Derivs Not Implemented Yet')
|
||||
_,S_eDeriv = src.getSourceDeriv(self.survey.prob, v, adjoint)
|
||||
if S_eDeriv is None:
|
||||
return None
|
||||
else:
|
||||
return - S_eDeriv
|
||||
def _jDeriv_u(self, src, v, adjoint=False):
|
||||
return self._jSecondaryDeriv_u(src,v,adjoint)
|
||||
|
||||
|
||||
# def calcFields(self, Solution, freq, fieldType, adjoint=False):
|
||||
# j = Solution
|
||||
# if fieldType == 'j':
|
||||
# return j
|
||||
# elif fieldType == 'h':
|
||||
# MeMuI = self._MeMuI
|
||||
# C = self.mesh.edgeCurl
|
||||
# MfRho = self._MfRho
|
||||
# if not adjoint:
|
||||
# h = -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( MfRho * j ) )
|
||||
# else:
|
||||
# h = -(1./(1j*omega(freq))) * MfRho.T * ( C * ( MeMuI.T * j ) )
|
||||
# return h
|
||||
# raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
# def calcFieldsDeriv(self, Solution, freq, fieldType, v, adjoint=False):
|
||||
# j = Solution
|
||||
# if fieldType == 'j':
|
||||
# return None
|
||||
# elif fieldType == 'h':
|
||||
# MeMuI = self._MeMuI
|
||||
# C = self.mesh.edgeCurl
|
||||
# sig = self._curModel.transform
|
||||
# sigi = 1/sig
|
||||
# dsig_dm = self._curModel.transformDeriv
|
||||
# dsigi_dsig = -Utils.sdiag(sigi)**2
|
||||
# dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j)
|
||||
# sigi = self._MfRho
|
||||
# if not adjoint:
|
||||
# return -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) )
|
||||
# else:
|
||||
# return -(1./(1j*omega(freq))) * dsig_dm.T * ( dsigi_dsig.T * ( dMf_dsigi.T * ( C * ( MeMuI.T * v ) ) ) )
|
||||
# raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
|
||||
# def calcFields(self, Solution, freq, fieldType, adjoint=False):
|
||||
# h = Solution
|
||||
# if fieldType == 'j':
|
||||
# C = self.mesh.edgeCurl
|
||||
# if adjoint:
|
||||
# return C.T*h
|
||||
# return C*h
|
||||
# elif fieldType == 'h':
|
||||
# return h
|
||||
# raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
# def calcFieldsDeriv(self, Solution, freq, fieldType, v, adjoint=False):
|
||||
# return None
|
||||
def _jDeriv_m(self, src, v, adjoint=False):
|
||||
# assuming the primary does not depend on the model
|
||||
return self._jSecondaryDeriv_m(src,v,adjoint)
|
||||
@@ -1,7 +1,7 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
import simpegEM as EM
|
||||
import sys
|
||||
import sys
|
||||
from scipy.constants import mu_0
|
||||
import copy
|
||||
|
||||
@@ -9,7 +9,7 @@ testDerivs = True
|
||||
testCrossCheck = True
|
||||
testAdjoint = True
|
||||
testEB = True
|
||||
testHJ = False
|
||||
testHJ = True
|
||||
|
||||
verbose = False
|
||||
|
||||
@@ -35,8 +35,8 @@ def getProblem(fdemType, comp):
|
||||
x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the source
|
||||
XYZ = Utils.ndgrid(x,x,np.r_[0.])
|
||||
Rx0 = EM.FDEM.RxFDEM(XYZ, comp)
|
||||
Src0 = EM.FDEM.SrcFDEM_MagDipole([Rx0], loc=np.r_[0.,0.,0.], freq=freq[0])
|
||||
Src1 = EM.FDEM.SrcFDEM_MagDipole([Rx0], loc=np.r_[0.,0.,0.], freq=freq[1])
|
||||
Src0 = EM.FDEM.SrcFDEM_MagDipole([Rx0], freq=freq[0], loc=np.r_[0.,0.,0.])
|
||||
Src1 = EM.FDEM.SrcFDEM_MagDipole([Rx0], freq=freq[1], loc=np.r_[0.,0.,0.])
|
||||
|
||||
survey = EM.FDEM.SurveyFDEM([Src0, Src1])
|
||||
|
||||
@@ -69,7 +69,7 @@ def adjointTest(fdemType, comp):
|
||||
print 'Adjoint %s formulation - %s' % (fdemType, comp)
|
||||
|
||||
m = np.log(np.ones(prb.mesh.nC)*CONDUCTIVITY)
|
||||
mu = np.log(np.ones(prb.mesh.nC)*MU)
|
||||
mu = np.log(np.ones(prb.mesh.nC))*MU
|
||||
|
||||
if addrandoms is True:
|
||||
m = m + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-1
|
||||
|
||||
Reference in New Issue
Block a user