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

derivatives for fields objects, UNTESTED. TODO: clean up how we get derivs of mass matrices wrt the physical properties

Browse Source
This commit is contained in:
Lindsey
2015-04-17 11:21:07 -07:00
parent 93778d13f1
commit 2864a976c6
3 changed files with 83 additions and 59 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegEM/Base.py
+48 -23
View File
@@ -48,59 +48,84 @@ class BaseEMProblem(Problem.BaseProblem):
# Mass Matrices
####################################################
@property
def MfMui(self):
if getattr(self, '_MfMui', None) is None:
self._MfMui = self.mesh.getFaceInnerProduct(1/self.mu)
return self._MfMui
@property
def MeMuI(self):
# TODO: Assuming isotropic mu
if getattr(self, '_MeMuI', None) is None:
self._MeMuI = self.mesh.getEdgeInnerProduct(self.mu, invMat=True)
return self._MeMuI
@property
def MeMu(self):
#TODO: Assuming isotropic mu
if getattr(self, '_MeMu', None) is None:
self._MeMu = self.mesh.getEdgeInnerProduct(self.mu)
return self._MeMu
@property
def Me(self):
if getattr(self, '_Me', None) is None:
self._Me = self.mesh.getEdgeInnerProduct()
return self._Me
@property
def Mf(self):
if getattr(self, '_Mf', None) is None:
self._Mf = self.mesh.getFaceInnerProduct()
return self._Mf
# ----- Magnetic Permeability ----- #
@property
def MfMui(self):
# TODO: hardcoded to assume diagonal mu
if getattr(self, '_MfMui', None) is None:
self._MfMui = self.mesh.getFaceInnerProduct(1/self.mu)
return self._MfMui
@property
def MeMuI(self):
if getattr(self, '_MeMuI', None) is None:
self._MeMuI = self.mesh.getEdgeInnerProduct(self.mu, invMat=True)
return self._MeMuI
@property
def MeMu(self):
if getattr(self, '_MeMu', None) is None:
self._MeMu = self.mesh.getEdgeInnerProduct(self.mu)
return self._MeMu
# ----- Electrical Conductivity ----- #
#TODO: hardcoded to sigma as the model
@property
def MeSigma(self):
#TODO: hardcoded to sigma as the model
if getattr(self, '_MeSigma', None) is None:
sigma = self.curModel.transform
self._MeSigma = self.mesh.getEdgeInnerProduct(sigma)
return self._MeSigma
@property
def MeSigmaI(self):
#TODO: hardcoded to sigma as the model
if getattr(self, '_MeSigmaI', None) is None:
sigma = self.curModel.transform
self._MeSigmaI = self.mesh.getEdgeInnerProduct(sigma, invMat=True)
return self._MeSigmaI
@property
def dMeSigmaI_dI(self):
# TODO: hardcoded that sigma is diagonal
if getattr(self, '_dMeSigmaI_dI', None) is None:
self._dMeSigmaI_dI = - self.MeSigmaI**2
return self._dMeSigmaI_dI
@property
def MfSigmai(self):
#TODO: hardcoded to sigma as the model
#TODO: hardcoded to sigma diagonal
if getattr(self, '_MfSigmai', None) is None:
sigma = self.curModel.transform
self._MfSigmai = self.mesh.getFaceInnerProduct(1/sigma)
return self._MfSigmai
@property
def dMfSigmai_dsig(self):
return self._dMfSigmai_dsig
deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI','_MfSigmai']
####################################################
# Fields
####################################################
def fields(self, m):
self.curModel = m
F = self.forward(m, self.getRHS)
simpegEM/FDEM/FDEM.py
+1 -2
View File
@@ -232,8 +232,7 @@ class ProblemFDEM_b(BaseFDEMProblem):
C = self.mesh.edgeCurl
sig = self.curModel.transform
dsig_dm = self.curModel.transformDeriv
#TODO: This only works if diagonal (no tensors)...
dMeSigmaI_dI = - self.MeSigmaI**2
dMeSigmaI_dI = self._dMeSigmaI_dI
vec = (C.T*(mui*u))
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(vec)
simpegEM/FDEM/FieldsFDEM.py
+34 -34
View File
@@ -6,7 +6,7 @@ class FieldsFDEM(Problem.Fields):
"""Fancy Field Storage for a FDEM survey."""
knownFields = {}
dtype = complex
class FieldsFDEM_e(FieldsFDEM):
knownFields = {'e':'E'}
@@ -26,7 +26,7 @@ class FieldsFDEM_e(FieldsFDEM):
def _b_sec(self, e, tx): #adjoint=False
return - 1./(1j*omega(tx.freq)) * (self.edgeCurl * e)
def _b_secDeriv(self, e, tx, adjoint=False):
def _b_secDeriv(self, e, tx, v, adjoint=False):
return None
def _b(self, e, tx): #adjoint=False
@@ -34,9 +34,9 @@ class FieldsFDEM_e(FieldsFDEM):
S_m,_ = self.getSource(tx.freq)
return b_sec + 1./(1j*omega(tx.freq)) * S_m
def _bDeriv(self, e, tx, adjoint=False):
S_mDeriv,_ = self.getSourceDeriv(tx.freq, adjoint)
b_secDeriv = self._b_secDeriv(e,tx.freq,adjoint)
def _bDeriv(self, e, tx, v, adjoint=False):
S_mDeriv,_ = self.getSourceDeriv(tx.freq, v, adjoint)
b_secDeriv = self._b_secDeriv(e, tx.freq, v, adjoint)
if S_mDeriv is None & b_secDeriv is None:
return None
elif b_secDeriv is None:
@@ -67,7 +67,7 @@ class FieldsFDEM_b(FieldsFDEM):
def _e_sec(self, b, tx):
return self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * b) )
def _e_secDeriv(self, b, tx, adjoint=False):
def _e_secDeriv(self, b, tx, v, adjoint=False):
return None
def _e(self, b, tx):
@@ -75,9 +75,9 @@ class FieldsFDEM_b(FieldsFDEM):
_, S_e = self.getSource(tx.freq)
return e_sec + S_e
def _eDeriv(self, b, tx, adjoint=False):
_,S_eDeriv = self.getSourceDeriv(tx.freq, adjoint)
e_secDeriv = self._e_secDeriv(b,tx,adjoint)
def _eDeriv(self, b, tx, v, adjoint=False):
_,S_eDeriv = self.getSourceDeriv(tx.freq, v, adjoint)
e_secDeriv = self._e_secDeriv(b, tx, v, adjoint)
if S_eDeriv is None & e_secDeriv is None:
return None
@@ -105,33 +105,33 @@ class FieldsFDEM_j(FieldsFDEM):
self.MfSigmai = self.survey.prob.MfSigmai
self.getSource = self.survey.prob.getSource
self.getSourceDeriv = self.survey.prob.getSourceDeriv
self.curModel = self.prob.curModel
def _h_sec(self, j, tx): #adjoint=False
def _h_sec(self, j, tx): #v, adjoint=False
return - 1./(1j*omega(tx.freq)) * self.MeMuI * (self.edgeCurl.T * (self.MfSigmai * j) )
def _h_secDeriv(self, j, tx, adjoint=False):
# 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.MfSigmai
# 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 _h_secDeriv(self, j, tx, v, adjoint=False):
MeMuI = self.MeMuI
C = self.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.MfSigmai
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 ) ) ) )
def _h(self, j, tx): #adjoint=False
def _h(self, j, tx): #v, adjoint=False
h_sec = self._h_sec(j,tx)
S_m,_ = self.getSource(tx.freq)
return h_sec + 1./(1j*omega(tx.freq)) * self.MeMuI * S_m
def _hDeriv(self, j, tx, adjoint=False):
S_mDeriv,_ = self.getSourceDeriv(tx.freq, adjoint)
h_secDeriv = self._h_secDeriv(j,tx.freq,adjoint)
def _hDeriv(self, j, tx, v, adjoint=False):
S_mDeriv,_ = self.getSourceDeriv(tx.freq, v, adjoint)
h_secDeriv = self._h_secDeriv(j,tx.freq, v, adjoint)
if S_mDeriv is None & h_secDeriv is None:
return None
elif h_secDeriv is None:
@@ -158,20 +158,20 @@ class FieldsFDEM_h(FieldsFDEM):
self.getSource = self.survey.prob.getSource
self.getSourceDeriv = self.survey.prob.getSourceDeriv
def _j_sec(self, h, tx): #adjoint=False
def _j_sec(self, h, tx): # adjoint=False
return self.edgeCurl*h
def _j_secDeriv(self, h, tx, adjoint=False):
def _j_secDeriv(self, h, tx, v, adjoint=False):
return None
def _j(self, h, tx): #adjoint=False
def _j(self, h, tx): # adjoint=False
j_sec = self._j_sec(h,tx)
_,S_e = self.getSource(tx.freq)
return j_sec - S_e
def _jDeriv(self, h, tx, adjoint=False):
_,S_eDeriv = self.getSourceDeriv(tx.freq, adjoint)
j_secDeriv = self._j_secDeriv(j,tx.freq,adjoint)
def _jDeriv(self, h, tx, v, adjoint=False):
_,S_eDeriv = self.getSourceDeriv(tx.freq, v, adjoint)
j_secDeriv = self._j_secDeriv(j,tx.freq, v, adjoint)
if S_eDeriv is None & j_secDeriv is None:
return None
elif j_secDeriv is None: