simpeg/notebooks/Derivative test MT1D.ipynb

Testing derivaties for 1D MT problem.

Especially the rx.projectFieldsDeriv

In [1]:

import SimPEG as simpeg
import simpegEM as simpegem, simpegMT as simpegmt
from SimPEG.Utils import meshTensor
import numpy as np

In [2]:

simpegmt.FieldsMT.FieldsMT_1D

Out [2]:

simpegMT.FieldsMT.FieldsMT_1D

In [3]:

# Setup the problem
sigmaHalf = 1e-2
# Frequency
nFreq = 33
# freqs = np.logspace(3,-3,nFreq)
freqs = np.array([100])
# Make the mesh
ct = 5
air = meshTensor([(ct,25,1.3)])
# coreT0 = meshTensor([(ct,15,1.2)])
# coreT1 = np.kron(meshTensor([(coreT0[-1],15,1.3)]),np.ones((7,)))
core = np.concatenate( (  np.kron(meshTensor([(ct,15,-1.2)]),np.ones((10,))) , meshTensor([(ct,20)]) ) )
bot = meshTensor([(core[0],10,-1.3)])
x0 = -np.array([np.sum(np.concatenate((core,bot)))])
# Change to use no air
m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core))], x0=x0)
# Make the model
sigma = np.zeros(m1d.nC) + sigmaHalf
sigma[ m1d.gridCC > 0 ] = 1e-8

rxList = []
for rxType in ['z1dr','z1di']:
    rxList.append(simpegmt.SurveyMT.RxMT(simpeg.mkvc(np.array([0.0]),2).T,rxType))
# Source list
srcList =[]
tD = False
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))
# Make the survey
survey = simpegmt.SurveyMT.SurveyMT(srcList)

# Set the problem
problem = simpegmt.ProblemMT1D.eForm_psField(m1d)
problem.pair(survey)

# Get the fields
fields = problem.fields(sigma)

# Project the data
data = survey.projectFields(fields)

In [ ]:

We need calculate this derivative. \begin{align} \underbrace{\frac{\partial P(f(u(m)),m^{fix})}{\partial f}}_{Rx} \end{align}

Use the rule \begin{align} \frac{d}{dx}\left( \frac{a(x)}{b(x)} \right) = \frac{\frac{d }{dx} a(x) b(x) - a(x)\frac{d }{dx} b(x) }{ b(x)^2 } \end{align}

In the case of the 1D MT problem the data is calculated as \begin{align} MT1Ddata = P(f(m)) &= \frac{P_{ex} f_e(src,m)}{P_{bx} f_b(src,m) \frac{1}{\mu_0}} = \frac{P_e u}{P_b f_b(u)} \ \frac{\partial P(f(m))}{\partial u} v &= \frac{P_e}{P_b \frac{1}{mu_0} f_b(u)}v - \frac{P_e u}{\left(P_b \frac{1}{mu_0} f_b(u)\right)^2} P_b \frac{1}{mu_0} \frac{d f_b}{du} v \end{align} where u is the fields that we solve for. \begin{align} \frac{d f_b}{du} = - \frac{1}{i \omega} \nabla \end{align}

In [4]:

# Unused code &= \frac{ P_{ex} P_{bx} \frac{1}{\mu_0} \left( f_b(src,m) - f_e(src,m) \right) } { \left(P_{bx}f_b(src,m) \frac{1}{\mu_0} \right)^2 }

As matrices the formulas above can be written as \begin{align} \left[ \frac{\partial P(f(m))}{\partial u} v \right] = \left[ diag \left[ \frac{1}{\left(P_b \frac{1}{mu_0} f_b(u)\right)} \right] [P_e v] , diag[P_e u] diag \left[ \frac{1}{\left(P_b \frac{1}{mu_0} f_b(u)\right)} \right]^T diag \left[ \frac{1}{\left(P_b \frac{1}{mu_0} f_b(u)\right)} \right] \left[ P_b \frac{1}{mu_0} \frac{d f_b}{du}(v) \right] \right] \end{align}

