mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-29 01:28:06 +08:00
Start of Jvec for problems j, h. failing at the moment
This commit is contained in:
Lindsey
+7
-2
@@ -94,6 +94,7 @@ class BaseEMProblem(Problem.BaseProblem):
|
||||
self._MeSigma = self.mesh.getEdgeInnerProduct(self.curModel.sigma)
|
||||
return self._MeSigma
|
||||
|
||||
# TODO: This should take a vector
|
||||
def MeSigmaDeriv(self, u):
|
||||
"""
|
||||
Deriv of MeSigma wrt sigma
|
||||
@@ -107,6 +108,7 @@ class BaseEMProblem(Problem.BaseProblem):
|
||||
self._MeSigmaI = self.mesh.getEdgeInnerProduct(self.curModel.sigma, invMat=True)
|
||||
return self._MeSigmaI
|
||||
|
||||
# TODO: This should take a vector
|
||||
def MeSigmaIDeriv(self, u):
|
||||
"""
|
||||
Deriv of MeSigma wrt sigma
|
||||
@@ -126,8 +128,9 @@ class BaseEMProblem(Problem.BaseProblem):
|
||||
self._MfRho = self.mesh.getFaceInnerProduct(self.curModel.rho)
|
||||
return self._MfRho
|
||||
|
||||
# TODO: This should take a vector
|
||||
def MfRhoDeriv(self,u):
|
||||
return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u)
|
||||
return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * self.curModel.rhoDeriv
|
||||
|
||||
@property
|
||||
def MfRhoI(self):
|
||||
@@ -135,5 +138,7 @@ class BaseEMProblem(Problem.BaseProblem):
|
||||
self._MfRhoI = self.mesh.getFaceInnerProduct(self.curModel.rho, invMat=True)
|
||||
return self._MfRhoI
|
||||
|
||||
# TODO: This isn't going to work yet
|
||||
# TODO: This should take a vector
|
||||
def dMfRhoIDeriv(self,u):
|
||||
return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)
|
||||
return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv
|
||||
|
||||
+83
-20
@@ -442,6 +442,7 @@ class ProblemFDEM_j(BaseFDEMProblem):
|
||||
|
||||
MeMuI = self.MeMuI
|
||||
MfRho = self.MfRho
|
||||
# print MfRho.max
|
||||
C = self.mesh.edgeCurl
|
||||
iomega = 1j * omega(freq) * sp.eye(self.mesh.nF)
|
||||
|
||||
@@ -463,19 +464,25 @@ class ProblemFDEM_j(BaseFDEMProblem):
|
||||
MeMuI = self.MeMuI
|
||||
MfRho = self.MfRho
|
||||
C = self.mesh.edgeCurl
|
||||
sigi = self.sigmai
|
||||
dsig_dm = self.curModel.transformDeriv
|
||||
dsigi_dsig = -Utils.sdiag(sigi)**2
|
||||
dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(u)
|
||||
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 dsig_dm.T * ( dsigi_dsig.T *( dMf_dsigi.T * ( C * ( MeMuI.T * ( C.T * v ) ) ) ) )
|
||||
return MfRhoDeriv_m.T * (C * (MeMuI.T * v))
|
||||
|
||||
if self._makeASymmetric is True:
|
||||
return MfRho.T * ( C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) )
|
||||
return C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) )
|
||||
# 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)))
|
||||
|
||||
|
||||
def getRHS(self, freq):
|
||||
@@ -497,9 +504,40 @@ class ProblemFDEM_j(BaseFDEMProblem):
|
||||
|
||||
return RHS
|
||||
|
||||
def getRHSDeriv(self, freq, v, adjoint=False):
|
||||
raise NotImplementedError('getRHSDeriv not implemented yet')
|
||||
return None
|
||||
def getRHSDeriv(self, src, v, adjoint=False):
|
||||
C = self.mesh.edgeCurl
|
||||
MeMuI = self.MeMuI
|
||||
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
|
||||
|
||||
if adjoint:
|
||||
if self._makeASymmetric:
|
||||
v = MfRho*v
|
||||
S_mDerivv = S_mDeriv(MeMuI.T * (C.T * v))
|
||||
S_eDerivv = S_eDeriv(v)
|
||||
if S_mDerivv is not None and S_eDerivv is not None:
|
||||
return S_mDerivv - 1j*omega(freq)*S_eDerivv
|
||||
elif S_mDerivv is not None:
|
||||
return S_mDerivv
|
||||
elif S_eDerivv is not None:
|
||||
return - 1j*omega(freq)*S_eDerivv
|
||||
else:
|
||||
return None
|
||||
else:
|
||||
S_mDerivv, S_eDerivv = S_mDeriv(v), S_eDeriv(v)
|
||||
|
||||
if S_mDerivv is not None and S_eDerivv is not None:
|
||||
RHSDeriv = C * (MeMuI * S_mDerivv) - 1j * omega(freq) * S_eDerivv
|
||||
elif S_mDerivv is not None:
|
||||
RHSDeriv = C * (MeMuI * S_mDerivv)
|
||||
elif S_eDerivv is not None:
|
||||
RHSDeriv = - 1j * omega(freq) * S_eDerivv
|
||||
else:
|
||||
return None
|
||||
|
||||
if self._makeASymmetric:
|
||||
return MfRho.T * RHSDeriv
|
||||
return RHSDeriv
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -552,16 +590,18 @@ 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
|
||||
# sigi = self.sigmai
|
||||
# dsig_dm = self.curModel.transformDeriv
|
||||
# dsigi_dsig = -Utils.sdiag(sigi)**2
|
||||
|
||||
dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(C*u)
|
||||
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 (C.T * (dMf_dsigi * (dsigi_dsig * (dsig_dm * v))))
|
||||
|
||||
# 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):
|
||||
"""
|
||||
@@ -578,7 +618,30 @@ class ProblemFDEM_h(BaseFDEMProblem):
|
||||
|
||||
return RHS
|
||||
|
||||
def getRHSDeriv(self, freq, v, adjoint=False):
|
||||
raise NotImplementedError('getRHSDeriv not implemented yet')
|
||||
return None
|
||||
def getRHSDeriv(self, src, v, adjoint=False):
|
||||
# raise NotImplementedError('getRHSDeriv not implemented yet')
|
||||
# return None
|
||||
_, S_e = self.getSourceTerm(src.freq)
|
||||
C = self.mesh.edgeCurl
|
||||
MfRho = self.MfRho
|
||||
MfRhoDeriv = self.MfRhoDeriv(S_e)
|
||||
|
||||
# if not adjoint:
|
||||
MfRhoDeriv_m = self.MfRhoDeriv(S_e)
|
||||
# elif adjoint:
|
||||
# MfRhoDeriv_m = self.MfRhoDeriv
|
||||
|
||||
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
|
||||
|
||||
RHSDeriv = C.T * (MfRhoDeriv * v)
|
||||
# = S_m + C.T * ( MfRho * S_e )
|
||||
|
||||
S_mDeriv = S_mDeriv(v)
|
||||
S_eDeriv = S_eDeriv(v)
|
||||
if S_mDeriv is not None:
|
||||
RHSDeriv += S_mDeriv(v)
|
||||
if S_eDeriv is not None:
|
||||
RHSDeriv += C.T * (MfRho * S_e)
|
||||
|
||||
return RHSDeriv
|
||||
|
||||
|
||||
+79
-24
@@ -157,10 +157,10 @@ class FieldsFDEM_b(FieldsFDEM):
|
||||
return self._MfMui.T * (self._edgeCurl * (self._MeSigmaI.T * v))
|
||||
|
||||
def _eSecondaryDeriv_m(self, src, v, adjoint=False):
|
||||
bsol = self[[src],'bSolution']
|
||||
bSolution = self[[src],'bSolution']
|
||||
_,S_e = src.eval(self.survey.prob)
|
||||
|
||||
w = self._edgeCurl.T * (self._MfMui * bsol)
|
||||
w = self._edgeCurl.T * (self._MfMui * bSolution)
|
||||
if S_e is not None:
|
||||
w += -S_e
|
||||
|
||||
@@ -207,6 +207,7 @@ class FieldsFDEM_j(FieldsFDEM):
|
||||
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)
|
||||
@@ -222,6 +223,13 @@ class FieldsFDEM_j(FieldsFDEM):
|
||||
def _j(self, jSolution, srcList):
|
||||
return self._jPrimary(jSolution, srcList) + self._jSecondary(jSolution, srcList)
|
||||
|
||||
def _jDeriv_u(self, src, v, adjoint=False):
|
||||
return None
|
||||
|
||||
def _jDeriv_m(self, src, v, adjoint=False):
|
||||
# assuming primary does not depend on the model
|
||||
return None
|
||||
|
||||
def _hPrimary(self, jSolution, srcList):
|
||||
hPrimary = np.zeros([self._edgeCurl.shape[1],jSolution.shape[1]],dtype = complex)
|
||||
for i, src in enumerate(srcList):
|
||||
@@ -242,27 +250,67 @@ class FieldsFDEM_j(FieldsFDEM):
|
||||
h[:,i] += 1./(1j*omega(src.freq)) * MeMuI * S_m
|
||||
return h
|
||||
|
||||
def _hSecondaryDeriv_u(self, src, v, adjoint=False):
|
||||
MeMuI = self._MeMuI
|
||||
C = self._edgeCurl
|
||||
MfRho = self._MfRho
|
||||
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))
|
||||
|
||||
def _hSecondaryDeriv_m(self, src, v, adjoint=False):
|
||||
jSolution = self[[src],'jSolution']
|
||||
MeMuI = self._MeMuI
|
||||
C = self._edgeCurl
|
||||
MfRho = self._MfRho
|
||||
MfRhoDeriv = self._MfRhoDeriv
|
||||
|
||||
if not adjoint:
|
||||
hDeriv_m = -1./(1j*omega(src.freq)) * MeMuI * (C.T * (MfRhoDeriv(jSolution)*v ) )
|
||||
elif adjoint:
|
||||
hDeriv_m = -1./(1j*omega(src.freq)) * MfRhoDeriv(jSolution).T * ( C * (MeMuI.T * v ) )
|
||||
|
||||
S_mDeriv,_ = src.evalDeriv(self.survey.prob, adjoint)
|
||||
|
||||
if not adjoint:
|
||||
S_mDeriv = S_mDeriv(v)
|
||||
if S_mDeriv is not None:
|
||||
hDeriv_m += 1./(1j*omega(src.freq)) * MeMuI * S_mDeriv
|
||||
elif adjoint:
|
||||
S_mDeriv = S_mDeriv(MeMuI.T * v)
|
||||
if S_mDeriv is not None:
|
||||
hDeriv_m += 1./(1j*omega(src.freq)) * S_mDeriv
|
||||
return h
|
||||
|
||||
|
||||
def _h(self, jSolution, srcList):
|
||||
return self._hPrimary(jSolution, srcList) + self._hSecondary(jSolution, srcList)
|
||||
|
||||
def _hDeriv(self, jSolution, srcList, v, adjoint=False):
|
||||
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
|
||||
# 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)
|
||||
# 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 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
|
||||
# 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)
|
||||
|
||||
def _hDeriv_m(self, src, v, adjoint=False):
|
||||
# assuming the primary doesn't depend on the model
|
||||
return _hSecondaryDeriv_u(self, src, v, adjoint)
|
||||
|
||||
|
||||
class FieldsFDEM_h(FieldsFDEM):
|
||||
@@ -298,6 +346,13 @@ class FieldsFDEM_h(FieldsFDEM):
|
||||
def _h(self, hSolution, srcList):
|
||||
return self._hPrimary(hSolution, srcList) + self._hSecondary(hSolution, srcList)
|
||||
|
||||
def _hDeriv_u(self, src, v, adjoint=False):
|
||||
return None
|
||||
|
||||
def _hDeriv_m(self, src, v, adjoint=False):
|
||||
# assuming primary does not depend on the model
|
||||
return None
|
||||
|
||||
def _jPrimary(self, hSolution, srcList):
|
||||
jPrimary = np.zeros([self._edgeCurl.shape[0], hSolution.shape[1]])
|
||||
for i, src in enumerate(srcList):
|
||||
@@ -326,8 +381,8 @@ class FieldsFDEM_h(FieldsFDEM):
|
||||
return - S_eDeriv
|
||||
|
||||
|
||||
# def calcFields(self, sol, freq, fieldType, adjoint=False):
|
||||
# j = sol
|
||||
# def calcFields(self, Solution, freq, fieldType, adjoint=False):
|
||||
# j = Solution
|
||||
# if fieldType == 'j':
|
||||
# return j
|
||||
# elif fieldType == 'h':
|
||||
@@ -341,8 +396,8 @@ class FieldsFDEM_h(FieldsFDEM):
|
||||
# return h
|
||||
# raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
# def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
|
||||
# j = sol
|
||||
# def calcFieldsDeriv(self, Solution, freq, fieldType, v, adjoint=False):
|
||||
# j = Solution
|
||||
# if fieldType == 'j':
|
||||
# return None
|
||||
# elif fieldType == 'h':
|
||||
@@ -361,8 +416,8 @@ class FieldsFDEM_h(FieldsFDEM):
|
||||
# raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
|
||||
# def calcFields(self, sol, freq, fieldType, adjoint=False):
|
||||
# h = sol
|
||||
# def calcFields(self, Solution, freq, fieldType, adjoint=False):
|
||||
# h = Solution
|
||||
# if fieldType == 'j':
|
||||
# C = self.mesh.edgeCurl
|
||||
# if adjoint:
|
||||
@@ -372,5 +427,5 @@ class FieldsFDEM_h(FieldsFDEM):
|
||||
# return h
|
||||
# raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
# def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
|
||||
# def calcFieldsDeriv(self, Solution, freq, fieldType, v, adjoint=False):
|
||||
# return None
|
||||
Reference in New Issue
Block a user