mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-15 11:26:09 +08:00
Mass matrix derivs, and replaced call of forward with fields
This commit is contained in:
+13
-23
@@ -8,11 +8,6 @@ class EMPropMap(Maps.PropMap):
|
||||
rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap))
|
||||
mui = Maps.Property("Inverse Magnetic Permeability", defaultVal = 1./mu_0, propertyLink=('mu', Maps.ReciprocalMap))
|
||||
|
||||
|
||||
|
||||
|
||||
# Do some error checking: only 1 of sigma, rho can be InvProp similar story with mu and mui
|
||||
# Also ensure that sigma and rho are reciprocals of one another ""
|
||||
|
||||
class BaseEMProblem(Problem.BaseProblem):
|
||||
|
||||
@@ -99,22 +94,25 @@ class BaseEMProblem(Problem.BaseProblem):
|
||||
self._MeSigma = self.mesh.getEdgeInnerProduct(self.curModel.sigma)
|
||||
return self._MeSigma
|
||||
|
||||
@property
|
||||
def MeSigmaDeriv(self, u):
|
||||
"""
|
||||
Deriv of MeSigma wrt sigma
|
||||
"""
|
||||
if getattr(self, 'MeSigmaDeriv', None) is None:
|
||||
self._MeSigmaDeriv = self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u)
|
||||
return self._MeSigmaDeriv
|
||||
"""
|
||||
return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u)
|
||||
|
||||
|
||||
@property
|
||||
def MeSigmaI(self):
|
||||
if getattr(self, '_MeSigmaI', None) is None:
|
||||
self._MeSigmaI = self.mesh.getEdgeInnerProduct(self.curModel.sigma, invMat=True)
|
||||
return self._MeSigmaI
|
||||
|
||||
# def dMeSigmaI_dsigma(self,u)
|
||||
def MeSigmaIDeriv(self, u):
|
||||
"""
|
||||
Deriv of MeSigma wrt sigma
|
||||
"""
|
||||
return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma, invMat=True)(u)
|
||||
|
||||
|
||||
@property
|
||||
def MfRho(self):
|
||||
@@ -122,7 +120,8 @@ class BaseEMProblem(Problem.BaseProblem):
|
||||
self._MfRho = self.mesh.getFaceInnerProduct(self.curModel.rho)
|
||||
return self._MfRho
|
||||
|
||||
# def dMfRho_dsigmai(self,u)
|
||||
def MfRhoDeriv(self,u):
|
||||
return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u)
|
||||
|
||||
@property
|
||||
def MfRhoI(self):
|
||||
@@ -130,14 +129,5 @@ class BaseEMProblem(Problem.BaseProblem):
|
||||
self._MfRhoI = self.mesh.getFaceInnerProduct(self.curModel.rho, invMat=True)
|
||||
return self._MfRhoI
|
||||
|
||||
# def dMfRhoI(self,u)
|
||||
|
||||
|
||||
####################################################
|
||||
# Fields
|
||||
####################################################
|
||||
|
||||
def fields(self, m):
|
||||
self.curModel = m
|
||||
F = self.forward(m)
|
||||
return F
|
||||
def dMfRhoIDeriv(self,u):
|
||||
return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)
|
||||
|
||||
@@ -19,8 +19,8 @@ class BaseFDEMProblem(BaseEMProblem):
|
||||
surveyPair = SurveyFDEM
|
||||
fieldsPair = FieldsFDEM
|
||||
|
||||
def forward(self, m):
|
||||
|
||||
def fields(self, m):
|
||||
self.curModel = m
|
||||
F = self.fieldsPair(self.mesh, self.survey)
|
||||
|
||||
for freq in self.survey.freqs:
|
||||
@@ -175,9 +175,11 @@ class ProblemFDEM_e(BaseFDEMProblem):
|
||||
|
||||
|
||||
def getADeriv(self, freq, u, v, adjoint=False):
|
||||
sig = self.sigma
|
||||
dsig_dm = self.curModel.transformDeriv
|
||||
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(u)
|
||||
# sig = self.sigma
|
||||
# dsig_dm = self.curModel.transformDeriv
|
||||
# dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(u)
|
||||
dsig_dm = self.curModel.sigmaDeriv
|
||||
dMe_dsig = self.MeSigmaDeriv(u)
|
||||
|
||||
if adjoint:
|
||||
return 1j * omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * v ) )
|
||||
|
||||
@@ -1,65 +0,0 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
import simpegEM as EM
|
||||
import sys
|
||||
from scipy.constants import mu_0
|
||||
import copy
|
||||
|
||||
testDerivs = False
|
||||
testCrossCheck = True
|
||||
testAdjoint = False
|
||||
testEB = True
|
||||
testHJ = True
|
||||
|
||||
verbose = False
|
||||
|
||||
TOL = 1e-4
|
||||
FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
|
||||
CONDUCTIVITY = 1e1
|
||||
MU = mu_0
|
||||
freq = [1e-1, 2e-1]
|
||||
addrandoms = True
|
||||
|
||||
def getProblem(fdemType):
|
||||
cs = 5.
|
||||
ncx, ncy, ncz = 6, 6, 6
|
||||
npad = 3
|
||||
hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)]
|
||||
hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
|
||||
hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
|
||||
mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C'])
|
||||
|
||||
mapping = Maps.ExpMap(mesh)
|
||||
|
||||
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.])
|
||||
Src0 = EM.FDEM.SrcFDEM_MagDipole([], loc=np.r_[0.,0.,0.], freq=freq[0])
|
||||
Src1 = EM.FDEM.SrcFDEM_MagDipole([], loc=np.r_[0.,0.,0.], freq=freq[1])
|
||||
|
||||
survey = EM.FDEM.SurveyFDEM([Src0])
|
||||
|
||||
|
||||
if verbose:
|
||||
print ' Fetching %s problem' % (fdemType)
|
||||
|
||||
if fdemType == 'e':
|
||||
prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping)
|
||||
elif fdemType == 'b':
|
||||
prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping)
|
||||
elif fdemType == 'j':
|
||||
prb = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping)
|
||||
elif fdemType == 'h':
|
||||
prb = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping)
|
||||
else:
|
||||
raise NotImplementedError()
|
||||
prb.pair(survey)
|
||||
|
||||
try:
|
||||
from pymatsolver import MumpsSolver
|
||||
prb.Solver = MumpsSolver
|
||||
except ImportError, e:
|
||||
pass
|
||||
|
||||
return prb
|
||||
|
||||
def test_MassMatDeriv()
|
||||
Reference in New Issue
Block a user