mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-08 01:23:08 +08:00
minor changes to tdem writeup and code (still buggy!)
This commit is contained in:
rowanc1
+38
-5
@@ -2,16 +2,49 @@
|
||||
|
||||
|
||||
.. math::
|
||||
|
||||
\newcommand{\dcurl}{{\mathbf C}}
|
||||
\renewcommand {\b} { {\vec b} }
|
||||
\newcommand {\e} { {\vec e} }
|
||||
\renewcommand {\j} { {\vec j} }
|
||||
|
||||
\renewcommand{\div}{\nabla\cdot\,}
|
||||
\newcommand{\grad}{\vec \nabla}
|
||||
\newcommand{\curl}{{\vec \nabla}\times\,}
|
||||
\newcommand {\J}{{\vec J}}
|
||||
\renewcommand{\H}{{\vec H}}
|
||||
\newcommand {\E}{{\vec E}}
|
||||
\newcommand{\dcurl}{{\mathbf C}}
|
||||
\newcommand{\dgrad}{{\mathbf G}}
|
||||
\newcommand{\Acf}{{\mathbf A_c^f}}
|
||||
\newcommand{\Ace}{{\mathbf A_c^e}}
|
||||
\renewcommand{\S}{{\mathbf \Sigma}}
|
||||
\newcommand{\St}{{\mathbf \Sigma_\tau}}
|
||||
\newcommand{\T}{{\mathbf T}}
|
||||
\newcommand{\Tt}{{\mathbf T_\tau}}
|
||||
\newcommand{\diag}[1]{\,{\sf diag}\left( #1 \right)}
|
||||
\newcommand{\M}{{\mathbf M}}
|
||||
\newcommand{\MfMui}{{\M^f_{\mu^{-1}}}}
|
||||
\newcommand{\MeSig}{{\M^e_\sigma}}
|
||||
\newcommand{\MeSigInf}{{\M^e_{\sigma_\infty}}}
|
||||
\newcommand{\MeSigO}{{\M^e_{\sigma_0}}}
|
||||
\newcommand{\Me}{{\M^e}}
|
||||
\newcommand{\Mes}[1]{{\M^e_{#1}}}
|
||||
\newcommand{\Mee}{{\M^e_e}}
|
||||
\newcommand{\Mej}{{\M^e_j}}
|
||||
\newcommand{\BigO}[1]{\mathcal{O}\bigl(#1\bigr)}
|
||||
\newcommand{\bE}{\mathbf{E}}
|
||||
\newcommand{\bH}{\mathbf{H}}
|
||||
\newcommand{\B}{\vec{B}}
|
||||
\newcommand{\D}{\vec{D}}
|
||||
\renewcommand{\H}{\vec{H}}
|
||||
\newcommand{\s}{\vec{s}}
|
||||
\newcommand{\bfJ}{\bf{J}}
|
||||
\newcommand{\vecm}{\vec m}
|
||||
\renewcommand{\Re}{\mathsf{Re}}
|
||||
\renewcommand{\Im}{\mathsf{Im}}
|
||||
\renewcommand {\j} { {\vec j} }
|
||||
\newcommand {\h} { {\vec h} }
|
||||
\renewcommand {\b} { {\vec b} }
|
||||
\newcommand {\e} { {\vec e} }
|
||||
\renewcommand {\d} { {\vec d} }
|
||||
\renewcommand {\u} { {\vec u} }
|
||||
\newcommand{\I}{\vec{I}}
|
||||
|
||||
|
||||
|
||||
|
||||
+44
-40
@@ -1,33 +1,33 @@
|
||||
\documentclass[]{article}
|
||||
|
||||
\renewcommand{\div}{\nabla\cdot\,}
|
||||
\newcommand{\grad}{\ensuremath {\vec \nabla}}
|
||||
\newcommand{\curl}{\ensuremath{{\vec \nabla}\times\,}}
|
||||
\newcommand {\J} { {\vec J} }
|
||||
\renewcommand {\H} { {\vec H} }
|
||||
\newcommand {\E} { {\vec E} }
|
||||
\newcommand{\dcurl}{\ensuremath{{\mathbf C}}}
|
||||
\newcommand{\dgrad}{\ensuremath{{\mathbf G}}}
|
||||
\newcommand{\Acf}{\ensuremath{{\mathbf A_c^f}}}
|
||||
\newcommand{\Ace}{\ensuremath{{\mathbf A_c^e}}}
|
||||
\renewcommand{\S}{\ensuremath{{\mathbf \Sigma}}}
|
||||
\newcommand{\St}{\ensuremath{{\mathbf \Sigma_\tau}}}
|
||||
\newcommand{\T}{\ensuremath{{\mathbf T}}}
|
||||
\newcommand{\Tt}{\ensuremath{{\mathbf T_\tau}}}
|
||||
\newcommand{\diag}[1]{\, {\sf diag}\left( #1 \right)}
|
||||
\newcommand{\grad}{\vec \nabla}
|
||||
\newcommand{\curl}{{\vec \nabla}\times\,}
|
||||
\newcommand {\J}{{\vec J}}
|
||||
\renewcommand{\H}{{\vec H}}
|
||||
\newcommand {\E}{{\vec E}}
|
||||
\newcommand{\dcurl}{{\mathbf C}}
|
||||
\newcommand{\dgrad}{{\mathbf G}}
|
||||
\newcommand{\Acf}{{\mathbf A_c^f}}
|
||||
\newcommand{\Ace}{{\mathbf A_c^e}}
|
||||
\renewcommand{\S}{{\mathbf \Sigma}}
|
||||
\newcommand{\St}{{\mathbf \Sigma_\tau}}
|
||||
\newcommand{\T}{{\mathbf T}}
|
||||
\newcommand{\Tt}{{\mathbf T_\tau}}
|
||||
\newcommand{\diag}[1]{\,{\sf diag}\left( #1 \right)}
|
||||
|
||||
%Common mass matricies
|
||||
\newcommand{\M}{\ensuremath{{\mathbf M}}}
|
||||
\newcommand{\MfMui}{\ensuremath{{\M^f_{\mu^{-1}}}}}
|
||||
\newcommand{\MeSig}{\ensuremath{{\M^e_\sigma}}}
|
||||
\newcommand{\MeSigInf}{\ensuremath{{\M^e_{\sigma_\infty}}}}
|
||||
\newcommand{\MeSigO}{\ensuremath{{\M^e_{\sigma_0}}}}
|
||||
\newcommand{\Me}{\ensuremath{{\M^e}}}
|
||||
\newcommand{\Mes}[1]{\ensuremath{{\M^e_{#1}}}}
|
||||
\newcommand{\Mee}{\ensuremath{{\M^e_e}}}
|
||||
\newcommand{\Mej}{\ensuremath{{\M^e_j}}}
|
||||
\newcommand{\M}{{\mathbf M}}
|
||||
\newcommand{\MfMui}{{\M^f_{\mu^{-1}}}}
|
||||
\newcommand{\MeSig}{{\M^e_\sigma}}
|
||||
\newcommand{\MeSigInf}{{\M^e_{\sigma_\infty}}}
|
||||
\newcommand{\MeSigO}{{\M^e_{\sigma_0}}}
|
||||
\newcommand{\Me}{{\M^e}}
|
||||
\newcommand{\Mes}[1]{{\M^e_{#1}}}
|
||||
\newcommand{\Mee}{{\M^e_e}}
|
||||
\newcommand{\Mej}{{\M^e_j}}
|
||||
|
||||
\newcommand{\BigO}[1]{\ensuremath{\mathcal{O}\bigl(#1\bigr)}}
|
||||
\newcommand{\BigO}[1]{\mathcal{O}\bigl(#1\bigr)}
|
||||
|
||||
% ********** TDIP paper
|
||||
|
||||
@@ -75,27 +75,27 @@ Using Gauss-Newton to solve the inverse problem requires the ability to calculat
|
||||
where
|
||||
\begin{subequations}
|
||||
\begin{align}
|
||||
\mathbf{A} =
|
||||
\mathbf{A} =
|
||||
\left[
|
||||
\begin{array}{cc}
|
||||
\frac{1}{\delta t} \mathbf{I} & \dcurl \\
|
||||
\dcurl^\top & -\MeSig
|
||||
\dcurl^\top \MfMui & -\MeSig
|
||||
\end{array}
|
||||
\right] \\
|
||||
\mathbf{B} =
|
||||
\mathbf{B} =
|
||||
\left[
|
||||
\begin{array}{cc}
|
||||
-\frac{1}{\delta t} \mathbf{I} & 0 \\
|
||||
0 & 0
|
||||
\end{array}
|
||||
\right] \\
|
||||
\u^{(k)} = \left[
|
||||
\u^{(k)} = \left[
|
||||
\begin{array}{c}
|
||||
\b^{(k)}\\
|
||||
\e^{(k)}
|
||||
\end{array}
|
||||
\right] \\
|
||||
\s^{(k)} = \left[
|
||||
\s^{(k)} = \left[
|
||||
\begin{array}{c}
|
||||
0\\
|
||||
\Me \j^{(k)}_s
|
||||
@@ -153,7 +153,7 @@ Defining the function $\vec{c}(m,\vec{u})$ to be
|
||||
\end{align}
|
||||
then
|
||||
\begin{align}
|
||||
\frac{\partial \vec{c}}{\partial m} \partial m
|
||||
\frac{\partial \vec{c}}{\partial m} \partial m
|
||||
+ \frac{\partial \vec{c}}{\partial \u} \partial \vec{u} = 0
|
||||
\end{align}
|
||||
or
|
||||
@@ -168,8 +168,8 @@ Differentiating, we find that
|
||||
\end{align}
|
||||
and
|
||||
\begin{align}
|
||||
\frac{\partial \vec{c}}{\partial \sigma} = \mathbf{G}_\sigma =
|
||||
\left[
|
||||
\frac{\partial \vec{c}}{\partial \sigma} = \mathbf{G}_\sigma =
|
||||
\left[
|
||||
\begin{array}{c}
|
||||
g_\sigma^{(1)}\\
|
||||
g_\sigma^{(2)}\\
|
||||
@@ -181,8 +181,8 @@ and
|
||||
with
|
||||
\begin{subequations}
|
||||
\begin{align}
|
||||
g_\sigma^{(n)} =
|
||||
\left[
|
||||
g_\sigma^{(n)} =
|
||||
\left[
|
||||
\begin{array}{c}
|
||||
\mathbf{0} \\
|
||||
- \diag{\e^{(n)}} \Ace \diag{\vec{V}}
|
||||
@@ -215,7 +215,7 @@ Multiplying $\mathbf{J}$ onto a vector can be broken into three steps
|
||||
|
||||
\begin{subequations}
|
||||
\begin{align}
|
||||
\frac{1}{\delta t} \vec{y}_{b}^{(1)} + \dcurl \vec{y}_{e}^{(1)} = 0 \\
|
||||
\dcurl \vec{y}_{e}^{(1)} + \frac{1}{\delta t} \vec{y}_{b}^{(1)} = 0 \\
|
||||
\dcurl^\top \MfMui \vec{y}_b^{(1)} - \MeSig \vec{y}_e^{(1)} = \vec{p}_e^{(1)}
|
||||
\end{align}
|
||||
\end{subequations}
|
||||
@@ -231,14 +231,18 @@ Multiplying $\mathbf{J}$ onto a vector can be broken into three steps
|
||||
|
||||
\begin{subequations}
|
||||
\begin{align}
|
||||
\dcurl \vec{y}_{e}^{(t+1)} + \frac{1}{\delta t} \vec{y}_{b}^{(t+1)} - \frac{1}{\delta t} \vec{y}_{b}^{(t)} = 0 \\
|
||||
\vec{y}_e^{(t+1)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t+1)} - \MeSig^{-1} \vec{p}_e^{(t+1)}
|
||||
\dcurl \vec{y}_{e}^{(t+1)} + \frac{1}{\delta t} \vec{y}_{b}^{(t+1)}
|
||||
{\color{red}- \frac{1}{\delta t} \vec{y}_{b}^{(t)} }
|
||||
= 0 \\
|
||||
\dcurl^\top \MfMui \vec{y}_b^{(t+1)} - \MeSig \vec{y}_e^{(t+1)} = \vec{p}_e^{(t+1)}
|
||||
\end{align}
|
||||
\end{subequations}
|
||||
|
||||
|
||||
\begin{subequations}
|
||||
\begin{align}
|
||||
\left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(1)} = \frac{1}{\delta t} \MfMui \vec{y}_b^{(t)} + \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(1)} \\
|
||||
\left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(t+1)} =
|
||||
{\color{red} \frac{1}{\delta t} \MfMui \vec{y}_b^{(t)} }
|
||||
+ \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(t+1)} \\
|
||||
\vec{y}_e^{(t+1)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t+1)} - \MeSig^{-1} \vec{p}_e^{(t+1)}
|
||||
\end{align}
|
||||
\end{subequations}
|
||||
@@ -247,4 +251,4 @@ Multiplying $\mathbf{J}$ onto a vector can be broken into three steps
|
||||
|
||||
|
||||
|
||||
\end{document}
|
||||
\end{document}
|
||||
|
||||
@@ -3,7 +3,7 @@ from SimPEG.Problem import BaseProblem
|
||||
from simpegEM.Utils import Sources
|
||||
from FieldsTDEM import FieldsTDEM
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.Utils import sdiag
|
||||
from SimPEG.Utils import sdiag, mkvc
|
||||
import numpy as np
|
||||
|
||||
|
||||
|
||||
+24
-57
@@ -16,7 +16,7 @@ class ProblemTDEM_b(ProblemBaseTDEM):
|
||||
ProblemBaseTDEM.__init__(self, mesh, model, **kwargs)
|
||||
|
||||
solType = 'b'
|
||||
|
||||
|
||||
####################################################
|
||||
# Internal Methods
|
||||
####################################################
|
||||
@@ -56,19 +56,19 @@ class ProblemTDEM_b(ProblemBaseTDEM):
|
||||
ei = u.get_e(i)
|
||||
pVal = np.empty_like(ei)
|
||||
for j in range(ei.shape[1]):
|
||||
pVal[:,j] = -ei[:,j]*c
|
||||
|
||||
pVal[:,j] = -ei[:,j]*c
|
||||
|
||||
p.set_e(pVal,i)
|
||||
p.set_b(np.zeros((self.mesh.nF,1)), i)
|
||||
return p
|
||||
|
||||
def solveAh(self, m, p):
|
||||
def AhRHS(tInd, u):
|
||||
rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p.get_e(tInd)
|
||||
if tInd == 0:
|
||||
return self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p.get_e(tInd)
|
||||
else:
|
||||
dt = self.getDt(tInd)
|
||||
return self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p.get_e(tInd) + 1./dt*self.MfMui*u.get_b(tInd-1)
|
||||
return rhs
|
||||
dt = self.getDt(tInd)
|
||||
return rhs + 1./dt*self.MfMui*u.get_b(tInd-1)
|
||||
|
||||
def AhCalcFields(sol, solType, tInd):
|
||||
b = sol
|
||||
@@ -126,61 +126,28 @@ if __name__ == '__main__':
|
||||
|
||||
prb = EM.TDEM.ProblemTDEM_b(mesh, model)
|
||||
# prb.setTimes([1e-5, 5e-5, 2.5e-4], [150, 150, 150])
|
||||
prb.setTimes([1e-5, 5e-5, 2.5e-4], [10, 10, 10])
|
||||
# prb.setTimes([1e-5], [10])
|
||||
# prb.setTimes([1e-5, 5e-5, 2.5e-4], [10, 10, 10])
|
||||
prb.setTimes([1e-5], [1])
|
||||
prb.pair(dat)
|
||||
|
||||
# sigma = np.ones(mesh.nCz)*1e-8
|
||||
# sigma[mesh.vectorCCz<0] = 0.1
|
||||
|
||||
|
||||
# u = prb.fields(sigma)
|
||||
# Ahu = prb.AhVec(sigma, u)
|
||||
|
||||
# Random fields
|
||||
sigma = np.random.rand(mesh.nCz)
|
||||
# f = FieldsTDEM(prb.mesh, 1, prb.times.size, 'b')
|
||||
# for i in range(f.nTimes):
|
||||
# f.set_b(np.random.rand(mesh.nF, 1), i)
|
||||
# f.set_e(np.random.rand(mesh.nE, 1), i)
|
||||
f = prb.fields(sigma)
|
||||
|
||||
dm = np.random.rand(mesh.nCz)
|
||||
|
||||
for h in np.logspace(0, -10, 10):
|
||||
# print h
|
||||
a = np.linalg.norm(prb.AhVec(sigma+h*dm, f).fieldVec() - prb.AhVec(sigma, f).fieldVec())
|
||||
b = np.linalg.norm(prb.AhVec(sigma+h*dm, f).fieldVec() - prb.AhVec(sigma, f).fieldVec() - h*prb.G(sigma, dm, u=f).fieldVec())
|
||||
print a, b, b/a
|
||||
# print
|
||||
# h = 1.
|
||||
plt.semilogy(np.abs(prb.AhVec(sigma+h*dm,f).fieldVec() - prb.AhVec(sigma, f).fieldVec()), 'ko')
|
||||
plt.semilogy(np.abs(h*prb.G(sigma, dm, u=f).fieldVec()), 'rx')
|
||||
# plt.semilogy(prb.AhVec(sigma+h*dm, f).fieldVec() - prb.AhVec(sigma, f).fieldVec() - h*prb.G(sigma, dm, u=f).fieldVec(),'ko')
|
||||
|
||||
f = FieldsTDEM(prb.mesh, 1, prb.times.size, 'b')
|
||||
for i in range(f.nTimes):
|
||||
f.set_b(np.zeros((mesh.nF, 1)), i)
|
||||
f.set_e(np.random.rand(mesh.nE, 1), i)
|
||||
|
||||
Ahf = prb.AhVec(sigma, f)
|
||||
f_test = prb.solveAh(sigma, Ahf)
|
||||
|
||||
e0 = f.get_e(0)
|
||||
e1 = f_test.get_e(0)
|
||||
b0 = f.get_b(0)
|
||||
b1 = f_test.get_b(0)
|
||||
plt.semilogy(np.abs(e0))
|
||||
plt.semilogy(np.abs(e1),'r')
|
||||
plt.show()
|
||||
|
||||
# plt.show()
|
||||
|
||||
# f = prb.fields(sigma)
|
||||
# print f.fieldVec()
|
||||
|
||||
# prb.AhVec(sigma,f)
|
||||
|
||||
# prb.G(prb.sigma, prb.sigma)
|
||||
# prb.solveAh(prb.sigma, f)
|
||||
# prb.J(prb.sigma, prb.sigma, f)
|
||||
|
||||
# from SimPEG.Tests import checkDerivative
|
||||
# m0 = sigma
|
||||
# dx = np.zeros_like(sigma)
|
||||
# dx[prb.mesh.vectorCCz<0] = 1e-4
|
||||
# derChk = lambda m: [dat.dpred(m), lambda mx: prb.J(m0, mx, u=f)]
|
||||
# passed = checkDerivative(derChk, m0, dx=dx, plotIt=False)
|
||||
|
||||
# bz_calc = dat.dpred(sigma)
|
||||
# bz_ana = mu_0*hzAnalyticDipoleT(dat.rxLoc[0], prb.times, sigma[0])
|
||||
|
||||
# plt.loglog(prb.times, np.abs(bz_calc.flatten()), label='TDEM_b')
|
||||
# plt.loglog(prb.times, np.abs(bz_ana), 'r', label='Analytic')
|
||||
# plt.legend()
|
||||
# plt.show()
|
||||
|
||||
Reference in New Issue
Block a user