mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-08 10:54:34 +08:00
cleaned up call of fields in Base.py, this will break the TDEM implementation
This commit is contained in:
Lindsey
+1
-3
@@ -6,8 +6,6 @@ class BaseEMProblem(Problem.BaseProblem):
|
||||
def __init__(self, mesh, **kwargs):
|
||||
Problem.BaseProblem.__init__(self, mesh, **kwargs)
|
||||
|
||||
solType = None
|
||||
storeTheseFields = ['e', 'b']
|
||||
|
||||
surveyPair = Survey.BaseSurvey
|
||||
dataPair = Survey.Data
|
||||
@@ -105,5 +103,5 @@ class BaseEMProblem(Problem.BaseProblem):
|
||||
|
||||
def fields(self, m):
|
||||
self.curModel = m
|
||||
F = self.forward(m, self.getRHS, self.calcFields)
|
||||
F = self.forward(m, self.getRHS)
|
||||
return F
|
||||
|
||||
+67
-110
@@ -1,56 +1,11 @@
|
||||
from SimPEG import Survey, Problem, Utils, np, sp, Solver as SimpegSolver
|
||||
from scipy.constants import mu_0
|
||||
from SurveyFDEM import SurveyFDEM, FieldsFDEM
|
||||
from SurveyFDEM import SurveyFDEM
|
||||
from FieldsFDEM import FieldsFDEM, FieldsFDEM_e, FieldsFDEM_b
|
||||
from simpegEM.Base import BaseEMProblem
|
||||
from simpegEM.Utils.EMUtils import omega
|
||||
|
||||
|
||||
|
||||
# class FieldsTDEM_e_from_b(FieldsFDEM):
|
||||
# """Fancy Field Storage for a TDEM survey."""
|
||||
# knownFields = {'b_sec': 'F'}
|
||||
# aliasFields = {
|
||||
# 'b': ['b_sec','F','b_from_bsec'],
|
||||
# 'e': ['b','E','e_from_b']
|
||||
# }
|
||||
|
||||
# def startup(self):
|
||||
# self.MeSigmaI = self.survey.prob.MeSigmaI
|
||||
# self.edgeCurlT = self.survey.prob.mesh.edgeCurl.T
|
||||
# self.MfMui = self.survey.prob.MfMui
|
||||
|
||||
# def e_from_b(self, b, txInd, timeInd):
|
||||
# # TODO: implement non-zero js
|
||||
# return self.MeSigmaI*(self.edgeCurlT*(self.MfMui*b))
|
||||
|
||||
# def e_from_bDeriv(self, b, txInd, timeInd):
|
||||
# # TODO: implement non-zero js
|
||||
# return self.MeSigmaI*(self.edgeCurlT*(self.MfMui*b))
|
||||
|
||||
|
||||
# def calcFields(self, sol, freq, fieldType, adjoint=False):
|
||||
# e = sol
|
||||
# if fieldType == 'e':
|
||||
# return e
|
||||
# elif fieldType == 'b':
|
||||
# if not adjoint:
|
||||
# b = - self.mesh.edgeCurl * e
|
||||
# b = 1./(1j*omega(freq)) * b
|
||||
# else:
|
||||
# b = -(1./(1j*omega(freq))) * ( self.mesh.edgeCurl.T * e )
|
||||
# return b
|
||||
# raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
# def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
|
||||
# e = sol
|
||||
# if fieldType == 'e':
|
||||
# return None
|
||||
# elif fieldType == 'b':
|
||||
# return None
|
||||
# raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
|
||||
|
||||
class BaseFDEMProblem(BaseEMProblem):
|
||||
"""
|
||||
We start by looking at Maxwell's equations in the electric field \\(\\vec{E}\\) and the magnetic flux density \\(\\vec{B}\\):
|
||||
@@ -62,26 +17,23 @@ class BaseFDEMProblem(BaseEMProblem):
|
||||
|
||||
"""
|
||||
surveyPair = SurveyFDEM
|
||||
# fieldsPair = FieldsFDEM
|
||||
fieldsPair = FieldsFDEM
|
||||
|
||||
def forward(self, m, RHS, CalcFields):
|
||||
def forward(self, m, RHS):
|
||||
|
||||
# F = self.fieldsPair(self.mesh, self.survey)
|
||||
|
||||
F = FieldsFDEM(self.mesh, self.survey)
|
||||
F = self.fieldsPair(self.mesh, self.survey)
|
||||
|
||||
for freq in self.survey.freqs:
|
||||
A = self.getA(freq)
|
||||
rhs = RHS(freq)
|
||||
Ainv = self.Solver(A, **self.solverOpts)
|
||||
sol = Ainv * rhs
|
||||
for fieldType in self.storeTheseFields:
|
||||
Txs = self.survey.getTransmitters(freq)
|
||||
F[Txs, fieldType] = CalcFields(sol, freq, fieldType)
|
||||
# for fieldType in self.storeTheseFields:
|
||||
# Txs = self.survey.getTransmitters(freq)
|
||||
# F[Txs, fieldType] = CalcFields(sol, freq, fieldType)
|
||||
|
||||
|
||||
# Txs = self.survey.getTransmitters(freq)
|
||||
# F[Txs, 'e_sec'] = sol
|
||||
Txs = self.survey.getTransmitters(freq)
|
||||
F[Txs, self._fieldType] = sol
|
||||
|
||||
return F
|
||||
|
||||
@@ -169,6 +121,10 @@ class BaseFDEMProblem(BaseEMProblem):
|
||||
return j_m, j_e
|
||||
# return np.concatenate(rhs).reshape((-1, len(Txs)), order='F') #, np.concatenate(j_e).reshape((-1, len(Txs)), order='F')
|
||||
|
||||
def getSourceDeriv(self,freq,adjoint=False):
|
||||
return None, None
|
||||
|
||||
|
||||
##########################################################################################
|
||||
################################ E-B Formulation #########################################
|
||||
##########################################################################################
|
||||
@@ -191,9 +147,9 @@ class ProblemFDEM_e(BaseFDEMProblem):
|
||||
|
||||
|
||||
"""
|
||||
solType = 'e'
|
||||
# _fieldType = 'e'
|
||||
# fieldsPair = FieldsFDEM_e
|
||||
|
||||
_fieldType = 'e'
|
||||
fieldsPair = FieldsFDEM_e
|
||||
|
||||
def __init__(self, model, **kwargs):
|
||||
BaseFDEMProblem.__init__(self, model, **kwargs)
|
||||
@@ -246,33 +202,34 @@ class ProblemFDEM_e(BaseFDEMProblem):
|
||||
return RHS
|
||||
|
||||
|
||||
def calcFields(self, sol, freq, fieldType, adjoint=False):
|
||||
e = sol
|
||||
if fieldType == 'e':
|
||||
return e
|
||||
elif fieldType == 'b':
|
||||
if not adjoint:
|
||||
b = - self.mesh.edgeCurl * e
|
||||
b = 1./(1j*omega(freq)) * b
|
||||
else:
|
||||
b = -(1./(1j*omega(freq))) * ( self.mesh.edgeCurl.T * e )
|
||||
return b
|
||||
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
# def calcFields(self, sol, freq, fieldType, adjoint=False):
|
||||
# e = sol
|
||||
# if fieldType == 'e':
|
||||
# return e
|
||||
# elif fieldType == 'b':
|
||||
# if not adjoint:
|
||||
# b = - self.mesh.edgeCurl * e
|
||||
# b = 1./(1j*omega(freq)) * b
|
||||
# else:
|
||||
# b = -(1./(1j*omega(freq))) * ( self.mesh.edgeCurl.T * e )
|
||||
# return b
|
||||
# raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
|
||||
e = sol
|
||||
if fieldType == 'e':
|
||||
return None
|
||||
elif fieldType == 'b':
|
||||
return None
|
||||
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
# def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
|
||||
# e = sol
|
||||
# if fieldType == 'e':
|
||||
# return None
|
||||
# elif fieldType == 'b':
|
||||
# return None
|
||||
# raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
|
||||
class ProblemFDEM_b(BaseFDEMProblem):
|
||||
"""
|
||||
Solving for b!
|
||||
"""
|
||||
solType = 'b'
|
||||
_fieldType = 'b'
|
||||
fieldsPair = FieldsFDEM_b
|
||||
|
||||
def __init__(self, model, **kwargs):
|
||||
BaseFDEMProblem.__init__(self, model, **kwargs)
|
||||
@@ -345,40 +302,40 @@ class ProblemFDEM_b(BaseFDEMProblem):
|
||||
return RHS
|
||||
|
||||
|
||||
def calcFields(self, sol, freq, fieldType, adjoint=False):
|
||||
b = sol
|
||||
if fieldType == 'e':
|
||||
if not adjoint:
|
||||
e = self.MeSigmaI * ( self.mesh.edgeCurl.T * ( self.MfMui * b ) )
|
||||
else:
|
||||
e = self.MfMui.T * ( self.mesh.edgeCurl * ( self.MeSigmaI.T * b ) )
|
||||
return e
|
||||
elif fieldType == 'b':
|
||||
return b
|
||||
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
# def calcFields(self, sol, freq, fieldType, adjoint=False):
|
||||
# b = sol
|
||||
# if fieldType == 'e':
|
||||
# if not adjoint:
|
||||
# e = self.MeSigmaI * ( self.mesh.edgeCurl.T * ( self.MfMui * b ) )
|
||||
# else:
|
||||
# e = self.MfMui.T * ( self.mesh.edgeCurl * ( self.MeSigmaI.T * b ) )
|
||||
# return e
|
||||
# elif fieldType == 'b':
|
||||
# return b
|
||||
# raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
|
||||
def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
|
||||
b = sol
|
||||
if fieldType == 'e':
|
||||
sig = self.curModel.transform
|
||||
dsig_dm = self.curModel.transformDeriv
|
||||
# def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
|
||||
# b = sol
|
||||
# if fieldType == 'e':
|
||||
# sig = self.curModel.transform
|
||||
# dsig_dm = self.curModel.transformDeriv
|
||||
|
||||
C = self.mesh.edgeCurl
|
||||
mui = self.MfMui
|
||||
# C = self.mesh.edgeCurl
|
||||
# mui = self.MfMui
|
||||
|
||||
#TODO: This only works if diagonal (no tensors)...
|
||||
dMeSigmaI_dI = - self.MeSigmaI**2
|
||||
# #TODO: This only works if diagonal (no tensors)...
|
||||
# dMeSigmaI_dI = - self.MeSigmaI**2
|
||||
|
||||
vec = C.T * ( mui * b )
|
||||
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(vec)
|
||||
if not adjoint:
|
||||
return dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) )
|
||||
else:
|
||||
return dsig_dm.T * ( dMe_dsig.T * ( dMeSigmaI_dI.T * v ) )
|
||||
elif fieldType == 'b':
|
||||
return None
|
||||
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
# vec = C.T * ( mui * b )
|
||||
# dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(vec)
|
||||
# if not adjoint:
|
||||
# return dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) )
|
||||
# else:
|
||||
# return dsig_dm.T * ( dMe_dsig.T * ( dMeSigmaI_dI.T * v ) )
|
||||
# elif fieldType == 'b':
|
||||
# return None
|
||||
# raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
|
||||
##########################################################################################
|
||||
|
||||
@@ -0,0 +1,120 @@
|
||||
from SimPEG import Survey, Problem, Utils, np, sp
|
||||
from simpegEM.Utils.EMUtils import omega
|
||||
|
||||
|
||||
class FieldsFDEM(Problem.Fields):
|
||||
"""Fancy Field Storage for a FDEM survey."""
|
||||
knownFields = {'b': 'F', 'e': 'E', 'j': 'F', 'h': 'E'} # TODO: a, phi
|
||||
dtype = complex
|
||||
|
||||
def calcFields(self,sol,txInd,freqInd,fieldType):
|
||||
if fieldType == 'e':
|
||||
return self._e(sol,txInd,freqInd)
|
||||
elif fieldType == 'e_sec':
|
||||
return self._e_sec(sol,txInd,freqInd)
|
||||
elif fieldType == 'b':
|
||||
return self._b(sol,txInd,freqInd)
|
||||
elif fieldType == 'b_sec':
|
||||
return self._b_sec(sol,txInd,freqInd)
|
||||
else:
|
||||
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
def calcFieldsDeriv(self,sol,txInd,freqInd,fieldType,adjoint=False):
|
||||
if fieldType == 'e':
|
||||
return self._eDeriv(sol,txInd,freqInd,adjoint)
|
||||
elif fieldType == 'e_sec':
|
||||
return self._e_secDeriv(sol,txInd,freqInd,adjoint)
|
||||
elif fieldType == 'b':
|
||||
return self._bDeriv(sol,txInd,freqInd,adjoint,adjoint)
|
||||
elif fieldType == 'b_sec':
|
||||
return self._b_secDeriv(sol,txInd,freqInd,adjoint,adjoint)
|
||||
else:
|
||||
raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
|
||||
|
||||
|
||||
|
||||
class FieldsFDEM_e(FieldsFDEM):
|
||||
knownFields = {'e':'E'}
|
||||
aliasFields = {
|
||||
'b_sec' : ['e','F','_b_sec'],
|
||||
'b' : ['e','F','_b']
|
||||
}
|
||||
|
||||
def __init__(self,mesh,survey,**kwargs):
|
||||
FieldsFDEM.__init__(self,mesh,survey,**kwargs)
|
||||
|
||||
def startup(self):
|
||||
self.edgeCurl = self.survey.prob.mesh.edgeCurl
|
||||
self.freqs = self.survey.freqs
|
||||
self.getSource = self.survey.prob.getSource
|
||||
self.getSourceDeriv = self.survey.prob.getSourceDeriv
|
||||
|
||||
def _e(self, e, txInd, freqInd):
|
||||
return e
|
||||
|
||||
def _eDeriv(self, e, txInd, freqInd, adjoint=False):
|
||||
return None
|
||||
|
||||
def _b_sec(self, e, txInd, freqInd): #adjoint=False
|
||||
iomegainv = 1./(1j*omega(self.freqs[freqInd]))
|
||||
return -iomegainv * (self.edgeCurl * e)
|
||||
|
||||
def _b_secDeriv(self, e, txInd, freqInd, adjoint=False):
|
||||
return None
|
||||
|
||||
def _b(self, e, txInd, freqInd): #adjoint=False
|
||||
freq = self.freqs[freqInd]
|
||||
b_sec = self._bsec(e,txInd,freqInd)
|
||||
j_m,_ = self.getSource(freq)
|
||||
return 1./(1j*omega(freq)) + b_sec
|
||||
|
||||
def _bDreiv(self, e, txInd, freqInd, adjoint=False):
|
||||
freq = self.freqs[freqInd]
|
||||
j_mDeriv,_ = self.getSourceDeriv(freq, adjoint)
|
||||
if j_mDeriv is None:
|
||||
return None
|
||||
else:
|
||||
return 1./(1j*omega(freq)) * j_mDeriv
|
||||
|
||||
|
||||
class FieldsFDEM_b(FieldsFDEM):
|
||||
knownFields = {'b':'F'}
|
||||
aliasFields = {
|
||||
'e_sec' : ['b','E','_e_sec'],
|
||||
'e' : ['b','E','_e']
|
||||
}
|
||||
|
||||
def __init__(self,mesh,survey,**kwargs):
|
||||
FieldsFDEM.__init__(self,mesh,survey,**kwargs)
|
||||
|
||||
def startup(self):
|
||||
self.edgeCurl = self.survey.prob.mesh.edgeCurl
|
||||
self.MeSigmaI = self.survey.prob.MeSigmaI
|
||||
self.MfMui = self.survey.prob.MfMui
|
||||
self.freqs = self.survey.freqs
|
||||
self.getSource = self.survey.prob.getSource
|
||||
self.getSourceDeriv = self.survey.prob.getSourceDeriv
|
||||
|
||||
def _b(self, b, txInd, freqInd):
|
||||
return b
|
||||
|
||||
def _bDeriv(self, b, txInd, freqInd, adjoint=False):
|
||||
return None
|
||||
|
||||
def _e_sec(self, b, txInd, freqInd):
|
||||
return self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * b) )
|
||||
|
||||
def _e_secDeriv(self, b, txInd, freqInd, adjoint=False):
|
||||
return None
|
||||
|
||||
def _e(self, b, txInd, freqInd):
|
||||
e_sec = _e_sec(self, b, txInd, freqInd)
|
||||
_, j_g = self.getSource(self.freqs[freqInd])
|
||||
return e_s - j_g
|
||||
|
||||
def _eDeriv(self, b, txInd, freqInd, adjoint=False):
|
||||
_,j_gDeriv = self.getSourceDeriv(self.freqs[freqInd], adjoint)
|
||||
if j_gDeriv is None:
|
||||
return None
|
||||
else:
|
||||
return -j_gDeriv
|
||||
@@ -2,9 +2,6 @@ from SimPEG import Survey, Problem, Utils, np, sp
|
||||
from simpegEM import Sources
|
||||
from simpegEM.Utils.EMUtils import omega
|
||||
|
||||
def omega(freq):
|
||||
"""Change frequency to angular frequency, omega"""
|
||||
return 2.*np.pi*freq
|
||||
|
||||
class RxFDEM(Survey.BaseRx):
|
||||
|
||||
@@ -105,7 +102,7 @@ class TxFDEM(Survey.BaseTx):
|
||||
|
||||
tx = self
|
||||
freq = tx.freq
|
||||
solType = prob.solType
|
||||
solType = prob._fieldType # Hack, should just ask whether j_m, j_g are defined on edges or faces
|
||||
|
||||
if solType == 'e' or solType == 'b':
|
||||
gridEJx = prob.mesh.gridEx
|
||||
@@ -212,12 +209,6 @@ class SimpleTxFDEM_m(TxFDEM):
|
||||
return self.vec, None
|
||||
|
||||
|
||||
class FieldsFDEM(Problem.Fields):
|
||||
"""Fancy Field Storage for a FDEM survey."""
|
||||
knownFields = {'b': 'F', 'e': 'E', 'j': 'F', 'h': 'E'}
|
||||
dtype = complex
|
||||
|
||||
|
||||
class SurveyFDEM(Survey.BaseSurvey):
|
||||
"""
|
||||
docstring for SurveyFDEM
|
||||
|
||||
@@ -1,2 +1,3 @@
|
||||
from SurveyFDEM import *
|
||||
from FDEM import BaseFDEMProblem, ProblemFDEM_e, ProblemFDEM_b, ProblemFDEM_j, ProblemFDEM_h
|
||||
from FieldsFDEM import *
|
||||
Reference in New Issue
Block a user