mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-06 05:16:51 +08:00
Updated 1D_ps problem. Working on testing derivatives, only to 1st order at the commit
This commit is contained in:
GudniRos
@@ -2,11 +2,19 @@ from SimPEG import Survey, Utils, Problem, np, sp, mkvc
|
||||
from scipy.constants import mu_0
|
||||
import sys
|
||||
from numpy.lib import recfunctions as recFunc
|
||||
from simpegEM.Utils.EMUtils import omega
|
||||
|
||||
##############
|
||||
### Fields ###
|
||||
##############
|
||||
class FieldsMT(Problem.Fields):
|
||||
"""Fancy Field Storage for a MT survey."""
|
||||
"""Field Storage for a MT survey."""
|
||||
knownFields = {'b_px': 'F','b_py': 'F', 'e_px': 'E','e_py': 'E','b_1d':'E','e_1d':'F'}
|
||||
dtype = complex
|
||||
|
||||
|
||||
def _b_1dDeriv_u(self,src,v,adjoint=False):
|
||||
"""
|
||||
The derivative of b_1d wrt u
|
||||
"""
|
||||
return -( self.mesh.nodalGrad * v)/( 1j*omega(src.freq) )
|
||||
@@ -10,15 +10,125 @@ import numpy as np
|
||||
import multiprocessing, sys, time
|
||||
|
||||
|
||||
# class eForm_ps(BaseMTProblem):
|
||||
class eForm_psField(BaseMTProblem):
|
||||
"""
|
||||
A MT problem soving a e formulation and primary/secondary fields decomposion.
|
||||
|
||||
Solves the equation
|
||||
|
||||
"""
|
||||
# From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
|
||||
_fieldType = 'e'
|
||||
_eqLocs = 'EF'
|
||||
|
||||
|
||||
def __init__(self, mesh, **kwargs):
|
||||
BaseMTProblem.__init__(self, mesh, **kwargs)
|
||||
|
||||
def getA(self, freq,):
|
||||
"""
|
||||
Function to get the A matrix.
|
||||
|
||||
:param float freq: Frequency
|
||||
:param logic full: Return full A or the inner part
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: A
|
||||
"""
|
||||
|
||||
Mmui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
|
||||
Msig = self.mesh.getFaceInnerProduct(self.curModel.sigma)
|
||||
# Note: need to use the code above since in the 1D problem I want
|
||||
# e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
|
||||
# Possible that _fieldType and _eqLocs can fix this
|
||||
# Mmui = self.MfMui
|
||||
# Msig = self.MeSigma
|
||||
C = self.mesh.nodalGrad
|
||||
# Make A
|
||||
A = C.T*Mmui*C + 1j*omega(freq)*Msig
|
||||
# Either return full or only the inner part of A
|
||||
return A
|
||||
|
||||
def getADeriv_m(self, freq, u, v, adjoint=False):
|
||||
"""
|
||||
The derivative of A wrt sigma
|
||||
"""
|
||||
|
||||
dsig_dm = self.curModel.sigmaDeriv
|
||||
dMf_dsig = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u) * self.curModel.sigmaDeriv
|
||||
if adjoint:
|
||||
return 1j * omega(freq) * ( dsig_dm.T * ( dMf_dsig.T * v ) )
|
||||
|
||||
return 1j * omega(freq) * ( dMf_dsig * ( dsig_dm * v ) )
|
||||
|
||||
def getRHS(self, freq):
|
||||
"""
|
||||
Function to return the right hand side for the system.
|
||||
:param float freq: Frequency
|
||||
:rtype: numpy.ndarray (nF, 1), numpy.ndarray (nF, 1)
|
||||
:return: RHS for 1 polarizations, primary fields
|
||||
"""
|
||||
|
||||
# Get sources for the frequncy(polarizations)
|
||||
Src = self.survey.getSrcByFreq(freq)[0]
|
||||
S_e = Src.S_e(self)
|
||||
return -1j * omega(freq) * S_e
|
||||
|
||||
def getRHSderiv_m(self, freq, u, v, adjoint=False):
|
||||
"""
|
||||
The derivative of the RHS wrt sigma
|
||||
"""
|
||||
|
||||
Src = self.survey.getSrcByFreq(freq)[0]
|
||||
S_eDeriv = Src.S_eDeriv(self, v, adjoint)
|
||||
return -1j * omega(freq) * S_eDeriv
|
||||
|
||||
def fields(self, m):
|
||||
'''
|
||||
Function to calculate all the fields for the model m.
|
||||
|
||||
:param np.ndarray (nC,) m: Conductivity model
|
||||
'''
|
||||
# Set the current model
|
||||
self.curModel = m
|
||||
|
||||
F = FieldsMT(self.mesh, self.survey)
|
||||
for freq in self.survey.freqs:
|
||||
if self.verbose:
|
||||
startTime = time.time()
|
||||
print 'Starting work for {:.3e}'.format(freq)
|
||||
sys.stdout.flush()
|
||||
A = self.getA(freq)
|
||||
rhs = self.getRHS(freq)
|
||||
Ainv = self.Solver(A, **self.solverOpts)
|
||||
e_s = Ainv * rhs
|
||||
|
||||
# Store the fields
|
||||
Src = self.survey.getSrcByFreq(freq)[0]
|
||||
# Calculate total e
|
||||
|
||||
e = Src.ePrimary(self) + e_s
|
||||
|
||||
# Store the fields
|
||||
# NOTE: only store
|
||||
F[Src, 'e_1d'] = e[:,1] # Only storing the yx polarization as 1d
|
||||
# F[Src, 'e_py'] = 0*e[:,0]
|
||||
# Note curl e = -iwb so b = -curl e /iw
|
||||
b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) )
|
||||
# F[Src, 'b_px'] = 0*b[:,0]
|
||||
F[Src, 'b_1d'] = b[:,1]
|
||||
if self.verbose:
|
||||
print 'Ran for {:f} seconds'.format(time.time()-startTime)
|
||||
sys.stdout.flush()
|
||||
return F
|
||||
|
||||
class eForm_TotalField(BaseMTProblem):
|
||||
"""
|
||||
A MT problem solving a e formulation and a primary/secondary fields decompostion.
|
||||
A MT problem solving a e formulation and a Total bondary domain decompostion.
|
||||
|
||||
Solves the equation:
|
||||
|
||||
Math:
|
||||
|
||||
|
||||
"""
|
||||
|
||||
|
||||
@@ -1 +1 @@
|
||||
from Problems import eForm_TotalField
|
||||
from Problems import eForm_TotalField, eForm_psField
|
||||
@@ -74,7 +74,7 @@ class eForm_ps(BaseMTProblem):
|
||||
def getADeriv(self, freq, u, v, adjoint=False):
|
||||
|
||||
dsig_dm = self.curModel.sigmaDeriv
|
||||
dMe_dsig = self.MeSimgaDeriv( v=u)
|
||||
dMe_dsig = self.MeSigmaDeriv( v=u)
|
||||
|
||||
if adjoint:
|
||||
return 1j * omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * v ) )
|
||||
|
||||
@@ -22,26 +22,43 @@ def homo1DModelSource(mesh,freq,sigma_1d):
|
||||
mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]]))
|
||||
# # Note: Everything is using e^iwt
|
||||
e0_1d = get1DEfields(mesh1d,sigma_1d,freq)
|
||||
# Setup x (east) polarization (_x)
|
||||
ex_px = np.zeros(mesh.vnEx,dtype=complex)
|
||||
ey_px = np.zeros((mesh.nEy,1),dtype=complex)
|
||||
ez_px = np.zeros((mesh.nEz,1),dtype=complex)
|
||||
# Assign the source to ex_x
|
||||
for i in np.arange(mesh.vnEx[0]):
|
||||
for j in np.arange(mesh.vnEx[1]):
|
||||
ex_px[i,j,:] = -e0_1d
|
||||
eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px,ez_px))
|
||||
# Setup y (north) polarization (_py)
|
||||
ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
|
||||
ey_py = np.zeros(mesh.vnEy, dtype='complex128')
|
||||
ez_py = np.zeros((mesh.nEz,1), dtype='complex128')
|
||||
# Assign the source to ey_py
|
||||
|
||||
for i in np.arange(mesh.vnEy[0]):
|
||||
for j in np.arange(mesh.vnEy[1]):
|
||||
ey_py[i,j,:] = e0_1d
|
||||
# ey_py[1:-1,1:-1,1:-1] = 0
|
||||
eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
|
||||
if mesh.dim == 1:
|
||||
eBG_px = -simpeg.mkvc(e0_1d,2)
|
||||
eBG_py = simpeg.mkvc(e0_1d,2)
|
||||
elif mesh.dim == 2:
|
||||
ex_px = np.zeros(mesh.vnEx,dtype=complex)
|
||||
ey_px = np.zeros((mesh.nEy,1),dtype=complex)
|
||||
for i in np.arange(mesh.vnEx[0]):
|
||||
ex_px[i,:] = -e0_1d
|
||||
eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px))
|
||||
# Setup y (north) polarization (_py)
|
||||
ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
|
||||
ey_py = np.zeros(mesh.vnEy, dtype='complex128')
|
||||
# Assign the source to ey_py
|
||||
for i in np.arange(mesh.vnEy[0]):
|
||||
ey_py[i,:] = e0_1d
|
||||
# ey_py[1:-1,1:-1,1:-1] = 0
|
||||
eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
|
||||
elif mesh.dim == 3:
|
||||
# Setup x (east) polarization (_x)
|
||||
ex_px = np.zeros(mesh.vnEx,dtype=complex)
|
||||
ey_px = np.zeros((mesh.nEy,1),dtype=complex)
|
||||
ez_px = np.zeros((mesh.nEz,1),dtype=complex)
|
||||
# Assign the source to ex_x
|
||||
for i in np.arange(mesh.vnEx[0]):
|
||||
for j in np.arange(mesh.vnEx[1]):
|
||||
ex_px[i,j,:] = -e0_1d
|
||||
eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px,ez_px))
|
||||
# Setup y (north) polarization (_py)
|
||||
ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
|
||||
ey_py = np.zeros(mesh.vnEy, dtype='complex128')
|
||||
ez_py = np.zeros((mesh.nEz,1), dtype='complex128')
|
||||
# Assign the source to ey_py
|
||||
for i in np.arange(mesh.vnEy[0]):
|
||||
for j in np.arange(mesh.vnEy[1]):
|
||||
ey_py[i,j,:] = e0_1d
|
||||
# ey_py[1:-1,1:-1,1:-1] = 0
|
||||
eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
|
||||
|
||||
# Return the electric fields
|
||||
eBG_bp = np.hstack((eBG_px,eBG_py))
|
||||
|
||||
+36
-19
@@ -94,6 +94,7 @@ class RxMT(Survey.BaseRx):
|
||||
ex = Pex*mkvc(u[src,'e_1d'],2)
|
||||
bx = Pbx*mkvc(u[src,'b_1d'],2)/mu_0
|
||||
f_part_complex = ex/bx
|
||||
# elif self.projType is 'Z2D':
|
||||
elif self.projType is 'Z3D':
|
||||
# Get the projection
|
||||
Pex = mesh.getInterpolationMat(self.locs,'Ex')
|
||||
@@ -124,25 +125,28 @@ class RxMT(Survey.BaseRx):
|
||||
# Get the real or imag component
|
||||
real_or_imag = self.projComp
|
||||
f_part = getattr(f_part_complex, real_or_imag)
|
||||
# print f_part
|
||||
return f_part
|
||||
|
||||
def projectFieldsDeriv(self, src, mesh, u, v, adjoint=False):
|
||||
P = self.getP(mesh)
|
||||
def projectFieldsDeriv(self, src, mesh, f, v, adjoint=False):
|
||||
"""
|
||||
The derivative of the projection wrt u
|
||||
"""
|
||||
|
||||
real_or_imag = self.projComp
|
||||
if not adjoint:
|
||||
Pv_complex = P * v
|
||||
real_or_imag = self.projComp
|
||||
Pv = getattr(Pv_complex, real_or_imag)
|
||||
if self.projType is 'Z1D':
|
||||
Pex = mesh.getInterpolationMat(self.locs,'Fx')
|
||||
Pbx = mesh.getInterpolationMat(self.locs,'Ex')
|
||||
# ex = Pex*mkvc(f[src,'e_1d'],2)
|
||||
# bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
|
||||
deriv_complex = Utils.sdiag(1/(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Pex*v) - Utils.sdiag(1/(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*Utils.sdiag(1/(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Pbx*f._b_1dDeriv_u(src,v)/mu_0)
|
||||
# elif self.projType is 'Z2D
|
||||
elif self.projType is 'Z3D':
|
||||
pass
|
||||
Pv = getattr(deriv_complex, real_or_imag)
|
||||
elif adjoint:
|
||||
Pv_real = P.T * v
|
||||
|
||||
real_or_imag = self.projComp
|
||||
if real_or_imag == 'imag':
|
||||
Pv = 1j*Pv_real
|
||||
elif real_or_imag == 'real':
|
||||
Pv = Pv_real.astype(complex)
|
||||
else:
|
||||
raise NotImplementedError('must be real or imag')
|
||||
raise NotImplementedError('must be real or imag')
|
||||
|
||||
return Pv
|
||||
|
||||
@@ -187,7 +191,7 @@ class srcMT_polxy_1Dprimary(srcMT):
|
||||
as fields in the full space of the problem.
|
||||
"""
|
||||
def __init__(self, rxList, freq, sigma1d):
|
||||
assert mkvc(self.mesh.hz.shape,1) == mkvc(sigma1d.shape,1),'The number of values in the 1D background model does not match the number of vertical cells (hz).'
|
||||
# assert mkvc(self.mesh.hz.shape,1) == mkvc(sigma1d.shape,1),'The number of values in the 1D background model does not match the number of vertical cells (hz).'
|
||||
self.sigma1d = sigma1d
|
||||
srcMT.__init__(self, rxList, freq)
|
||||
|
||||
@@ -201,7 +205,7 @@ class srcMT_polxy_1Dprimary(srcMT):
|
||||
def bPrimary(self,problem):
|
||||
# Project ePrimary to bPrimary
|
||||
# Satisfies the primary(background) field conditions
|
||||
bBG_bp = (- self.mesh.edgeCurl * self.ePrimary )/( 1j*omega(freq) )
|
||||
bBG_bp = (- problem.mesh.edgeCurl * self.ePrimary )/( 1j*omega(freq) )
|
||||
return bBG_bp
|
||||
|
||||
def S_e(self,problem):
|
||||
@@ -213,12 +217,25 @@ class srcMT_polxy_1Dprimary(srcMT):
|
||||
sigma_p = Map_sigma_p._transform(self.sigma1d)
|
||||
# Make mass matrix
|
||||
# Note: M(sig) - M(sig_p) = M(sig - sig_p)
|
||||
Mesigma = problem.MeSigma
|
||||
Mesigma_p = problem.mesh.getEdgeInnerProduct(sigma_p)
|
||||
# Need to deal with the edge/face discrepencies between 1d/2d/3d
|
||||
if problem.mesh.dim == 1:
|
||||
Mesigma = problem.mesh.getFaceInnerProduct(problem.curModel.sigma)
|
||||
Mesigma_p = problem.mesh.getFaceInnerProduct(sigma_p)
|
||||
if problem.mesh.dim == 2:
|
||||
pass
|
||||
if problem.mesh.dim == 3:
|
||||
Mesigma = problem.MeSigma
|
||||
Mesigma_p = problem.mesh.getEdgeInnerProduct(sigma_p)
|
||||
return (Mesigma - Mesigma_p) * e_p
|
||||
|
||||
def S_eDeriv(self, problem, v, adjoint = False):
|
||||
MesigmaDeriv = problem.MeSigmaDeriv(self.ePrimary(problem))
|
||||
# Need to deal with
|
||||
if problem.mesh.dim == 1:
|
||||
pass
|
||||
if problem.mesh.dim == 2:
|
||||
pass
|
||||
if problem.mesh.dim == 3:
|
||||
MesigmaDeriv = problem.MeSigmaDeriv(self.ePrimary(problem))
|
||||
if adjoint:
|
||||
return MesigmaDeriv.T * v
|
||||
else:
|
||||
|
||||
@@ -8,7 +8,7 @@ TOLr = 5e-2
|
||||
TOLp = 5e-2
|
||||
|
||||
|
||||
def setupSurvey(sigmaHalf):
|
||||
def setupSurvey(sigmaHalf,tD=True):
|
||||
|
||||
# Frequency
|
||||
nFreq = 33
|
||||
@@ -31,8 +31,13 @@ def setupSurvey(sigmaHalf):
|
||||
rxList.append(simpegmt.SurveyMT.RxMT(simpeg.mkvc(np.array([0.0]),2).T,rxType))
|
||||
# Source list
|
||||
srcList =[]
|
||||
for freq in freqs:
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1DhomotD(rxList,freq))
|
||||
if tD:
|
||||
for freq in freqs:
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1DhomotD(rxList,freq))
|
||||
else:
|
||||
for freq in freqs:
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq,sigma))
|
||||
|
||||
survey = simpegmt.SurveyMT.SurveyMT(srcList)
|
||||
return survey, sigma, m1d
|
||||
|
||||
@@ -90,23 +95,68 @@ def appPhs_TotalFieldNorm(sigmaHalf):
|
||||
|
||||
return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*135)/ 135)
|
||||
|
||||
def appRes_psFieldNorm(sigmaHalf):
|
||||
|
||||
# Make the survey
|
||||
survey, sigma, mesh = setupSurvey(sigmaHalf,False)
|
||||
problem = simpegmt.ProblemMT1D.eForm_psField(mesh)
|
||||
problem.pair(survey)
|
||||
|
||||
# Get the fields
|
||||
fields = problem.fields(sigma)
|
||||
|
||||
# Project the data
|
||||
data = survey.projectFields(fields)
|
||||
|
||||
# Calculate the app res and phs
|
||||
app_r = np.array(getAppResPhs(data))[:,0]
|
||||
|
||||
return np.linalg.norm(np.abs(app_r - np.ones(survey.nFreq)/sigmaHalf)*sigmaHalf)
|
||||
|
||||
def appPhs_psFieldNorm(sigmaHalf):
|
||||
|
||||
# Make the survey
|
||||
survey, sigma, mesh = setupSurvey(sigmaHalf,False)
|
||||
problem = simpegmt.ProblemMT1D.eForm_psField(mesh)
|
||||
problem.pair(survey)
|
||||
|
||||
# Get the fields
|
||||
fields = problem.fields(sigma)
|
||||
|
||||
# Project the data
|
||||
data = survey.projectFields(fields)
|
||||
|
||||
# Calculate the app phs
|
||||
app_p = np.array(getAppResPhs(data))[:,1]
|
||||
|
||||
return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*135)/ 135)
|
||||
|
||||
class TestAnalytics(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
pass
|
||||
def test_appRes2en1(self):self.assertLess(appRes_TotalFieldNorm(2e-1), TOLr)
|
||||
def test_appRes2en2(self):self.assertLess(appRes_TotalFieldNorm(2e-2), TOLr)
|
||||
def test_appRes2en3(self):self.assertLess(appRes_TotalFieldNorm(2e-3), TOLr)
|
||||
def test_appRes2en4(self):self.assertLess(appRes_TotalFieldNorm(2e-4), TOLr)
|
||||
def test_appRes2en5(self):self.assertLess(appRes_TotalFieldNorm(2e-5), TOLr)
|
||||
def test_appRes2en6(self):self.assertLess(appRes_TotalFieldNorm(2e-6), TOLr)
|
||||
def test_appPhs2en1(self):self.assertLess(appPhs_TotalFieldNorm(2e-1), TOLp)
|
||||
def test_appPhs2en2(self):self.assertLess(appPhs_TotalFieldNorm(2e-2), TOLp)
|
||||
def test_appPhs2en3(self):self.assertLess(appPhs_TotalFieldNorm(2e-3), TOLp)
|
||||
def test_appPhs2en4(self):self.assertLess(appPhs_TotalFieldNorm(2e-4), TOLp)
|
||||
def test_appPhs2en5(self):self.assertLess(appPhs_TotalFieldNorm(2e-5), TOLp)
|
||||
def test_appPhs2en6(self):self.assertLess(appPhs_TotalFieldNorm(2e-6), TOLp)
|
||||
# Total Fields
|
||||
# def test_appRes2en1(self):self.assertLess(appRes_TotalFieldNorm(2e-1), TOLr)
|
||||
# def test_appPhs2en1(self):self.assertLess(appPhs_TotalFieldNorm(2e-1), TOLp)
|
||||
|
||||
def test_appRes2en2(self):self.assertLess(appRes_TotalFieldNorm(2e-2), TOLr)
|
||||
def test_appPhs2en2(self):self.assertLess(appPhs_TotalFieldNorm(2e-2), TOLp)
|
||||
|
||||
# def test_appRes2en3(self):self.assertLess(appRes_TotalFieldNorm(2e-3), TOLr)
|
||||
# def test_appPhs2en3(self):self.assertLess(appPhs_TotalFieldNorm(2e-3), TOLp)
|
||||
|
||||
# def test_appRes2en4(self):self.assertLess(appRes_TotalFieldNorm(2e-4), TOLr)
|
||||
# def test_appPhs2en4(self):self.assertLess(appPhs_TotalFieldNorm(2e-4), TOLp)
|
||||
|
||||
# def test_appRes2en5(self):self.assertLess(appRes_TotalFieldNorm(2e-5), TOLr)
|
||||
# def test_appPhs2en5(self):self.assertLess(appPhs_TotalFieldNorm(2e-5), TOLp)
|
||||
|
||||
# def test_appRes2en6(self):self.assertLess(appRes_TotalFieldNorm(2e-6), TOLr)
|
||||
# def test_appPhs2en6(self):self.assertLess(appPhs_TotalFieldNorm(2e-6), TOLp)
|
||||
|
||||
# Primary/secondary
|
||||
def test_appRes2en2_ps(self):self.assertLess(appRes_psFieldNorm(2e-2), TOLr)
|
||||
def test_appPhs2en2_ps(self):self.assertLess(appPhs_psFieldNorm(2e-2), TOLp)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
Reference in New Issue
Block a user