diff --git a/notebooks/MT 1D code_testLayered.ipynb b/notebooks/MT 1D code_testLayered.ipynb index 0170b5a9..7cf83993 100644 --- a/notebooks/MT 1D code_testLayered.ipynb +++ b/notebooks/MT 1D code_testLayered.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 16, + "execution_count": 6, "metadata": { "collapsed": false }, @@ -22,7 +22,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 7, "metadata": { "collapsed": false }, @@ -39,7 +39,7 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 8, "metadata": { "collapsed": false }, @@ -62,7 +62,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 9, "metadata": { "collapsed": false }, @@ -70,10 +70,10 @@ { "data": { "text/plain": [ - "[]" + "[]" ] }, - "execution_count": 19, + "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, @@ -186,7 +186,7 @@ "ZSzcklQZC7ckVcbCLUmVsXBLUmX+H3zfa/Qjde8mAAAAAElFTkSuQmCC\n" ], "text/plain": [ - "" + "" ] }, "metadata": {}, @@ -199,7 +199,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 10, "metadata": { "collapsed": false }, @@ -217,7 +217,7 @@ " 7888.671875 , 11733.0078125 , 17499.51171875])" ] }, - "execution_count": 20, + "execution_count": 10, "metadata": {}, "output_type": "execute_result" } @@ -228,7 +228,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 28, "metadata": { "collapsed": false }, @@ -249,7 +249,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 29, "metadata": { "collapsed": false }, @@ -260,7 +260,7 @@ "(12836609.712174654+24128916.329733729j)" ] }, - "execution_count": 22, + "execution_count": 29, "metadata": {}, "output_type": "execute_result" } @@ -271,7 +271,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 30, "metadata": { "collapsed": false }, @@ -343,7 +343,7 @@ " 1.00000000e+00 -0.00000000e+00j]])" ] }, - "execution_count": 23, + "execution_count": 30, "metadata": {}, "output_type": "execute_result" } @@ -354,7 +354,7 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 31, "metadata": { "collapsed": false }, @@ -362,11 +362,11 @@ { "data": { "text/plain": [ - "[,\n", - " ]" + "[,\n", + " ]" ] }, - "execution_count": 24, + "execution_count": 31, "metadata": {}, "output_type": "execute_result" }, @@ -562,7 +562,7 @@ "ThxmZlYrThxmZlYrThxmZlYr/x/lWzp+eZujggAAAABJRU5ErkJggg==\n" ], "text/plain": [ - "" + "" ] }, "metadata": {}, @@ -576,7 +576,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 32, "metadata": { "collapsed": false }, @@ -584,10 +584,10 @@ { "data": { "text/plain": [ - "[]" + "[]" ] }, - "execution_count": 25, + "execution_count": 32, "metadata": {}, "output_type": "execute_result" }, @@ -788,7 +788,7 @@ "ZlYRFw4zM6uIC4eZmVXEhcPMzCriwmFmZhVx4TAzs4r8f6jZeSiIHLknAAAAAElFTkSuQmCC\n" ], "text/plain": [ - "" + "" ] }, "metadata": {}, @@ -801,7 +801,7 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 33, "metadata": { "collapsed": false }, @@ -809,11 +809,11 @@ { "data": { "text/plain": [ - "[,\n", - " ]" + "[,\n", + " ]" ] }, - "execution_count": 26, + "execution_count": 33, "metadata": {}, "output_type": "execute_result" }, @@ -1048,7 +1048,7 @@ "c865pPGk4pxzLmk8qTjnnEsaTyrOOeeS5v8Bn8/HW3Alj/EAAAAASUVORK5CYII=\n" ], "text/plain": [ - "" + "" ] }, "metadata": {}, @@ -1061,7 +1061,7 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 34, "metadata": { "collapsed": false }, @@ -1077,10 +1077,10 @@ { "data": { "text/plain": [ - "[]" + "[]" ] }, - "execution_count": 27, + "execution_count": 34, "metadata": {}, "output_type": "execute_result" }, @@ -1252,7 +1252,7 @@ "rkJggg==\n" ], "text/plain": [ - "" + "" ] }, "metadata": {}, @@ -1265,7 +1265,7 @@ }, { "cell_type": "code", - "execution_count": 28, + "execution_count": 35, "metadata": { "collapsed": false }, @@ -1273,11 +1273,11 @@ { "data": { "text/plain": [ - "[,\n", - " ]" + "[,\n", + " ]" ] }, - "execution_count": 28, + "execution_count": 35, "metadata": {}, "output_type": "execute_result" }, @@ -1485,7 +1485,7 @@ "5P8DEIvS62qo2M0AAAAASUVORK5CYII=\n" ], "text/plain": [ - "" + "" ] }, "metadata": {}, @@ -1498,7 +1498,7 @@ }, { "cell_type": "code", - "execution_count": 29, + "execution_count": 36, "metadata": { "collapsed": false }, @@ -1506,10 +1506,10 @@ { "data": { "text/plain": [ - "[]" + "[]" ] }, - "execution_count": 29, + "execution_count": 36, "metadata": {}, "output_type": "execute_result" }, @@ -1689,7 +1689,7 @@ "rkJggg==\n" ], "text/plain": [ - "" + "" ] }, "metadata": {}, @@ -1702,7 +1702,7 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 37, "metadata": { "collapsed": false }, @@ -1729,7 +1729,7 @@ }, { "cell_type": "code", - "execution_count": 31, + "execution_count": 38, "metadata": { "collapsed": false }, @@ -1740,7 +1740,7 @@ "(30,)" ] }, - "execution_count": 31, + "execution_count": 38, "metadata": {}, "output_type": "execute_result" } @@ -1751,7 +1751,7 @@ }, { "cell_type": "code", - "execution_count": 32, + "execution_count": 25, "metadata": { "collapsed": false }, @@ -1762,7 +1762,7 @@ "31" ] }, - "execution_count": 32, + "execution_count": 25, "metadata": {}, "output_type": "execute_result" } @@ -1773,7 +1773,7 @@ }, { "cell_type": "code", - "execution_count": 33, + "execution_count": 26, "metadata": { "collapsed": false }, @@ -1781,11 +1781,11 @@ { "data": { "text/plain": [ - "[,\n", - " ]" + "[,\n", + " ]" ] }, - "execution_count": 33, + "execution_count": 26, "metadata": {}, "output_type": "execute_result" }, @@ -1985,7 +1985,7 @@ "AABJRU5ErkJggg==\n" ], "text/plain": [ - "" + "" ] }, "metadata": {}, @@ -1998,7 +1998,7 @@ }, { "cell_type": "code", - "execution_count": 34, + "execution_count": 27, "metadata": { "collapsed": false }, @@ -2016,7 +2016,7 @@ " 9810.83984375, 14616.25976562])" ] }, - "execution_count": 34, + "execution_count": 27, "metadata": {}, "output_type": "execute_result" } @@ -2049,18 +2049,6 @@ "display_name": "Python 2", "language": "python", "name": "python2" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 2 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.9" } }, "nbformat": 4, diff --git a/notebooks/MT Script-3D_twoHalfspace-simpegMT.ipynb b/notebooks/MT Script-3D_twoHalfspace-simpegMT.ipynb index 4230295b..36847d04 100644 --- a/notebooks/MT Script-3D_twoHalfspace-simpegMT.ipynb +++ b/notebooks/MT Script-3D_twoHalfspace-simpegMT.ipynb @@ -2813,10 +2813,7 @@ ], "source": [ "ind = np.where(np.sum(np.power(rx_loc - np.array([-3000,0,elev]),2),axis=1)< 5)\n", - "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", + "m\n", "print appResPhs(freq,zyx[ind])\n", "print appResPhs(freq,zxy[ind])" ] diff --git a/notebooks/Test ProblemMT1D.eFrom_TotalField.ipynb b/notebooks/Test ProblemMT1D.eFrom_TotalField.ipynb new file mode 100644 index 00000000..d4061ba1 --- /dev/null +++ b/notebooks/Test ProblemMT1D.eFrom_TotalField.ipynb @@ -0,0 +1,263 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# Notebook to test 1D code" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "import SimPEG as simpeg\n", + "import simpegMT as simpegmt" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "ct = 10\n", + "m1d = simpeg.Mesh.TensorMesh([[(ct,20,-1.5),(ct,100),(ct,20,1.5)]], x0=['C'])\n", + "sigma = np.zeros(m1d.nC) + 2e-3\n", + "sigma[m1d.gridCC[:]>200] = 1e-8" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Make the rx and src\n", + "freqs = np.logspace(3,-3,31)\n", + "rxList = []\n", + "for rxType in ['zxyr','zxyi']:\n", + " rxList.append(simpegmt.SurveyMT.RxMT(simpeg.mkvc(np.array([405]),2).T,rxType))\n", + "# Source list\n", + "srcList =[]\n", + "for freq in freqs: \n", + " srcList.append(simpegmt.SurveyMT.srcMT(freq,rxList)) \n", + "survey = simpegmt.SurveyMT.SurveyMT(srcList)\n", + "problem = simpegmt.ProblemMT1D.eForm_TotalField(m1d)\n", + "problem.pair(survey)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false, + "scrolled": true + }, + "outputs": [], + "source": [ + "fields = problem.fields(sigma)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 1.52442595e-125 -7.20764170e-126j,\n", + " 7.74672164e-101 +1.76864736e-100j,\n", + " -2.07557099e-080 -3.74608600e-081j, ...,\n", + " 3.10393307e-001 -2.86962945e-001j,\n", + " 4.10341179e-001 -2.87349774e-001j,\n", + " 5.05105758e-001 -2.74745005e-001j],\n", + " [ 8.67440346e-030 +6.35706181e-030j,\n", + " -1.53155881e-027 +3.75444030e-027j,\n", + " -1.08396611e-024 -4.98500246e-026j, ...,\n", + " -4.05147339e-001 +2.76576534e-001j,\n", + " -4.96380502e-001 +2.67133906e-001j,\n", + " -5.80363237e-001 +2.48768467e-001j],\n", + " [ -6.16527128e-026 +8.41573172e-026j,\n", + " -2.29821697e-023 -9.36787887e-024j,\n", + " 1.90707675e-022 -4.18621866e-021j, ...,\n", + " -4.75055964e-001 +2.59780137e-001j,\n", + " -5.57847068e-001 +2.46025119e-001j,\n", + " -6.32948088e-001 +2.25820916e-001j],\n", + " ..., \n", + " [ -4.42815267e-001 +4.54129730e-002j,\n", + " -4.45450945e-001 +2.94352343e-002j,\n", + " -4.46746292e-001 +1.94083017e-002j, ...,\n", + " -7.96037136e-001 +1.08182287e-001j,\n", + " -8.29064328e-001 +1.00291076e-001j,\n", + " -8.58609530e-001 +9.06236881e-002j],\n", + " [ -6.64333571e-001 +4.65810415e-002j,\n", + " -6.66721482e-001 +2.99034482e-002j,\n", + " -6.67823222e-001 +1.93821193e-002j, ...,\n", + " -8.77622277e-001 +6.49094362e-002j,\n", + " -8.97438594e-001 +6.01746869e-002j,\n", + " -9.15165716e-001 +5.43742395e-002j],\n", + " [ 1.00000000e+000 +0.00000000e+000j,\n", + " 1.00000000e+000 +0.00000000e+000j,\n", + " 1.00000000e+000 +0.00000000e+000j, ...,\n", + " 1.00000000e+000 +0.00000000e+000j,\n", + " 1.00000000e+000 +0.00000000e+000j,\n", + " 1.00000000e+000 +0.00000000e+000j]])" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "fields[:,'e_1d']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "m1d.nN" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "self = problem\n", + "Mmui = self.MfMui\n", + "Msig = self.mesh.getFaceInnerProduct(self.curModel)\n", + "C = self.mesh.nodalGrad" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "self" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "Mmui" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "self.Me[1,1]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "self.mesh.vol" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "m1d.h" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "m1d.gridCC" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.9" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/simpegMT/BaseMT.py b/simpegMT/BaseMT.py new file mode 100644 index 00000000..fbc87187 --- /dev/null +++ b/simpegMT/BaseMT.py @@ -0,0 +1,25 @@ +from simpegEM.FDEM import BaseFDEMProblem +from SurveyMT import SurveyMT +from DataMT import DataMT +from FieldsMT import FieldsMT +from SimPEG import SolverLU as SimpegSolver + +class BaseMTProblem(BaseFDEMProblem): + + def __init__(self, mesh, **kwargs): + BaseFDEMProblem.__init__(self, mesh, **kwargs) + + # Set the default pairs of the problem + surveyPair = SurveyMT + dataPair = DataMT + fieldsPair = FieldsMT + + # Set the solver + Solver = SimpegSolver + solverOpts = {} + + verbose = False + # Notes: + # Use the forward and devs from BaseFDEMProblem + # Might need to add more stuff here. + diff --git a/simpegMT/DataMT.py b/simpegMT/DataMT.py new file mode 100644 index 00000000..6ef20939 --- /dev/null +++ b/simpegMT/DataMT.py @@ -0,0 +1,70 @@ +from SimPEG import Survey, Utils, Problem, np, sp, mkvc +from scipy.constants import mu_0 +import sys +from numpy.lib import recfunctions as recFunc + +############ +### Data ### +############ +class DataMT(Survey.Data): + ''' + Data class for MTdata + + :param SimPEG survey object survey: + :param v vector with data + + ''' + def __init__(self, survey, v=None): + # Pass the variables to the "parent" method + Survey.Data.__init__(self, survey, v) + + def toRecArray(self,returnType='RealImag'): + ''' + Function that returns a numpy.recarray for a SimpegMT impedance data object. + + :param str returnType: Switches between returning a rec array where the impedance is split to real and imaginary ('RealImag') or is a complex ('Complex') + + ''' + + def rec2ndarr(x,dt=float): + return x.view((dt, len(x.dtype.names))) + # Define the record fields + dtRI = [('freq',float),('x',float),('y',float),('z',float),('zxxr',float),('zxxi',float),('zxyr',float),('zxyi',float),('zyxr',float),('zyxi',float),('zyyr',float),('zyyi',float)] + dtCP = [('freq',float),('x',float),('y',float),('z',float),('zxx',complex),('zxy',complex),('zyx',complex),('zyy',complex)] + impList = ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi'] + for src in self.survey.srcList: + # Temp array for all the receivers of the source. + # Note: needs to be written more generally, using diffterent rxTypes and not all the data at the locaitons + # Assume the same locs for all RX + locs = src.rxList[0].locs + tArrRec = np.concatenate((src.freq*np.ones((locs.shape[0],1)),locs,np.nan*np.ones((locs.shape[0],8))),axis=1).view(dtRI) + # np.array([(src.freq,rx.locs[0,0],rx.locs[0,1],rx.locs[0,2],np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ) for rx in src.rxList],dtype=dtRI) + # Get the type and the value for the DataMT object as a list + typeList = [[rx.rxType,self[src,rx]] for rx in src.rxList] + # Insert the values to the temp array + for nr,(key,val) in enumerate(typeList): + tArrRec[key] = mkvc(val,2) + # Masked array + mArrRec = np.ma.MaskedArray(rec2ndarr(tArrRec),mask=np.isnan(rec2ndarr(tArrRec))).view(dtype=tArrRec.dtype) + # Unique freq and loc of the masked array + uniFLmarr = np.unique(mArrRec[['freq','x','y','z']]) + + try: + outTemp = recFunc.stack_arrays((outTemp,mArrRec)) + #outTemp = np.concatenate((outTemp,dataBlock),axis=0) + except NameError as e: + outTemp = mArrRec + + if 'RealImag' in returnType: + outArr = outTemp + if 'Complex' in returnType: + # Add the real and imaginary to a complex number + + outArr = np.empty(outTemp.shape,dtype=dtCP) + for comp in ['freq','x','y','z']: + outArr[comp] = outTemp[comp].copy() + for comp in ['zxx','zxy','zyx','zyy']: + outArr[comp] = outTemp[comp+'r'].copy() + 1j*outTemp[comp+'i'].copy() + + # Return + return outArr diff --git a/simpegMT/FieldsMT.py b/simpegMT/FieldsMT.py new file mode 100644 index 00000000..9487d122 --- /dev/null +++ b/simpegMT/FieldsMT.py @@ -0,0 +1,12 @@ +from SimPEG import Survey, Utils, Problem, np, sp, mkvc +from scipy.constants import mu_0 +import sys +from numpy.lib import recfunctions as recFunc + +############## +### Fields ### +############## +class FieldsMT(Problem.Fields): + """Fancy Field Storage for a MT survey.""" + knownFields = {'b_px': 'F','b_py': 'F', 'e_px': 'E','e_py': 'E','b_1d':'E','e_1d':'F'} + dtype = complex diff --git a/simpegMT/ProblemMT1D/Problems.py b/simpegMT/ProblemMT1D/Problems.py new file mode 100644 index 00000000..35fee655 --- /dev/null +++ b/simpegMT/ProblemMT1D/Problems.py @@ -0,0 +1,123 @@ +from simpegEM.Utils.EMUtils import omega +from SimPEG import mkvc +from scipy.constants import mu_0 +from simpegMT.BaseMT import BaseMTProblem +from simpegMT.SurveyMT import SurveyMT +from simpegMT.FieldsMT import FieldsMT +from simpegMT.DataMT import DataMT +from simpegMT.Utils.MT1Danalytic import getEHfields +import numpy as np +import multiprocessing, sys, time + + +class eForm_TotalField(BaseMTProblem): + """ + A MT problem solving a e formulation and a primary/secondary fields decompostion. + + Solves the equation: + + + """ + + # From FDEMproblem: Used to project the fields. Currently not used for MTproblem. + _fieldType = 'e' + _eqLocs = 'FE' + + + def __init__(self, mesh, **kwargs): + BaseMTProblem.__init__(self, mesh, **kwargs) + + def getA(self, freq, full=False): + """ + Function to get the A matrix. + + :param float freq: Frequency + :param logic full: Return full A or the inner part + :rtype: scipy.sparse.csr_matrix + :return: A + """ + Mmui = self.mesh.getEdgeInnerProduct(1.0/mu_0) + Msig = self.mesh.getFaceInnerProduct(self.curModel) + C = self.mesh.nodalGrad + # Make A + A = C.T*Mmui*C + 1j*omega(freq)*Msig + # Either return full or only the inner part of A + if full: + return A + else: + return A[1:-1,1:-1] + + def getADeriv(self, freq, u, v, adjoint=False): + sig = self.curTModel + dsig_dm = self.curTModelDeriv + dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=u) + + if adjoint: + return 1j * omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * v ) ) + + return 1j * omega(freq) * ( dMe_dsig * ( dsig_dm * v ) ) + + def getRHS(self, freq): + """ + Function to return the right hand side for the system. + :param float freq: Frequency + :rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2) + :return: RHS for both polarizations, primary fields + """ + # Get sources for the frequency + # NOTE: Need to use the source information, doesn't really apply in 1D + src = self.survey.getSources(freq) + # Get the full A + A = self.getA(freq,full=True) + # Define the outer part of the solution matrix + Aio = A[1:-1,[0,-1]] + Ed, Eu, Hd, Hu = getEHfields(self.mesh,self.curModel,freq,self.mesh.vectorNx) + Etot = (Ed + Eu) + sourceAmp = 1.0 + Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top + ## Note: The analytic solution is derived with e^iwt + eBC = np.r_[Etot[0],Etot[-1]] + # The right hand side + + return Aio*eBC, eBC + + def getRHSderiv(self, freq, backSigma, u, v, adjoint=False): + raise NotImplementedError('getRHSDeriv not implemented yet') + return None + + def fields(self, m): + ''' + Function to calculate all the fields for the model m. + + :param np.ndarray (nC,) m: Conductivity model + :param np.ndarray (nC,) m_back: Background conductivity model + ''' + self.curModel = m + # RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields + + F = FieldsMT(self.mesh, self.survey) + for freq in self.survey.freqs: + if self.verbose: + startTime = time.time() + print 'Starting work for {:.3e}'.format(freq) + sys.stdout.flush() + A = self.getA(freq) + rhs, e_o = self.getRHS(freq) + Ainv = self.Solver(A, **self.solverOpts) + e_i = Ainv * rhs + e = mkvc(np.r_[e_o[0], e_i, e_o[1]],2) + # Store the fields + Src = self.survey.getSources(freq) + # Store the fields + # NOTE: only store + F[Src, 'e_1d'] = e + # F[Src, 'e_py'] = 0*e[:,0] + # Note curl e = -iwb so b = -curl e /iw + b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) ) + # F[Src, 'b_px'] = 0*b[:,0] + F[Src, 'b_1d'] = b + if self.verbose: + print 'Ran for {:f} seconds'.format(time.time()-startTime) + sys.stdout.flush() + return F + \ No newline at end of file diff --git a/simpegMT/ProblemMT1D/__init__.py b/simpegMT/ProblemMT1D/__init__.py new file mode 100644 index 00000000..767c81df --- /dev/null +++ b/simpegMT/ProblemMT1D/__init__.py @@ -0,0 +1 @@ +from Problems import eForm_TotalField \ No newline at end of file diff --git a/simpegMT/ProblemMT2D/Problems.py b/simpegMT/ProblemMT2D/Problems.py new file mode 100644 index 00000000..e69de29b diff --git a/simpegMT/ProblemMT2D/__init__.py b/simpegMT/ProblemMT2D/__init__.py new file mode 100644 index 00000000..fc80254b --- /dev/null +++ b/simpegMT/ProblemMT2D/__init__.py @@ -0,0 +1 @@ +pass \ No newline at end of file diff --git a/simpegMT/ProblemMT.py b/simpegMT/ProblemMT3D/Problems.py similarity index 91% rename from simpegMT/ProblemMT.py rename to simpegMT/ProblemMT3D/Problems.py index 34ed5fc9..efec3d64 100644 --- a/simpegMT/ProblemMT.py +++ b/simpegMT/ProblemMT3D/Problems.py @@ -1,41 +1,25 @@ from SimPEG import Survey, Problem, Utils, Models, np, sp, SolverLU as SimpegSolver -from simpegEM.FDEM import BaseFDEMProblem from simpegEM.Utils.EMUtils import omega from scipy.constants import mu_0 -from SurveyMT import SurveyMT, FieldsMT, DataMT +from simpegMT.BaseMT import BaseMTProblem +from simpegMT.SurveyMT import SurveyMT +from simpegMT.FieldsMT import FieldsMT +from simpegMT.DataMT import DataMT import multiprocessing, sys, time -class BaseMTProblem(BaseFDEMProblem): - - def __init__(self, mesh, **kwargs): - BaseFDEMProblem.__init__(self, mesh, **kwargs) - - - surveyPair = SurveyMT - dataPair = DataMT - - Solver = SimpegSolver - solverOpts = {} - - verbose = False - # Notes: - # Use the forward and devs from BaseFDEMProblem - # Might need to add more stuff here. - - - -class ProblemMT_eForm_ps(BaseMTProblem): +class eForm_ps(BaseMTProblem): """ A MT problem solving a e formulation and a primary/secondary fields decompostion. Solves the equation """ + # From FDEMproblem: Used to project the fields. Currently not used for MTproblem. _fieldType = 'e' _eqLocs = 'FE' - fieldsPair = FieldsMT + # Set new properties # Background model @@ -96,8 +80,8 @@ class ProblemMT_eForm_ps(BaseMTProblem): Function to return the right hand side for the system. :param float freq: Frequency :param numpy.ndarray (nC,) backSigma: Background conductivity model - :rtype: numpy.ndarray (nE, 2) - :return: one RHS for both polarizations + :rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2) + :return: RHS for both polarizations, primary fields """ # Get sources for the frequency src = self.survey.getSources(freq) @@ -151,9 +135,9 @@ class ProblemMT_eForm_ps(BaseMTProblem): sys.stdout.flush() return F -class ProblemMT_eForm_Tp(BaseMTProblem): +class eForm_Tp(BaseMTProblem): """ - A MT problem solving a e formulation and a primary/secondary fields decompostion. + A MT problem solving a e formulation and a total/primary fields decompostion. Solves the equation """ @@ -263,7 +247,7 @@ class ProblemMT_eForm_Tp(BaseMTProblem): rhs, e_p = self.getRHS(freq,m_back) Ainv = self.Solver(A, **self.solverOpts) e_s = Ainv * rhs - e = e_p + e_s + e = e_s # Store the fields Src = self.survey.getSources(freq) # Store the fieldss diff --git a/simpegMT/ProblemMT3D/__init__.py b/simpegMT/ProblemMT3D/__init__.py new file mode 100644 index 00000000..4a98814c --- /dev/null +++ b/simpegMT/ProblemMT3D/__init__.py @@ -0,0 +1 @@ +from Problems import eForm_ps, eForm_Tp \ No newline at end of file diff --git a/simpegMT/Sources/backgroundModelSources.py b/simpegMT/Sources/backgroundModelSources.py index c878baf5..74d21964 100644 --- a/simpegMT/Sources/backgroundModelSources.py +++ b/simpegMT/Sources/backgroundModelSources.py @@ -16,8 +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]])) - # 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() + # Note: Everything is using e^iwt + e0_1d = get1DEfields(mesh1d,mesh.r(m_back,'CC','CC','M')[0,0,:],freq) # Setup x (east) polarization (_x) ex_px = np.zeros(mesh.vnEx,dtype=complex) ey_px = np.zeros((mesh.nEy,1),dtype=complex) @@ -26,7 +26,6 @@ def homo1DModelSource(mesh,freq,m_back): 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(ex_px,2),ey_px,ez_px)) # Setup y (north) polarization (_py) ex_py = np.zeros((mesh.nEx,1), dtype='complex128') diff --git a/simpegMT/SurveyMT.py b/simpegMT/SurveyMT.py index 47d6103c..96bc7723 100644 --- a/simpegMT/SurveyMT.py +++ b/simpegMT/SurveyMT.py @@ -2,19 +2,28 @@ from SimPEG import Survey, Utils, Problem, np, sp, mkvc from scipy.constants import mu_0 import sys from numpy.lib import recfunctions as recFunc +from DataMT import DataMT +################# +### Receivers ### +################# class RxMT(Survey.BaseRx): knownRxTypes = { - 'zxxr':[['e', 'Ex'],['b','Fx'], 'real'], - 'zxyr':[['e', 'Ex'],['b','Fy'], 'real'], - 'zyxr':[['e', 'Ey'],['b','Fx'], 'real'], - 'zyyr':[['e', 'Ey'],['b','Fy'], 'real'], - 'zxxi':[['e', 'Ex'],['b','Fx'], 'imag'], - 'zxyi':[['e', 'Ex'],['b','Fy'], 'imag'], - 'zyxi':[['e', 'Ey'],['b','Fx'], 'imag'], - 'zyyi':[['e', 'Ey'],['b','Fy'], 'imag'], - + # 3D impedance + 'zxxr':['Z3D', 'real'], + 'zxyr':['Z3D', 'real'], + 'zyxr':['Z3D', 'real'], + 'zyyr':['Z3D', 'real'], + 'zxxi':['Z3D', 'imag'], + 'zxyi':['Z3D', 'imag'], + 'zyxi':['Z3D', 'imag'], + 'zyyi':['Z3D', 'imag'], + # 2D impedance + # TODO: + # 1D impedance + 'z1dr':['Z1D', 'real'], + 'z1di':['Z1D', 'imag'] #TODO: Add tipper fractions as well. Bz/B(x|y) # 'exi':['e', 'Ex', 'imag'], # 'eyi':['e', 'Ey', 'imag'], @@ -32,7 +41,7 @@ class RxMT(Survey.BaseRx): Survey.BaseRx.__init__(self, locs, rxType) @property - def projField(self,fracPos): + def projField(self): """ Field Type projection (e.g. e b ...) :param str fracPos: Position of the field in the data ratio @@ -46,7 +55,7 @@ class RxMT(Survey.BaseRx): raise Exception('{s} is an unknown option. Use numerator or denominator.') @property - def projGLoc(self,fracPos): + def projGLoc(self): """ Grid Location projection (e.g. Ex Fy ...) :param str fracPos: Position of the field in the data ratio @@ -58,52 +67,58 @@ class RxMT(Survey.BaseRx): return self.knownRxTypes[self.rxType][0][1] else: raise Exception('{s} is an unknown option. Use numerator or denominator.') + @property + def projType(self): + """ + Receiver type for projection. + """ + return self.knownRxTypes[self.rxType][0] + @property def projComp(self): """Component projection (real/imag)""" - return self.knownRxTypes[self.rxType][2] + return self.knownRxTypes[self.rxType][1] def projectFields(self, src, mesh, u): ''' Project the fields and return the ''' - # Get the projection - # Pex = self.getP(mesh,'Ex') - # Pey = self.getP(mesh,'Ey') - # Pbx = self.getP(mesh,'Fx') - # Pby = self.getP(mesh,'Fy') - Pex = mesh.getInterpolationMat(self.locs,'Ex') - Pey = mesh.getInterpolationMat(self.locs,'Ey') - Pbx = mesh.getInterpolationMat(self.locs,'Fx') - Pby = mesh.getInterpolationMat(self.locs,'Fy') - # Get the fields at location - ex_px = Pex*u[src,'e_px'] - ey_px = Pey*u[src,'e_px'] - ex_py = Pex*u[src,'e_py'] - ey_py = Pey*u[src,'e_py'] - hx_px = Pbx*u[src,'b_px']/mu_0 - hy_px = Pby*u[src,'b_px']/mu_0 - hx_py = Pbx*u[src,'b_py']/mu_0 - hy_py = Pby*u[src,'b_py']/mu_0 - if 'zxx' in self.rxType: - f_part_complex = (ex_px*hy_py - ex_py*hy_px)/(hx_px*hy_py - hx_py*hy_px) - elif 'zxy' in self.rxType: - f_part_complex = (-ex_px*hx_py + ex_py*hx_px)/(hx_px*hy_py - hx_py*hy_px) - elif 'zyx' in self.rxType: - f_part_complex = (ey_px*hy_py - ey_py*hy_px)/(hx_px*hy_py - hx_py*hy_px) - elif 'zyy' in self.rxType: - f_part_complex = (-ey_px*hx_py + ey_py*hx_px)/(hx_px*hy_py - hx_py*hy_px) - - # P_num = self.getP(mesh,self.projGLoc('numerator')) - # u_num_complex = u[src, self.projField('numerator')] - # # Get the denominator information - # P_den = self.getP(mesh,self.projGLoc('denominator')) - # u_den_complex = u[src, self.projField('denominator')] - # # Calculate the fraction - # f_part_complex = (P_num*u_num_complex)/(P_den*u_den_complex) - # get the real or imag component + 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 + elif self.projType is 'Z3D': + # Get the projection + Pex = mesh.getInterpolationMat(self.locs,'Ex') + Pey = mesh.getInterpolationMat(self.locs,'Ey') + Pbx = mesh.getInterpolationMat(self.locs,'Fx') + Pby = mesh.getInterpolationMat(self.locs,'Fy') + # Get the fields at location + # px: x-polaration and py: y-polaration. + ex_px = Pex*u[src,'e_px'] + ey_px = Pey*u[src,'e_px'] + ex_py = Pex*u[src,'e_py'] + ey_py = Pey*u[src,'e_py'] + hx_px = Pbx*u[src,'b_px']/mu_0 + hy_px = Pby*u[src,'b_px']/mu_0 + hx_py = Pbx*u[src,'b_py']/mu_0 + hy_py = Pby*u[src,'b_py']/mu_0 + # Make the complex data + if 'zxx' in self.rxType: + f_part_complex = (ex_px*hy_py - ex_py*hy_px)/(hx_px*hy_py - hx_py*hy_px) + elif 'zxy' in self.rxType: + f_part_complex = (-ex_px*hx_py + ex_py*hx_px)/(hx_px*hy_py - hx_py*hy_px) + elif 'zyx' in self.rxType: + f_part_complex = (ey_px*hy_py - ey_py*hy_px)/(hx_px*hy_py - hx_py*hy_px) + elif 'zyy' in self.rxType: + f_part_complex = (-ey_px*hx_py + ey_py*hx_px)/(hx_px*hy_py - hx_py*hy_px) + else: + NotImplementedError('Projection of {:s} receiver type is not implemented.'.format(self.rxType)) + # Get the real or imag component real_or_imag = self.projComp f_part = getattr(f_part_complex, real_or_imag) return f_part @@ -130,6 +145,10 @@ class RxMT(Survey.BaseRx): # Note: Might need to add tests to make sure that both polarization have the same rxList. + +############### +### Sources ### +############### class srcMT(Survey.BaseSrc): ''' Sources for the MT problem. @@ -151,13 +170,9 @@ class srcMT(Survey.BaseSrc): Survey.BaseSrc.__init__(self, None, srcPol, rxList) - -class FieldsMT(Problem.Fields): - """Fancy Field Storage for a MT survey.""" - knownFields = {'b_px': 'F','b_py': 'F', 'e_px': 'E','e_py': 'E'} - dtype = complex - - +############## +### Survey ### +############## class SurveyMT(Survey.BaseSurvey): """ Survey class for MT. Contains all the sources associated with the survey. @@ -210,65 +225,3 @@ class SurveyMT(Survey.BaseSurvey): def projectFieldsDeriv(self, u): raise Exception('Use Transmitters to project fields deriv.') -class DataMT(Survey.Data): - ''' - Data class for MTdata - - :param SimPEG survey object survey: - :param v vector with data - - ''' - def __init__(self, survey, v=None): - # Pass the variables to the "parent" method - Survey.Data.__init__(self, survey, v) - - def toRecArray(self,returnType='RealImag'): - ''' - Function that returns a numpy.recarray for a SimpegMT impedance data object. - - :param str returnType: Switches between returning a rec array where the impedance is split to real and imaginary ('RealImag') or is a complex ('Complex') - - ''' - - def rec2ndarr(x,dt=float): - return x.view((dt, len(x.dtype.names))) - # Define the record fields - dtRI = [('freq',float),('x',float),('y',float),('z',float),('zxxr',float),('zxxi',float),('zxyr',float),('zxyi',float),('zyxr',float),('zyxi',float),('zyyr',float),('zyyi',float)] - dtCP = [('freq',float),('x',float),('y',float),('z',float),('zxx',complex),('zxy',complex),('zyx',complex),('zyy',complex)] - impList = ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi'] - for src in self.survey.srcList: - # Temp array for all the receivers of the source. - # Note: needs to be written more generally, using diffterent rxTypes and not all the data at the locaitons - # Assume the same locs for all RX - locs = src.rxList[0].locs - tArrRec = np.concatenate((src.freq*np.ones((locs.shape[0],1)),locs,np.nan*np.ones((locs.shape[0],8))),axis=1).view(dtRI) - # np.array([(src.freq,rx.locs[0,0],rx.locs[0,1],rx.locs[0,2],np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ) for rx in src.rxList],dtype=dtRI) - # Get the type and the value for the DataMT object as a list - typeList = [[rx.rxType,self[src,rx]] for rx in src.rxList] - # Insert the values to the temp array - for nr,(key,val) in enumerate(typeList): - tArrRec[key] = mkvc(val,2) - # Masked array - mArrRec = np.ma.MaskedArray(rec2ndarr(tArrRec),mask=np.isnan(rec2ndarr(tArrRec))).view(dtype=tArrRec.dtype) - # Unique freq and loc of the masked array - uniFLmarr = np.unique(mArrRec[['freq','x','y','z']]) - - try: - outTemp = recFunc.stack_arrays((outTemp,mArrRec)) - #outTemp = np.concatenate((outTemp,dataBlock),axis=0) - except NameError as e: - outTemp = mArrRec - - if 'RealImag' in returnType: - outArr = outTemp - if 'Complex' in returnType: - # Add the real and imaginary to a complex number - - outArr = np.empty(outTemp.shape,dtype=dtCP) - for comp in ['freq','x','y','z']: - outArr[comp] = outTemp[comp].copy() - for comp in ['zxx','zxy','zyx','zyy']: - outArr[comp] = outTemp[comp+'r'].copy() + 1j*outTemp[comp+'i'].copy() - - # Return - return outArr \ No newline at end of file diff --git a/simpegMT/Tests/test_ApparentResistivityLayer.py b/simpegMT/Tests/test_ApparentResistivityLayer.py deleted file mode 100644 index 8fbd3202..00000000 --- a/simpegMT/Tests/test_ApparentResistivityLayer.py +++ /dev/null @@ -1,48 +0,0 @@ -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() \ No newline at end of file diff --git a/simpegMT/Tests/test_Problem1D_againstAnalyticHalfspace.py b/simpegMT/Tests/test_Problem1D_againstAnalyticHalfspace.py new file mode 100644 index 00000000..04d491ee --- /dev/null +++ b/simpegMT/Tests/test_Problem1D_againstAnalyticHalfspace.py @@ -0,0 +1,112 @@ +import unittest +import SimPEG as simpeg +import simpegMT as simpegmt +from SimPEG.Utils import meshTensor +import numpy as np +# Define the tolerances +TOLr = 5e-2 +TOLp = 5e-2 + + +def setupSurvey(sigmaHalf): + + # Frequency + nFreq = 33 + freqs = np.logspace(3,-3,nFreq) + # 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)))]) + m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], 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 =[] + for freq in freqs: + srcList.append(simpegmt.SurveyMT.srcMT(freq,rxList)) + survey = simpegmt.SurveyMT.SurveyMT(srcList) + return survey, sigma, m1d + +def getAppResPhs(MTdata): + # Make impedance + 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 + zList = [] + for src in MTdata.survey.srcList: + zc = [src.freq] + for rx in src.rxList: + if 'i' in rx.rxType: + m=1j + else: + m = 1 + zc.append(m*MTdata[src,rx]) + zList.append(zc) + return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))] + +def appRes_TotalFieldNorm(sigmaHalf): + + # Make the survey + survey, sigma, mesh = setupSurvey(sigmaHalf) + problem = simpegmt.ProblemMT1D.eForm_TotalField(mesh) + problem.pair(survey) + + # Get the fields + fields = problem.fields(sigma) + + # Project the data + data = survey.projectFields(fields) + + # Calculate the app res and phs + app_r = np.array(getAppResPhs(data))[:,0] + + return np.linalg.norm(np.abs(app_r - np.ones(survey.nFreq)/sigmaHalf)*sigmaHalf) + +def appPhs_TotalFieldNorm(sigmaHalf): + + # Make the survey + survey, sigma, mesh = setupSurvey(sigmaHalf) + problem = simpegmt.ProblemMT1D.eForm_TotalField(mesh) + problem.pair(survey) + + # Get the fields + fields = problem.fields(sigma) + + # Project the data + data = survey.projectFields(fields) + + # Calculate the app phs + app_p = np.array(getAppResPhs(data))[:,1] + + return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*135)/ 135) + +class TestAnalytics(unittest.TestCase): + + def setUp(self): + pass + def test_appRes2en1(self):self.assertLess(appRes_TotalFieldNorm(2e-1), TOLr) + def test_appRes2en2(self):self.assertLess(appRes_TotalFieldNorm(2e-2), TOLr) + def test_appRes2en3(self):self.assertLess(appRes_TotalFieldNorm(2e-3), TOLr) + def test_appRes2en4(self):self.assertLess(appRes_TotalFieldNorm(2e-4), TOLr) + def test_appRes2en5(self):self.assertLess(appRes_TotalFieldNorm(2e-5), TOLr) + def test_appRes2en6(self):self.assertLess(appRes_TotalFieldNorm(2e-6), TOLr) + def test_appPhs2en1(self):self.assertLess(appPhs_TotalFieldNorm(2e-1), TOLp) + def test_appPhs2en2(self):self.assertLess(appPhs_TotalFieldNorm(2e-2), TOLp) + def test_appPhs2en3(self):self.assertLess(appPhs_TotalFieldNorm(2e-3), TOLp) + def test_appPhs2en4(self):self.assertLess(appPhs_TotalFieldNorm(2e-4), TOLp) + def test_appPhs2en5(self):self.assertLess(appPhs_TotalFieldNorm(2e-5), TOLp) + def test_appPhs2en6(self):self.assertLess(appPhs_TotalFieldNorm(2e-6), TOLp) + + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/simpegMT/Tests/test_Problem3D_againstAnalytic.py b/simpegMT/Tests/test_Problem3D_againstAnalytic.py new file mode 100644 index 00000000..8b459ad7 --- /dev/null +++ b/simpegMT/Tests/test_Problem3D_againstAnalytic.py @@ -0,0 +1,123 @@ +# Test functions +from glob import glob +import numpy as np, sys, os, time, scipy, subprocess +import simpegMT as simpegmt, SimPEG as simpeg +import unittest +import SimPEG as simpeg +import simpegMT as simpegmt +from SimPEG.Utils import meshTensor + + +TOLr = 5e-2 + +def getInputs(): + """ + Function that returns Mesh, freqs, rx_loc, elev. + """ + # Make a mesh + # M = simpeg.Mesh.TensorMesh([[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,1.6),(100.,10),(100,3,2)]], x0=['C','C',-3529.5360]) + M = simpeg.Mesh.TensorMesh([[(1000,6,-1.5),(1000.,6),(1000,6,1.5)],[(1000,6,-1.5),(1000.,2),(1000,6,1.5)],[(1000,10,-1.3),(1000.,2),(1000,10,1.3)]], x0=['C','C','C'])# Setup the model + # Set the frequencies + freqs = np.logspace(3,-3,7) + elev = 0 + + ## Setup the the survey object + # Receiver locations + rx_x, rx_y = np.meshgrid(np.arange(-3000,3001,500),np.arange(-1000,1001,500)) + rx_loc = np.array([[0, 0, 0]]) #np.hstack((simpeg.Utils.mkvc(rx_x,2),simpeg.Utils.mkvc(rx_y,2),elev+np.zeros((np.prod(rx_x.shape),1)))) + + return M, freqs, rx_loc, elev + + +def halfSpace(conds): + + M, freqs, rx_loc, elev = getInputs() + + # Model + ccM = M.gridCC + # conds = [1e-2] + groundInd = ccM[:,2] < elev + sig = np.zeros(M.nC) + 1e-8 + sig[groundInd] = conds + # Set the background, not the same as the model + sigBG = np.zeros(M.nC) + 1e-8 + sigBG[groundInd] = conds + + return (M, freqs, sig, sigBG, rx_loc) + +def twoLayer(): + M, freqs, rx_loc, elev = getInputs() + + # Model + ccM = M.gridCC + conds = [1e-2,1] + groundInd = ccM[:,2] < elev + botInd = ccM[:,2] < -3000 + sig = np.zeros(M.nC) + 1e-8 + sig[groundInd] = conds[1] + sig[botInd] = conds[0] + # Set the background, not the same as the model + sigBG = np.zeros(M.nC) + 1e-8 + sigBG[groundInd] = conds[1] + + + return (M, freqs, sig, sigBG, rx_loc) + +def runSimpegMTfwd_eForm_ps(inputsProblem): + M,freqs,sig,sigBG,rx_loc = inputsProblem + # Make a receiver list + rxList = [] + for rxType in ['zxyr','zxyi','zyxr','zyxi']: + rxList.append(simpegmt.SurveyMT.RxMT(rx_loc,rxType)) + # Source list + srcList =[] + for freq in freqs: + srcList.append(simpegmt.SurveyMT.srcMT(freq,rxList)) + # Survey MT + survey = simpegmt.SurveyMT.SurveyMT(srcList) + + ## Setup the problem object + problem = simpegmt.ProblemMT3D.eForm_ps(M) + problem.verbose = False + from pymatsolver import MumpsSolver + problem.Solver = MumpsSolver + problem.pair(survey) + + fields = problem.fields(sig,sigBG) + mtData = survey.projectFields(fields) + + return (survey, problem, fields, mtData) + + +def getAppResPhs(MTdata): + # Make impedance + 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 + recData = MTdata.toRecArray('Complex') + return appResPhs(recData['freq'],recData['zxy']), appResPhs(recData['freq'],recData['zyx']) + +def appResPhsHalfspace_eFrom_ps_Norm(sigmaHalf,appR=True): + + # Make the survey + survey, problem, fields, data = runSimpegMTfwd_eForm_ps(halfSpace(sigmaHalf)) + # Calculate the app phs + app_rpxy, app_rpyx = np.array(getAppResPhs(data)) + if appR: + return np.linalg.norm(np.abs(app_rpxy[0,:] - np.ones(survey.nFreq)/sigmaHalf) * sigmaHalf) + else: + return np.linalg.norm(np.abs(app_rpxy[1,:] - np.ones(survey.nFreq)/135) * 135) + +class TestAnalytics(unittest.TestCase): + + def setUp(self): + pass + def test_appRes2en1(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(2e-1), TOLr) + def test_appRes2en2(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(2e-2), TOLr) + def test_appRes2en3(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(2e-3), TOLr) + def test_appRes2en1(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(2e-1,False), TOLr) + def test_appRes2en2(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(2e-2,False), TOLr) + def test_appRes2en3(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(2e-3,False), TOLr) +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/simpegMT/Utils/MT1Dsolutions.py b/simpegMT/Utils/MT1Dsolutions.py index 33af58f5..1de508d9 100644 --- a/simpegMT/Utils/MT1Dsolutions.py +++ b/simpegMT/Utils/MT1Dsolutions.py @@ -24,7 +24,7 @@ def get1DEfields(m1d,sigma,freq,sourceAmp=1.0): 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 + ## Note: The analytic solution is derived with e^iwt bc = np.r_[Etot[0],Etot[-1]] # The right hand side rhs = -Aio*bc diff --git a/simpegMT/Utils/__init__.py b/simpegMT/Utils/__init__.py index 73005269..0b98a3c3 100644 --- a/simpegMT/Utils/__init__.py +++ b/simpegMT/Utils/__init__.py @@ -1,2 +1,3 @@ from MT1Dsolutions import * # Add the names of the functions from MT1Danalytic import * +from dataUtils import * diff --git a/simpegMT/Utils/srcUtils.py b/simpegMT/Utils/srcUtils.py new file mode 100644 index 00000000..f0163564 --- /dev/null +++ b/simpegMT/Utils/srcUtils.py @@ -0,0 +1,46 @@ +import SimPEG as simpeg, numpy as np + +def homo1DModelSource(mesh,freq,m_back): + ''' + Function that calculates and return background fields for a 3D mesh and model. + The calculuations use 1D field solution for a vertical slice throught model (south-western most column), + which is assigned at the fields everywhere for the respective polarizations.2 + + :param Simpeg mesh object mesh: Holds information on the discretization + :param float freq: The frequency to solve at + :param np.array m_back: Background model of conductivity to base the calculations on. + :rtype: numpy.ndarray (mesh.nE,2) + :return: eBG_bp, E fields for the background model at both polarizations. + + ''' + + # import + from simpegMT.Utils import get1DEfields + # Get a 1d solution for a halfspace background + mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]])) + # Note: Everything is using e^iwt + e0_1d = get1DEfields(mesh1d,mesh.r(m_back,'CC','CC','M')[0,0,:],freq) + # Setup x (east) polarization (_x) + ex_px = np.zeros(mesh.vnEx,dtype=complex) + ey_px = np.zeros((mesh.nEy,1),dtype=complex) + ez_px = np.zeros((mesh.nEz,1),dtype=complex) + # 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 + 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 + eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py)) + + # Return the electric fields + eBG_bp = np.hstack((eBG_px,eBG_py)) + return eBG_bp \ No newline at end of file diff --git a/simpegMT/__init__.py b/simpegMT/__init__.py index fe5cf927..a578bbe8 100644 --- a/simpegMT/__init__.py +++ b/simpegMT/__init__.py @@ -2,5 +2,7 @@ import Utils # import Tests import Sources -import ProblemMT -import SurveyMT +# from BaseMT import SurveyMT, RxMT, srcMT, DataMT, FieldsMT +import BaseMT +import SurveyMT, DataMT, FieldsMT +import ProblemMT1D, ProblemMT2D, ProblemMT3D