mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-01 04:51:05 +08:00
Adding notebooks on the MT1D problem
This commit is contained in:
File diff suppressed because one or more lines are too long
@@ -0,0 +1,447 @@
|
||||
{
|
||||
"metadata": {
|
||||
"name": "",
|
||||
"signature": "sha256:74ba4ac0804b189d65dce7f7af6f88a605b83052cf38c5acb492600297a13398"
|
||||
},
|
||||
"nbformat": 3,
|
||||
"nbformat_minor": 0,
|
||||
"worksheets": [
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"%pylab inline"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"output_type": "stream",
|
||||
"stream": "stdout",
|
||||
"text": [
|
||||
"Populating the interactive namespace from numpy and matplotlib\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 1
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Simple notebook on performing a 1D MT problem.\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"#from IPython.display import Latex"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 1
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Maxwell's equations in 1D are as follows\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"$i \\omega b = - \\partial_z e $ \n",
|
||||
"\n",
|
||||
"$ s = \\partial_z(\\mu^{-1} b) - \\sigma(z) e$\n",
|
||||
"\n",
|
||||
"$b(0) = 1 \\hspace{1cm} ;\\hspace{1cm} b(-\\infty) = 0$\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"where $e = \\widehat{\\overrightarrow{E}}_x $ and $b = \\widehat{\\overrightarrow{B}}_y$\n",
|
||||
"\n",
|
||||
"In weak form the equations become\n",
|
||||
"\n",
|
||||
"$i \\omega(b,f) = - (\\partial_ze,f)$ \n",
|
||||
"\n",
|
||||
"$(s,w) = -(\\mu^{-1} b, \\partial_z w) - (\\sigma(z) e, w)$\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"where f and w are abritrary functions living in same discritizational space as b and e, respectivily.\n",
|
||||
"\n",
|
||||
"We consider e on nodes and b on cell centers. This way the derivative of any nodal function becomes\n",
|
||||
"\n",
|
||||
"$ (\\partial_z u)_k \\approx $\n",
|
||||
"\n",
|
||||
"$h_{k}^{-1} ( u_{k+\\frac{1}{2}} - u_{k-\\frac{1}{2}}) + O(h^2) $\n",
|
||||
"\n",
|
||||
"Matrix form\n",
|
||||
"\n",
|
||||
"$ e_z \\approx \\textbf{L}^{-1} \\textbf{G} e = \\begin{bmatrix} h_1^{-1} & & & \\\\\\\\ & h_2^{-1} & & \\\\\\\\ & & \\ddots & \\\\\\\\ & & & h_n^{-1} \\end{bmatrix}^{(n,n)}\n",
|
||||
"\\begin{bmatrix} -1 & 1 & & & & \\\\\\\\ & -1 & 1 & & & \\\\\\\\ & & \\ddots & \\ddots & \\\\\\\\ & & & -1 & 1 \\end{bmatrix}^{(n,n+1)} \n",
|
||||
"\\begin{bmatrix} e_1 \\\\\\\\ \\\\\\\\ \\vdots \\\\\\\\ \\\\\\\\ e_{n+1} \\end{bmatrix}^{(n+1,1)} $\n",
|
||||
"\n",
|
||||
"where $ \\textbf{L} = diag(h) $ is the cell size and $ \\textbf{G}$ is the gradient operator with -1,1 representing the topology of the mesh, taking the difference between adjoint cells.\n",
|
||||
"\n",
|
||||
"We need to compute 2 inner products, on cell centers and from nodes to cell centers.\n",
|
||||
"\n",
|
||||
"Cell centers inner product is\n",
|
||||
"\n",
|
||||
"$ (b,f) \\approx \\sum\\limits_k h_k \\textbf{b}_k \\textbf{f}_k + O(h^2) $\n",
|
||||
"\n",
|
||||
"and in matrix from\n",
|
||||
"\n",
|
||||
"$ (b,) \\approx \\textbf{b}^T \\textbf{M}^f \\textbf{f}$ and $ (\\mu^{-1} b,f) \\approx \\textbf{b}^T \\textbf{M}_{\\mu}^f \\textbf{f}$ \n",
|
||||
"\n",
|
||||
"where $ \\textbf{M}_{\\mu}^f = diag(\\textbf{h} \\odot \\mu^{-1}) $ and $ \\textbf{M}^f = diag(\\textbf{h}) $ are the matrices.\n",
|
||||
"Nodes to cell centers inner product is\n",
|
||||
"\n",
|
||||
"$ (\\sigma e, w) \\approx \\sum\\limits_k \\frac{h_k \\sigma_k}{4} ( e_{k+\\frac{1}{2}} w_{k+\\frac{1}{2}} + e_{k+\\frac{1}{2}} w_{k+\\frac{1}{2}} )$\n",
|
||||
"\n",
|
||||
"and in matrix from\n",
|
||||
"\n",
|
||||
"$ (\\sigma e, w ) \\approx (\\textbf{h} \\odot \\sigma )^T ( \\textbf{A}_v (\\textbf{e} \\odot \\textbf{w} = \\textbf{w}^T diag(\\textbf{A}_v^T (\\textbf{h} \\odot \\sigma)) \\textbf{e} $\n",
|
||||
"\n",
|
||||
"Here $\\odot$ is a point wise Hadamard product and $\\textbf{A}_v$ is the averaging operator/matrix from nodes to cell centers\n",
|
||||
"\n",
|
||||
"$ \\textbf{A}_v = \\begin{bmatrix} \\frac{1}{2} & \\frac{1}{2} & & \\\\\\\\ & \\ddots & \\ddots & \\\\\\\\ & & \\frac{1}{2} & \\frac{1}{2} \\end{bmatrix} ^{(n+1 , n)} $\n",
|
||||
"\n",
|
||||
"The sigma mass matrix is defined as \n",
|
||||
"\n",
|
||||
"$ \\textbf{M}_{\\sigma}^{e} = diag(\\textbf{A}_v^T (\\textbf{h} \\odot \\sigma)) $\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"In the MT problem there is no source in the domain, so $ s = 0 $. How ever the boundary conditions provide the right hand side where \n",
|
||||
"\n",
|
||||
"$ (\\partial_z \\mu^{-1} b, w ) = - (\\mu^{-1} b, \\partial_z w ) + (\\mu^{-1} b w )|_0^{end} $\n",
|
||||
"\n",
|
||||
"where \n",
|
||||
"\n",
|
||||
"$ (\\mu^{-1} b w )|_0^{end} = \\textbf{bc}^T (\\textbf{BC w}) $\n",
|
||||
"\n",
|
||||
"here $\\textbf{BC}$ is an matrix operator that extracts the boundary elements from $ \\textbf{w}$ and $\\textbf{bc} $ are the known boundary condintions. For the 1D case with homogenous boundary conditions we have \n",
|
||||
"\n",
|
||||
"$ \\textbf{B} = \\begin{bmatrix} -1 & 0 \\\\\\\\ \\vdots & \\vdots \\\\\\\\ 0 & 1 \\end{bmatrix} $ \n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"The weak form is \n",
|
||||
"\n",
|
||||
"$ (\\mu^{-1} b, \\partial_z w) + (\\sigma e, w) = (\\mu^{-1} b w )|_0^{end} $\n",
|
||||
"\n",
|
||||
"$ (i \\omega b,f) + (\\partial_z e , f) = 0 $\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Using the above matrix represntation we get the Maxwells equations in following form\n",
|
||||
"\n",
|
||||
"$ i \\omega \\textbf{f}^T \\textbf{M}^f \\textbf{b} + \\textbf{f}^T \\textbf{M}^f \\textbf{L}^{-1} \\textbf{G} \\textbf{e} = 0 $ \n",
|
||||
"\n",
|
||||
"$ \\textbf{w}^T \\textbf{G}^T \\textbf{L}^{-1} \\textbf{M}^f_{\\mu} \\textbf{b} + \\textbf{w}^T \\textbf{M}_{\\sigma}^e \\textbf{e} = \\textbf{w}^T \\textbf{bc}^T \\textbf{BC} $ \n",
|
||||
"\n",
|
||||
"Here we use that \n",
|
||||
"\n",
|
||||
"$ (\\textbf{b},\\textbf{f}) \\approx \\textbf{b}^T \\textbf{M}^f \\textbf{f} = \\textbf{f}^T \\textbf{M}^f \\textbf{b} $ \n",
|
||||
"\n",
|
||||
"since $\\textbf{M}^f$ is a symmetric diaoganal matrix of size n by n and $\\textbf{b}$ and $\\textbf{f}$ are vectors of length n."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We eliminate testing vectors and get system of equations to solve\n",
|
||||
"\n",
|
||||
"$ i \\omega \\textbf{b} + \\textbf{L}^{-1} \\textbf{G} \\textbf{e} = 0 $ \n",
|
||||
"\n",
|
||||
"$ \\textbf{G}^T \\textbf{L}^{-1} \\textbf{M}^f_{\\mu} \\textbf{b} + \\textbf{M}_{\\sigma}^e \\textbf{e} = \\textbf{bc}^T \\textbf{BC} $ \n",
|
||||
"\n",
|
||||
"and as $ \\textbf{A} \\textbf{x} = \\textbf{bc} $ system that we will solve, where \n",
|
||||
"\n",
|
||||
"$ \\textbf{A} = \\begin{bmatrix} \\textbf{G}^T \\textbf{L}^{-1} \\textbf{M}^f_{\\mu} & \\textbf{M}_{\\sigma}^e \\\\\\\\ i \\omega & \\textbf{L}^{-1} \\textbf{G} \\end{bmatrix} $\n",
|
||||
"\n",
|
||||
"$ \\textbf{x} = \\begin{bmatrix} \\textbf{b} \\\\\\\\ \\textbf{e} \\end{bmatrix} $\n",
|
||||
"\n",
|
||||
"$ \\textbf{bc} = \\begin{bmatrix} \\textbf{bc}^T \\textbf{BC} \\\\\\\\ \\textbf{0} \\end{bmatrix} $ "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"import sys\n",
|
||||
"sys.path.append('C:/GudniWork/Codes/python/simpeg')\n",
|
||||
"import SimPEG as simpeg, numpy as np, scipy, scipy.sparse as sp"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"output_type": "stream",
|
||||
"stream": "stdout",
|
||||
"text": [
|
||||
"Efficiency Warning: Interpolation will be slow, use setup.py!\n",
|
||||
"\n",
|
||||
" python setup.py build_ext --inplace\n",
|
||||
" \n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 2
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We have\n",
|
||||
"\n",
|
||||
"$ i \\omega b = - \\partial_z e \\hspace{1cm} ; \\hspace{1cm} \\partial_z(\\mu^{-1} b) - \\sigma(z) e \\hspace{1cm} ;\\hspace{1cm} b(0) = 1 \\hspace{1cm} ;\\hspace{1cm} b(-\\infty) = 0 $\n",
|
||||
"\n",
|
||||
" \n",
|
||||
"To deal with boundary: we assume that below depth L both $ \\sigma $ and $ \\mu $ are constants ($ z < - L $). At the boundary we have that \n",
|
||||
"\n",
|
||||
"$ e = c \\exp(ikz) \\hspace{0.2cm} where \\hspace{0.2cm} k = \\sqrt{i\\omega\\mu\\sigma} $. \n",
|
||||
"\n",
|
||||
"Therefore for $ z < - L $ we have that\n",
|
||||
"\n",
|
||||
"$\\omega b - k e = 0 $.\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We discretize the e field on the nodes and b field at the cell centers. The system we want to solve is \n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"$\\begin{bmatrix} i \\omega & \\frac{\\partial}{\\partial z} \\\\\\\\ \\frac{1}{\\mu} \\frac{\\partial}{\\partial z} & -\\sigma \\end{bmatrix} \n",
|
||||
"\\begin{bmatrix} b \\\\\\\\ e \\end{bmatrix} = \\begin{bmatrix} s1 \\\\\\\\ s2 \\end{bmatrix}$\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 2
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
" # Set up the problem\n",
|
||||
"mu = 4*np.pi*1e-7\n",
|
||||
"eps0 = 8.85e-12\n",
|
||||
"# Frequency\n",
|
||||
"fr = np.array([1e1]) #np.logspace(0,5,200) #np.array([2000]) #np.logspace(-4,5,82)\n",
|
||||
"omega = 2*np.pi*fr\n",
|
||||
"# Mesh\n",
|
||||
"sig0 = 1e-2\n",
|
||||
"#L = 3*np.sqrt(2/(mu*omega[0]*sig0))\n",
|
||||
"#nn=np.ceil(np.log(0.3*L + 1)/np.log(1.3))\n",
|
||||
"#h = 5*(1.3**(np.arange(nn+1)))\n",
|
||||
"\n",
|
||||
"h = np.ones(18)\n",
|
||||
"x0 = np.array([0])\n",
|
||||
"#sig = sig0*np.ones((len(h),1)) \n",
|
||||
"#sig[0:50] = 0.1\n",
|
||||
"#sig[50:100] = 1\n",
|
||||
"# Make the mesh\n",
|
||||
"mesh = simpeg.Mesh.TensorMesh([h],x0)\n",
|
||||
"sig = np.zeros(mesh.nC) + 1e-8\n",
|
||||
"sig[mesh.vectorCCx<=0] = 1e-2"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 3
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"fr,omega"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "pyout",
|
||||
"prompt_number": 4,
|
||||
"text": [
|
||||
"(array([ 10.]), array([ 62.83185307]))"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 4
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"# Make the operators\n",
|
||||
"G = mesh.nodalGrad\n",
|
||||
"Av = mesh.aveN2CC\n",
|
||||
"Li = scipy.sparse.spdiags(1/mesh.hx,0,mesh.nNx,mesh.nNx)\n",
|
||||
"Mmu = scipy.sparse.spdiags(mesh.hx/mu,0,mesh.nCx,mesh.nCx)\n",
|
||||
"Msig = scipy.sparse.spdiags(Av.T.dot(mesh.hx*sig.ravel()),0,mesh.nNx,mesh.nNx)\n",
|
||||
"# The boundaries\n",
|
||||
"bc_b = np.zeros((mesh.nCx,1))\n",
|
||||
"bc_b[0] = -1 # Set the top b field to 1\n",
|
||||
"bc_e = np.zeros((mesh.nNx,1))\n",
|
||||
"# Make the sparse matrix\n",
|
||||
"bc = sp.vstack((bc_b,bc_e))\n"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 5
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"b = np.empty((mesh.nCx,len(omega)),dtype=np.complex64)\n",
|
||||
"e = np.empty((mesh.nNx,len(omega)),dtype=np.complex64)\n",
|
||||
"# Loop all the frequencies\n",
|
||||
"for nrOm, om in enumerate(omega):\n",
|
||||
" # Left hand side\n",
|
||||
" A = sp.vstack((sp.hstack(( -G.conj().T.dot(Mmu), - Msig)), sp.hstack((1j*om*scipy.sparse.identity(mesh.nCx) , G))))\n",
|
||||
" #A = A.tocsr\n",
|
||||
" # Solve the system\n",
|
||||
" bef = scipy.sparse.linalg.spsolve(A,bc)\n",
|
||||
" # Sort the output\n",
|
||||
" b[:,nrOm] = bef[0:mesh.nCx]\n",
|
||||
" e[:,nrOm] = bef[mesh.nCx::]\n"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"output_type": "stream",
|
||||
"stream": "stderr",
|
||||
"text": [
|
||||
"/home/Gudni/anaconda/lib/python2.7/site-packages/scipy/sparse/linalg/dsolve/linsolve.py:90: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format\n",
|
||||
" SparseEfficiencyWarning)\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 6
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"# Plot the solution\n",
|
||||
"z=e[0,:]/(b[0,:]/mu)\n",
|
||||
"app_res = ((1./(8e-7*np.pi**2))/fr)*np.abs(z)**2\n",
|
||||
"app_phs = np.arctan(z.imag/z.real)*(180/np.pi)\n",
|
||||
"ax_res = plt.subplot(2,1,1)\n",
|
||||
"ax_res.loglog(fr,app_res)\n",
|
||||
"ax_phs = plt.subplot(2,1,2)\n",
|
||||
"ax_phs.semilogx(fr,app_phs)"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "pyout",
|
||||
"prompt_number": 7,
|
||||
"text": [
|
||||
"[<matplotlib.lines.Line2D at 0x7fa8671e1fd0>]"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 7
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"plt.show()"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 8
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"# Calculate the impedance\n",
|
||||
"z = e[0,:]/(b[0,:]/mu)\n",
|
||||
"z"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "pyout",
|
||||
"prompt_number": 9,
|
||||
"text": [
|
||||
"array([-5714285.5+0.00050081j], dtype=complex64)"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 9
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"app_res"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "pyout",
|
||||
"prompt_number": 10,
|
||||
"text": [
|
||||
"array([ 4.13555827e+17])"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 10
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 10
|
||||
}
|
||||
],
|
||||
"metadata": {}
|
||||
}
|
||||
]
|
||||
}
|
||||
@@ -0,0 +1,191 @@
|
||||
{
|
||||
"metadata": {
|
||||
"name": "",
|
||||
"signature": "sha256:4f51688cd2ee8a11dad3df1928925d3c9cad0da43a3f6a3c3c840024caae5fe1"
|
||||
},
|
||||
"nbformat": 3,
|
||||
"nbformat_minor": 0,
|
||||
"worksheets": [
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Issues with padding cells and high frequencies in the analytic MT layered earth."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"import SimPEG as simpeg\n",
|
||||
"elev = 300\n",
|
||||
"# 3D mesh and model\n",
|
||||
"M = simpeg.Mesh.TensorMesh([[(100,5,-1.5),(100.,5),(100,5,1.5)],[(100,5,-1.5),(100.,5),(100,5,1.5)],[(100,10,-1.5),(100.,10),(100,10,1.5)]], x0=['C','C','C'])\n",
|
||||
"conds = [1,1e-2]\n",
|
||||
"sig = simpeg.Utils.ModelBuilder.defineBlock(M.gridCC,[-10000,-10000,-200],[10000,10000,0],conds)\n",
|
||||
"sig[M.gridCC[:,2]>elev] = 1e-8\n",
|
||||
"sig[M.gridCC[:,2]<-600] = 1e-1\n",
|
||||
"# Make the 1D mesh and model\n",
|
||||
"mesh1d = simpeg.Mesh.TensorMesh([M.hz],np.array([M.x0[2]]))\n",
|
||||
"sig1D = M.r(sig,'CC','CC','M')[0,0,:]"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"output_type": "stream",
|
||||
"stream": "stdout",
|
||||
"text": [
|
||||
"Efficiency Warning: Interpolation will be slow, use setup.py!\n",
|
||||
"\n",
|
||||
" python setup.py build_ext --inplace\n",
|
||||
" \n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 1
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"# Run for high frequency\n",
|
||||
"freq = 1e4"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 3
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"# Run the analytic problem\n",
|
||||
"import simpegMT as simpegmt\n",
|
||||
"anaEd, anaEu, anaHd, anaHu = simpegmt.Utils.MT1Danalytic.getEHfields(mesh1d,sig1D,freq,np.array([300]))\n",
|
||||
"anaE = anaEd+anaEu\n",
|
||||
"anaH = anaHd+anaHu\n",
|
||||
"anaZ = anaE/anaH\n"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 4
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"anaZ"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "pyout",
|
||||
"prompt_number": 5,
|
||||
"text": [
|
||||
"array([ nan+nanj])"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 5
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Returns nan because in the analytic solution the propagation of the fields in the layer \"blows\" up."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"sig = 10\n",
|
||||
"sig0 = 20\n",
|
||||
"mu = 4*np.pi*1e-7\n",
|
||||
"eps = 8.85*1e-12\n",
|
||||
"for h in [10000,5000,1000,500,100,50,10]:\n",
|
||||
" w = 2*np.pi*freq\n",
|
||||
" k0 = np.sqrt(eps*mu*w**2-1j*mu*sig0*w)\n",
|
||||
" k = np.sqrt(eps*mu*w**2-1j*mu*sig*w)\n",
|
||||
" zp = (w*mu)/k\n",
|
||||
" yp1 = k0/(w*mu)\n",
|
||||
" # Convert fields to down/up going components in layer below current layer\n",
|
||||
" Pj1 = np.array([[1,1],[yp1,-yp1]])\n",
|
||||
" # Convert fields to down/up going components in current layer\n",
|
||||
" Pjinv = 1./2*np.array([[1,zp],[1,-zp]])\n",
|
||||
" # Propagate down and up components through the current layer\n",
|
||||
" elamh = np.array([[np.exp(-1j*k*h),0],[0,np.exp(1j*k*h)]])\n",
|
||||
" UD = elamh.dot(Pjinv.dot(Pj1)).dot([1,0])\n",
|
||||
" print h, w, k \n",
|
||||
" print elamh\n",
|
||||
" #print Pj1, Pjinv, elamh\n",
|
||||
" print UD"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"output_type": "stream",
|
||||
"stream": "stdout",
|
||||
"text": [
|
||||
"10000 62831.8530718 (0.628318548187-0.628318513249j)\n",
|
||||
"[[ 0. -0.j 0. +0.j]\n",
|
||||
" [ 0. +0.j inf+infj]]\n",
|
||||
"[ 0. +0.j nan+nanj]\n",
|
||||
"5000 62831.8530718 (0.628318548187-0.628318513249j)\n",
|
||||
"[[ 0. -0.j 0. +0.j]\n",
|
||||
" [ 0. +0.j inf+infj]]\n",
|
||||
"[ 0. +0.j nan+nanj]\n",
|
||||
"1000 62831.8530718 (0.628318548187-0.628318513249j)\n",
|
||||
"[[ 1.33271357e-273 -2.32814399e-278j 0.00000000e+000 +0.00000000e+000j]\n",
|
||||
" [ 0.00000000e+000 +0.00000000e+000j 7.50348781e+272 +1.31079930e+268j]]\n",
|
||||
"[ 1.60872758e-273 -2.81162844e-278j -1.55402321e+272 -2.70737840e+267j]\n",
|
||||
"500 62831.8530718 (0.628318548187-0.628318513249j)\n",
|
||||
"[[ 3.65063497e-137 -3.18868363e-142j 0.00000000e+000 +0.00000000e+000j]\n",
|
||||
" [ 0.00000000e+000 +0.00000000e+000j 2.73924950e+136 +2.39262488e+131j]]\n",
|
||||
"[ 4.40670622e-137 -3.85267017e-142j -5.67317146e+135 -4.92836188e+130j]\n",
|
||||
"100 62831.8530718 (0.628318548187-0.628318513249j)\n",
|
||||
"[[ 5.15790907e-28 -9.01045454e-34j 0.00000000e+00 +0.00000000e+00j]\n",
|
||||
" [ 0.00000000e+00 +0.00000000e+00j 1.93877012e+27 +3.38687631e+21j]]\n",
|
||||
"[ 6.22614702e-28 -1.09272824e-33j -4.01532439e+26 -6.82387175e+20j]\n",
|
||||
"50 62831.8530718 (0.628318548187-0.628318513249j)\n",
|
||||
"[[ 2.27110305e-14 -1.98371768e-20j 0.00000000e+00 +0.00000000e+00j]\n",
|
||||
" [ 0.00000000e+00 +0.00000000e+00j 4.40314674e+13 +3.84597256e+07j]]\n",
|
||||
"[ 2.74146389e-14 -2.41688373e-20j -9.11921548e+12 -7.53244599e+06j]\n",
|
||||
"10 62831.8530718 (0.628318548187-0.628318513249j)\n",
|
||||
"[[ 1.86744306e-03 -3.26227365e-10j 0.00000000e+00 +0.00000000e+00j]\n",
|
||||
" [ 0.00000000e+00 +0.00000000e+00j 5.35491562e+02 +9.35460925e-05j]]\n",
|
||||
"[ 2.25420318e-03 -4.12148003e-10j -1.10903934e+02 -1.41102131e-05j]\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 11
|
||||
},
|
||||
{
|
||||
"cell_type": "heading",
|
||||
"level": 3,
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Is there a smart way to \"fix\" this so that the 1D layering can be used for the analytic solution?"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": []
|
||||
}
|
||||
],
|
||||
"metadata": {}
|
||||
}
|
||||
]
|
||||
}
|
||||
Reference in New Issue
Block a user