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

Fixed 1D test and current code to work, where the src in the 1D problem is partly implemented

Browse Source
This commit is contained in:
GudniRos
2015-06-11 16:26:11 -07:00
parent 21d788edf4
commit 422911a95f
7 changed files with 174 additions and 101 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegMT/ProblemMT1D/Problems.py
+18 -10
View File
@@ -10,8 +10,11 @@ import numpy as np
import multiprocessing, sys, time
# class eForm_ps(BaseMTProblem):
class eForm_TotalField(BaseMTProblem):
"""
"""
A MT problem solving a e formulation and a primary/secondary fields decompostion.
Solves the equation:
@@ -21,8 +24,8 @@ class eForm_TotalField(BaseMTProblem):
# From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
_fieldType = 'e'
_eqLocs = 'FE'
_eqLocs = 'EF'
def __init__(self, mesh, **kwargs):
BaseMTProblem.__init__(self, mesh, **kwargs)
@@ -36,8 +39,14 @@ class eForm_TotalField(BaseMTProblem):
:rtype: scipy.sparse.csr_matrix
:return: A
"""
Mmui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
Msig = self.mesh.getFaceInnerProduct(self.curModel)
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
@@ -66,12 +75,12 @@ class eForm_TotalField(BaseMTProblem):
"""
# Get sources for the frequency
# NOTE: Need to use the source information, doesn't really apply in 1D
src = self.survey.getSources(freq)
src = self.survey.getSrcByFreq(freq)
# Get the full A
A = self.getA(freq,full=True)
# Define the outer part of the solution matrix
Aio = A[1:-1,[0,-1]]
Ed, Eu, Hd, Hu = getEHfields(self.mesh,self.curModel,freq,self.mesh.vectorNx)
Ed, Eu, Hd, Hu = getEHfields(self.mesh,self.curModel.sigma,freq,self.mesh.vectorNx)
Etot = (Ed + Eu)
sourceAmp = 1.0
Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
@@ -104,12 +113,12 @@ class eForm_TotalField(BaseMTProblem):
A = self.getA(freq)
rhs, e_o = self.getRHS(freq)
Ainv = self.Solver(A, **self.solverOpts)
e_i = Ainv * rhs
e_i = Ainv * rhs
e = mkvc(np.r_[e_o[0], e_i, e_o[1]],2)
# Store the fields
Src = self.survey.getSources(freq)
Src = self.survey.getSrcByFreq(freq)
# Store the fields
# NOTE: only store
# NOTE: only store
F[Src, 'e_1d'] = e
# F[Src, 'e_py'] = 0*e[:,0]
# Note curl e = -iwb so b = -curl e /iw
@@ -120,4 +129,3 @@ class eForm_TotalField(BaseMTProblem):
print 'Ran for {:f} seconds'.format(time.time()-startTime)
sys.stdout.flush()
return F
simpegMT/ProblemMT3D/Problems.py
+67 -62
View File
@@ -4,49 +4,55 @@ from scipy.constants import mu_0
from simpegMT.BaseMT import BaseMTProblem
from simpegMT.SurveyMT import SurveyMT
from simpegMT.FieldsMT import FieldsMT
from simpegMT.DataMT import DataMT
from simpegMT.DataMT import DataMT
import multiprocessing, sys, time
class eForm_ps(BaseMTProblem):
"""
"""
A MT problem solving a e formulation and a primary/secondary fields decompostion.
Solves the equation
Solves the equation:
"""
# From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
_fieldType = 'e'
_eqLocs = 'FE'
# Set new properties
# Need to add the src ....
# Set new properties
# Background model
@property
def backModel(self):
"""
Sets the model, and removes dependent mass matrices.
"""
return getattr(self, '_backModel', None)
# Shouldn't need the commented block.
# @property
# def backModel(self):
# """
# Sets the model, and removes dependent mass matrices.
# """
# return getattr(self, '_backModel', None)
@backModel.setter
def backModel(self, value):
if value is self.backModel:
return # it is the same!
self._backModel = Models.Model(value, self.mapping)
for prop in self.deleteTheseOnModelUpdate:
if hasattr(self, prop):
delattr(self, prop)
# @backModel.setter
# def backModel(self, value):
# if value is self.backModel:
# return # it is the same!
# self._backModel = Models.Model(value, self.mapping)
# for prop in self.deleteTheseOnModelUpdate:
# if hasattr(self, prop):
# delattr(self, prop)
@property
def MeDeltaSigma(self):
#TODO: hardcoded to sigma as the model
if getattr(self, '_MeDeltaSigma', None) is None:
sigma = self.curModel
sigmaBG = self.backModel
self._MeDeltaSigma = self.mesh.getEdgeInnerProduct(sigma - sigmaBG)
return self._MeDeltaSigma
# @property
# def MeDeltaSigma(self):
# #TODO: hardcoded to sigma as the model
# if getattr(self, '_MeDeltaSigma', None) is None:
# sigma = self.curModel
# sigmaBG = self.backModel
# self._MeDeltaSigma = self.mesh.getEdgeInnerProduct(sigma - sigmaBG)
# return self._MeDeltaSigma
def __init__(self, mesh, **kwargs):
BaseMTProblem.__init__(self, mesh, **kwargs)
@@ -59,56 +65,52 @@ class eForm_ps(BaseMTProblem):
:rtype: scipy.sparse.csr_matrix
:return: A
"""
mui = self.MfMui
sig = self.MeSigma
Mmui = self.MfMui
Msig = self.MeSigma
C = self.mesh.edgeCurl
return C.T*mui*C + 1j*omega(freq)*sig
return C.T*Mmui*C + 1j*omega(freq)*Msig
def getADeriv(self, freq, u, v, adjoint=False):
sig = self.curTModel
dsig_dm = self.curTModelDeriv
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=u)
dsig_dm = self.curModel.sigmaDeriv
dMe_dsig = self.MeSimgaDeriv( v=u)
if adjoint:
return 1j * omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * v ) )
return 1j * omega(freq) * ( dMe_dsig * ( dsig_dm * v ) )
def getRHS(self, freq, backSigma):
def getRHS(self, freq):
"""
Function to return the right hand side for the system.
:param float freq: Frequency
:param numpy.ndarray (nC,) backSigma: Background conductivity model
:rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
:return: RHS for both polarizations, primary fields
"""
# Get sources for the frequency
src = self.survey.getSources(freq)
# Make sure that there is 2 polarizations.
# assert len()
# Get the background electric fields
from simpegMT.Sources import homo1DModelSource
eBG_bp = homo1DModelSource(self.mesh,freq,backSigma)
deltM = self.MeDeltaSigma
Abg = -1j*omega(freq)*deltM
return Abg*eBG_bp, eBG_bp
# 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(self, freq, backSigma, u, v, adjoint=False):
raise NotImplementedError('getRHSDeriv not implemented yet')
return None
def getRHSderiv(self, freq, u, v, adjoint=False):
"""
The derivative of the RHS with respect to sigma
"""
def fields(self, m, m_back):
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
:param np.ndarray (nC,) m_back: Background conductivity model
'''
# Set the current model
self.curModel = m
self.backModel = m_back
# RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields
F = FieldsMT(self.mesh, self.survey)
for freq in self.survey.freqs:
@@ -117,12 +119,15 @@ class eForm_ps(BaseMTProblem):
print 'Starting work for {:.3e}'.format(freq)
sys.stdout.flush()
A = self.getA(freq)
rhs, e_p = self.getRHS(freq,m_back)
rhs = self.getRHS(freq)
Ainv = self.Solver(A, **self.solverOpts)
e_s = Ainv * rhs
e = e_p + e_s
e_s = Ainv * rhs
# Store the fields
Src = self.survey.getSources(freq)
Src = self.survey.getSrcByFreq(freq)[0]
# Calculate total e
e = Src.ePrimary(self) + e_s
# Store the fieldss
F[Src, 'e_px'] = e[:,0]
F[Src, 'e_py'] = e[:,1]
@@ -134,9 +139,9 @@ class eForm_ps(BaseMTProblem):
print 'Ran for {:f} seconds'.format(time.time()-startTime)
sys.stdout.flush()
return F
class eForm_Tp(BaseMTProblem):
"""
"""
A MT problem solving a e formulation and a total/primary fields decompostion.
Solves the equation
@@ -146,7 +151,7 @@ class eForm_Tp(BaseMTProblem):
_eqLocs = 'FE'
fieldsPair = FieldsMT
# Set new properties
# Set new properties
# Background model
@property
def backModel(self):
@@ -210,7 +215,7 @@ class eForm_Tp(BaseMTProblem):
"""
# Get sources for the frequency
src = self.survey.getSources(freq)
# Make sure that there is 2 polarizations.
# Make sure that there is 2 polarizations.
# assert len()
# Get the background electric fields
from simpegMT.Sources import homo1DModelSource
@@ -246,7 +251,7 @@ class eForm_Tp(BaseMTProblem):
A = self.getA(freq)
rhs, e_p = self.getRHS(freq,m_back)
Ainv = self.Solver(A, **self.solverOpts)
e_s = Ainv * rhs
e_s = Ainv * rhs
e = e_s
# Store the fields
Src = self.survey.getSources(freq)
@@ -261,4 +266,4 @@ class eForm_Tp(BaseMTProblem):
print 'Ran for {:f} seconds'.format(time.time()-startTime)
sys.stdout.flush()
return F
simpegMT/Sources/backgroundModelSources.py
+11 -7
View File
@@ -1,23 +1,27 @@
import SimPEG as simpeg, numpy as np
def homo1DModelSource(mesh,freq,m_back):
def homo1DModelSource(mesh,freq,sigma_1d):
'''
Function that calculates and return background fields
:param Simpeg mesh object mesh: Holds information on the discretization
:param float freq: The frequency to solve at
:param np.array m_back: Background model of conductivity to base the calculations on.
:param np.array sigma_1d: Background model of conductivity to base the calculations on, 1d model.
:rtype: numpy.ndarray (mesh.nE,2)
:return: eBG_bp, E fields for the background model at both polarizations.
'''
# import
from simpegMT.Utils import get1DEfields
# Get a 1d solution for a halfspace background
mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]]))
# Note: Everything is using e^iwt
e0_1d = get1DEfields(mesh1d,mesh.r(m_back,'CC','CC','M')[0,0,:],freq)
if mesh.dim == 1:
mesh1d = mesh
elif mesh.dim == 2:
mesh1d = simpeg.Mesh.TensorMesh([mesh.hy],np.array([mesh.x0[1]]))
elif mesh.dim == 3:
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)
@@ -32,7 +36,7 @@ def homo1DModelSource(mesh,freq,m_back):
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
simpegMT/SurveyMT.py
+70 -15
View File
@@ -1,8 +1,11 @@
from SimPEG import Survey, Utils, Problem, np, sp, mkvc
from SimPEG import Survey, Utils, Problem, Maps, np, sp, mkvc
from simpegEM.FDEM.SurveyFDEM import SrcFDEM
from simpegEM.Utils.EMUtils import omega
from scipy.constants import mu_0
import sys
from numpy.lib import recfunctions as recFunc
from DataMT import DataMT
from simpegMT.Sources import homo1DModelSource
#################
### Receivers ###
#################
@@ -45,7 +48,7 @@ class RxMT(Survey.BaseRx):
"""
Field Type projection (e.g. e b ...)
:param str fracPos: Position of the field in the data ratio
"""
if 'numerator' in fracPos:
return self.knownRxTypes[self.rxType][0][0]
@@ -59,7 +62,7 @@ class RxMT(Survey.BaseRx):
"""
Grid Location projection (e.g. Ex Fy ...)
:param str fracPos: Position of the field in the data ratio
"""
if 'numerator' in fracPos:
return self.knownRxTypes[self.rxType][0][1]
@@ -74,7 +77,7 @@ class RxMT(Survey.BaseRx):
"""
return self.knownRxTypes[self.rxType][0]
@property
def projComp(self):
"""Component projection (real/imag)"""
@@ -82,12 +85,12 @@ class RxMT(Survey.BaseRx):
def projectFields(self, src, mesh, u):
'''
Project the fields and return the
Project the fields and return the
'''
if self.projType is 'Z1D':
Pex = mesh.getInterpolationMat(self.locs,'Fx')
Pbx = mesh.getInterpolationMat(self.locs,'Ex')
Pbx = mesh.getInterpolationMat(self.locs,'Ex')
ex = Pex*mkvc(u[src,'e_1d'],2)
bx = Pbx*mkvc(u[src,'b_1d'],2)/mu_0
f_part_complex = ex/bx
@@ -144,30 +147,82 @@ class RxMT(Survey.BaseRx):
return Pv
# Note: Might need to add tests to make sure that both polarization have the same rxList.
# Note: Might need to add tests to make sure that both polarization have the same rxList.
###############
### Sources ###
###############
class srcMT(Survey.BaseSrc):
'''
Sources for the MT problem.
Sources for the MT problem.
Use the SimPEG BaseSrc, since the source fields share properties with the transmitters.
:param float freq: The frequency of the source
:param list rxList: A list of receivers associated with the source
:param str srcPol: The polarization of the source
'''
freq = None #: Frequency (float)
rxPair = RxMT
knownSrcTypes = ['pol_xy','pol_x','pol_y'] # ORThogonal POLarization
def __init__(self, freq, rxList, srcPol = 'pol_xy'): # remove rxType? hardcode to one thing. always polarizations
def __init__(self, rxList, freq):
self.freq = float(freq)
Survey.BaseSrc.__init__(self, None, srcPol, rxList)
Survey.BaseSrc.__init__(self, rxList)
# 1D sources
class srcMT_polxy_1DhomotD(srcMT):
"""
MT source for both polarizations (x and y) for the total Domain. It calculates fields calculated based on conditions on the boundary of the domain.
"""
def __init__(self, rxList, freq):
srcMT.__init__(self, rxList, freq)
# TODO: need to add the primary fields calc and source terms into the problem.
# Need to implement such that it works for all dims.
class srcMT_polxy_1Dprimary(srcMT):
"""
MT source for both polarizations (x and y) given a 1D primary models. It assigns fields calculated from the 1D model
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).'
self.sigma1d = sigma1d
srcMT.__init__(self, rxList, freq)
def ePrimary(self,problem):
# Get primary fields for both polarizations
eBG_bp = homo1DModelSource(problem.mesh,self.freq,self.sigma1d)
return eBG_bp
def bPrimary(self,problem):
# Project ePrimary to bPrimary
# Satisfies the primary(background) field conditions
bBG_bp = (- self.mesh.edgeCurl * self.ePrimary )/( 1j*omega(freq) )
return bBG_bp
def S_e(self,problem):
"""
Get the electrical field source
"""
e_p = self.ePrimary(problem)
Map_sigma_p = Maps.Vertical1DMap(problem.mesh)
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)
return (Mesigma - Mesigma_p) * e_p
def S_eDeriv(self, problem, v, adjoint = False):
MesigmaDeriv = problem.MeSigmaDeriv(self.ePrimary(problem))
if adjoint:
return MesigmaDeriv.T * v
else:
return MesigmaDeriv * v
##############
@@ -208,7 +263,7 @@ class SurveyMT(Survey.BaseSurvey):
return len(self._freqDict)
# TODO: Rename to getSources
def getSources(self, freq):
def getSrcByFreq(self, freq):
"""Returns the sources associated with a specific frequency."""
assert freq in self._freqDict, "The requested frequency is not in this survey."
return self._freqDict[freq]
simpegMT/Tests/test_Problem1D_againstAnalyticHalfspace.py
+1 -1
View File
@@ -32,7 +32,7 @@ def setupSurvey(sigmaHalf):
# Source list
srcList =[]
for freq in freqs:
srcList.append(simpegmt.SurveyMT.srcMT(freq,rxList))
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1DhomotD(rxList,freq))
survey = simpegmt.SurveyMT.SurveyMT(srcList)
return survey, sigma, m1d
simpegMT/Tests/test_Problem3D_againstAnalytic.py
+5 -4
View File
@@ -71,8 +71,9 @@ def runSimpegMTfwd_eForm_ps(inputsProblem):
rxList.append(simpegmt.SurveyMT.RxMT(rx_loc,rxType))
# Source list
srcList =[]
sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:]
for freq in freqs:
srcList.append(simpegmt.SurveyMT.srcMT(freq,rxList))
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq,sigma1d))
# Survey MT
survey = simpegmt.SurveyMT.SurveyMT(srcList)
@@ -83,7 +84,7 @@ def runSimpegMTfwd_eForm_ps(inputsProblem):
problem.Solver = MumpsSolver
problem.pair(survey)
fields = problem.fields(sig,sigBG)
fields = problem.fields(sig)
mtData = survey.projectFields(fields)
return (survey, problem, fields, mtData)
@@ -93,7 +94,7 @@ def getAppResPhs(MTdata):
# Make impedance
def appResPhs(freq,z):
app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
app_phs = np.arctan2(-z.imag,z.real)*(180/np.pi)
app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
return app_res, app_phs
recData = MTdata.toRecArray('Complex')
return appResPhs(recData['freq'],recData['zxy']), appResPhs(recData['freq'],recData['zyx'])
@@ -107,7 +108,7 @@ def appResPhsHalfspace_eFrom_ps_Norm(sigmaHalf,appR=True):
if appR:
return np.linalg.norm(np.abs(app_rpxy[0,:] - np.ones(survey.nFreq)/sigmaHalf) * sigmaHalf)
else:
return np.linalg.norm(np.abs(app_rpxy[1,:] - np.ones(survey.nFreq)/135) * 135)
return np.linalg.norm(np.abs(app_rpxy[1,:] + np.ones(survey.nFreq)*135) / 135)
class TestAnalytics(unittest.TestCase):
simpegMT/Utils/MT1Dsolutions.py
+2 -2
View File
@@ -13,7 +13,7 @@ def get1DEfields(m1d,sigma,freq,sourceAmp=1.0):
# Conductivity
Msig = m1d.getFaceInnerProduct(sigma)
# Set up the solution matrix
A = G.T*Mmu*G - 1j*2.*np.pi*freq*Msig
A = G.T*Mmu*G + 1j*2.*np.pi*freq*Msig
# Define the inner part of the solution matrix
Aii = A[1:-1,1:-1]
# Define the outer part of the solution matrix
@@ -27,7 +27,7 @@ def get1DEfields(m1d,sigma,freq,sourceAmp=1.0):
## Note: The analytic solution is derived with e^iwt
bc = np.r_[Etot[0],Etot[-1]]
# The right hand side
rhs = -Aio*bc
rhs = Aio*bc
# Solve the system
Aii_inv = simpeg.Solver(Aii)
eii = Aii_inv*rhs