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

Working Jvec for 2.5D DC code

Browse Source
This commit is contained in:
seogi_macbook
2016-04-28 11:18:37 -07:00
parent d14cd444ac
commit 0610289fdf
6 changed files with 180 additions and 39 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/EM/Static/DC/FieldsDC_2D.py
+2 -1
View File
@@ -1,5 +1,6 @@
import SimPEG
import Utils, numpy as np, scipy.sparse as sp
from SimPEG.Utils import Identity, Zero
import numpy as np
class Fields_ky(SimPEG.Problem.TimeFields):
SimPEG/EM/Static/DC/ProblemDC_2D.py
+40 -22
View File
@@ -12,7 +12,7 @@ class BaseDCProblem_2D(BaseEMProblem):
surveyPair = Survey_ky
fieldsPair = Fields_ky
nky = 15
ky = np.logspace(-4, 1, nky)
kys = np.logspace(-4, 1, nky)
Ainv = [None for i in range(nky)]
nT = nky # Only for using TimeFields
@@ -26,7 +26,7 @@ class BaseDCProblem_2D(BaseEMProblem):
f = self.fieldsPair(self.mesh, self.survey)
Srcs = self.survey.srcList
for iky in range(self.nky):
ky = self.ky[iky]
ky = self.kys[iky]
A = self.getA(ky)
self.Ainv[iky] = self.Solver(A, **self.solverOpts)
RHS = self.getRHS(ky)
@@ -34,28 +34,44 @@ class BaseDCProblem_2D(BaseEMProblem):
f[Srcs, self._solutionType, iky] = u
return f
# def Jvec(self, m, v, f=None):
def Jvec(self, m, v, f=None):
# if f is None:
# f = self.fields(m)
if f is None:
f = self.fields(m)
# self.curModel = m
self.curModel = m
# Jv = self.dataPair(self.survey) #same size as the data
Jv = self.dataPair(self.survey) #same size as the data
Jv0 = self.dataPair(self.survey)
# A = self.getA()
# Assume y=0.
# This needs some thoughts to implement in general when src is dipole
dky = np.diff(self.kys)
dky = np.r_[dky[0], dky]
y = 0.
# for src in self.survey.srcList:
# u_src = f[src, self._solutionType] # solution vector
# dA_dm_v = self.getADeriv(u_src, v)
# dRHS_dm_v = self.getRHSDeriv(src, v)
# du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v )
# for rx in src.rxList:
# df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
# df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
# Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
# return Utils.mkvc(Jv)
for iky in range(self.nky):
ky = self.kys[iky]
A = self.getA(ky)
for src in self.survey.srcList:
u_src = f[src, self._solutionType, iky] # solution vector
dA_dm_v = self.getADeriv(ky, u_src, v)
dRHS_dm_v = self.getRHSDeriv(ky, src, v)
du_dm_v = self.Ainv[iky] * ( - dA_dm_v + dRHS_dm_v )
for rx in src.rxList:
df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
df_dm_v = df_dmFun(iky, src, du_dm_v, v, adjoint=False)
# Trapezoidal intergration
Jv1_temp = 1./np.pi*rx.evalDeriv(ky, src, self.mesh, f, df_dm_v)
if iky==0:
#First assigment
Jv[src, rx] = Jv1_temp*dky[iky]*np.cos(ky*y)
else:
Jv[src, rx] += Jv1_temp*dky[iky] /2.*np.cos(ky*y)
Jv[src, rx] += Jv0[src, rx]*dky[iky]/2.*np.cos(ky*y)
Jv0[src, rx] = Jv1_temp.copy()
JV[iky,isrc,:] = Jv1_temp.copy()
return Utils.mkvc(Jv)
# def Jtvec(self, m, v, f=None):
# if f is None:
@@ -146,11 +162,13 @@ class Problem2D_CC(BaseDCProblem_2D):
D = self.Div
G = self.Grad
vol = self.mesh.vol
MfRhoIDeriv = self.MfRhoIDeriv
rho = self.curModel.rho
if adjoint:
return(MfRhoIDeriv( G * u ).T) * ( D.T * v) + Utils.sdiag(ky**2*mesh.vol)*v
return D * ((MfRhoIDeriv( G * u )) * v) + Utils.sdiag(ky**2*mesh.vol)*v
return(MfRhoIDeriv( G * u ).T) * ( D.T * v) + ky**2*Utils.sdiag(u.flatten()*vol*(-1./rho**2))*v
return D * ((MfRhoIDeriv( G * u )) * v) + ky**2*Utils.sdiag(u.flatten()*vol*(-1./rho**2))*v
def getRHS(self, ky):
"""
SimPEG/EM/Static/DC/RxDC.py
+9 -9
View File
@@ -102,10 +102,10 @@ class Dipole_ky(BaseRx):
self._Ps[mesh] = P
return P
def eval(self, ky, src, mesh, f):
def eval(self, kys, src, mesh, f):
P = self.getP(mesh, self.projGLoc(f))
Pf = P*f[src, self.projField,:]
return self.IntTrapezoidal(ky, Pf, y=0.)
return self.IntTrapezoidal(kys, Pf, y=0.)
def evalDeriv(self, ky, src, mesh, f, v, adjoint=False):
P = self.getP(mesh, self.projGLoc(f))
@@ -114,16 +114,16 @@ class Dipole_ky(BaseRx):
elif adjoint:
return P.T*v
def IntTrapezoidal(self, ky, Pf, y=0.):
def IntTrapezoidal(self, kys, Pf, y=0.):
phi = np.zeros(Pf.shape[0])
nky = ky.size
dky = np.diff(ky)
nky = kys.size
dky = np.diff(kys)
dky = np.r_[dky[0], dky]
phi0 = Pf[:,0]
phi0 = 1./np.pi*Pf[:,0]
for iky in range(nky):
phi1 = 2./np.pi*Pf[:,iky]/2.
phi += phi1*dky[iky]/2.*np.cos(ky[iky]*y)
phi += phi0*dky[iky]/2.*np.cos(ky[iky]*y)
phi1 = 1./np.pi*Pf[:,iky]
phi += phi1*dky[iky]/2.*np.cos(kys[iky]*y)
phi += phi0*dky[iky]/2.*np.cos(kys[iky]*y)
phi0 = phi1.copy()
return phi
SimPEG/EM/Static/DC/SrcDC.py
-5
View File
@@ -36,11 +36,6 @@ class Dipole(BaseSrc):
q = self.current * mkvc(qa+qb)
return q
# def bc_contribution
# How to treat boundary conditions here
class Pole(BaseSrc):
def __init__(self, rxList, loc, **kwargs):
SimPEG/EM/Static/DC/SurveyDC.py
+2 -2
View File
@@ -29,8 +29,8 @@ class Survey_ky(BaseEMSurvey):
:return: data
"""
data = SimPEG.Survey.Data(self)
ky = self.prob.ky
kys = self.prob.kys
for src in self.srcList:
for rx in src.rxList:
data[src, rx] = rx.eval(ky, src, self.mesh, f)
data[src, rx] = rx.eval(kys, src, self.mesh, f)
return data
tests/em/static/test_DC_2D_jvecjtvecadj.py
+127
View File
@@ -0,0 +1,127 @@
import unittest
from SimPEG import *
import SimPEG.EM.Static.DC as DC
class DCProblem_2DTestsCC(unittest.TestCase):
def setUp(self):
cs = 12.5
hx = [(cs,7, -1.3),(cs,61),(cs,7, 1.3)]
hy = [(cs,7, -1.3),(cs,20)]
mesh = Mesh.TensorMesh([hx, hy],x0="CN")
x = np.linspace(-135, 250., 20)
M = Utils.ndgrid(x-12.5, np.r_[0.])
N = Utils.ndgrid(x+12.5, np.r_[0.])
A0loc = np.r_[-150, 0.]
A1loc = np.r_[-130, 0.]
rxloc = [np.c_[M, np.zeros(20)], np.c_[N, np.zeros(20)]]
rx = DC.Rx.Dipole_ky(M, N)
src0 = DC.Src.Pole([rx], A0loc)
src1 = DC.Src.Pole([rx], A1loc)
survey = DC.Survey_ky([src0, src1])
problem = DC.Problem2D_CC(mesh, mapping=[('rho', Maps.IdentityMap(mesh))])
problem.pair(survey)
mSynth = np.ones(mesh.nC)
survey.makeSyntheticData(mSynth)
# Now set up the problem to do some minimization
dmis = DataMisfit.l2_DataMisfit(survey)
reg = Regularization.Tikhonov(mesh)
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
inv = Inversion.BaseInversion(invProb)
self.inv = inv
self.reg = reg
self.p = problem
self.mesh = mesh
self.m0 = mSynth
self.survey = survey
self.dmis = dmis
def test_misfit(self):
derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
self.assertTrue(passed)
# def test_adjoint(self):
# # Adjoint Test
# u = np.random.rand(self.mesh.nC*self.survey.nSrc)
# v = np.random.rand(self.mesh.nC)
# w = np.random.rand(self.survey.dobs.shape[0])
# wtJv = w.dot(self.p.Jvec(self.m0, v))
# vtJtw = v.dot(self.p.Jtvec(self.m0, w))
# passed = np.abs(wtJv - vtJtw) < 1e-10
# print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
# self.assertTrue(passed)
# def test_dataObj(self):
# derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
# passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
# self.assertTrue(passed)
# class DCProblemTestsN(unittest.TestCase):
# def setUp(self):
# aSpacing=2.5
# nElecs=10
# surveySize = nElecs*aSpacing - aSpacing
# cs = surveySize/nElecs/4
# mesh = Mesh.TensorMesh([
# [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)],
# [(cs,3, -1.3),(cs,3,1.3)],
# # [(cs,5, -1.3),(cs,10)]
# ],'CN')
# srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True)
# survey = DC.Survey(srcList)
# problem = DC.Problem3D_N(mesh, mapping=[('rho', Maps.IdentityMap(mesh))])
# problem.pair(survey)
# mSynth = np.ones(mesh.nC)
# survey.makeSyntheticData(mSynth)
# # Now set up the problem to do some minimization
# dmis = DataMisfit.l2_DataMisfit(survey)
# reg = Regularization.Tikhonov(mesh)
# opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
# invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
# inv = Inversion.BaseInversion(invProb)
# self.inv = inv
# self.reg = reg
# self.p = problem
# self.mesh = mesh
# self.m0 = mSynth
# self.survey = survey
# self.dmis = dmis
# def test_misfit(self):
# derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
# passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
# self.assertTrue(passed)
# def test_adjoint(self):
# # Adjoint Test
# u = np.random.rand(self.mesh.nC*self.survey.nSrc)
# v = np.random.rand(self.mesh.nC)
# w = np.random.rand(self.survey.dobs.shape[0])
# wtJv = w.dot(self.p.Jvec(self.m0, v))
# vtJtw = v.dot(self.p.Jtvec(self.m0, w))
# passed = np.abs(wtJv - vtJtw) < 1e-8
# print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
# self.assertTrue(passed)
# def test_dataObj(self):
# derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
# passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
# self.assertTrue(passed)
if __name__ == '__main__':
unittest.main()