From 732ff882e0c1fd79d4b00497ec564f6e67db0d4c Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 12 Feb 2014 15:02:55 -0800 Subject: [PATCH] minor changes to tdem writeup and code (still buggy!) --- docs/api_TDEM.rst | 43 +++++++++++++++++--- notes/tem/tem.tex | 84 ++++++++++++++++++++------------------- simpegEM/TDEM/BaseTDEM.py | 2 +- simpegEM/TDEM/TDEM_b.py | 81 +++++++++++-------------------------- 4 files changed, 107 insertions(+), 103 deletions(-) diff --git a/docs/api_TDEM.rst b/docs/api_TDEM.rst index 5d7f9a65..3c940b10 100644 --- a/docs/api_TDEM.rst +++ b/docs/api_TDEM.rst @@ -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}} diff --git a/notes/tem/tem.tex b/notes/tem/tem.tex index 4d7058c8..e0de4c1e 100644 --- a/notes/tem/tem.tex +++ b/notes/tem/tem.tex @@ -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} \ No newline at end of file +\end{document} diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py index c3663c3c..180f7b27 100644 --- a/simpegEM/TDEM/BaseTDEM.py +++ b/simpegEM/TDEM/BaseTDEM.py @@ -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 diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 154c8c45..1ad15b49 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -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()