diff --git a/SimPEG/Examples/DC_Analytic_Dipole.py b/SimPEG/Examples/DC_Analytic_Dipole.py index 239f88d9..aed5eed4 100644 --- a/SimPEG/Examples/DC_Analytic_Dipole.py +++ b/SimPEG/Examples/DC_Analytic_Dipole.py @@ -1,7 +1,5 @@ from SimPEG import * import SimPEG.DCIP as DC -import matplotlib.pyplot as plt - def run(plotIt=False): cs = 25. @@ -51,6 +49,7 @@ def run(plotIt=False): Y = xyz_rxM[:,1].reshape((21, 21), order = 'F') if plotIt: + import matplotlib.pyplot as plt fig, ax = plt.subplots(1,2, figsize = (12, 5)) vmin = np.r_[data, data_ana].min() vmax = np.r_[data, data_ana].max() diff --git a/SimPEG/Examples/DC_Forward_PseudoSection.py b/SimPEG/Examples/DC_Forward_PseudoSection.py index f7fdcbaa..da69fe9a 100644 --- a/SimPEG/Examples/DC_Forward_PseudoSection.py +++ b/SimPEG/Examples/DC_Forward_PseudoSection.py @@ -1,21 +1,31 @@ -from SimPEG import * +from SimPEG import Mesh, Utils, np, sp import SimPEG.DCIP as DC -import scipy.interpolate as interpolation -import matplotlib.pyplot as plt import time -import re -def run(loc=np.c_[[-50.,0.,-50.],[50.,0.,-50.]], sig=np.r_[1e-2,1e-1,1e-3], radi=np.r_[25.,25.], param = np.r_[30.,30.,5], stype = 'dpdp', plotIt=True): +def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True): """ - DC Forward Simulation + DC Forward Simulation + ===================== - Forward model conductive spheres in a half-space and plot a pseudo-section + Forward model conductive spheres in a half-space and plot a pseudo-section - Created on Mon Feb 01 19:28:06 2016 + Created by @fourndo on Mon Feb 01 19:28:06 2016 - @fourndo """ + assert stype in ['pdp', 'dpdp'], "Source type (stype) must be pdp or dpdp (pole dipole or dipole dipole)" + + + if loc is None: + loc = np.c_[[-50.,0.,-50.],[50.,0.,-50.]] + if sig is None: + sig = np.r_[1e-2,1e-1,1e-3] + if radi is None: + radi = np.r_[25.,25.] + if param is None: + param = np.r_[30.,30.,5] + + # First we need to create a mesh and a model. # This is our mesh @@ -104,24 +114,21 @@ def run(loc=np.c_[[-50.,0.,-50.],[50.,0.,-50.]], sig=np.r_[1e-2,1e-1,1e-3], radi # For usual cases "dpdp" or "gradient" - if not re.match(stype,'pdp'): - inds = Utils.closestPoints(mesh, np.asarray(Tx[ii]).T ) - RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1,1] / mesh.vol[inds] ) - - else: - + if stype == 'pdp': # Create an "inifinity" pole tx = np.squeeze(Tx[ii][:,0:1]) tinf = tx + np.array([dl_x,dl_y,0])*dl_len*2 inds = Utils.closestPoints(mesh, np.c_[tx,tinf].T) RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1] / mesh.vol[inds] ) - + else: + inds = Utils.closestPoints(mesh, np.asarray(Tx[ii]).T ) + RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1,1] / mesh.vol[inds] ) # Iterative Solve Ainvb = sp.linalg.bicgstab(P*A,P*RHS, tol=1e-5) # We now have the potential everywhere - phi = mkvc(Ainvb[0]) + phi = Utils.mkvc(Ainvb[0]) # Solve for phi on pole locations P1 = mesh.getInterpolationMat(rxloc_M, 'CC') @@ -143,6 +150,7 @@ def run(loc=np.c_[[-50.,0.,-50.],[50.,0.,-50.]], sig=np.r_[1e-2,1e-1,1e-3], radi # Here is an example for the first tx-rx array if plotIt: + import matplotlib.pyplot as plt fig = plt.figure() ax = plt.subplot(2,1,1, aspect='equal') mesh.plotSlice(np.log10(model), ax =ax, normal = 'Y', ind = indy,grid=True) diff --git a/SimPEG/Examples/DC_PseudoSection_Simulation.py b/SimPEG/Examples/DC_PseudoSection_Simulation.py deleted file mode 100644 index 6a27de96..00000000 --- a/SimPEG/Examples/DC_PseudoSection_Simulation.py +++ /dev/null @@ -1,187 +0,0 @@ -from SimPEG import Mesh, Utils, np, sp -import SimPEG.DCIP as DC -import time - -def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True): - """ - DC Forward Simulation - - Forward model conductive spheres in a half-space and plot a pseudo-section - - Created on Mon Feb 01 19:28:06 2016 - - @fourndo - """ - - assert stype in ['pdp', 'dpdp'], "Source type (stype) must be pdp or dpdp (pole dipole or dipole dipole)" - - - if loc is None: - loc = np.c_[[-50.,0.,-50.],[50.,0.,-50.]] - if sig is None: - sig = np.r_[1e-2,1e-1,1e-3] - if radi is None: - radi = np.r_[25.,25.] - if param is None: - param = np.r_[30.,30.,5] - - - # First we need to create a mesh and a model. - - # This is our mesh - dx = 5. - - hxind = [(dx,15,-1.3), (dx, 75), (dx,15,1.3)] - hyind = [(dx,15,-1.3), (dx, 10), (dx,15,1.3)] - hzind = [(dx,15,-1.3),(dx, 15)] - - mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCN') - - - # Set background conductivity - model = np.ones(mesh.nC) * sig[0] - - # First anomaly - ind = Utils.ModelBuilder.getIndicesSphere(loc[:,0],radi[0],mesh.gridCC) - model[ind] = sig[1] - - # Second anomaly - ind = Utils.ModelBuilder.getIndicesSphere(loc[:,1],radi[1],mesh.gridCC) - model[ind] = sig[2] - - # Get index of the center - indy = int(mesh.nCy/2) - - - # Plot the model for reference - # Define core mesh extent - xlim = 200 - zlim = 125 - - # Specify the survey type: "pdp" | "dpdp" - - - # Then specify the end points of the survey. Let's keep it simple for now and survey above the anomalies, top of the mesh - ends = [(-175,0),(175,0)] - ends = np.c_[np.asarray(ends),np.ones(2).T*mesh.vectorNz[-1]] - - # Snap the endpoints to the grid. Easier to create 2D section. - indx = Utils.closestPoints(mesh, ends ) - locs = np.c_[mesh.gridCC[indx,0],mesh.gridCC[indx,1],np.ones(2).T*mesh.vectorNz[-1]] - - # We will handle the geometry of the survey for you and create all the combination of tx-rx along line - [Tx, Rx] = DC.gen_DCIPsurvey(locs, mesh, stype, param[0], param[1], param[2]) - - # Define some global geometry - dl_len = np.sqrt( np.sum((locs[0,:] - locs[1,:])**2) ) - dl_x = ( Tx[-1][0,1] - Tx[0][0,0] ) / dl_len - dl_y = ( Tx[-1][1,1] - Tx[0][1,0] ) / dl_len - azm = np.arctan(dl_y/dl_x) - - #Set boundary conditions - mesh.setCellGradBC('neumann') - - # Define the differential operators needed for the DC problem - Div = mesh.faceDiv - Grad = mesh.cellGrad - Msig = Utils.sdiag(1./(mesh.aveF2CC.T*(1./model))) - - A = Div*Msig*Grad - - # Change one corner to deal with nullspace - A[0,0] = 1 - A = sp.csc_matrix(A) - - # We will solve the system iteratively, so a pre-conditioner is helpful - # This is simply a Jacobi preconditioner (inverse of the main diagonal) - dA = A.diagonal() - P = sp.spdiags(1/dA,0,A.shape[0],A.shape[0]) - - # Now we can solve the system for all the transmitters - # We want to store the data - data = [] - - # There is probably a more elegant way to do this, but we can just for-loop through the transmitters - for ii in range(len(Tx)): - - start_time = time.time() # Let's time the calculations - - #print("Transmitter %i / %i\r" % (ii+1,len(Tx))) - - # Select dipole locations for receiver - rxloc_M = np.asarray(Rx[ii][:,0:3]) - rxloc_N = np.asarray(Rx[ii][:,3:]) - - - # For usual cases "dpdp" or "gradient" - if stype == 'pdp': - # Create an "inifinity" pole - tx = np.squeeze(Tx[ii][:,0:1]) - tinf = tx + np.array([dl_x,dl_y,0])*dl_len*2 - inds = Utils.closestPoints(mesh, np.c_[tx,tinf].T) - RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1] / mesh.vol[inds] ) - else: - inds = Utils.closestPoints(mesh, np.asarray(Tx[ii]).T ) - RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1,1] / mesh.vol[inds] ) - - # Iterative Solve - Ainvb = sp.linalg.bicgstab(P*A,P*RHS, tol=1e-5) - - # We now have the potential everywhere - phi = Utils.mkvc(Ainvb[0]) - - # Solve for phi on pole locations - P1 = mesh.getInterpolationMat(rxloc_M, 'CC') - P2 = mesh.getInterpolationMat(rxloc_N, 'CC') - - # Compute the potential difference - dtemp = (P1*phi - P2*phi)*np.pi - - data.append( dtemp ) - print '\rTransmitter {0} of {1} -> Time:{2} sec'.format(ii,len(Tx),time.time()- start_time), - - print 'Transmitter {0} of {1}'.format(ii,len(Tx)) - print 'Forward completed' - - - # Let's just convert the 3D format into 2D (distance along line) and plot - [Tx2d, Rx2d] = DC.convertObs_DC3D_to_2D(Tx,Rx) - - - # Here is an example for the first tx-rx array - if plotIt: - import matplotlib.pyplot as plt - fig = plt.figure() - ax = plt.subplot(2,1,1, aspect='equal') - mesh.plotSlice(np.log10(model), ax =ax, normal = 'Y', ind = indy,grid=True) - ax.set_title('E-W section at '+str(mesh.vectorCCy[indy])+' m') - plt.gca().set_aspect('equal', adjustable='box') - - plt.scatter(Tx[0][0,:],Tx[0][2,:],s=40,c='g', marker='v') - plt.scatter(Rx[0][:,0::3],Rx[0][:,2::3],s=40,c='y') - plt.xlim([-xlim,xlim]) - plt.ylim([-zlim,mesh.vectorNz[-1]+dx]) - - - ax = plt.subplot(2,1,2, aspect='equal') - - # Plot the location of the spheres for reference - circle1=plt.Circle((loc[0,0]-Tx[0][0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3) - circle2=plt.Circle((loc[0,1]-Tx[0][0,0],loc[2,1]),radi[1],color='k',fill=False, lw=3) - ax.add_artist(circle1) - ax.add_artist(circle2) - - # Add the speudo section - DC.plot_pseudoSection(Tx2d,Rx2d,data,mesh.vectorNz[-1],stype) - - plt.scatter(Tx2d[0][:],Tx[0][2,:],s=40,c='g', marker='v') - plt.scatter(Rx2d[0][:],Rx[0][:,2::3],s=40,c='y') - plt.plot(np.r_[Tx2d[0][0],Rx2d[-1][-1,-1]],np.ones(2)*mesh.vectorNz[-1], color='k') - plt.ylim([-zlim,mesh.vectorNz[-1]+dx]) - - plt.show() - - return fig, ax - -if __name__ == '__main__': - run() diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index 6dadba2e..9c231068 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_PseudoSection_Simulation import EM_FDEM_1D_Inversion import EM_FDEM_Analytic_MagDipoleWholespace import EM_TDEM_1D_Inversion @@ -18,7 +17,7 @@ import Mesh_QuadTree_FaceDiv import Mesh_QuadTree_HangingNodes import 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"] +__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "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/docs/examples/DC_Forward_PseudoSection.rst b/docs/examples/DC_Forward_PseudoSection.rst index 87c9fd5a..1a500cae 100644 --- a/docs/examples/DC_Forward_PseudoSection.rst +++ b/docs/examples/DC_Forward_PseudoSection.rst @@ -9,13 +9,13 @@ .. --------------------------------- .. -orward Simulation +DC Forward Simulation +===================== -ard model conductive spheres in a half-space and plot a pseudo-section +Forward model conductive spheres in a half-space and plot a pseudo-section -ted on Mon Feb 01 19:28:06 2016 +Created by @fourndo on Mon Feb 01 19:28:06 2016 -rndo .. plot:: diff --git a/docs/examples/DC_Forward_WennerArra.rst b/docs/examples/DC_Forward_WennerArra.rst deleted file mode 100644 index 04ed685a..00000000 --- a/docs/examples/DC_Forward_WennerArra.rst +++ /dev/null @@ -1,21 +0,0 @@ -.. _examples_DC_Forward_WennerArra: - -.. --------------------------------- .. -.. .. -.. THIS FILE IS AUTO GENEREATED .. -.. .. -.. SimPEG/Examples/__init__.py .. -.. .. -.. --------------------------------- .. - -DC Forward WennerArra -===================== - -.. plot:: - - from SimPEG import Examples - Examples.DC_Forward_WennerArra.run() - -.. literalinclude:: ../../SimPEG/Examples/DC_Forward_WennerArra.py - :language: python - :linenos: diff --git a/docs/examples/DC_PseudoSection_Simulation.rst b/docs/examples/DC_PseudoSection_Simulation.rst deleted file mode 100644 index 8f1c3ce1..00000000 --- a/docs/examples/DC_PseudoSection_Simulation.rst +++ /dev/null @@ -1,28 +0,0 @@ -.. _examples_DC_PseudoSection_Simulation: - -.. --------------------------------- .. -.. .. -.. THIS FILE IS AUTO GENEREATED .. -.. .. -.. SimPEG/Examples/__init__.py .. -.. .. -.. --------------------------------- .. - - -orward Simulation - -ard model conductive spheres in a half-space and plot a pseudo-section - -ted on Mon Feb 01 19:28:06 2016 - -rndo - - -.. plot:: - - from SimPEG import Examples - Examples.DC_PseudoSection_Simulation.run() - -.. literalinclude:: ../../SimPEG/Examples/DC_PseudoSection_Simulation.py - :language: python - :linenos: