mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-28 11:26:12 +08:00
Merge branch 'master' of https://github.com/simpeg/simpegmt
This commit is contained in:
Gudni Karl
+1
-1
@@ -30,7 +30,7 @@ install:
|
||||
|
||||
# Run test
|
||||
script:
|
||||
- nosetests --with-cov --cov SimPEG --cov-config .coveragerc -v -s
|
||||
- nosetests --with-cov --cov simpegMT --cov-config .coveragerc -v -s
|
||||
|
||||
# Calculate coverage
|
||||
after_success:
|
||||
|
||||
+378
-37
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
+961
-260
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": {}
|
||||
}
|
||||
]
|
||||
}
|
||||
@@ -0,0 +1,128 @@
|
||||
{
|
||||
"metadata": {
|
||||
"name": "",
|
||||
"signature": "sha256:260d0d5c93203805d78b8fa9e20a9079ce0e65e4b6a3453fedfcda305fffaa80"
|
||||
},
|
||||
"nbformat": 3,
|
||||
"nbformat_minor": 0,
|
||||
"worksheets": [
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"import SimPEG as simpeg\n",
|
||||
"import simpegMT as simpegmt\n"
|
||||
],
|
||||
"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": [
|
||||
"#Define the mesh\n",
|
||||
"m1d = simpeg.Mesh.TensorMesh([[(100,5,1.5),(100.,10),(100,5,1.5)]], x0=['C'])\n",
|
||||
"sigma = np.zeros(m1d.nC) + 2e-3\n",
|
||||
"sigma[m1d.gridCC[:]>200] = 1e-8\n"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 6
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"# Calculate the analytic fields\n",
|
||||
"freqs = np.logspace(4,-1,26)\n",
|
||||
"Z = []\n",
|
||||
"for freq in freqs:\n",
|
||||
" Ed, Eu, Hd, Hu = simpegmt.Utils.getEHfields(m1d,sigma,freq,np.array([200]))\n",
|
||||
" Z.append((Ed + Eu)/(Hd + Hu))"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 17
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"Zarr = np.concatenate(Z)"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 19
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"def appResPhs(freq,z):\n",
|
||||
" app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2\n",
|
||||
" app_phs = np.arctan2(-z.imag,z.real)*(180/np.pi)\n",
|
||||
" return app_res, app_phs\n",
|
||||
"app_r, app_p = appResPhs(freqs,Zarr)"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 20
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"app_r"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "pyout",
|
||||
"prompt_number": 21,
|
||||
"text": [
|
||||
"array([ 499.99998067, 499.99999231, 499.99999694, 499.99999878,\n",
|
||||
" 499.99999951, 499.99999981, 499.99999992, 499.99999997,\n",
|
||||
" 499.99999999, 500. , 500. , 500. ,\n",
|
||||
" 500. , 500. , 500. , 500. ,\n",
|
||||
" 500. , 500. , 500. , 500. ,\n",
|
||||
" 500. , 500. , 500. , 500. ,\n",
|
||||
" 500. , 500. ])"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 21
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": []
|
||||
}
|
||||
],
|
||||
"metadata": {}
|
||||
}
|
||||
]
|
||||
}
|
||||
@@ -26,12 +26,12 @@ for loc in rx_loc:
|
||||
rxList.append(simpegmt.SurveyMT.RxMT(simpeg.mkvc(loc,2).T,rxType))
|
||||
# Source list
|
||||
srcList =[]
|
||||
for freq in np.logspace(3,-1,5):
|
||||
for freq in np.logspace(3,-3,7):
|
||||
srcList.append(simpegmt.SurveyMT.srcMT(freq,rxList))
|
||||
# Survey MT
|
||||
survey = simpegmt.SurveyMT.SurveyMT(srcList)
|
||||
|
||||
## Setup the problem objec
|
||||
## Setup the problem object
|
||||
problem = simpegmt.ProblemMT.MTProblem(M)
|
||||
problem.pair(survey)
|
||||
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
from SimPEG import Survey, Problem, Utils, Models, np, sp, Solver as SimpegSolver
|
||||
from SimPEG import Survey, Problem, Utils, Models, np, sp, SolverLU as SimpegSolver
|
||||
from scipy.constants import mu_0
|
||||
from SurveyMT import SurveyMT, FieldsMT
|
||||
import multiprocessing
|
||||
@@ -121,7 +121,8 @@ class MTProblem(Problem.BaseProblem):
|
||||
# Store the fieldss
|
||||
F[Src, 'e_px'] = e[:,0]
|
||||
F[Src, 'e_py'] = e[:,1]
|
||||
b = self.mesh.edgeCurl * e
|
||||
# Note curl e = -iwb so b = -curl/iw
|
||||
b = -( self.mesh.edgeCurl * e )/( 1j*omega(freq) )
|
||||
F[Src, 'b_px'] = b[:,0]
|
||||
F[Src, 'b_py'] = b[:,1]
|
||||
return F
|
||||
@@ -150,7 +151,7 @@ class MTProblem(Problem.BaseProblem):
|
||||
sig = self.MeSigma
|
||||
C = self.mesh.edgeCurl
|
||||
|
||||
return C.T*mui*C + 1j*omega(freq)*sig
|
||||
return C.T*mui*C - 1j*omega(freq)*sig
|
||||
|
||||
def getAbg(self, freq):
|
||||
"""
|
||||
@@ -164,7 +165,7 @@ class MTProblem(Problem.BaseProblem):
|
||||
sigBG = self.MeSigmaBG
|
||||
C = self.mesh.edgeCurl
|
||||
|
||||
return C.T*mui*C + 1j*omega(freq)*sigBG
|
||||
return C.T*mui*C - 1j*omega(freq)*sigBG
|
||||
|
||||
def getADeriv(self, freq, u, v, adjoint=False):
|
||||
sig = self.curTModel
|
||||
@@ -193,7 +194,7 @@ class MTProblem(Problem.BaseProblem):
|
||||
eBG_bp = homo1DModelSource(self.mesh,freq,backSigma)
|
||||
Abg = self.getAbg(freq)
|
||||
|
||||
return Abg*eBG_bp
|
||||
return -Abg*eBG_bp
|
||||
|
||||
##################################################################
|
||||
# Inversion stuff
|
||||
|
||||
@@ -16,7 +16,8 @@ def homo1DModelSource(mesh,freq,m_back):
|
||||
from simpegMT.Utils import get1DEfields
|
||||
# Get a 1d solution for a halfspace background
|
||||
mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]]))
|
||||
e0_1d = get1DEfields(mesh1d,mesh.r(m_back,'CC','CC','M')[0,0,:],freq)
|
||||
# Note: Need to conjugate the source field to comply with orientations
|
||||
e0_1d = get1DEfields(mesh1d,mesh.r(m_back,'CC','CC','M')[0,0,:],freq).conj()
|
||||
# Setup x (east) polarization (_x)
|
||||
ex_px = np.zeros(mesh.vnEx,dtype=complex)
|
||||
ey_px = np.zeros((mesh.nEy,1),dtype=complex)
|
||||
@@ -24,18 +25,19 @@ def homo1DModelSource(mesh,freq,m_back):
|
||||
# Assign the source to ex_x
|
||||
for i in np.arange(mesh.vnEx[0]):
|
||||
for j in np.arange(mesh.vnEx[1]):
|
||||
ex_px[i,j,:] = e0_1d
|
||||
ex_px[1:-1,1:-1,1:-1] = 0
|
||||
eBG_px = np.vstack((simpeg.Utils.mkvc(mesh.r(ex_px,'Ex','Ex','V'),2),ey_px,ez_px))
|
||||
ex_px[i,j,:] = -e0_1d
|
||||
# ex_px[1:-1,1:-1,1:-1] = 0
|
||||
eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px,ez_px))
|
||||
# Setup y (north) polarization (_py)
|
||||
ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
|
||||
ey_py = np.zeros(mesh.vnEy, dtype='complex128')
|
||||
ez_py = np.zeros((mesh.nEz,1), dtype='complex128')
|
||||
# Assign the source to ey_py
|
||||
|
||||
for i in np.arange(mesh.vnEy[0]):
|
||||
for j in np.arange(mesh.vnEy[1]):
|
||||
ey_py[i,j,:] = e0_1d
|
||||
ey_py[1:-1,1:-1,1:-1] = 0
|
||||
ey_py[i,j,:] = e0_1d
|
||||
# ey_py[1:-1,1:-1,1:-1] = 0
|
||||
eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
|
||||
|
||||
# Return the electric fields
|
||||
|
||||
@@ -269,7 +269,7 @@ class DataMT(Survey.Data):
|
||||
from numpy.lib import recfunctions as recFunc
|
||||
dt = dtCP
|
||||
for uniFL in uniFLmarr:
|
||||
mTemp = simpeg.mkvc(rec2ndarr(mArrRec[np.ma.where(mArrRec[['freq','x','y','z']].data == np.array(uniFL))][impList]).sum(axis=0),2).T
|
||||
mTemp = mkvc(rec2ndarr(mArrRec[np.ma.where(mArrRec[['freq','x','y','z']].data == np.array(uniFL))][impList]).sum(axis=0),2).T
|
||||
compBlock = np.sum(mTemp.data.reshape((4,2))*np.array([[1,1j],[1,1j],[1,1j],[1,1j]]),axis=1).copy().view(dt[4::])
|
||||
dataBlock = mkvc(recFunc.merge_arrays((np.array(uniFL),compBlock),flatten=True),2).T
|
||||
try:
|
||||
|
||||
@@ -0,0 +1,48 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
import simpegMT as MT
|
||||
|
||||
TOL = 1e-6
|
||||
|
||||
def appResPhs(freq,z):
|
||||
app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
|
||||
app_phs = np.arctan2(-z.imag,z.real)*(180/np.pi)
|
||||
return app_res, app_phs
|
||||
|
||||
def appResNorm(sigmaHalf):
|
||||
nFreq = 26
|
||||
|
||||
m1d = Mesh.TensorMesh([[(100,5,1.5),(100.,10),(100,5,1.5)]], x0=['C'])
|
||||
sigma = np.zeros(m1d.nC) + sigmaHalf
|
||||
sigma[m1d.gridCC[:]>200] = 1e-8
|
||||
|
||||
# Calculate the analytic fields
|
||||
freqs = np.logspace(4,-4,nFreq)
|
||||
Z = []
|
||||
for freq in freqs:
|
||||
Ed, Eu, Hd, Hu = MT.Utils.getEHfields(m1d,sigma,freq,np.array([200]))
|
||||
Z.append((Ed + Eu)/(Hd + Hu))
|
||||
|
||||
Zarr = np.concatenate(Z)
|
||||
|
||||
app_r, app_p = appResPhs(freqs,Zarr)
|
||||
|
||||
return np.linalg.norm(np.abs(app_r - np.ones(nFreq)/sigmaHalf)) / np.log10(sigmaHalf)
|
||||
|
||||
|
||||
class TestAnalytics(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
pass
|
||||
def test_appRes2en1(self):self.assertLess(appResNorm(2e-1), TOL)
|
||||
def test_appRes2en2(self):self.assertLess(appResNorm(2e-2), TOL)
|
||||
def test_appRes2en3(self):self.assertLess(appResNorm(2e-3), TOL)
|
||||
def test_appRes2en4(self):self.assertLess(appResNorm(2e-4), TOL)
|
||||
def test_appRes2en5(self):self.assertLess(appResNorm(2e-5), TOL)
|
||||
def test_appRes2en6(self):self.assertLess(appResNorm(2e-6), TOL)
|
||||
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -1,103 +0,0 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
import simpegMT as MT
|
||||
|
||||
class FieldsTest(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30])
|
||||
x = np.linspace(5,10,3)
|
||||
XYZ = Utils.ndgrid(x,x,np.r_[0.])
|
||||
txLoc = np.r_[0,0,0.]
|
||||
rxList0 = MT.FDEM.RxFDEM(XYZ, 'exi')
|
||||
Tx0 = MT.FDEM.TxFDEM(txLoc, 'VMD', 3., [rxList0])
|
||||
rxList1 = MT.FDEM.RxFDEM(XYZ, 'bxi')
|
||||
Tx1 = MT.FDEM.TxFDEM(txLoc, 'VMD', 3., [rxList1])
|
||||
rxList2 = MT.FDEM.RxFDEM(XYZ, 'bxi')
|
||||
Tx2 = MT.FDEM.TxFDEM(txLoc, 'VMD', 2., [rxList2])
|
||||
rxList3 = MT.FDEM.RxFDEM(XYZ, 'bxi')
|
||||
Tx3 = MT.FDEM.TxFDEM(txLoc, 'VMD', 2., [rxList3])
|
||||
Tx4 = MT.FDEM.TxFDEM(txLoc, 'VMD', 1., [rxList0, rxList1, rxList2, rxList3])
|
||||
txList = [Tx0,Tx1,Tx2,Tx3,Tx4]
|
||||
survey = MT.FDEM.SurveyFDEM(txList)
|
||||
self.F = MT.FDEM.FieldsFDEM(mesh, survey)
|
||||
self.Tx0 = Tx0
|
||||
self.Tx1 = Tx1
|
||||
self.mesh = mesh
|
||||
self.XYZ = XYZ
|
||||
|
||||
def test_SetGet(self):
|
||||
F = self.F
|
||||
for freq in F.survey.freqs:
|
||||
nFreq = F.survey.nTxByFreq[freq]
|
||||
Txs = F.survey.getTransmitters(freq)
|
||||
e = np.random.rand(F.mesh.nE, nFreq)
|
||||
F[Txs, 'e'] = e
|
||||
b = np.random.rand(F.mesh.nF, nFreq)
|
||||
F[Txs, 'b'] = b
|
||||
if nFreq == 1:
|
||||
F[Txs, 'b'] = Utils.mkvc(b)
|
||||
if e.shape[1] == 1:
|
||||
e, b = Utils.mkvc(e), Utils.mkvc(b)
|
||||
self.assertTrue(np.all(F[Txs, 'e'] == e))
|
||||
self.assertTrue(np.all(F[Txs, 'b'] == b))
|
||||
F[Txs] = {'b':b,'e':e}
|
||||
self.assertTrue(np.all(F[Txs, 'e'] == e))
|
||||
self.assertTrue(np.all(F[Txs, 'b'] == b))
|
||||
|
||||
lastFreq = F[Txs]
|
||||
self.assertTrue(type(lastFreq) is dict)
|
||||
self.assertTrue(sorted([k for k in lastFreq]) == ['b','e'])
|
||||
self.assertTrue(np.all(lastFreq['b'] == b))
|
||||
self.assertTrue(np.all(lastFreq['e'] == e))
|
||||
|
||||
Tx_f3 = F.survey.getTransmitters(3.)
|
||||
self.assertTrue(F[Tx_f3,'b'].shape == (F.mesh.nF, 2))
|
||||
|
||||
b = np.random.rand(F.mesh.nF, 2)
|
||||
Tx_f0 = F.survey.getTransmitters(self.Tx0.freq)
|
||||
F[Tx_f0,'b'] = b
|
||||
self.assertTrue(F[self.Tx0]['b'].shape == (F.mesh.nF,))
|
||||
self.assertTrue(F[self.Tx0,'b'].shape == (F.mesh.nF,))
|
||||
self.assertTrue(np.all(F[self.Tx0,'b'] == b[:,0]))
|
||||
self.assertTrue(np.all(F[self.Tx1,'b'] == b[:,1]))
|
||||
|
||||
def test_assertions(self):
|
||||
freq = self.F.survey.freqs[0]
|
||||
Txs = self.F.survey.getTransmitters(freq)
|
||||
bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nTxByFreq[freq])
|
||||
def fun(): self.F[Txs, 'b'] = bWrongSize
|
||||
self.assertRaises(ValueError, fun)
|
||||
def fun(): self.F[-999.]
|
||||
self.assertRaises(KeyError, fun)
|
||||
def fun(): self.F['notRight']
|
||||
self.assertRaises(KeyError, fun)
|
||||
def fun(): self.F[Txs,'notThere']
|
||||
self.assertRaises(KeyError, fun)
|
||||
|
||||
def test_FieldProjections(self):
|
||||
F = self.F
|
||||
for freq in F.survey.freqs:
|
||||
nFreq = F.survey.nTxByFreq[freq]
|
||||
Txs = F.survey.getTransmitters(freq)
|
||||
e = np.random.rand(F.mesh.nE, nFreq)
|
||||
b = np.random.rand(F.mesh.nF, nFreq)
|
||||
F[Txs] = {'b':b,'e':e}
|
||||
|
||||
Txs = F.survey.getTransmitters(freq)
|
||||
for ii, tx in enumerate(Txs):
|
||||
for jj, rx in enumerate(tx.rxList):
|
||||
dat = rx.projectFields(tx, self.mesh, F)
|
||||
self.assertTrue(dat.dtype == float)
|
||||
fieldType = rx.projField
|
||||
u = {'b':b[:,ii], 'e': e[:,ii]}[fieldType]
|
||||
real_or_imag = rx.projComp
|
||||
u = getattr(u, real_or_imag)
|
||||
gloc = rx.projGLoc
|
||||
d = self.mesh.getInterpolationMat(self.XYZ, gloc)*u
|
||||
self.assertTrue(np.all(dat == d))
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -47,12 +47,19 @@ def getEHfields(m1d,sigma,freq,zd):
|
||||
UDp[:,lnr+1] = elamh.dot(Pjinv.dot(Pj1)).dot(UDp[:,lnr])
|
||||
|
||||
# Calculate the fields
|
||||
Ed = np.zeros((zd.size,),dtype=complex)
|
||||
Eu = np.zeros((zd.size,),dtype=complex)
|
||||
Hd = np.zeros((zd.size,),dtype=complex)
|
||||
Hu = np.zeros((zd.size,),dtype=complex)
|
||||
Ed = np.empty((zd.size,),dtype=complex)
|
||||
Eu = np.empty((zd.size,),dtype=complex)
|
||||
Hd = np.empty((zd.size,),dtype=complex)
|
||||
Hu = np.empty((zd.size,),dtype=complex)
|
||||
|
||||
# Loop over the layers and calculate the fields
|
||||
# In the halfspace below the mesh
|
||||
dup = m1d.vectorNx[0]
|
||||
dind = dup >= zd
|
||||
Ed[dind] = UDp[1,0]*np.exp(-1j*k[0]*(dup-zd[dind]))
|
||||
Eu[dind] = UDp[0,0]*np.exp(1j*k[0]*(dup-zd[dind]))
|
||||
Hd[dind] = (k[0]/(w*mu[0]))*UDp[1,0]*np.exp(-1j*k[0]*(dup-zd[dind]))
|
||||
Hu[dind] = -(k[0]/(w*mu[0]))*UDp[0,0]*np.exp(1j*k[0]*(dup-zd[dind]))
|
||||
for ki,mui,epsi,dlow,dup,Up,Dp in zip(k[1::],mu[1::],eps[1::],m1d.vectorNx[:-1],m1d.vectorNx[1::],UDp[0,1::],UDp[1,1::]):
|
||||
dind = np.logical_and(dup >= zd, zd > dlow)
|
||||
Ed[dind] = Dp*np.exp(-1j*ki*(dup-zd[dind]))
|
||||
|
||||
@@ -20,10 +20,12 @@ def get1DEfields(m1d,sigma,freq,sourceAmp=1.0):
|
||||
Aio = A[1:-1,[0,-1]]
|
||||
|
||||
# Set the boundary conditions
|
||||
Ed_low, Eu_low, Hd_low, Hu_low = getEHfields(m1d,sigma,freq,np.array([m1d.vectorNx[0]]))
|
||||
Etot_low = Ed_low + Eu_low
|
||||
Ed, Eu, Hd, Hu = getEHfields(m1d,sigma,freq,m1d.vectorNx)
|
||||
Etot = Ed + Eu
|
||||
if sourceAmp is not None:
|
||||
Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
|
||||
## Note: need to use conjugate of the analytic solution. It is derived with e^iwt
|
||||
bc = np.r_[Etot_low.conj(),sourceAmp]
|
||||
bc = np.r_[Etot[0],Etot[-1]]
|
||||
# The right hand side
|
||||
rhs = -Aio*bc
|
||||
# Solve the system
|
||||
|
||||
Reference in New Issue
Block a user