Zero, Identity for E formulation

This commit is contained in:
Lindsey Heagy
2015-11-06 10:00:30 -08:00
parent 6fb93ee6e4
commit d102ec9c00
3 changed files with 44 additions and 40 deletions
+33 -31
View File
@@ -73,26 +73,29 @@ class BaseFDEMProblem(BaseEMProblem):
ftype = self._fieldType + 'Solution'
u_src = f[src, ftype]
dA_dm = self.getADeriv_m(freq, u_src, v)
dRHS_dm = self.getRHSDeriv_m(src, v)
if dRHS_dm is None:
du_dm = dA_duI * ( - dA_dm )
else:
du_dm = dA_duI * ( - dA_dm + dRHS_dm )
dRHS_dm = self.getRHSDeriv_m(freq, src, v)
# if dRHS_dm is None:
# du_dm = dA_duI * ( - dA_dm )
# else:
du_dm = dA_duI * ( - dA_dm + dRHS_dm )
for rx in src.rxList:
# df_duFun = u.deriv_u(rx.fieldsUsed, m)
df_duFun = getattr(f, '_%sDeriv_u'%rx.projField, None)
df_du = df_duFun(src, du_dm, adjoint=False)
if df_du is not None:
du_dm = df_du
# if df_du is not None:
du_dm = df_du
df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None)
df_dm = df_dmFun(src, v, adjoint=False)
if df_dm is not None:
du_dm += df_dm
# if df_dm is not None:
du_dm += df_dm
P = lambda v: rx.projectFieldsDeriv(src, self.mesh, f, v) # wrt u, also have wrt m
Jv[src, rx] = P(du_dm)
return Utils.mkvc(Jv)
@@ -249,12 +252,11 @@ class ProblemFDEM_e(BaseFDEMProblem):
C = self.mesh.edgeCurl
MfMui = self.MfMui
# RHS = C.T * (MfMui * S_m) -1j * omega(freq) * Me * S_e
RHS = C.T * (MfMui * S_m) -1j * omega(freq) * S_e
return RHS
def getRHSDeriv_m(self, src, v, adjoint=False):
def getRHSDeriv_m(self, freq, src, v, adjoint=False):
C = self.mesh.edgeCurl
MfMui = self.MfMui
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
@@ -263,25 +265,25 @@ class ProblemFDEM_e(BaseFDEMProblem):
dRHS = MfMui * (C * v)
S_mDerivv = S_mDeriv(dRHS)
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
# 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:
return C.T * (MfMui * S_mDerivv) -1j * omega(freq) * S_eDerivv
elif S_mDerivv is not None:
return C.T * (MfMui * S_mDerivv)
elif S_eDerivv is not None:
return -1j * omega(freq) * S_eDerivv
else:
return None
# if S_mDerivv is not None and S_eDerivv is not None:
return C.T * (MfMui * S_mDerivv) -1j * omega(freq) * S_eDerivv
# elif S_mDerivv is not None:
# return C.T * (MfMui * S_mDerivv)
# elif S_eDerivv is not None:
# return -1j * omega(freq) * S_eDerivv
# else:
# return None
class ProblemFDEM_b(BaseFDEMProblem):
@@ -372,7 +374,7 @@ class ProblemFDEM_b(BaseFDEMProblem):
return RHS
def getRHSDeriv_m(self, src, v, adjoint=False):
def getRHSDeriv_m(self, freq, src, v, adjoint=False):
C = self.mesh.edgeCurl
S_m, S_e = src.eval(self)
MfMui = self.MfMui
@@ -519,7 +521,7 @@ class ProblemFDEM_j(BaseFDEMProblem):
return RHS
def getRHSDeriv_m(self, src, v, adjoint=False):
def getRHSDeriv_m(self, freq, src, v, adjoint=False):
C = self.mesh.edgeCurl
MeMuI = self.MeMuI
S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint)
@@ -627,7 +629,7 @@ class ProblemFDEM_h(BaseFDEMProblem):
return RHS
def getRHSDeriv_m(self, src, v, adjoint=False):
def getRHSDeriv_m(self, freq, src, v, adjoint=False):
_, S_e = src.eval(self)
C = self.mesh.edgeCurl
MfRho = self.MfRho
+8 -7
View File
@@ -1,5 +1,6 @@
from SimPEG import Survey, Problem, Utils, np, sp
from SimPEG.EM.Utils.EMUtils import omega
from SimPEG.Utils import Zero, Identity
class FieldsFDEM(Problem.Fields):
@@ -40,11 +41,11 @@ class FieldsFDEM_e(FieldsFDEM):
return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList)
def _eDeriv_u(self, src, v, adjoint = False):
return None
return Identity()*v
def _eDeriv_m(self, src, v, adjoint = False):
# assuming primary does not depend on the model
return None
return Zero()
def _bPrimary(self, eSolution, srcList):
bPrimary = np.zeros([self._edgeCurl.shape[0],eSolution.shape[1]],dtype = complex)
@@ -73,9 +74,9 @@ class FieldsFDEM_e(FieldsFDEM):
def _bSecondaryDeriv_m(self, src, v, adjoint = False):
S_mDeriv, _ = src.evalDeriv(self.prob, adjoint)
S_mDeriv = S_mDeriv(v)
if S_mDeriv is not None:
return 1./(1j * omega(src.freq)) * S_mDeriv
return None
# if S_mDeriv is not None:
return 1./(1j * omega(src.freq)) * S_mDeriv
# return None
def _b(self, eSolution, srcList):
return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList)
@@ -126,11 +127,11 @@ class FieldsFDEM_b(FieldsFDEM):
return self._bPrimary(bSolution, srcList) + self._bSecondary(bSolution, srcList)
def _bDeriv_u(self, src, v, adjoint=False):
return None
return Identity()*v
def _bDeriv_m(self, src, v, adjoint=False):
# assuming primary does not depend on the model
return None
return Zero()
def _ePrimary(self, bSolution, srcList):
ePrimary = np.zeros([self._edgeCurl.shape[1],bSolution.shape[1]],dtype = complex)
+3 -2
View File
@@ -2,6 +2,7 @@ from SimPEG import Survey, Problem, Utils, np, sp
from SimPEG.EM.Utils import SrcUtils
from SimPEG.EM.Utils.EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b
from scipy.constants import mu_0
from SimPEG.Utils import Zero, Identity
####################################################
# Receivers
@@ -123,10 +124,10 @@ class SrcFDEM(Survey.BaseSrc):
return None
def S_mDeriv(self, prob, v, adjoint = False):
return None
return Zero()
def S_eDeriv(self, prob, v, adjoint = False):
return None
return Zero()
class SrcFDEM_RawVec_e(SrcFDEM):