diff --git a/SimPEG/DCIP/BaseDC.py b/SimPEG/DCIP/BaseDC.py index 33800188..c9aeb264 100644 --- a/SimPEG/DCIP/BaseDC.py +++ b/SimPEG/DCIP/BaseDC.py @@ -1,7 +1,5 @@ from SimPEG import * - - class FieldsDC_CC(Problem.Fields): knownFields = {'phi_sol':'CC'} aliasFields = { @@ -22,7 +20,7 @@ class FieldsDC_CC(Problem.Fields): # for i, src in enumerate(srcList): # phi_p = src.phi_p(self.survey.prob) # if phi_p is not None: - # phi[:,i] += phi_p + # phi[:,i] += phi_p return phi def _e(self, phi_sol, srcList): @@ -30,9 +28,9 @@ class FieldsDC_CC(Problem.Fields): # for i, src in enumerate(srcList): # e_p = src.e_p(self.survey.prob) # if e_p is not None: - # e[:,i] += e_p + # e[:,i] += e_p return e - + def _j(self, phi_sol, srcList): j = -self._Mfinv*self.survey.prob.Msig*self._cellGrad*phi_sol @@ -56,7 +54,7 @@ class SrcDipole(Survey.BaseSrc): super(SrcDipole, self).__init__(rxList, **kwargs) def eval(self, prob): - # Recompute rhs + # Recompute rhs # if getattr(self, '_rhsDict', None) is None: # self._rhsDict = {} # if mesh not in self._rhsDict: @@ -191,12 +189,12 @@ class ProblemDC_CC(Problem.BaseProblem): def fields(self, m): - F = self.fieldsPair(self.mesh, self.survey) + F = self.fieldsPair(self.mesh, self.survey) self.curModel = m A = self.A self.Ainv = self.Solver(A, **self.solverOpts) RHS = self.getRHS() - Phi = self.Ainv * RHS + Phi = self.Ainv * RHS Srcs = self.survey.srcList F[Srcs, 'phi_sol'] = Phi @@ -302,51 +300,51 @@ def readUBC_DC2DModel(fileName): Output: :param SimPEG TensorMesh 2D object :return - + Created on Thu Nov 12 13:14:10 2015 @author: dominiquef """ - + # Open fileand skip header... assume that we know the mesh already obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!') - + dim = np.array(obsfile[0].split(),dtype=float) - + temp = np.array(obsfile[1].split(),dtype=float) - + if len(temp) > 1: model = np.zeros(dim) - + for ii in range(len(obsfile)-1): mm = np.array(obsfile[ii+1].split(),dtype=float) model[:,ii] = mm - + model = model[:,::-1] - + else: - + if len(obsfile[1:])==1: mm = np.array(obsfile[1:].split(),dtype=float) - + else: mm = np.array(obsfile[1:],dtype=float) - + # Permute the second dimension to flip the order model = mm.reshape(dim[1],dim[0]) - + model = model[::-1,:] model = np.transpose(model, (1, 0)) - + model = mkvc(model) return model def plot_pseudoSection(Tx,Rx,data,z0, stype): - + from SimPEG import np, mkvc from scipy.interpolate import griddata from matplotlib.colors import LogNorm @@ -355,265 +353,265 @@ def plot_pseudoSection(Tx,Rx,data,z0, stype): """ Read list of 2D tx-rx location and plot a speudo-section of apparent resistivity. - + Assumes flat topo for now... - + Input: :param d2D, z0 :switch stype -> Either 'pdp' (pole-dipole) | 'dpdp' (dipole-dipole) - + Output: :figure scatter plot overlayed on image - + Created on Mon December 7th, 2015 - + @author: dominiquef - + """ #d2D = np.asarray(d2D) - + midl = [] midz = [] rho = [] - + for ii in range(len(Tx)): # Get distances between each poles - rC1P1 = np.abs(Tx[ii][0] - Rx[ii][:,0]) + rC1P1 = np.abs(Tx[ii][0] - Rx[ii][:,0]) rC2P1 = np.abs(Tx[ii][1] - Rx[ii][:,0]) rC1P2 = np.abs(Tx[ii][1] - Rx[ii][:,1]) rC2P2 = np.abs(Tx[ii][0] - Rx[ii][:,1]) - rP1P2 = np.abs(Rx[ii][:,1] - Rx[ii][:,0]) - + rP1P2 = np.abs(Rx[ii][:,1] - Rx[ii][:,0]) + # Compute apparent resistivity if re.match(stype,'pdp'): rho = np.hstack([rho, data[ii] * 2*np.pi * rC1P1 * ( rC1P1 + rP1P2 ) / rP1P2] ) - + elif re.match(stype,'dpdp'): rho = np.hstack([rho, data[ii] * 2*np.pi / ( 1/rC1P1 - 1/rC2P1 - 1/rC1P2 + 1/rC2P2 ) ]) - + Cmid = (Tx[ii][0] + Tx[ii][1])/2 Pmid = (Rx[ii][:,0] + Rx[ii][:,1])/2 - + midl = np.hstack([midl, ( Cmid + Pmid )/2 ]) midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + z0 ]) - - + + # Grid points grid_x, grid_z = np.mgrid[np.min(midl):np.max(midl), np.min(midz):np.max(midz)] grid_rho = griddata(np.c_[midl,midz], np.log10(abs(1/rho.T)), (grid_x, grid_z), method='linear') - - + + #plt.subplot(2,1,2) plt.imshow(grid_rho.T, extent = (np.min(midl),np.max(midl),np.min(midz),np.max(midz)), origin='lower', alpha=0.8) cbar = plt.colorbar(format = '%.2f',fraction=0.02) cmin,cmax = cbar.get_clim() ticks = np.linspace(cmin,cmax,3) cbar.set_ticks(ticks) - + # Plot apparent resistivity plt.scatter(midl,midz,s=50,c=np.log10(abs(1/rho.T))) def gen_DCIPsurvey(endl, mesh, stype, a, b, n): - + from SimPEG import np import re """ Load in endpoints and survey specifications to generate Tx, Rx location stations. - + Assumes flat topo for now... - + Input: :param endl -> input endpoints [x1, y1, z1, x2, y2, z2] :object mesh -> SimPEG mesh object :switch stype -> "dpdp" (dipole-dipole) | "pdp" (pole-dipole) | 'gradient' : param a, n -> pole seperation, number of rx dipoles per tx - + Output: :param Tx, Rx -> List objects for each tx location Lines: P1x, P1y, P1z, P2x, P2y, P2z - + Created on Wed December 9th, 2015 - + @author: dominiquef - + """ def xy_2_r(x1,x2,y1,y2): r = np.sqrt( np.sum((x2 - x1)**2 + (y2 - y1)**2) ) - return r - + return r + ## Evenly distribute electrodes and put on surface # Mesure survey length and direction dl_len = xy_2_r(endl[0,0],endl[1,0],endl[0,1],endl[1,1]) - + dl_x = ( endl[1,0] - endl[0,0] ) / dl_len dl_y = ( endl[1,1] - endl[0,1] ) / dl_len - + nstn = np.floor( dl_len / a ) - + # Compute discrete pole location along line stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*a stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*a - + # Create line of P1 locations M = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]] - + # Create line of P2 locations N = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]] - + ## Build list of Tx-Rx locations depending on survey type # Dipole-dipole: Moving tx with [a] spacing -> [AB a MN1 a MN2 ... a MNn] # Pole-dipole: Moving pole on one end -> [A a MN1 a MN2 ... MNn a B] Tx = [] Rx = [] - + if not re.match(stype,'gradient'): - - for ii in range(0, int(nstn)-1): - - + + for ii in range(0, int(nstn)-1): + + if re.match(stype,'dpdp'): tx = np.c_[M[ii,:],N[ii,:]] elif re.match(stype,'pdp'): tx = np.c_[M[ii,:],M[ii,:]] - + #Rx.append(np.c_[M[ii+1:indx,:],N[ii+1:indx,:]]) - + # Current elctrode seperation AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1]) - + # Number of receivers to fit nstn = np.min([np.floor( (AB - b) / a ) , n]) - + # Check if there is enough space, else break the loop if nstn <= 0: continue - + # Compute discrete pole location along line stn_x = N[ii,0] + dl_x*b + np.array(range(int(nstn)))*dl_x*a stn_y = N[ii,1] + dl_y*b + np.array(range(int(nstn)))*dl_y*a - + # Create receiver poles # Create line of P1 locations P1 = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]] - + # Create line of P2 locations P2 = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]] - + Rx.append(np.c_[P1,P2]) - Tx.append(tx) - + Tx.append(tx) + #============================================================================== # elif re.match(stype,'dpdp'): -# -# for ii in range(0, int(nstn)-2): -# +# +# for ii in range(0, int(nstn)-2): +# # indx = np.min([ii+n+1,nstn]) # Tx.append(np.c_[M[ii,:],N[ii,:]]) # Rx.append(np.c_[M[ii+2:indx,:],N[ii+2:indx,:]]) #============================================================================== - + elif re.match(stype,'gradient'): - + # Gradient survey only requires Tx at end of line and creates a square # grid of receivers at in the middle at a pre-set minimum distance Tx.append(np.c_[M[0,:],N[-1,:]]) - + # Get the edge limit of survey area min_x = endl[0,0] + dl_x * b min_y = endl[0,1] + dl_y * b - + max_x = endl[1,0] - dl_x * b max_y = endl[1,1] - dl_y * b - + box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 ) box_w = box_l/2. - + nstn = np.floor( box_l / a ) - + # Compute discrete pole location along line stn_x = min_x + np.array(range(int(nstn)))*dl_x*a stn_y = min_y + np.array(range(int(nstn)))*dl_y*a - + # Define number of cross lines nlin = int(np.floor( box_w / a )) - lind = range(-nlin,nlin+1) - + lind = range(-nlin,nlin+1) + ngrad = nstn * len(lind) - + rx = np.zeros([ngrad,6]) for ii in range( len(lind) ): - + # Move line in perpendicular direction by dipole spacing lxx = stn_x - lind[ii]*a*dl_y lyy = stn_y + lind[ii]*a*dl_x - - + + M = np.c_[ lxx, lyy , np.ones(nstn).T*mesh.vectorNz[-1]] N = np.c_[ lxx+a*dl_x, lyy+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]] - + rx[(ii*nstn):((ii+1)*nstn),:] = np.c_[M,N] - + Rx.append(rx) - + else: print """stype must be either 'pdp', 'dpdp' or 'gradient'. """ - - - return Tx, Rx + + + return Tx, Rx def writeUBC_DCobs(fileName,Tx,Rx,d,wd, dtype): - + from SimPEG import np, mkvc import re """ Read UBC GIF DCIP 3D observation file and generate arrays for tx-rx location - + Input: :param fileName, path to the UBC GIF 3D obs file - + Output: :param rx, tx, d, wd :return - + Created on Mon December 7th, 2015 - + @author: dominiquef - + """ fid = open(fileName,'w') - fid.write('! GENERAL FORMAT\n') - + fid.write('! GENERAL FORMAT\n') + for ii in range(len(Tx)): - + tx = np.asarray(Tx[ii]) rx = np.asarray(Rx[ii]) nrx = rx.shape[0] - + fid.write('\n') - + if re.match(dtype,'2D'): - + for jj in range(nrx): - + fid.writelines("%e " % ii for ii in mkvc(tx)) fid.writelines("%e " % ii for ii in mkvc(rx[jj])) fid.write('%e %e\n'% (d[ii][jj],wd[ii][jj])) - #np.savetxt(fid, np.c_[ rx ,np.asarray(d[ii]), np.asarray(wd[ii]) ], fmt='%e',delimiter=' ',newline='\n') - + #np.savetxt(fid, np.c_[ rx ,np.asarray(d[ii]), np.asarray(wd[ii]) ], fmt='%e',delimiter=' ',newline='\n') + elif re.match(dtype,'3D'): - + fid.write('\n') fid.writelines("%e " % ii for ii in mkvc(tx)) fid.write('%i\n'% nrx) - np.savetxt(fid, np.c_[ rx ,np.asarray(d[ii]), np.asarray(wd[ii]) ], fmt='%e',delimiter=' ',newline='\n') - - - fid.close() + np.savetxt(fid, np.c_[ rx ,np.asarray(d[ii]), np.asarray(wd[ii]) ], fmt='%e',delimiter=' ',newline='\n') + + + fid.close() def convertObs_DC3D_to_2D(Tx,Rx): - + from SimPEG import np import numpy.matlib as npm """ @@ -622,35 +620,35 @@ def convertObs_DC3D_to_2D(Tx,Rx): First transmitter pole is assumed to be at the origin Assumes flat topo for now... - + Input: :param Tx, Rx - + Output: :figure Tx2d, Rx2d - + Created on Mon December 7th, 2015 - + @author: dominiquef - + """ - - + + Tx2d = [] Rx2d = [] for ii in range(len(Tx)): - + if ii == 0: endp = Tx[0][0:2,0] - + nrx = Rx[ii].shape[0] - + rP1 = np.sqrt( np.sum( ( endp - Tx[ii][0:2,0] )**2 , axis=0)) rP2 = np.sqrt( np.sum( ( endp - Tx[ii][0:2,1] )**2 , axis=0)) rC1 = np.sqrt( np.sum( ( npm.repmat(endp.T,nrx,1) - Rx[ii][:,0:2] )**2 , axis=1)) rC2 = np.sqrt( np.sum( ( npm.repmat(endp.T,nrx,1) - Rx[ii][:,3:5] )**2 , axis=1)) - + Tx2d.append( np.r_[rP1, rP2] ) Rx2d.append( np.c_[rC1, rC2] ) #np.savetxt(fid, data, fmt='%e',delimiter=' ',newline='\n') @@ -658,72 +656,72 @@ def convertObs_DC3D_to_2D(Tx,Rx): return Tx2d, Rx2d def readUBC_DC3Dobs(fileName): - + from SimPEG import np """ Read UBC GIF DCIP 3D observation file and generate arrays for tx-rx location - + Input: :param fileName, path to the UBC GIF 3D obs file - + Output: :param rx, tx, d, wd :return - + Created on Mon December 7th, 2015 - + @author: dominiquef - + """ - + # Load file obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!') - + # Pre-allocate Tx = [] Rx = [] d = [] wd = [] - + # Countdown for number of obs/tx count = 0 for ii in range(obsfile.shape[0]): - + if not obsfile[ii]: continue - + # First line is transmitter with number of receivers if count==0: - + temp = (np.fromstring(obsfile[ii], dtype=float,sep=' ').T) count = int(temp[-1]) temp = np.reshape(temp[0:-1],[2,3]).T - + Tx.append(temp) rx = [] continue - + temp = np.fromstring(obsfile[ii], dtype=float,sep=' ') - - - rx.append(temp) - - count = count -1 - - # Reach the end of + + + rx.append(temp) + + count = count -1 + + # Reach the end of if count == 0: temp = np.asarray(rx) Rx.append(temp[:,0:6]) - + # Check for data + uncertainties if temp.shape[1]==8: d.append(temp[:,6]) wd.append(temp[:,7]) - - # Check for data only + + # Check for data only elif temp.shape[1]==7: d.append(temp[:,6]) - + return Tx, Rx, d, wd def readUBC_DC2DLoc(fileName): @@ -738,13 +736,13 @@ def readUBC_DC2DLoc(fileName): Output: :param rx, tx :return - + Created on Thu Nov 12 13:14:10 2015 @author: dominiquef """ - + # Open fileand skip header... assume that we know the mesh already #============================================================================== # fopen = open(fileName,'r') @@ -754,32 +752,32 @@ def readUBC_DC2DLoc(fileName): # Load file obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!') - + # Check first line and figure out if 2D or 3D file format - line = np.array(obsfile[0].split(),dtype=float) - + line = np.array(obsfile[0].split(),dtype=float) + tx_A = [] tx_B = [] rx_M = [] rx_N = [] d = [] wd = [] - + for ii in range(obsfile.shape[0]): - + # If len==3, then simple format where tx-rx is listed on each line if len(line) == 4: - + temp = np.fromstring(obsfile[ii], dtype=float,sep=' ') tx_A = np.hstack((tx_A,temp[0])) tx_B = np.hstack((tx_B,temp[1])) rx_M = np.hstack((rx_M,temp[2])) rx_N = np.hstack((rx_N,temp[3])) - - - rx = np.transpose(np.array((rx_M,rx_N))) + + + rx = np.transpose(np.array((rx_M,rx_N))) tx = np.transpose(np.array((tx_A,tx_B))) - + return tx, rx, d, wd def readUBC_DC2DMesh(fileName): @@ -794,59 +792,59 @@ def readUBC_DC2DMesh(fileName): Output: :param SimPEG TensorMesh 2D object :return - + Created on Thu Nov 12 13:14:10 2015 @author: dominiquef """ - - # Open file + + # Open file fopen = open(fileName,'r') - + # Read down the file and unpack dx vector def unpackdx(fid,nrows): for ii in range(nrows): - + line = fid.readline() var = np.array(line.split(),dtype=float) - + if ii==0: x0= var[0] xvec = np.ones(int(var[2])) * (var[1] - var[0]) / int(var[2]) - xend = var[1] - + xend = var[1] + else: xvec = np.hstack((xvec,np.ones(int(var[1])) * (var[0] - xend) / int(var[1]))) - xend = var[0] - + xend = var[0] + return x0, xvec - + #%% Start with dx block # First line specifies the number of rows for x-cells line = fopen.readline() nl = np.array(line.split(),dtype=float) - [x0, dx] = unpackdx(fopen,nl) - + [x0, dx] = unpackdx(fopen,nl) + #%% Move down the file until reaching the z-block line = fopen.readline() if not line: line = fopen.readline() - + #%% End with dz block # First line specifies the number of rows for z-cells line = fopen.readline() nl = np.array(line.split(),dtype=float) - [z0, dz] = unpackdx(fopen,nl) - + [z0, dz] = unpackdx(fopen,nl) + # Flip z0 to be the bottom of the mesh for SimPEG z0 = z0 - sum(dz) dz = dz[::-1] #%% Make the mesh using SimPEG - + from SimPEG import Mesh tensMsh = Mesh.TensorMesh([dx,dz],(x0, z0)) return tensMsh diff --git a/SimPEG/Examples/DC_Forward_WennerArray.py b/SimPEG/DCIP/Utils.py similarity index 56% rename from SimPEG/Examples/DC_Forward_WennerArray.py rename to SimPEG/DCIP/Utils.py index 8a481937..702a2333 100644 --- a/SimPEG/Examples/DC_Forward_WennerArray.py +++ b/SimPEG/DCIP/Utils.py @@ -1,9 +1,8 @@ -from SimPEG import * -import SimPEG.DCIP as DC -import matplotlib.pyplot as plt +import numpy as np +def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False): -def getSrcList(nElecs, aSpacing, in2D=False, plotIt=False): + import SimPEG.DCIP as DC elocs = np.arange(0,aSpacing*nElecs,aSpacing) elocs -= (nElecs*aSpacing - aSpacing)/2 @@ -37,29 +36,3 @@ def getSrcList(nElecs, aSpacing, in2D=False, plotIt=False): srcList += [src] return srcList - - - -def run(plotIt=False,aSpacing=2.5, nElecs=10): - - surveySize = nElecs*aSpacing - aSpacing - cs = surveySize/nElecs/4 - - mesh = Mesh.TensorMesh([ - [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)], - [(cs,3, -1.3),(cs,3,1.3)], - # [(cs,5, -1.3),(cs,10)] - ],'CN') - if plotIt: - mesh.plotGrid(showIt=True) - - srcList = getSrcList(nElecs, aSpacing, in2D=True) - survey = DC.SurveyDC(srcList) - problem = DC.ProblemDC_CC(mesh) - problem.pair(survey) - - return mesh, survey, problem - - -if __name__ == '__main__': - run(plotIt=True) diff --git a/SimPEG/DCIP/__init__.py b/SimPEG/DCIP/__init__.py index c9a67453..6913124c 100644 --- a/SimPEG/DCIP/__init__.py +++ b/SimPEG/DCIP/__init__.py @@ -1,2 +1,3 @@ from BaseDC import * from BaseIP import * +import Utils diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index 1101ae83..6dadba2e 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -3,7 +3,6 @@ ##### AUTOIMPORTS ##### import DC_Analytic_Dipole import DC_Forward_PseudoSection -import DC_Forward_WennerArray import DC_PseudoSection_Simulation import EM_FDEM_1D_Inversion import EM_FDEM_Analytic_MagDipoleWholespace @@ -19,7 +18,7 @@ import Mesh_QuadTree_FaceDiv import Mesh_QuadTree_HangingNodes import Mesh_Tensor_Creation -__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "DC_Forward_WennerArray", "DC_PseudoSection_Simulation", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation"] +__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "DC_PseudoSection_Simulation", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation"] ##### AUTOIMPORTS ##### diff --git a/tests/dcip/test_forward_DCproblem.py b/tests/dcip/test_forward_DCproblem.py index 97e92383..a2708936 100644 --- a/tests/dcip/test_forward_DCproblem.py +++ b/tests/dcip/test_forward_DCproblem.py @@ -1,14 +1,28 @@ import unittest from SimPEG import * -import simpegDCIP as DC +import SimPEG.DCIP as DC class DCProblemTests(unittest.TestCase): def setUp(self): - mesh, survey, problem = DC.Examples.WennerArray.example() + aSpacing=2.5 + nElecs=10 + surveySize = nElecs*aSpacing - aSpacing + cs = surveySize/nElecs/4 + + mesh = Mesh.TensorMesh([ + [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)], + [(cs,3, -1.3),(cs,3,1.3)], + # [(cs,5, -1.3),(cs,10)] + ],'CN') + + srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True) + survey = DC.SurveyDC(srcList) + problem = DC.ProblemDC_CC(mesh) + problem.pair(survey) mSynth = np.ones(mesh.nC) survey.makeSyntheticData(mSynth) diff --git a/tests/dcip/test_forward_DCproblem_analytic.py b/tests/dcip/test_forward_DCproblem_analytic.py deleted file mode 100644 index 15493115..00000000 --- a/tests/dcip/test_forward_DCproblem_analytic.py +++ /dev/null @@ -1,12 +0,0 @@ -import unittest -import simpegDCIP as DC - - -class DCAnalyticTests(unittest.TestCase): - - def test_forwardAnalytic(self): - self.assertTrue(DC.Examples.Verification.run() < 0.1) - - -if __name__ == '__main__': - unittest.main() diff --git a/tests/dcip/test_forward_IPproblem.py b/tests/dcip/test_forward_IPproblem.py index 2a64ee7b..261becd2 100644 --- a/tests/dcip/test_forward_IPproblem.py +++ b/tests/dcip/test_forward_IPproblem.py @@ -1,56 +1,64 @@ import unittest -import simpegDCIP as DC +import SimPEG.DCIP as DC from SimPEG import * -from pymatsolver import MumpsSolver class IPforwardTests(unittest.TestCase): def test_IPforward(self): - cs = 12.5 - nc = 500/cs+1 - hx = [(cs,7, -1.3),(cs,nc),(cs,7, 1.3)] - hy = [(cs,7, -1.3),(cs,int(nc/2+1)),(cs,7, 1.3)] - hz = [(cs,7, -1.3),(cs,int(nc/2+1))] - mesh = Mesh.TensorMesh([hx, hy, hz], 'CCN') - sighalf = 1e-2 - sigma = np.ones(mesh.nC)*sighalf - p0 = np.r_[-50., 50., -50.] - p1 = np.r_[ 50.,-50., -150.] - blk_ind = Utils.ModelBuilder.getIndicesBlock(p0, p1, mesh.gridCC) - sigma[blk_ind] = 1e-3 - eta = np.zeros_like(sigma) - eta[blk_ind] = 0.1 - sigmaInf = sigma.copy() - sigma0 = sigma*(1-eta) + cs = 12.5 + nc = 500/cs+1 + hx = [(cs,7, -1.3),(cs,nc),(cs,7, 1.3)] + hy = [(cs,7, -1.3),(cs,int(nc/2+1)),(cs,7, 1.3)] + hz = [(cs,7, -1.3),(cs,int(nc/2+1))] + mesh = Mesh.TensorMesh([hx, hy, hz], 'CCN') + sighalf = 1e-2 + sigma = np.ones(mesh.nC)*sighalf + p0 = np.r_[-50., 50., -50.] + p1 = np.r_[ 50.,-50., -150.] + blk_ind = Utils.ModelBuilder.getIndicesBlock(p0, p1, mesh.gridCC) + sigma[blk_ind] = 1e-3 + eta = np.zeros_like(sigma) + eta[blk_ind] = 0.1 + sigmaInf = sigma.copy() + sigma0 = sigma*(1-eta) - nElecs = 11 - x_temp = np.linspace(-250, 250, nElecs) - aSpacing = x_temp[1]-x_temp[0] - y_temp = 0. - xyz = Utils.ndgrid(x_temp, np.r_[y_temp], np.r_[0.]) - srcList = DC.Examples.WennerArray.getSrcList(nElecs,aSpacing) - survey = DC.SurveyDC(srcList) + nElecs = 11 + x_temp = np.linspace(-250, 250, nElecs) + aSpacing = x_temp[1]-x_temp[0] + y_temp = 0. + xyz = Utils.ndgrid(x_temp, np.r_[y_temp], np.r_[0.]) + srcList = DC.Utils.WennerSrcList(nElecs,aSpacing) + survey = DC.SurveyDC(srcList) - imap = Maps.IdentityMap(mesh) - problem = DC.ProblemDC_CC(mesh, mapping= imap ) - problem.Solver = MumpsSolver - problem.pair(survey) + imap = Maps.IdentityMap(mesh) + problem = DC.ProblemDC_CC(mesh, mapping=imap) + try: + from pymatsolver import MumpsSolver + problem.Solver = MumpsSolver + except Exception, e: + pass + problem.pair(survey) - phi0 = survey.dpred(sigma0) - phiInf = survey.dpred(sigmaInf) + phi0 = survey.dpred(sigma0) + phiInf = survey.dpred(sigmaInf) - phiIP_true = phi0-phiInf + phiIP_true = phi0-phiInf - surveyIP = DC.SurveyIP(srcList) - problemIP = DC.ProblemIP(mesh, sigma=sigma) - problemIP.pair(surveyIP) - problemIP.Solver = MumpsSolver - phiIP_approx = surveyIP.dpred(eta) + surveyIP = DC.SurveyIP(srcList) + problemIP = DC.ProblemIP(mesh, sigma=sigma) + problemIP.pair(surveyIP) - err = np.linalg.norm(phiIP_true-phiIP_approx) / np.linalg.norm(phiIP_true) + try: + from pymatsolver import MumpsSolver + problemIP.Solver = MumpsSolver + except Exception, e: + pass + phiIP_approx = surveyIP.dpred(eta) - self.assertTrue(err < 0.02) + err = np.linalg.norm(phiIP_true-phiIP_approx) / np.linalg.norm(phiIP_true) + + self.assertTrue(err < 0.02) if __name__ == '__main__': diff --git a/tests/dcip/test_sens_IPproblem.py b/tests/dcip/test_sens_IPproblem.py index 2145f477..5f158ed3 100644 --- a/tests/dcip/test_sens_IPproblem.py +++ b/tests/dcip/test_sens_IPproblem.py @@ -1,7 +1,7 @@ import unittest from SimPEG import * -import simpegDCIP as DC -from pymatsolver import MumpsSolver +import SimPEG.DCIP as DC + @@ -29,12 +29,17 @@ class IPProblemTests(unittest.TestCase): aSpacing = x_temp[1]-x_temp[0] y_temp = 0. xyz = Utils.ndgrid(x_temp, np.r_[y_temp], np.r_[0.]) - srcList = DC.Examples.WennerArray.getSrcList(nElecs,aSpacing) + srcList = DC.Utils.WennerSrcList(nElecs,aSpacing) survey = DC.SurveyIP(srcList) imap = Maps.IdentityMap(mesh) problem = DC.ProblemIP(mesh, sigma=sigma, mapping= imap) problem.pair(survey) - problem.Solver = MumpsSolver + + try: + from pymatsolver import MumpsSolver + problem.Solver = MumpsSolver + except Exception, e: + pass mSynth = eta survey.makeSyntheticData(mSynth)