The adjoint problem is done simliarly \begin{align} \left[ \frac{\partial P(f(m))}{\partial u} v \right]^T = [P_e]^T diag \left[ \frac{1}{\left(P_b \frac{1}{mu_0} f_b(u)\right)} \right]^T v - \left[ P_b \frac{d f_b}{du} \frac{1}{mu_0} \right]^T diag \left[ \frac{1}{\left(P_b \frac{1}{mu_0} f_b(u)\right)} \right] diag \left[ \frac{1}{\left(P_b \frac{1}{mu_0} f_b(u)\right)} \right]^T diag \left[ P_e u \right]^T v \end{align}

In [5]:

# def projectFields(self, src, mesh, u):
#         '''
#         Project the fields and return the
#         '''

#         if self.projType is 'Z1D':
#             Pex = mesh.getInterpolationMat(self.locs,'Fx')
#             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
#         real_or_imag = self.projComp
#         f_part = getattr(f_part_complex, real_or_imag)
#         return f_part

In [6]:

# Initate things for the derivs Test
src = survey.srcList[0]
rx = src.rxList[0]
v = np.random.randn(m1d.nN)
v0 = np.random.randn(m1d.nF+m1d.nE)
u0 = np.random.randn(m1d.nN)+np.random.randn(m1d.nN)*1j
f0 = problem.fieldsPair(m1d,survey)
f0[src,'e_1dSolution'] = u0
# f0[src,'b_1d'] = -1/(1j*simpegem.Utils.EMUtils.omega(src.freq))*m1d.nodalGrad*u0

Cell:

[Cell type raw - unsupported, skipped]

In [7]:

# Run a test
def fun(u):
    f = problem.fieldsPair(m1d,survey)
    f[src,'e_1dSolution'] = u
    return rx.projectFields(src,m1d,f), lambda t: rx.projectFieldsDeriv(src,m1d,f0,t)
simpeg.Tests.checkDerivative(fun,u0,num=5,plotIt=False)

Out [7]:

==================== checkDerivative ====================
iter    h         |ft-f0|   |ft-f0-h*J0*dx|  Order
---------------------------------------------------------
 0   1.00e-01    1.884e-05     1.227e-07      nan
 1   1.00e-02    1.873e-06     1.265e-09      1.987
 2   1.00e-03    1.872e-07     1.269e-11      1.999
 3   1.00e-04    1.872e-08     1.269e-13      2.000
 4   1.00e-05    1.872e-09     1.269e-15      2.000
========================= PASS! =========================
That was easy!

True

In [8]:

rx.projectFieldsDeriv(src,m1d,f0,u0)

Out [8]:

array([ 0.00052762])

In [9]:

rx.projectFields(src,m1d,f0)

Out [9]:

array([[ 0.00124017]])

In [10]:

# Test the Jvec derivative.

In [11]:

# print '%s formulation - %s' % (fdemType, comp)
CONDUCTIVITY = 0.01
m0 = np.log(np.ones(problem.mesh.nC)*CONDUCTIVITY)
# mu = np.log(np.ones(problem.mesh.nC)*MU)

if True:
    m0  = m0 + np.random.randn(problem.mesh.nC)*CONDUCTIVITY*1e-1 
#     mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1

# prb.mu = mu
# survey = prb.survey
def fun(x):
    
    return survey.dpred(x), lambda x: problem.Jvec(m0, x)
simpeg.Tests.checkDerivative(fun, m0, num=4, plotIt=False)

Out [11]:

==================== checkDerivative ====================
iter    h         |ft-f0|   |ft-f0-h*J0*dx|  Order
---------------------------------------------------------
 0   1.00e-01    4.417e-08     4.873e-09      nan
 1   1.00e-02    4.132e-09     4.832e-11      2.004
 2   1.00e-03    4.105e-10     4.828e-13      2.000
 3   1.00e-04    4.103e-11     4.827e-15      2.000
========================= PASS! =========================
You get a gold star!

True

In [12]:

### Adjoint test

We have \begin{align} Jvec =&

In [13]:

# Run a test
TOL = 1e-4
FLR = 1e-20

