diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index c67652a2..7b2e1f70 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -1,27 +1,28 @@ # Run this file to add imports. ##### AUTOIMPORTS ##### -import DC_Analytic_Dipole -import DC_Forward_PseudoSection import EM_FDEM_1D_Inversion -import EM_FDEM_Analytic_MagDipoleWholespace -import EM_Schenkel_Morrison_Casing +import Mesh_QuadTree_Creation import EM_TDEM_1D_Inversion +import Mesh_QuadTree_FaceDiv +import Mesh_Tensor_Creation import FLOW_Richards_1D_Celia1990 -import Forward_BasicDirectCurrent +import DC_Forward_PseudoSection +import Mesh_Operators_CahnHilliard +import Mesh_Basic_Types import Inversion_IRLS import Inversion_Linear -import Mesh_Basic_PlotImage -import Mesh_Basic_Types -import Mesh_Operators_CahnHilliard -import Mesh_QuadTree_Creation -import Mesh_QuadTree_FaceDiv -import Mesh_QuadTree_HangingNodes -import Mesh_Tensor_Creation -import MT_1D_ForwardAndInversion +import EM_Schenkel_Morrison_Casing import MT_3D_Foward +import MT_1D_ForwardAndInversion +import MT_1D_analytic_nlayer_Earth +import Forward_BasicDirectCurrent +import EM_FDEM_Analytic_MagDipoleWholespace +import Mesh_Basic_PlotImage +import DC_Analytic_Dipole +import Mesh_QuadTree_HangingNodes -__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward"] +__examples__ = ["EM_FDEM_1D_Inversion", "Mesh_QuadTree_Creation", "EM_TDEM_1D_Inversion", "Mesh_QuadTree_FaceDiv", "Mesh_Tensor_Creation", "FLOW_Richards_1D_Celia1990", "DC_Forward_PseudoSection", "Mesh_Operators_CahnHilliard", "Mesh_Basic_Types", "Inversion_IRLS", "Inversion_Linear", "EM_Schenkel_Morrison_Casing", "MT_3D_Foward", "MT_1D_ForwardAndInversion", "MT_1D_analytic_nlayer_Earth", "Forward_BasicDirectCurrent", "EM_FDEM_Analytic_MagDipoleWholespace", "Mesh_Basic_PlotImage", "DC_Analytic_Dipole", "Mesh_QuadTree_HangingNodes"] ##### AUTOIMPORTS ##### diff --git a/SimPEG/NSEM/Utils/__init__.py b/SimPEG/NSEM/Utils/__init__.py index b683f8b4..acfc1a65 100644 --- a/SimPEG/NSEM/Utils/__init__.py +++ b/SimPEG/NSEM/Utils/__init__.py @@ -1,4 +1,5 @@ -from MT1Dsolutions import * # Add the names of the functions -from MT1Danalytic import * +from MT1Dsolutions import get1DEfields # Add the names of the functions +from MT1Danalytic import getEHfields, getImpedance from dataUtils import * from ediFilesUtils import * +from testUtils import * diff --git a/SimPEG/NSEM/Utils/srcUtils.py b/SimPEG/NSEM/Utils/srcUtils.py deleted file mode 100644 index 9bf141d8..00000000 --- a/SimPEG/NSEM/Utils/srcUtils.py +++ /dev/null @@ -1,46 +0,0 @@ -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 SimPEG.NSEM.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 diff --git a/SimPEG/NSEM/Utils/testUtils.py b/SimPEG/NSEM/Utils/testUtils.py new file mode 100644 index 00000000..8ac2b940 --- /dev/null +++ b/SimPEG/NSEM/Utils/testUtils.py @@ -0,0 +1,198 @@ +import unittest +import sys +from scipy.constants import mu_0 +import SimPEG as simpeg + +from SimPEG.Utils import meshTensor +import numpy as np + +np.random.seed(1100) +# Define the tolerances +TOLr = 5e-2 +TOLp = 5e-2 + + +def getAppResPhs(NSEMdata): + # Make impedance + from SimPEG.NSEM.Utils import appResPhs + zList = [] + for src in NSEMdata.survey.srcList: + zc = [src.freq] + for rx in src.rxList: + if 'i' in rx.rxType: + m=1j + else: + m = 1 + zc.append(m*NSEMdata[src,rx]) + zList.append(zc) + return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))] + + +def setup1DSurvey(sigmaHalf,tD=True,structure=False): + from SimPEG import NSEM + # 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],20,-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 + sigmaBack = sigma.copy() + # Add structure + if structure: + shallow = (m1d.gridCC < -200) * (m1d.gridCC > -600) + deep = (m1d.gridCC < -3000) * (m1d.gridCC > -5000) + sigma[shallow] = 1 + sigma[deep] = 0.1 + + rxList = [] + for rxType in ['z1dr','z1di']: + rxList.append(NSEM.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType)) + # Source list + srcList =[] + if tD: + for freq in freqs: + srcList.append(NSEM.SrcNSEM.polxy_1DhomotD(rxList,freq)) + else: + for freq in freqs: + srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq)) + + survey = NSEM.Survey(srcList) + return survey, sigma, m1d + + +def setupSimpegNSEM_ePrimSec(inputSetup,comp='Imp',singleFreq=False,expMap=True): + from SimPEG import NSEM + + M,freqs,sig,sigBG,rx_loc = inputSetup + # Make a receiver list + rxList = [] + if comp == 'All': + for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']: + rxList.append(NSEM.Rx(rx_loc,rxType)) + elif comp == 'Imp': + for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']: + rxList.append(NSEM.Rx(rx_loc,rxType)) + elif comp == 'Tip': + for rxType in ['tzxr','tzxi','tzyr','tzyi']: + rxList.append(NSEM.Rx(rx_loc,rxType)) + else: + rxList.append(NSEM.Rx(rx_loc,comp)) + # Source list + srcList =[] + + if singleFreq: + srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,singleFreq)) + else: + for freq in freqs: + srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq)) + # Survey NSEM + survey = NSEM.Survey(srcList) + + ## Setup the problem object + sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:] + if expMap: + problem = NSEM.Problem3D_ePrimSec(M,sigmaPrimary= np.log(sigma1d) ) + problem.mapping = simpeg.Maps.ExpMap(problem.mesh) + problem.curModel = np.log(sig) + else: + problem = NSEM.Problem3D_ePrimSec(M,sigmaPrimary= sigma1d) + problem.curModel = sig + problem.pair(survey) + problem.verbose = False + try: + from pymatsolver import MumpsSolver + problem.Solver = MumpsSolver + except: + pass + + return (survey, problem) + +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,6,-1.3),(1000.,6),(1000,6,1.3)]], x0=['C','C','C'])# Setup the model + M = simpeg.Mesh.TensorMesh([[(200,6,-1.5),(200.,4),(200,6,1.5)],[(200,6,-1.5),(200.,4),(200,6,1.5)],[(200,8,-1.5),(200.,8),(200,8,1.5)]], x0=['C','C','C'])# Setup the model + # Set the frequencies + freqs = np.logspace(1,-3,5) + elev = 0 + + ## Setup the the survey object + # Receiver locations + rx_x, rx_y = np.meshgrid(np.arange(-350,350,200),np.arange(-350,350,200)) + rx_loc = 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 random(conds): + ''' Returns a halfspace model based on the inputs''' + M, freqs, rx_loc, elev = getInputs() + + # Backround + sigBG = np.ones(M.nC)*conds + # Add randomness to the model (10% of the value). + sig = np.exp( np.log(sigBG) + np.random.randn(M.nC)*(conds)*1e-1 ) + + return (M, freqs, sig, sigBG, rx_loc) + +def halfSpace(conds): + ''' Returns a halfspace model based on the inputs''' + 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 blockInhalfSpace(conds): + ''' Returns a halfspace model based on the inputs''' + M, freqs, rx_loc, elev = getInputs() + + # Model + ccM = M.gridCC + # conds = [1e-2] + groundInd = ccM[:,2] < elev + sig = simpeg.Utils.ModelBuilder.defineBlock(M.gridCC,np.array([-1000,-1000,-1500]),np.array([1000,1000,-1000]),conds) + sig[~groundInd] = 1e-8 + # 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 twoLayer(conds): + ''' Returns a 2 layer model based on the conductivity values given''' + M, freqs, rx_loc, elev = getInputs() + + # Model + ccM = M.gridCC + 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) + diff --git a/tests/mt/__init__.py b/tests/mt/forward/__init__.py similarity index 100% rename from tests/mt/__init__.py rename to tests/mt/forward/__init__.py diff --git a/tests/mt/test_ApparentResistivityAnalytic.py b/tests/mt/forward/test_AnalyticFunctionVsAppResPhs.py similarity index 85% rename from tests/mt/test_ApparentResistivityAnalytic.py rename to tests/mt/forward/test_AnalyticFunctionVsAppResPhs.py index 44bc3c23..9a838b4b 100644 --- a/tests/mt/test_ApparentResistivityAnalytic.py +++ b/tests/mt/forward/test_AnalyticFunctionVsAppResPhs.py @@ -4,10 +4,7 @@ from SimPEG import NSEM 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 @@ -25,7 +22,7 @@ def appResNorm(sigmaHalf): Zarr = np.concatenate(Z) - app_r, app_p = appResPhs(freqs,Zarr) + app_r, app_p = NSEM.Utils.appResPhs(freqs,Zarr) return np.linalg.norm(np.abs(app_r - np.ones(nFreq)/sigmaHalf)) / np.log10(sigmaHalf) diff --git a/tests/mt/test_Problem1D_totalDvsPSvsAnalytic.py b/tests/mt/forward/test_Problem1D_AnalyticVsNumeric.py similarity index 57% rename from tests/mt/test_Problem1D_totalDvsPSvsAnalytic.py rename to tests/mt/forward/test_Problem1D_AnalyticVsNumeric.py index 09c68400..90a4b5ca 100644 --- a/tests/mt/test_Problem1D_totalDvsPSvsAnalytic.py +++ b/tests/mt/forward/test_Problem1D_AnalyticVsNumeric.py @@ -8,45 +8,6 @@ TOLr = 5e-2 TOLp = 5e-2 -def setupSurvey(sigmaHalf,tD=True): - - # 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],15,-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 - sigmaBack = sigma.copy() - # Add structure - shallow = (m1d.gridCC < -200) * (m1d.gridCC > -600) - deep = (m1d.gridCC < -3000) * (m1d.gridCC > -5000) - sigma[shallow] = 1 - sigma[deep] = 0.1 - - rxList = [] - for rxType in ['z1dr','z1di']: - rxList.append(NSEM.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType)) - # Source list - srcList =[] - if tD: - for freq in freqs: - srcList.append(NSEM.SrcNSEM.polxy_1DhomotD(rxList,freq)) - else: - for freq in freqs: - srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq)) - - survey = NSEM.Survey(srcList) - return survey, sigma, m1d - def getAppResPhs(NSEMdata): # Make impedance def appResPhs(freq,z): @@ -70,7 +31,7 @@ def calculateAnalyticSolution(srcList,mesh,model): data1D = NSEM.Data(surveyAna) for src in surveyAna.srcList: elev = src.rxList[0].locs[0] - anaEd, anaEu, anaHd, anaHu = NSEM.Utils.NSEM1Danalytic.getEHfields(mesh,model,src.freq,elev) + anaEd, anaEu, anaHd, anaHu = NSEM.Utils.MT1Danalytic.getEHfields(mesh,model,src.freq,elev) anaE = anaEd+anaEu anaH = anaHd+anaHu # Scale the solution @@ -86,8 +47,8 @@ def dataMis_AnalyticTotalDomain(sigmaHalf): # Make the survey # Total domain solution - surveyTD, sigma, mesh = setupSurvey(sigmaHalf) - problemTD = NSEM.Problem1D.eForm_TotalField(mesh) + surveyTD, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf) + problemTD = NSEM.Problem1D_eTotal(mesh) # This not fully implemented problemTD.pair(surveyTD) # Analytic data dataAnaObj = calculateAnalyticSolution(surveyTD.srcList,mesh,sigma) @@ -108,16 +69,16 @@ def dataMis_AnalyticPrimarySecondary(sigmaHalf): # Make the survey # Primary secondary - surveyPS, sigmaPS, mesh = setupSurvey(sigmaHalf,tD=False) - problemPS = NSEM.Problem1D.eForm_psField(mesh) - problemPS.sigmaPrimary = sigmaPS - problemPS.pair(surveyPS) + survey, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf,False,structure=True) # Analytic data - dataAnaObj = calculateAnalyticSolution(surveyPS.srcList,mesh,sigmaPS) + problem = NSEM.Problem1D_ePrimSec(mesh, sigmaPrimary = sigma) + problem.pair(survey) - dataPS = surveyPS.dpred(sigmaPS) + dataAnaObj = calculateAnalyticSolution(survey.srcList,mesh,sigma) + + data = survey.dpred(sigma) dataAna = simpeg.mkvc(dataAnaObj) - return np.all((dataPS - dataAna)/dataAna < 2.) + return np.all((data - dataAna)/dataAna < 2.) diff --git a/tests/mt/forward/test_Problem1D_VsAnalyticHalfspace.py b/tests/mt/forward/test_Problem1D_VsAnalyticHalfspace.py new file mode 100644 index 00000000..098ca0ef --- /dev/null +++ b/tests/mt/forward/test_Problem1D_VsAnalyticHalfspace.py @@ -0,0 +1,102 @@ +import unittest +import SimPEG as simpeg +from SimPEG import NSEM +from SimPEG.Utils import meshTensor +import numpy as np +# Define the tolerances +TOLr = 5e-1 +TOLp = 5e-1 + + + +def appRes_TotalFieldNorm(sigmaHalf): + + # Make the survey + survey, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf) + problem = NSEM.Problem1D_eTotal(mesh) + problem.pair(survey) + + # Get the fields + fields = problem.fields(sigma) + + # Project the data + data = survey.eval(fields) + + # Calculate the app res and phs + app_r = np.array(NSEM.Utils.testUtils.getAppResPhs(data))[:,0] + + return np.linalg.norm(np.abs(np.log(app_r) - np.log(np.ones(survey.nFreq)/sigmaHalf))*np.log(sigmaHalf)) + +def appPhs_TotalFieldNorm(sigmaHalf): + + # Make the survey + survey, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf) + problem = NSEM.Problem1D_eTotal(mesh) + problem.pair(survey) + + # Get the fields + fields = problem.fields(sigma) + + # Project the data + data = survey.eval(fields) + + # Calculate the app phs + app_p = np.array(NSEM.Utils.testUtils.getAppResPhs(data))[:,1] + + return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45) + +def appRes_psFieldNorm(sigmaHalf): + + # Make the survey + survey, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf,False) + problem = NSEM.Problem1D_ePrimSec(mesh, sigmaPrimary = sigma) + problem.pair(survey) + + # Get the fields + fields = problem.fields(sigma) + + # Project the data + data = survey.eval(fields) + + # Calculate the app res and phs + app_r = np.array(NSEM.Utils.testUtils.getAppResPhs(data))[:,0] + + return np.linalg.norm(np.abs(np.log(app_r) - np.log(np.ones(survey.nFreq)/sigmaHalf))*np.log(sigmaHalf)) + +def appPhs_psFieldNorm(sigmaHalf): + + # Make the survey + survey, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf,False) + problem = NSEM.Problem1D_ePrimSec(mesh, sigmaPrimary = sigma) + problem.pair(survey) + + # Get the fields + fields = problem.fields(sigma) + + # Project the data + data = survey.eval(fields) + + # Calculate the app phs + app_p = np.array(NSEM.Utils.testUtils.getAppResPhs(data))[:,1] + + return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45) + +class TestAnalytics(unittest.TestCase): + + def setUp(self): + pass + # Total Fields + # def test_appRes2en1(self):self.assertLess(appRes_TotalFieldNorm(2e-1), TOLr) + # def test_appPhs2en1(self):self.assertLess(appPhs_TotalFieldNorm(2e-1), TOLp) + + # Primary/secondary + def test_appRes1en0_ps(self):self.assertLess(appRes_psFieldNorm(1e-0), TOLr) + def test_appPhs1en0_ps(self):self.assertLess(appPhs_psFieldNorm(1e-0), TOLp) + def test_appRes2en1_ps(self):self.assertLess(appRes_psFieldNorm(2e-1), TOLr) + def test_appPhs2en1_ps(self):self.assertLess(appPhs_psFieldNorm(2e-1), TOLp) + def test_appRes2en3_ps(self):self.assertLess(appRes_psFieldNorm(2e-3), TOLr) + def test_appPhs2en3_ps(self):self.assertLess(appPhs_psFieldNorm(2e-3), TOLp) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/mt/forward/test_Problem3D_VsAnalyticSolution.py b/tests/mt/forward/test_Problem3D_VsAnalyticSolution.py new file mode 100644 index 00000000..b736bc17 --- /dev/null +++ b/tests/mt/forward/test_Problem3D_VsAnalyticSolution.py @@ -0,0 +1,54 @@ +# Test functions +from glob import glob +import numpy as np, sys, os, time, scipy, subprocess +import SimPEG as simpeg +import unittest +from SimPEG import NSEM +from SimPEG.Utils import meshTensor +from scipy.constants import mu_0 + +np.random.seed(1100) + +TOLr = 1 +TOLp = 2 +FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order +CONDUCTIVITY = 1e1 +MU = mu_0 +freq = [1e-1, 2e-1] +addrandoms = True + +def appResPhsHalfspace_eFrom_ps_Norm(sigmaHalf,appR=True,expMap=False): + if appR: + label = 'resistivity' + else: + label = 'phase' + print 'Apperent {:s} test of eFormulation primary/secondary at {:g}\n\n'.format(label,sigmaHalf) + + # Calculate the app phs + survey, problem = NSEM.Utils.testUtils.setupSimpegNSEM_ePrimSec(NSEM.Utils.testUtils.halfSpace(sigmaHalf),expMap=expMap) + data = problem.dataPair(survey,survey.dpred(problem.curModel)) + recData = data.toRecArray('Complex') + app_rpxy, app_rpyx = NSEM.Utils.appResPhs(recData['freq'],recData['zxy'])[0], NSEM.Utils.appResPhs(recData['freq'],recData['zyx'])[0] + if appR: + return np.linalg.norm( np.abs(np.log10(app_rpxy[0]) - np.log10(1./sigmaHalf)) * np.log10(sigmaHalf )) + else: + return np.linalg.norm( np.abs(app_rpxy[1] + 135) / 135 ) + + +class TestAnalytics(unittest.TestCase): + + def setUp(self): + # Make the survey and the problem + pass + + # # Test apparent resistivity and phase + def test_appRes1en2(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(1e-2),TOLr) + def test_appPhs1en2(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(1e-2,False),TOLp) + + def test_appRes1en1(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(1e-1),TOLr) + def test_appPhs1en1(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(1e-1,False),TOLp) + + + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/tests/mt/inversion/__init__.py b/tests/mt/inversion/__init__.py new file mode 100644 index 00000000..420388ef --- /dev/null +++ b/tests/mt/inversion/__init__.py @@ -0,0 +1,12 @@ +import os +import glob +import unittest + +if __name__ == '__main__': + test_file_strings = glob.glob('test_*.py') + module_strings = [str[0:len(str)-3] for str in test_file_strings] + suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str + in module_strings] + testSuite = unittest.TestSuite(suites) + + unittest.TextTestRunner(verbosity=2).run(testSuite) diff --git a/tests/mt/inversion/test_Problem3D_Adjoint.py b/tests/mt/inversion/test_Problem3D_Adjoint.py new file mode 100644 index 00000000..2061e1ff --- /dev/null +++ b/tests/mt/inversion/test_Problem3D_Adjoint.py @@ -0,0 +1,58 @@ +# Test functions +from glob import glob +import numpy as np, sys, os, time, scipy, subprocess +import SimPEG as simpeg +import unittest +from SimPEG import NSEM +from SimPEG.Utils import meshTensor +from scipy.constants import mu_0 + + +TOLr = 5e-2 +TOL = 1e-4 +FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order +CONDUCTIVITY = 1e1 +MU = mu_0 +freq = [1e-1, 2e-1] +addrandoms = True + + + +def JvecAdjointTest(inputSetup,comp='All',freq=False): + (M, freqs, sig, sigBG, rx_loc) = inputSetup + survey, problem = NSEM.Utils.testUtils.setupSimpegNSEM_ePrimSec(inputSetup,comp='All',singleFreq=freq) + print 'Adjoint test of eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,str(survey.freqs)) + + m = sig + u = problem.fields(m) + + v = np.random.rand(survey.nD,) + # print problem.PropMap.PropModel.nP + w = np.random.rand(problem.mesh.nC,) + + vJw = v.ravel().dot(problem.Jvec(m, w, u)) + wJtv = w.ravel().dot(problem.Jtvec(m, v, u)) + tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) + print ' vJw wJtv vJw - wJtv tol abs(vJw - wJtv) < tol' + print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol + return np.abs(vJw - wJtv) < tol + + +class NSEM_AdjointTests(unittest.TestCase): + + def setUp(self): + pass + + # Test the adjoint of Jvec and Jtvec + # def test_JvecAdjoint_zxxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxr',.1)) + # def test_JvecAdjoint_zxxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxi',.1)) + # def test_JvecAdjoint_zxyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyr',.1)) + # def test_JvecAdjoint_zxyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyi',.1)) + # def test_JvecAdjoint_zyxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxr',.1)) + # def test_JvecAdjoint_zyxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxi',.1)) + # def test_JvecAdjoint_zyyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyr',.1)) + # def test_JvecAdjoint_zyyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyi',.1)) + def test_JvecAdjoint_All(self):self.assertTrue(JvecAdjointTest(NSEM.Utils.testUtils.random(1e-2),'All',.1)) + +if __name__ == '__main__': + unittest.main() diff --git a/tests/mt/inversion/test_Problem3D_Derivs.py b/tests/mt/inversion/test_Problem3D_Derivs.py new file mode 100644 index 00000000..9e365ec9 --- /dev/null +++ b/tests/mt/inversion/test_Problem3D_Derivs.py @@ -0,0 +1,83 @@ +# Test functions +from glob import glob +import numpy as np, sys, os, time, scipy, subprocess +import SimPEG as simpeg +import unittest +from SimPEG import NSEM +from SimPEG.Utils import meshTensor +from scipy.constants import mu_0 + +np.random.seed(1100) + +TOLr = 5e-2 +TOL = 1e-4 +FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order +CONDUCTIVITY = 1e1 +MU = mu_0 +freq = [1e-1, 2e-1] +addrandoms = True + + +# Test the Jvec derivative +def DerivJvecTest(inputSetup,comp='All',freq=False,expMap=True): + (M, freqs, sig, sigBG, rx_loc) = inputSetup + survey, problem = NSEM.Utils.testUtils.setupSimpegNSEM_ePrimSec(inputSetup,comp=comp,singleFreq=freq,expMap=expMap) + print 'Derivative test of Jvec for eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,survey.freqs) + # problem.mapping = simpeg.Maps.ExpMap(problem.mesh) + # problem.sigmaPrimary = np.log(sigBG) + x0 = np.log(sigBG) + # cond = sig[0] + # x0 = np.log(np.ones(problem.mesh.nC)*cond) + # problem.sigmaPrimary = x0 + # if True: + # x0 = x0 + np.random.randn(problem.mesh.nC)*cond*1e-1 + survey = problem.survey + def fun(x): + return survey.dpred(x), lambda x: problem.Jvec(x0, x) + return simpeg.Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR) + +def DerivProjfieldsTest(inputSetup,comp='All',freq=False): + + survey, problem = NSEM.Utils.testUtils.setupSimpegNSEM_ePrimSec(inputSetup,comp,freq) + print 'Derivative test of data projection for eFormulation primary/secondary\n\n' + # problem.mapping = simpeg.Maps.ExpMap(problem.mesh) + # Initate things for the derivs Test + src = survey.srcList[0] + rx = src.rxList[0] + + u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j + u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j + u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2))) + f0 = problem.fieldsPair(survey.mesh,survey) + # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2))) + f0[src,'e_pxSolution'] = u0[:len(u0)/2]#u0x + f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y + + def fun(u): + f = problem.fieldsPair(survey.mesh,survey) + f[src,'e_pxSolution'] = u[:len(u)/2] + f[src,'e_pySolution'] = u[len(u)/2::] + return rx.eval(src,survey.mesh,f), lambda t: rx.evalDeriv(src,survey.mesh,f0,simpeg.mkvc(t,2)) + + return simpeg.Tests.checkDerivative(fun, u0, num=3, plotIt=False, eps=FLR) + + + +class NSEM_DerivTests(unittest.TestCase): + + def setUp(self): + pass + + # Do a derivative test of Jvec + # def test_derivJvec_zxxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxr',.1)) + # def test_derivJvec_zxxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxi',.1)) + # def test_derivJvec_zxyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyr',.1)) + # def test_derivJvec_zxyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyi',.1)) + # def test_derivJvec_zyxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxr',.1)) + # def test_derivJvec_zyxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxi',.1)) + # def test_derivJvec_zyyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyr',.1)) + # def test_derivJvec_zyyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyi',.1)) + def test_derivJvec_All(self):self.assertTrue(DerivJvecTest(NSEM.Utils.testUtils.random(1e-2),'All',.1)) + +if __name__ == '__main__': + unittest.main() diff --git a/tests/mt/test_Problem1D_againstAnalyticHalfspace.py b/tests/mt/test_Problem1D_againstAnalyticHalfspace.py deleted file mode 100644 index d848c1ba..00000000 --- a/tests/mt/test_Problem1D_againstAnalyticHalfspace.py +++ /dev/null @@ -1,147 +0,0 @@ -import unittest -import SimPEG as simpeg -from SimPEG import NSEM -from SimPEG.Utils import meshTensor -import numpy as np -# Define the tolerances -TOLr = 5e-2 -TOLp = 5e-2 - - -def setupSurvey(sigmaHalf,tD=True): - - # 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(NSEM.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType)) - # Source list - srcList =[] - if tD: - for freq in freqs: - srcList.append(NSEM.SrcNSEM.polxy_1DhomotD(rxList,freq)) - else: - for freq in freqs: - srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq)) - - survey = NSEM.Survey(srcList) - return survey, sigma, m1d - -def getAppResPhs(NSEMdata): - # 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 NSEMdata.survey.srcList: - zc = [src.freq] - for rx in src.rxList: - if 'i' in rx.rxType: - m=1j - else: - m = 1 - zc.append(m*NSEMdata[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 = NSEM.Problem1D.eForm_TotalField(mesh) - problem.pair(survey) - - # Get the fields - fields = problem.fields(sigma) - - # Project the data - data = survey.eval(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 = NSEM.Problem1D.eForm_TotalField(mesh) - problem.pair(survey) - - # Get the fields - fields = problem.fields(sigma) - - # Project the data - data = survey.eval(fields) - - # Calculate the app phs - app_p = np.array(getAppResPhs(data))[:,1] - - return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45) - -def appRes_psFieldNorm(sigmaHalf): - - # Make the survey - survey, sigma, mesh = setupSurvey(sigmaHalf,False) - problem = NSEM.Problem1D_ePrimSec(mesh, sigmaPrimary = sigma) - problem.pair(survey) - - # Get the fields - fields = problem.fields(sigma) - - # Project the data - data = survey.eval(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_psFieldNorm(sigmaHalf): - - # Make the survey - survey, sigma, mesh = setupSurvey(sigmaHalf,False) - problem = NSEM.Problem1D_ePrimSec(mesh, sigmaPrimary = sigma) - problem.pair(survey) - - # Get the fields - fields = problem.fields(sigma) - - # Project the data - data = survey.eval(fields) - - # Calculate the app phs - app_p = np.array(getAppResPhs(data))[:,1] - - return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45) - -class TestAnalytics(unittest.TestCase): - - def setUp(self): - pass - # Total Fields - # def test_appRes2en1(self):self.assertLess(appRes_TotalFieldNorm(2e-1), TOLr) - # def test_appPhs2en1(self):self.assertLess(appPhs_TotalFieldNorm(2e-1), TOLp) - - # Primary/secondary - def test_appRes2en2_ps(self):self.assertLess(appRes_psFieldNorm(2e-2), TOLr) - def test_appPhs2en2_ps(self):self.assertLess(appPhs_psFieldNorm(2e-2), TOLp) - -if __name__ == '__main__': - unittest.main() diff --git a/tests/mt/test_Problem3D_againstAnalytic.py b/tests/mt/test_Problem3D_againstAnalytic.py deleted file mode 100644 index a267c9d7..00000000 --- a/tests/mt/test_Problem3D_againstAnalytic.py +++ /dev/null @@ -1,270 +0,0 @@ -# Test functions -from glob import glob -import numpy as np, sys, os, time, scipy, subprocess -import SimPEG as simpeg -import unittest -from SimPEG import NSEM -from SimPEG.Utils import meshTensor -from scipy.constants import mu_0 - -np.random.seed(1100) - -TOLr = 5e-2 -TOL = 1e-4 -FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order -CONDUCTIVITY = 1e1 -MU = mu_0 -freq = [1e-1, 2e-1] -addrandoms = True - - -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,6,-1.3),(1000.,6),(1000,6,1.3)]], x0=['C','C','C'])# Setup the model - M = simpeg.Mesh.TensorMesh([[(1000,6,-1.5),(1000.,4),(1000,6,1.5)],[(1000,6,-1.5),(1000.,4),(1000,6,1.5)],[(500,8,-1.3),(500.,8),(500,8,1.3)]], x0=['C','C','C'])# Setup the model - # Set the frequencies - freqs = np.logspace(1,-3,5) - elev = 0 - - ## Setup the the survey object - # Receiver locations - rx_x, rx_y = np.meshgrid(np.arange(-1000,1001,500),np.arange(-1000,1001,500)) - rx_loc = 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 random(conds): - ''' Returns a halfspace model based on the inputs''' - M, freqs, rx_loc, elev = getInputs() - - # Backround - sigBG = np.ones(M.nC)*conds - # Add randomness to the model (10% of the value). - sig = np.exp( np.log(sigBG) + np.random.randn(M.nC)*(conds)*1e-1 ) - - return (M, freqs, sig, sigBG, rx_loc) - -def halfSpace(conds): - ''' Returns a halfspace model based on the inputs''' - 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 blockInhalfSpace(conds): - ''' Returns a halfspace model based on the inputs''' - M, freqs, rx_loc, elev = getInputs() - - # Model - ccM = M.gridCC - # conds = [1e-2] - groundInd = ccM[:,2] < elev - sig = simpeg.Utils.ModelBuilder.defineBlock(M.gridCC,np.array([-1000,-1000,-1500]),np.array([1000,1000,-1000]),conds) - sig[~groundInd] = 1e-8 - # 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 twoLayer(conds): - ''' Returns a 2 layer model based on the conductivity values given''' - M, freqs, rx_loc, elev = getInputs() - - # Model - ccM = M.gridCC - 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 setupSimpegNSEMfwd_eForm_ps(inputSetup,comp='Imp',singleFreq=False,expMap=True): - M,freqs,sig,sigBG,rx_loc = inputSetup - # Make a receiver list - rxList = [] - if comp == 'All': - for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']: - rxList.append(NSEM.Rx(rx_loc,rxType)) - elif comp == 'Imp': - for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']: - rxList.append(NSEM.Rx(rx_loc,rxType)) - elif comp == 'Tip': - for rxType in ['tzxr','tzxi','tzyr','tzyi']: - rxList.append(NSEM.Rx(rx_loc,rxType)) - else: - rxList.append(NSEM.Rx(rx_loc,comp)) - # Source list - srcList =[] - - if singleFreq: - srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,singleFreq)) - else: - for freq in freqs: - srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq)) - # Survey NSEM - survey = NSEM.Survey(srcList) - - ## Setup the problem object - sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:] - if expMap: - problem = NSEM.Problem3D_ePrimSec(M,sigmaPrimary= np.log(sigma1d) ) - problem.mapping = simpeg.Maps.ExpMap(problem.mesh) - problem.curModel = np.log(sig) - else: - problem = NSEM.Problem3D_ePrimSec(M,sigmaPrimary= sigma1d) - problem.curModel = sig - problem.pair(survey) - problem.verbose = False - try: - from pymatsolver import MumpsSolver - problem.Solver = MumpsSolver - except: - pass - - return (survey, problem) - -def getAppResPhs(NSEMdata): - # 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 = NSEMdata.toRecArray('Complex') - return appResPhs(recData['freq'],recData['zxy']), appResPhs(recData['freq'],recData['zyx']) - -def JvecAdjointTest(inputSetup,comp='All',freq=False): - (M, freqs, sig, sigBG, rx_loc) = inputSetup - survey, problem = setupSimpegNSEMfwd_eForm_ps(inputSetup,comp='All',singleFreq=freq) - print 'Adjoint test of eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,str(survey.freqs)) - - m = sig - u = problem.fields(m) - - v = np.random.rand(survey.nD,) - # print problem.PropMap.PropModel.nP - w = np.random.rand(problem.mesh.nC,) - - vJw = v.ravel().dot(problem.Jvec(m, w, u)) - wJtv = w.ravel().dot(problem.Jtvec(m, v, u)) - tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) - print ' vJw wJtv vJw - wJtv tol abs(vJw - wJtv) < tol' - print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol - return np.abs(vJw - wJtv) < tol - -# Test the Jvec derivative -def DerivJvecTest(inputSetup,comp='All',freq=False,expMap=True): - (M, freqs, sig, sigBG, rx_loc) = inputSetup - survey, problem = setupSimpegNSEMfwd_eForm_ps(inputSetup,comp=comp,singleFreq=freq,expMap=expMap) - print 'Derivative test of Jvec for eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,survey.freqs) - # problem.mapping = simpeg.Maps.ExpMap(problem.mesh) - # problem.sigmaPrimary = np.log(sigBG) - x0 = np.log(sigBG) - # cond = sig[0] - # x0 = np.log(np.ones(problem.mesh.nC)*cond) - # problem.sigmaPrimary = x0 - # if True: - # x0 = x0 + np.random.randn(problem.mesh.nC)*cond*1e-1 - survey = problem.survey - def fun(x): - return survey.dpred(x), lambda x: problem.Jvec(x0, x) - return simpeg.Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR) - -def DerivProjfieldsTest(inputSetup,comp='All',freq=False): - - survey, problem = setupSimpegNSEMfwd_eForm_ps(inputSetup,comp,freq) - print 'Derivative test of data projection for eFormulation primary/secondary\n\n' - # problem.mapping = simpeg.Maps.ExpMap(problem.mesh) - # Initate things for the derivs Test - src = survey.srcList[0] - rx = src.rxList[0] - - u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j - u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j - u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2))) - f0 = problem.fieldsPair(survey.mesh,survey) - # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2))) - f0[src,'e_pxSolution'] = u0[:len(u0)/2]#u0x - f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y - - def fun(u): - f = problem.fieldsPair(survey.mesh,survey) - f[src,'e_pxSolution'] = u[:len(u)/2] - f[src,'e_pySolution'] = u[len(u)/2::] - return rx.eval(src,survey.mesh,f), lambda t: rx.evalDeriv(src,survey.mesh,f0,simpeg.mkvc(t,2)) - - return simpeg.Tests.checkDerivative(fun, u0, num=3, plotIt=False, eps=FLR) - -def appResPhsHalfspace_eFrom_ps_Norm(sigmaHalf,appR=True,expMap=False): - if appR: - label = 'resistivity' - else: - label = 'phase' - # Make the survey and the problem - survey, problem = setupSimpegNSEMfwd_eForm_ps(halfSpace(sigmaHalf),expMap=expMap) - print 'Apperent {:s} test of eFormulation primary/secondary at {:g}\n\n'.format(label,sigmaHalf) - - data = problem.dataPair(survey,survey.dpred(problem.curModel)) - # Calculate the app phs - app_rpxy, app_rpyx = np.array(getAppResPhs(data)) - if appR: - return np.all(np.abs(app_rpxy[0,:] - 1./sigmaHalf) * sigmaHalf < .4) - else: - return np.all(np.abs(app_rpxy[1,:] + 135) / 135 < .4) - -class TestAnalytics(unittest.TestCase): - - def setUp(self): - pass - # # Test apparent resistivity and phase - def test_appRes1en2(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-2)) - def test_appPhs1en2(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-2,False)) - - def test_appRes1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3)) - def test_appPhs1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3,False)) - - # Do a derivative test of Jvec - # def test_derivJvec_zxxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxr',.1)) - # def test_derivJvec_zxxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxi',.1)) - # def test_derivJvec_zxyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyr',.1)) - # def test_derivJvec_zxyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyi',.1)) - # def test_derivJvec_zyxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxr',.1)) - # def test_derivJvec_zyxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxi',.1)) - # def test_derivJvec_zyyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyr',.1)) - # def test_derivJvec_zyyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyi',.1)) - def test_derivJvec_All(self):self.assertTrue(DerivJvecTest(random(1e-2),'All',.1)) - - # Test the adjoint of Jvec and Jtvec - # def test_JvecAdjoint_zxxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxr',.1)) - # def test_JvecAdjoint_zxxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxi',.1)) - # def test_JvecAdjoint_zxyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyr',.1)) - # def test_JvecAdjoint_zxyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyi',.1)) - # def test_JvecAdjoint_zyxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxr',.1)) - # def test_JvecAdjoint_zyxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxi',.1)) - # def test_JvecAdjoint_zyyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyr',.1)) - # def test_JvecAdjoint_zyyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyi',.1)) - def test_JvecAdjoint_All(self):self.assertTrue(JvecAdjointTest(random(1e-2),'All',.1)) - -if __name__ == '__main__': - unittest.main()