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

Analytics, b and e formulations, and better solver access

Browse Source
This commit is contained in:
rowanc1
2014-03-18 10:53:20 -07:00
parent d177960abb
commit 49eddb4b83
5 changed files with 157 additions and 18 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegEM/FDEM/FDEM.py
+136 -15
View File
@@ -1,4 +1,4 @@
from SimPEG import Problem, Solver, Utils, np, sp
from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver
from scipy.constants import mu_0
from SurveyFDEM import SurveyFDEM, DataFDEM, FieldsFDEM
from simpegEM.Utils import Sources
@@ -7,7 +7,7 @@ def omega(freq):
"""Change frequency to angular frequency, omega"""
return 2.*np.pi*freq
class ProblemFDEM_e(Problem.BaseProblem):
class BaseProblemFDEM(Problem.BaseProblem):
"""
Frequency-Domain EM problem - E-formulation
@@ -26,8 +26,8 @@ class ProblemFDEM_e(Problem.BaseProblem):
surveyPair = SurveyFDEM
dataPair = DataFDEM
solveOpts = {'factorize':False, 'backend':'scipy'}
Solver = SimpegSolver
solverOpts = {'doDirect':True, 'options':{'factorize':False, 'backend':'scipy'}}
####################################################
# Mass Matrices
@@ -55,9 +55,14 @@ class ProblemFDEM_e(Problem.BaseProblem):
#TODO: assuming constant mu
self._MfMui = self.mesh.getFaceInnerProduct(1/mu_0)
####################################################
# Internal Methods
####################################################
class ProblemFDEM_e(BaseProblemFDEM):
"""
Solving for e!
"""
def __init__(self, model, **kwargs):
BaseProblemFDEM.__init__(self, model, **kwargs)
def getA(self, freq):
"""
@@ -85,9 +90,11 @@ class ProblemFDEM_e(Problem.BaseProblem):
SRCz = src(tx.loc, self.mesh.gridEz, 'z')
rhs[i] = np.concatenate((SRCx, SRCy, SRCz))
#TODO: this is completely wrong. b0 = self.mesh.edgeCurl*rhs but we are doing an e formulation...
j_s = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F')
return -1j*omega(freq)*self.Me*j_s
a = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F')
C = self.mesh.edgeCurl
j_s = C.T*self.MfMui*C*a
#TODO: self.Me* ??
return -1j*omega(freq)*j_s
def fields(self, m, useThisRhs=None):
@@ -99,11 +106,13 @@ class ProblemFDEM_e(Problem.BaseProblem):
for freq in self.survey.freqs:
A = self.getA(freq)
b = self.getRHS(freq)
e = Solver(A, options=self.solveOpts).solve(b)
rhs = self.getRHS(freq)
solver = self.Solver(A, **self.solverOpts)
e = solver.solve(rhs)
print np.linalg.norm(A*Utils.mkvc(e) - Utils.mkvc(rhs)) / np.linalg.norm(Utils.mkvc(rhs))
F[freq, 'e'] = e
#TODO: check if mass matrices needed:
b = -1./(1j*omega(freq))*self.mesh.edgeCurl*e
F[freq, 'b'] = b
@@ -121,7 +130,7 @@ class ProblemFDEM_e(Problem.BaseProblem):
for i, freq in enumerate(self.survey.freqs):
e = u[freq, 'e']
A = self.getA(freq)
solver = Solver(A, options=self.solveOpts)
solver = self.Solver(A, **self.solverOpts)
for txi, tx in enumerate(self.survey.getTransmitters(freq)):
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=e[:,txi])
@@ -150,7 +159,7 @@ class ProblemFDEM_e(Problem.BaseProblem):
for i, freq in enumerate(self.survey.freqs):
e = u[freq, 'e']
AT = self.getA(freq).T
solver = Solver(AT, options=self.solveOpts)
solver = self.Solver(AT, **self.solverOpts)
for txi, tx in enumerate(self.survey.getTransmitters(freq)):
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=e[:,txi])
@@ -161,7 +170,119 @@ class ProblemFDEM_e(Problem.BaseProblem):
return Jtv
class ProblemFDEM_b(BaseProblemFDEM):
"""
Solving for b!
"""
def __init__(self, model, **kwargs):
BaseProblemFDEM.__init__(self, model, **kwargs)
def getA(self, freq):
"""
:param float freq: Frequency
:rtype: scipy.sparse.csr_matrix
:return: A
"""
return self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui + 1j*omega(freq)*self.MfMui
def getRHS(self, freq):
"""
:param float freq: Frequency
:rtype: numpy.ndarray (nE, nTx)
:return: RHS
"""
Txs = self.survey.getTransmitters(freq)
rhs = range(len(Txs))
for i, tx in enumerate(Txs):
if tx.txType == 'VMD':
src = Sources.MagneticDipoleVectorPotential
else:
raise NotImplemented('%s txType is not implemented' % tx.txType)
SRCx = src(tx.loc, self.mesh.gridEx, 'x')
SRCy = src(tx.loc, self.mesh.gridEy, 'y')
SRCz = src(tx.loc, self.mesh.gridEz, 'z')
rhs[i] = np.concatenate((SRCx, SRCy, SRCz))
a = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F')
C = self.mesh.edgeCurl
b_0 = C*a
return -1j*omega(freq)*self.MfMui*b_0
def fields(self, m, useThisRhs=None):
RHS = useThisRhs or self.getRHS
self.makeMassMatrices(m)
F = FieldsFDEM(self.mesh, self.survey)
for freq in self.survey.freqs:
A = self.getA(freq)
rhs = self.getRHS(freq)
solver = self.Solver(A, **self.solverOpts)
#Note that we are solving for b_s
b = solver.solve(rhs)
print np.linalg.norm(A*Utils.mkvc(b) - Utils.mkvc(rhs)) / np.linalg.norm(Utils.mkvc(rhs))
F[freq, 'b'] = b
e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b
F[freq, 'e'] = e
return F
def Jvec(self, m, v, u=None):
if u is None:
u = self.fields(m)
raise NotImplemented('')
# Jv = self.dataPair(self.survey)
# sig = self.model.transform(m)
# dsig_dm = self.model.transformDeriv(m)
# for i, freq in enumerate(self.survey.freqs):
# e = u[freq, 'e']
# A = self.getA(freq)
# solver = self.Solver(A, **self.solverOpts)
# for txi, tx in enumerate(self.survey.getTransmitters(freq)):
# dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=e[:,txi])
# P = tx.projectFieldsDeriv(self.mesh, u)
# b = 1j*omega(freq) * ( dMe_dsig * ( dsig_dm * v ) )
# Ainvb = solver.solve(b)
# Jv[tx] = -P*Ainvb
# return Utils.mkvc(Jv)
def Jtvec(self, m, v, u=None):
if u is None:
u = self.fields(m)
# Ensure v is a data object.
if not isinstance(v, self.dataPair):
v = self.dataPair(self.survey, v)
raise NotImplemented('')
# Jtv = np.zeros(self.model.nP, dtype=complex)
# sig = self.model.transform(m)
# dsig_dm = self.model.transformDeriv(m)
# for i, freq in enumerate(self.survey.freqs):
# e = u[freq, 'e']
# AT = self.getA(freq).T
# solver = self.Solver(AT, **self.solverOpts)
# for txi, tx in enumerate(self.survey.getTransmitters(freq)):
# dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=e[:,txi])
# P = tx.projectFieldsDeriv(self.mesh, u)
# w = solver.solve(P.T * v[tx])
# Jtv += - 1j*omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * w ) )
# return Jtv
simpegEM/FDEM/__init__.py
+1 -1
View File
@@ -1,2 +1,2 @@
from SurveyFDEM import *
from FDEM import ProblemFDEM_e
from FDEM import ProblemFDEM_e, ProblemFDEM_b
simpegEM/Tests/test_FDEM.py
+1 -1
View File
@@ -41,7 +41,7 @@ class FDEM_bDerivTests(unittest.TestCase):
x0 = self.sigma
def fun(x):
return self.survey.dpred(x), lambda x: self.prb.Jvec(x0, x)
passed = Tests.checkDerivative(fun, x0, num=3, plotIt=False)
passed = Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=1e-18)
self.assertTrue(passed)
def test_Jtvec_adjointTest(self):
simpegEM/Utils/Ana/FEM.py
+17
View File
@@ -0,0 +1,17 @@
import numpy as np
from scipy.constants import mu_0, pi
from scipy.special import erf
def hzAnalyticDipoleF(r, freq, sigma):
"""
4.56 in Ward and Hohmann
"""
r = np.abs(r)
k = np.sqrt(-1j*2.*np.pi*freq*mu_0*sigma)
m = 1
front = m / (2. * np.pi * (k**2) * (r**5) )
back = 9 - ( 9 + 9j * k * r - 4 * (k**2) * (r**2) - 1j * (k**3) * (r**3)) * np.exp(-1j*k*r)
hz = front*back
hp =-1/(4*np.pi*r**3)
return hz-hp
simpegEM/Utils/Ana/__init__.py
+2 -1
View File
@@ -1 +1,2 @@
from TEM import hzAnalyticDipoleT
from TEM import hzAnalyticDipoleT
from FEM import hzAnalyticDipoleF