def projectFieldsAdjointTest(fdemType, comp):
    print 'Adjoint %s formulation - %s' % (fdemType, comp)

    m  = np.log(np.ones(problem.mesh.nC)*0.01)
    if True:
        m  = m + np.random.randn(problem.mesh.nC)*0.01*1e-1 

    u = problem.fields(m)
    v = np.random.randn(1)#+np.random.randn(1)*1j
    # print prb.PropMap.PropModel.nP
    w = np.random.randn(m1d.nN)+np.random.randn(m1d.nN)*1j

    vJw = v.dot(rx.projectFieldsDeriv(src,m1d,f0,w))
    wJtv = w.dot(rx.projectFieldsDeriv(src,m1d,f0,v,adjoint=True)).real
    tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) 
    print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
    return np.abs(vJw - wJtv) < tol
projectFieldsAdjointTest('e','projectFieldsDeriv')

Out [13]:

Adjoint e formulation - projectFieldsDeriv
-2.26989762698e-05 -2.26989762698e-05 0.0 1e-08 True

True

In [14]:

# Run a test
TOL = 1e-4
FLR = 1e-20

def getADeriv_mAdjointTest():
    print 'Adjoint test e formulation - getADeriv_m' 

    m  = np.log(np.ones(problem.mesh.nC)*0.01)
    if True:
        m  = m + np.random.randn(problem.mesh.nC)*0.01*1e-1 

    u = problem.fields(m)
    v = np.random.randn(m1d.nN)#+np.random.randn(1)*1j
    # print prb.PropMap.PropModel.nP
    w = np.random.randn(m1d.nC)#+np.random.randn(m1d.nN)*1j

    vJw = v.dot(problem.getADeriv_m(freq,u0,w))
    wJtv = w.dot(problem.getADeriv_m(freq,u0,v,adjoint=True))
    tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) 
    print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
    return np.abs(vJw - wJtv) < tol
getADeriv_mAdjointTest()

Out [14]:

Adjoint test e formulation - getADeriv_m
(-1977540.36505+2093781.70221j) (-1977540.36505+2093781.70221j) (-1.86264514923e-09+2.79396772385e-09j) 100.0 True

True

In [ ]:

In [15]:

# Run a test
TOL = 1e-4
FLR = 1e-20

def getRHSDeriv_mAdjointTest():
    print 'Adjoint test e formulation - getRHSDeriv_m'

    m  = np.log(np.ones(problem.mesh.nC)*0.01)
    if True:
        m  = m + np.random.randn(problem.mesh.nC)*0.01*1e-1 

    u = problem.fields(m)
    v = np.random.randn(m1d.nN)#+np.random.randn(1)*1j
    # print prb.PropMap.PropModel.nP
    w = np.random.randn(m1d.nC)#+np.random.randn(m1d.nN)*1j

    vJw = v.dot(problem.getRHSDeriv_m(freq,w))
    wJtv = w.dot(problem.getRHSDeriv_m(freq,v,adjoint=True))
    tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) 
    print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
    return np.abs(vJw - wJtv) < tol
getRHSDeriv_mAdjointTest( )

Out [15]:

Adjoint test e formulation - getRHSDeriv_m
(13201.2196403+13827.5790776j) (13201.2196403+13827.5790776j) (-5.45696821064e-12+3.63797880709e-12j) 1.0 True

True

In [16]:

simpeg.mkvc(np.random.randn(survey.nD)+np.random.randn(survey.nD)*1j,2)

print survey.nD

2

In [17]:

TOL = 1e-4
FLR = 1e-20

def JvecAdjointTest():
    print 'Adjoint e formulation - Jvec' 

    m  = np.log(np.ones(problem.mesh.nC)*0.01)
    if True:
        m  = m + np.random.randn(problem.mesh.nC)*0.01*1e-1 

    u = problem.fields(m)

    v = np.random.rand(survey.nD)
    # print prb.PropMap.PropModel.nP
    w = np.random.rand(problem.mesh.nC)

    vJw = v.dot(problem.Jvec(m, w, u))
    wJtv = w.dot(problem.Jtvec(m, v, u))
    tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) 
    print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
    return np.abs(vJw - wJtv) < tol

In [18]:

JvecAdjointTest()

Out [18]:

Adjoint e formulation - Jvec
1.96695386678e-05 1.96695386678e-05 3.38813178902e-21 1e-08 True

True