From df78f7b33aca3a0c77d202d1bd9c11b99f9aa64b Mon Sep 17 00:00:00 2001 From: D Fournier Date: Wed, 3 Feb 2016 14:35:03 -0800 Subject: [PATCH 01/18] Branch off master Add DC_Pseudo_Section example. --- .../Examples/DC_PseudoSection_Simulation.py | 179 ++++++++++++++++++ 1 file changed, 179 insertions(+) create mode 100644 SimPEG/Examples/DC_PseudoSection_Simulation.py diff --git a/SimPEG/Examples/DC_PseudoSection_Simulation.py b/SimPEG/Examples/DC_PseudoSection_Simulation.py new file mode 100644 index 00000000..c3517360 --- /dev/null +++ b/SimPEG/Examples/DC_PseudoSection_Simulation.py @@ -0,0 +1,179 @@ +from SimPEG import * +import simpegDCIP 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): + """ + 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 + """ + + # 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 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: + + # 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] ) + + + # Iterative Solve + Ainvb = sp.linalg.bicgstab(P*A,P*RHS, tol=1e-5) + + # We now have the potential everywhere + phi = 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: + 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() \ No newline at end of file From 3eeb5dbd3cbd33c941b8821e814e3f7f78f7af67 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Wed, 3 Feb 2016 20:03:05 -0800 Subject: [PATCH 02/18] Remove dependency from Utils. Replace by internal function. --- .../Examples/DC_PseudoSection_Simulation.py | 41 ++++++++++++++++++- 1 file changed, 39 insertions(+), 2 deletions(-) diff --git a/SimPEG/Examples/DC_PseudoSection_Simulation.py b/SimPEG/Examples/DC_PseudoSection_Simulation.py index c3517360..e7c73474 100644 --- a/SimPEG/Examples/DC_PseudoSection_Simulation.py +++ b/SimPEG/Examples/DC_PseudoSection_Simulation.py @@ -16,6 +16,43 @@ def run(loc=np.c_[[-50.,0.,-50.],[50.,0.,-50.]], sig=np.r_[1e-2,1e-1,1e-3], radi @fourndo """ + def getIndicesSphere(center,radius,ccMesh): + """ + Creates a vector containing the sphere indices in the cell centers mesh. + Returns a tuple + + The sphere is defined by the points + + p0, describe the position of the center of the cell + + r, describe the radius of the sphere. + + ccMesh represents the cell-centered mesh + + The points p0 must live in the the same dimensional space as the mesh. + + """ + + # Validation: mesh and point (p0) live in the same dimensional space + dimMesh = np.size(ccMesh[0,:]) + assert len(center) == dimMesh, "Dimension mismatch. len(p0) != dimMesh" + + if dimMesh == 1: + # Define the reference points + + ind = np.abs(center[0] - ccMesh[:,0]) < radius + + elif dimMesh == 2: + # Define the reference points + + ind = np.sqrt( ( center[0] - ccMesh[:,0] )**2 + ( center[1] - ccMesh[:,1] )**2 ) < radius + + elif dimMesh == 3: + # Define the points + ind = np.sqrt( ( center[0] - ccMesh[:,0] )**2 + ( center[1] - ccMesh[:,1] )**2 + ( center[2] - ccMesh[:,2] )**2 ) < radius + + # Return a tuple + return ind # First we need to create a mesh and a model. # This is our mesh @@ -32,11 +69,11 @@ def run(loc=np.c_[[-50.,0.,-50.],[50.,0.,-50.]], sig=np.r_[1e-2,1e-1,1e-3], radi model = np.ones(mesh.nC) * sig[0] # First anomaly - ind = Utils.ModelBuilder.getIndicesSphere(loc[:,0],radi[0],mesh.gridCC) + ind = getIndicesSphere(loc[:,0],radi[0],mesh.gridCC) model[ind] = sig[1] # Second anomaly - ind = Utils.ModelBuilder.getIndicesSphere(loc[:,1],radi[1],mesh.gridCC) + ind = getIndicesSphere(loc[:,1],radi[1],mesh.gridCC) model[ind] = sig[2] # Get index of the center From 433457f6497cbf926828d5afd80609ed3793bd9e Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 3 Feb 2016 20:12:55 -0800 Subject: [PATCH 03/18] docs and import for DC_PseudoSection_Simulation.rst --- SimPEG/Examples/__init__.py | 3 +- docs/examples/DC_PseudoSection_Simulation.rst | 28 +++++++++++++++++++ 2 files changed, 30 insertions(+), 1 deletion(-) create mode 100644 docs/examples/DC_PseudoSection_Simulation.rst diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index 7ef15451..4fe604ab 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -1,6 +1,7 @@ # Run this file to add imports. ##### AUTOIMPORTS ##### +import DC_PseudoSection_Simulation import EM_FDEM_Analytic_MagDipoleWholespace import EM_TDEM_1D_Inversion import FLOW_Richards_1D_Celia1990 @@ -14,7 +15,7 @@ import Mesh_QuadTree_FaceDiv import Mesh_QuadTree_HangingNodes import Mesh_Tensor_Creation -__examples__ = ["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_PseudoSection_Simulation", "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_PseudoSection_Simulation.rst b/docs/examples/DC_PseudoSection_Simulation.rst new file mode 100644 index 00000000..8f1c3ce1 --- /dev/null +++ b/docs/examples/DC_PseudoSection_Simulation.rst @@ -0,0 +1,28 @@ +.. _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: From 429d8b119117e17daa6c63592c87ea7ec26ccba3 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Wed, 10 Feb 2016 23:45:19 -0800 Subject: [PATCH 04/18] Add EM_FDEM_SusEffects example --- SimPEG/Examples/__init__.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index 4fe604ab..b821e10b 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -1,6 +1,7 @@ # Run this file to add imports. ##### AUTOIMPORTS ##### +import EM_FDEM_SusEffects import DC_PseudoSection_Simulation import EM_FDEM_Analytic_MagDipoleWholespace import EM_TDEM_1D_Inversion @@ -15,7 +16,8 @@ import Mesh_QuadTree_FaceDiv import Mesh_QuadTree_HangingNodes import Mesh_Tensor_Creation -__examples__ = ["DC_PseudoSection_Simulation", "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__ = ["EM_FDEM_SusEffects","DC_PseudoSection_Simulation", "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 ##### From 47895ef2701ebfbae59c344b2addcaa846c2ea97 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Wed, 10 Feb 2016 23:58:48 -0800 Subject: [PATCH 05/18] consistent file name --- SimPEG/Examples/EM_FDEM_SusEffects.py | 7 +++-- SimPEG/Examples/__init__.py | 6 ++-- docs/examples/EM_FDEM_1D_Inversion.rst | 26 ---------------- docs/examples/EM_FDEM_SusEffects.rst | 41 ++++++++++++++++++++++++++ 4 files changed, 48 insertions(+), 32 deletions(-) delete mode 100644 docs/examples/EM_FDEM_1D_Inversion.rst create mode 100644 docs/examples/EM_FDEM_SusEffects.rst diff --git a/SimPEG/Examples/EM_FDEM_SusEffects.py b/SimPEG/Examples/EM_FDEM_SusEffects.py index 5021e87f..b7de07bc 100644 --- a/SimPEG/Examples/EM_FDEM_SusEffects.py +++ b/SimPEG/Examples/EM_FDEM_SusEffects.py @@ -5,7 +5,7 @@ from scipy.constants import mu_0 def run(plotIt=True): """ - FDEM: Effects of susceptibility + EM: FDEM: Effects of susceptibility =============================== When airborne freqeuncy domain EM (AFEM) survey is flown over @@ -138,8 +138,9 @@ def run(plotIt=True): axtemp.set_xlim(10, 100.) axtemp.set_title(titles[i]) plt.show() - vizfields(1, primsec="primary", realimag="real") - vizfields(1, primsec="secondary", realimag="real") + return fig, ax + fig1, ax1 = vizfields(1, primsec="primary", realimag="real") + fig2, ax2 = vizfields(1, primsec="secondary", realimag="real") if __name__ == '__main__': run() diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index b821e10b..df65e16b 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -1,9 +1,10 @@ # Run this file to add imports. ##### AUTOIMPORTS ##### -import EM_FDEM_SusEffects import DC_PseudoSection_Simulation +import EM_FDEM_1D_Inversion import EM_FDEM_Analytic_MagDipoleWholespace +import EM_FDEM_SusEffects import EM_TDEM_1D_Inversion import FLOW_Richards_1D_Celia1990 import Forward_BasicDirectCurrent @@ -16,8 +17,7 @@ import Mesh_QuadTree_FaceDiv import Mesh_QuadTree_HangingNodes import Mesh_Tensor_Creation - -__examples__ = ["EM_FDEM_SusEffects","DC_PseudoSection_Simulation", "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_PseudoSection_Simulation", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_FDEM_SusEffects", "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/EM_FDEM_1D_Inversion.rst b/docs/examples/EM_FDEM_1D_Inversion.rst deleted file mode 100644 index acbc8cdc..00000000 --- a/docs/examples/EM_FDEM_1D_Inversion.rst +++ /dev/null @@ -1,26 +0,0 @@ -.. _examples_EM_FDEM_1D_Inversion: - -.. --------------------------------- .. -.. .. -.. THIS FILE IS AUTO GENEREATED .. -.. .. -.. SimPEG/Examples/__init__.py .. -.. .. -.. --------------------------------- .. - - -EM: FDEM: 1D: Inversion -======================= - -Here we will create and run a FDEM 1D inversion. - - - -.. plot:: - - from SimPEG import Examples - Examples.EM_FDEM_1D_Inversion.run() - -.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_1D_Inversion.py - :language: python - :linenos: diff --git a/docs/examples/EM_FDEM_SusEffects.rst b/docs/examples/EM_FDEM_SusEffects.rst new file mode 100644 index 00000000..47cd3644 --- /dev/null +++ b/docs/examples/EM_FDEM_SusEffects.rst @@ -0,0 +1,41 @@ +.. _examples_EM_FDEM_SusEffects: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +FDEM: Effects of susceptibility +=============================== + +When airborne freqeuncy domain EM (AFEM) survey is flown over +the earth including significantly susceptible bodies (magnetite-rich rocks), +negative data is often observed in the real part of the lowest frequency +(e.g. Dighem system 900 Hz). This phenomenon mostly based upon magnetization +occurs due to a susceptible body when the magnetic field applied. + +To clarify what is happening in the earth when we are exciting the earth with +a loop source in the frequency domain we run three forward modelling: + + - F[:math:`\sigma`, :math:`\mu`]: Anomalous conductivity and susceptibility + - F[:math:`\sigma`, :math:`\mu_0`]: Anomalous conductivity + - F[:math:`\sigma_{air}`, :math:`\mu_0`]: primary field + +We plot vector magnetic fields in the earth. For secondary fields we provide +F[:math:`\sigma`, :math:`\mu`]-F[:math:`\sigma`, :math:`\mu_0`]. Following +figure show only real part, since that is our interest. + + + +.. plot:: + + from SimPEG import Examples + Examples.EM_FDEM_SusEffects.run() + +.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_SusEffects.py + :language: python + :linenos: From fa6bcd3ffa263f62c9e7eebd2f00afba0c041bbd Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Thu, 11 Feb 2016 00:10:25 -0800 Subject: [PATCH 06/18] Minor changes --- SimPEG/Examples/EM_FDEM_SusEffects.py | 6 ++++-- docs/examples/EM_FDEM_1D_Inversion.rst | 26 ++++++++++++++++++++++++++ docs/examples/EM_FDEM_SusEffects.rst | 6 +++--- 3 files changed, 33 insertions(+), 5 deletions(-) create mode 100644 docs/examples/EM_FDEM_1D_Inversion.rst diff --git a/SimPEG/Examples/EM_FDEM_SusEffects.py b/SimPEG/Examples/EM_FDEM_SusEffects.py index b7de07bc..3e50bf91 100644 --- a/SimPEG/Examples/EM_FDEM_SusEffects.py +++ b/SimPEG/Examples/EM_FDEM_SusEffects.py @@ -12,7 +12,7 @@ def run(plotIt=True): the earth including significantly susceptible bodies (magnetite-rich rocks), negative data is often observed in the real part of the lowest frequency (e.g. Dighem system 900 Hz). This phenomenon mostly based upon magnetization - occurs due to a susceptible body when the magnetic field applied. + occurs due to a susceptible body when the magnetic field is applied. To clarify what is happening in the earth when we are exciting the earth with a loop source in the frequency domain we run three forward modelling: @@ -23,7 +23,7 @@ def run(plotIt=True): We plot vector magnetic fields in the earth. For secondary fields we provide F[:math:`\sigma`, :math:`\mu`]-F[:math:`\sigma`, :math:`\mu_0`]. Following - figure show only real part, since that is our interest. + figure show both real and parts. """ # Generate Cylindrical mesh @@ -128,6 +128,7 @@ def run(plotIt=True): ax = [ax1, ax2, ax3] ax3.text(30, 50, ("Frequency=%5.2f Hz")%(frequency[ifreq]), color="k", fontsize=18) ax2.text(30, 50, primsec, color="k", fontsize=18) + ax1.text(30, 50, realimag, color="k", fontsize=18) for i, axtemp in enumerate(ax): axtemp.plot(np.r_[0, 29.75], np.r_[-50, -50], 'w', lw=3) @@ -141,6 +142,7 @@ def run(plotIt=True): return fig, ax fig1, ax1 = vizfields(1, primsec="primary", realimag="real") fig2, ax2 = vizfields(1, primsec="secondary", realimag="real") + fig4, ax4 = vizfields(1, primsec="secondary", realimag="imag") if __name__ == '__main__': run() diff --git a/docs/examples/EM_FDEM_1D_Inversion.rst b/docs/examples/EM_FDEM_1D_Inversion.rst new file mode 100644 index 00000000..acbc8cdc --- /dev/null +++ b/docs/examples/EM_FDEM_1D_Inversion.rst @@ -0,0 +1,26 @@ +.. _examples_EM_FDEM_1D_Inversion: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +EM: FDEM: 1D: Inversion +======================= + +Here we will create and run a FDEM 1D inversion. + + + +.. plot:: + + from SimPEG import Examples + Examples.EM_FDEM_1D_Inversion.run() + +.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_1D_Inversion.py + :language: python + :linenos: diff --git a/docs/examples/EM_FDEM_SusEffects.rst b/docs/examples/EM_FDEM_SusEffects.rst index 47cd3644..da4b00a0 100644 --- a/docs/examples/EM_FDEM_SusEffects.rst +++ b/docs/examples/EM_FDEM_SusEffects.rst @@ -9,14 +9,14 @@ .. --------------------------------- .. -FDEM: Effects of susceptibility +EM: FDEM: Effects of susceptibility =============================== When airborne freqeuncy domain EM (AFEM) survey is flown over the earth including significantly susceptible bodies (magnetite-rich rocks), negative data is often observed in the real part of the lowest frequency (e.g. Dighem system 900 Hz). This phenomenon mostly based upon magnetization -occurs due to a susceptible body when the magnetic field applied. +occurs due to a susceptible body when the magnetic field is applied. To clarify what is happening in the earth when we are exciting the earth with a loop source in the frequency domain we run three forward modelling: @@ -27,7 +27,7 @@ a loop source in the frequency domain we run three forward modelling: We plot vector magnetic fields in the earth. For secondary fields we provide F[:math:`\sigma`, :math:`\mu`]-F[:math:`\sigma`, :math:`\mu_0`]. Following -figure show only real part, since that is our interest. +figure show both real and parts. From 3a506b9051cdfcb449f05a725a3c1505568e6059 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Thu, 11 Feb 2016 08:54:44 -0800 Subject: [PATCH 07/18] Fix DC examples ... --- .../Examples/DC_PseudoSection_Simulation.py | 154 +++++++++--------- SimPEG/Examples/EM_FDEM_SusEffects.py | 2 +- docs/examples/DC_PseudoSection_Simulation.rst | 11 +- docs/examples/EM_FDEM_SusEffects.rst | 2 +- 4 files changed, 88 insertions(+), 81 deletions(-) diff --git a/SimPEG/Examples/DC_PseudoSection_Simulation.py b/SimPEG/Examples/DC_PseudoSection_Simulation.py index e7c73474..3b9c4c05 100644 --- a/SimPEG/Examples/DC_PseudoSection_Simulation.py +++ b/SimPEG/Examples/DC_PseudoSection_Simulation.py @@ -1,5 +1,6 @@ from SimPEG import * -import simpegDCIP as DC +# import simpegDCIP as DC +from SimPEG import DCIP as DC import scipy.interpolate as interpolation import matplotlib.pyplot as plt import time @@ -7,118 +8,121 @@ 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): """ - 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 + + 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 + """ - + def getIndicesSphere(center,radius,ccMesh): """ Creates a vector containing the sphere indices in the cell centers mesh. Returns a tuple - + The sphere is defined by the points - + p0, describe the position of the center of the cell - + r, describe the radius of the sphere. - + ccMesh represents the cell-centered mesh - + The points p0 must live in the the same dimensional space as the mesh. - + """ - + # Validation: mesh and point (p0) live in the same dimensional space dimMesh = np.size(ccMesh[0,:]) assert len(center) == dimMesh, "Dimension mismatch. len(p0) != dimMesh" - + if dimMesh == 1: # Define the reference points - + ind = np.abs(center[0] - ccMesh[:,0]) < radius - + elif dimMesh == 2: # Define the reference points - + ind = np.sqrt( ( center[0] - ccMesh[:,0] )**2 + ( center[1] - ccMesh[:,1] )**2 ) < radius - + elif dimMesh == 3: # Define the points ind = np.sqrt( ( center[0] - ccMesh[:,0] )**2 + ( center[1] - ccMesh[:,1] )**2 + ( center[2] - ccMesh[:,2] )**2 ) < radius - + # Return a tuple return ind # 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 = getIndicesSphere(loc[:,0],radi[0],mesh.gridCC) model[ind] = sig[1] - + # Second anomaly ind = 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_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() @@ -127,56 +131,56 @@ def run(loc=np.c_[[-50.,0.,-50.],[50.,0.,-50.]], sig=np.r_[1e-2,1e-1,1e-3], radi # 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 not re.match(stype,'pdp'): + 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: - + RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1,1] / mesh.vol[inds] ) + + else: + # 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] ) - - + RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-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]) - + # 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 ) + + 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: @@ -190,27 +194,27 @@ def run(loc=np.c_[[-50.,0.,-50.],[50.,0.,-50.]], sig=np.r_[1e-2,1e-1,1e-3], radi 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') - + + + 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) - + 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() \ No newline at end of file + run() diff --git a/SimPEG/Examples/EM_FDEM_SusEffects.py b/SimPEG/Examples/EM_FDEM_SusEffects.py index 3e50bf91..1abbb16f 100644 --- a/SimPEG/Examples/EM_FDEM_SusEffects.py +++ b/SimPEG/Examples/EM_FDEM_SusEffects.py @@ -6,7 +6,7 @@ from scipy.constants import mu_0 def run(plotIt=True): """ EM: FDEM: Effects of susceptibility - =============================== + =================================== When airborne freqeuncy domain EM (AFEM) survey is flown over the earth including significantly susceptible bodies (magnetite-rich rocks), diff --git a/docs/examples/DC_PseudoSection_Simulation.rst b/docs/examples/DC_PseudoSection_Simulation.rst index 8f1c3ce1..1d4330a1 100644 --- a/docs/examples/DC_PseudoSection_Simulation.rst +++ b/docs/examples/DC_PseudoSection_Simulation.rst @@ -9,13 +9,16 @@ .. --------------------------------- .. -orward Simulation -ard model conductive spheres in a half-space and plot a pseudo-section +DC Forward Simulation +===================== -ted on Mon Feb 01 19:28:06 2016 +Forward model conductive spheres in a half-space and plot a pseudo-section + +Created on Mon Feb 01 19:28:06 2016 + +@fourndo -rndo .. plot:: diff --git a/docs/examples/EM_FDEM_SusEffects.rst b/docs/examples/EM_FDEM_SusEffects.rst index da4b00a0..3e7dea80 100644 --- a/docs/examples/EM_FDEM_SusEffects.rst +++ b/docs/examples/EM_FDEM_SusEffects.rst @@ -10,7 +10,7 @@ EM: FDEM: Effects of susceptibility -=============================== +=================================== When airborne freqeuncy domain EM (AFEM) survey is flown over the earth including significantly susceptible bodies (magnetite-rich rocks), From 88ef74ac38ed0e29bf95f907912a201c212467b7 Mon Sep 17 00:00:00 2001 From: Thibaut Astic Date: Fri, 12 Feb 2016 15:40:48 -0800 Subject: [PATCH 08/18] MT_1D_analytic example --- .../Examples/MT_1D_analytic_nlayer_Earth.py | 443 ++++++++++++++++++ 1 file changed, 443 insertions(+) create mode 100644 SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py diff --git a/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py b/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py new file mode 100644 index 00000000..0cccc5f1 --- /dev/null +++ b/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py @@ -0,0 +1,443 @@ +from scipy.constants import epsilon_0, mu_0 +import matplotlib.pyplot as plt +import numpy as np +from ipywidgets import * +#from SimPEG.EM.Utils import k, omega + +""" +MT1D: n layered earth problem +***************************** + +Author: Thibaut Astic +Contact: thast@eos.ubc.ca +Date: January 2016 + +This code compute the analytic response of a n-layered Earth to a plane wave (Magneto-Tellurics). + +We start by looking at Maxwell's equations in the electric +field \\\(\\\mathbf{E}\\) and the magnetic flux +\\\(\\\mathbf{H}\\) to write the wave equations +\\(\\ \nabla ^2 \mathbf{E_x} + k^2 \mathbf{E_x} = 0 \\) & +\\(\\ \nabla ^2 \mathbf{H_y} + k^2 \mathbf{H_y} = 0 \\) + +Then solving the equations in each layer "j" between z_{j-1} and z_j in the form of +\\(\\ E_{x,j} (z) = U_j e^{i k (z-z_{j-1})} + D_j e^{-i k (z-z_{j-1})} \\) +\\(\\ H_{y,j} (z) = \frac{1}{Z_j} (D_j e^{-i k (z-z_{j-1})} - U_j e^{i k (z-z_{j-1})}) \\) + +With U and D the Up and Down components of the E-field. + +The iteration from one layer to another is ensure by: + +\\(\\ \left(\begin{matrix} E_{x,j} \\ H_{y,j} \end{matrix} \right) = + P_j T_j P^{-1}_J \left(\begin{matrix} E_{x,j+1} \\ H_{y,j+1} \end{matrix} \right) \\) + +And the Boundary Condition is set for the E-field in the last layer, with no Up component (=0) +and only a down component (=1 then normalized by the highest amplitude to ensure numeric stability) + +The layer 0 is assumed to be the air layer. + +""" + +#Frequency conversion +omega = lambda f: 2.*np.pi*f + +#Evaluate k wavenumber +k = lambda mu,sig,eps,f: np.sqrt(mu*mu_0*eps*epsilon_0*(2.*np.pi*f)**2.-1.j*mu*mu_0*sig*omega(f)) + +#Define a frquency range for a survey +frange = lambda minfreq, maxfreq, step: np.logspace(minfreq,maxfreq,num = step, base = 10.) + +#Functions to create random physical Perties for a n-layered earth +thick = lambda minthick, maxthick, nlayer: np.append(np.array([1.2*10.**5]), + np.ndarray.round(minthick + (maxthick-minthick)* np.random.rand(nlayer-1,1) + ,decimals =1)) + +sig = lambda minsig, maxsig, nlayer: np.append(np.array([0.]), + np.ndarray.round(10.**minsig + (10.**maxsig-10.**minsig)* np.random.rand(nlayer,1) + ,decimals=3)) + +mu = lambda minmu, maxmu, nlayer: np.append(np.array([1.]), + np.ndarray.round(minmu + (maxmu-minmu)* np.random.rand(nlayer,1) + ,decimals=1)) + +eps = lambda mineps, maxeps, nlayer: np.append(np.array([1.]), + np.ndarray.round(mineps + (maxeps-mineps)* np.random.rand(nlayer,1) + ,decimals=1)) + +#Evaluate Impedance Z of a layer +ImpZ = lambda f, mu, k: omega(f)*mu*mu_0/k + +#Complex Cole-Cole Conductivity - EM utils +PCC= lambda siginf,m,t,c,f: siginf*(1.-(m/(1.+(1j*omega(f)*t)**c))) + + +#Converted thickness array into top of layer array +def top(thick): + topv= np.zeros(len(thick)+1) + + topv[0]=-thick[0] + + for i in range(1,len(topv),1): + topv[i] = topv[i-1] + thick[i-1] + + return topv + +#Propagation Matrix and theirs inverses + +#matrix T for transition of Up and Down components accross a layer +T = lambda h,k: np.matrix([[np.exp(1j*k*h),0.],[0.,np.exp(-1j*k*h)]],dtype='complex_') + +Tinv = lambda h,k: np.matrix([[np.exp(-1j*k*h),0.],[0.,np.exp(1j*k*h)]],dtype='complex_') + +#transition of Up and Down components accross a layer +UD_Z = lambda UD,z,zj,k : T((z-zj),k)*UD + + +#matrix P relating Up and Down components with E and H fields +P = lambda z: np.matrix([[1.,1,],[-1./z,1./z]],dtype='complex_') + +Pinv = lambda z: np.matrix([[1.,-z],[1.,z]],dtype='complex_')/2. + + +#Time Variation of E and H +E_ZT = lambda U,D,f,t : np.exp(1j*omega(f)*t)*(U+D) +H_ZT = lambda U,D,Z,f,t : (1./Z)*np.exp(1j*omega(f)*t)*(D-U) + +#Plot the configuration of the problem +def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z): + + topn = top(thick) + widthn = np.arange(-widthg,widthg+widthg/10.,widthg/10.) + + ax.set_ylim([z.min(),z.max()]) + ax.set_xlim([-widthg,widthg]) + + ax.set_ylabel("Depth (m)", fontsize=16.) + ax.yaxis.tick_right() + ax.yaxis.set_label_position("right") + + #define filling for the different layers + hatches=['/' , '+', 'x', '|' , '\\', '-' , 'o' , 'O' , '.' , '*' ] + + #Write the physical properties of air + ax.annotate(("Air, $\sigma$ =%1.0f mS/m")%(sig[0]*10**(3)), + xy=(-widthg/2., -np.abs(z.max())/2.), xycoords='data', + xytext=(-widthg/2., -np.abs(z.max())/2.), textcoords='data', + fontsize=14.) + + ax.annotate(("$\epsilon_r$= %1i")%(eps[0]), + xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data', + xytext=(-widthg/2., -np.abs(z.max())/3.), textcoords='data', + fontsize=14.) + + ax.annotate(("$\mu_r$= %1i")%(mu[0]), + xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data', + xytext=(0, -np.abs(z.max())/3.), textcoords='data', + fontsize=14.) + + #Write the physical properties of the differents layers up to the (n-1)-th and fill it with pattern + for i in range(1,len(topn)-1,1): + if topn[i] == topn[i+1]: + pass + else: + ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[i]*10**(3)), + xy=(0., (2.*topn[i]+topn[i+1])/3), xycoords='data', + xytext=(0., (2.*topn[i]+topn[i+1])/3), textcoords='data', + fontsize=14.) + + ax.annotate(("$\epsilon_r$= %1i")%(eps[i]), + xy=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), xycoords='data', + xytext=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), textcoords='data', + fontsize=14.) + + ax.annotate(("$\mu_r$= %1.2f")%(mu[i]), + xy=(-widthg/2., (2.*topn[i]+topn[i+1])/3), xycoords='data', + xytext=(-widthg/2., (2.*topn[i]+topn[i+1])/3), textcoords='data', + fontsize=14.) + + ax.plot(widthn,topn[i]*np.ones_like(widthn),color='black') + ax.fill_between(widthn,topn[i],topn[i+1],alpha=0.3,color="none",edgecolor='black', hatch=hatches[(i-1)%10]) + + #Write the physical properties of the n-th layer and fill it with pattern + ax.plot(widthn,topn[-1]*np.ones_like(widthn),color='black') + ax.fill_between(widthn,topn[-1],z.max(),alpha=0.3,color="none",edgecolor='black', hatch=hatches[(len(topn)-2)%10]) + + ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[-1]*10**(3)), + xy=(0., (2.*topn[-1]+z.max())/3), xycoords='data', + xytext=(0., (2.*topn[-1]+z.max())/3), textcoords='data', + fontsize=14.) + + ax.annotate(("$\epsilon_r$= %1i")%(eps[-1]), + xy=(-widthg/1.1, (2.*topn[-1]+z.max())/3), xycoords='data', + xytext=(-widthg/1.1, (2.*topn[-1]+z.max())/3), textcoords='data', + fontsize=14.) + + ax.annotate(("$\mu_r$= %1.2f")%(mu[-1]), + xy=(-widthg/2., (2.*topn[-1]+z.max())/3), xycoords='data', + xytext=(-widthg/2., (2.*topn[-1]+z.max())/3), textcoords='data', + fontsize=14.) + + #plot Trees! + ax.annotate("", + xy=(widthg/2., -1.*z.max()/5.), xycoords='data', + xytext=(widthg/2., 0.), textcoords='data', + arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.) + ) + + ax.annotate("", + xy=(widthg/2., -3./4.*z.max()/5.), xycoords='data', + xytext=(widthg/2., 0.), textcoords='data', + arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.) + ) + + ax.annotate("", + xy=(widthg/2., -1./2.*z.max()/5.), xycoords='data', + xytext=(widthg/2., 0.), textcoords='data', + arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.) + ) + + ax.annotate("", + xy=(1.2*widthg/2., -1.*z.max()/5.), xycoords='data', + xytext=(1.2*widthg/2., 0.), textcoords='data', + arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.) + ) + + ax.annotate("", + xy=(1.2*widthg/2., -3./4.*z.max()/5.), xycoords='data', + xytext=(1.2*widthg/2., 0.), textcoords='data', + arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.) + ) + + ax.annotate("", + xy=(1.2*widthg/2., -1./2.*z.max()/5.), xycoords='data', + xytext=(1.2*widthg/2., 0.), textcoords='data', + arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.) + ) + + ax.annotate("", + xy=(1.5*widthg/2., -1.*z.max()/5.), xycoords='data', + xytext=(1.5*widthg/2., 0.), textcoords='data', + arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.) + ) + + ax.annotate("", + xy=(1.5*widthg/2., -3./4.*z.max()/5.), xycoords='data', + xytext=(1.5*widthg/2., 0.), textcoords='data', + arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.) + ) + + ax.annotate("", + xy=(1.5*widthg/2., -1./2.*z.max()/5.), xycoords='data', + xytext=(1.5*widthg/2., 0.), textcoords='data', + arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.) + ) + + + ax.invert_yaxis() + + return ax + +#Propagate Up and Down component for a certain frequency & evaluate E and H field + +def Propagate(f,H,sig,chg,taux,c,mu,eps,n): + + sigcm = np.zeros_like(sig,dtype='complex_') + + for j in range(1,len(sig)): + sigcm[j]=PCC(sig[j],chg[j],taux[j],c[j],f) + + K = k(mu,sigcm,eps,f) + Z = ImpZ(f,mu,K) + + EH = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_') + UD = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_') + + UD[1,-1] = 1. + + for i in range(-2,-(n+2),-1): + + UD[:,i] = Tinv(H[i+1],K[i])*Pinv(Z[i])*P(Z[i+1])*UD[:,i+1] + UD = UD/((np.abs(UD[0,:]+UD[1,:])).max()) + + for j in range(0,n+1): + EH[:,j] = np.matrix([[1.,1,],[-1./Z[j],1./Z[j]]])*UD[:,j] + + return UD, EH, Z ,K + + +#Evaluate the apparent resistivity and phase for a frequency range +def appres(F,H,sig,chg,taux,c,mu,eps,n): + + Res = np.zeros_like(F) + Phase = np.zeros_like(F) + App_ImpZ= np.zeros_like(F,dtype='complex_') + + for i in range(0,len(F)): + + UD,EH,Z ,K = Propagate(F[i],H,sig,chg,taux,c,mu,eps,n) + + App_ImpZ[i] = EH[0,1]/EH[1,1] + + Res[i] = np.abs(App_ImpZ[i])**2./(mu_0*omega(F[i])) + Phase[i] = np.angle(App_ImpZ[i], deg = True) + + return Res,Phase + +#Evaluate Up, Down components, E and H field, for a frequency range, +#a discretized depth range and a time range (use to calculate envelope) +def calculateEHzt(F,H,sig,chg,taux,c,mu,eps,n,zsample,tsample): + + topc = top(H) + + layer = np.zeros(len(zsample),dtype=np.int)-1 + + Exzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_') + Hyzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_') + Uz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_') + Dz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_') + UDaux = np.matrix(np.zeros((2,len(zsample)),dtype = 'complex_'),dtype = 'complex_') + + for i in range(0,n+1,1): + layer = layer+(zsample>=topc[i])*1 + + for j in range(0,len(F)): + + UD,EH,Z ,K = Propagate(F[j],H,sig,chg,taux,c,mu,eps,n) + + for p in range(0,len(zsample)): + + UDaux[:,p] = UD_Z(UD[:,layer[p]],zsample[p],topc[layer[p]],K[layer[p]]) + + for q in range(0,len(tsample)): + + Exzt[p,q] = Exzt[p,q] + E_ZT(UDaux[0,p],UDaux[1,p],F[j],tsample[q])/len(F) + Hyzt[p,q] = Hyzt[p,q] + H_ZT(UDaux[0,p],UDaux[1,p],Z[layer[p]],F[j],tsample[q])/len(F) + Uz[p,q] = Uz[p,q] + UDaux[0,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F) + Dz[p,q] = Dz[p,q] + UDaux[1,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F) + + return Exzt,Hyzt,Uz,Dz,UDaux,layer + + +#Function to Plot Apparent Resistivity and Phase +def PlotAppRes(F,H,sig,chg,taux,c,mu,eps,n,fenvelope,PlotEnvelope): + + Res, Phase = appres(F,H,sig,chg,taux,c,mu,eps,n) + + fig,ax = plt.subplots(1,2,figsize=(16,10)) + + ax[0].scatter(Res,F,color='black') + ax[0].set_xscale('Log') + ax[0].set_yscale('Log') + ax[0].set_xlim([10.**(np.log10(Res.min())-1.),10.**(np.log10(Res.max())+1.)]) + ax[0].set_ylim([F.min(),F.max()]) + ax[0].set_xlabel('Apparent Resistivity (Ohm*m)',fontsize=16.,color="black") + ax[0].set_ylabel('Frequency (Hz)',fontsize=16.) + ax[0].grid(which='major') + + ax0 = ax[0].twiny() + + ax0.set_xlim([0.,90.]) + ax0.set_ylim([F.min(),F.max()]) + ax0.scatter(Phase,F,color='purple') + ax0.set_xlabel('Phase (Degrees)',fontsize=16.,color="purple") + + zc=np.arange(-(H[1:].max()+10)*n,(H[1:].max()+10)*n,10.) + + ax[0].tick_params(labelsize=16) + ax[1].tick_params(labelsize=16) + ax0.tick_params(labelsize=16) + + if PlotEnvelope: + + widthn=np.logspace(np.log10(Res.min())-1., np.log10(Res.max())+1., num=100, endpoint=True, base=10.0) + fenvelope1n=np.ones(100)*fenvelope + ax[0].plot(widthn,fenvelope1n,linestyle='dashed',color='black') + + tc=np.arange(0.,1./fenvelope,0.01/(fenvelope)) + Exzt,Hyzt,Uz,Dz,UDaux,layer = calculateEHzt(np.array([fenvelope]),H,sig,chg,taux,c,mu,eps,n,zc,tc) + + ax1=ax[1].twiny() + + ax[1].tick_params(labelsize=16) + ax1.tick_params(labelsize=16) + + ax[1].set_xlabel('Amplitude Electric Field E (V/m)',color='blue',fontsize=16) + + ax1.set_xlabel('Amplitude Magnetic Field H (A/m)',color='red',fontsize=16) + + ax[1].fill_betweenx(zc,np.squeeze(np.asarray(np.real(Exzt.min(axis=1)))), + np.squeeze(np.asarray(np.real(Exzt.max(axis=1)))), + color='blue', alpha=0.1) + + ax1.fill_betweenx(zc,np.squeeze(np.asarray(np.real(Hyzt.min(axis=1)))), + np.squeeze(np.asarray(np.real(Hyzt.max(axis=1)))), + color='red', alpha=0.1) + + ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],(1.5*np.abs(Exzt).max()),zc) + ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()]) + ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()]) + else: + print 'No envelop (if True, might be slow)' + ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],1.,zc) + ax[1].get_xaxis().set_ticks([]) + + plt.show() + +#Interactive MT for Notebook +def PlotAppRes3LayersInteract(h1,h2,sigl1,sigl2,sigl3,mul1,mul2,mul3,epsl1,epsl2,epsl3,PlotEnvelope,F_Envelope): + + frangn=frange(-5,5,100.) + sig3= np.array([0.,0.001,0.1, 0.001]) + thick3 = np.array([120000.,50.,50.]) + eps3=np.array([1.,1.,1.,1]) + mu3=np.array([1.,1.,1.,1]) + chg3=np.array([0.,0.1,0.,0.2]) + chg3_0=np.array([0.,0.1,0.,0.]) + taux3=np.array([0.,0.1,0.,0.1]) + c3=np.array([1.,1.,1.,1.]) + + sig3[1]=sigl1 + sig3[1]=10.**sig3[1] + sig3[2]=sigl2 + sig3[2]=10.**sig3[2] + sig3[3]=sigl3 + sig3[3]=10.**sig3[3] + mu3[1]=mul1 + mu3[2]=mul2 + mu3[3]=mul3 + eps3[1]=epsl1 + eps3[2]=epsl2 + eps3[3]=epsl3 + thick3[1]=h1 + thick3[2]=h2 + + PlotAppRes(frangn,thick3,sig3,chg3_0,taux3,c3,mu3,eps3,3,F_Envelope,PlotEnvelope) + + +def run(n,plotIt=True): + # something to make a plot + + F = frange(-5.,5.,20) + H = thick(50.,100.,n) + sign = sig(-5.,0.,n) + mun = mu(1.,2.,n) + epsn = eps(1.,9.,n) + chg = np.zeros_like(sign) + taux = np.zeros_like(sign) + c = np.zeros_like(sign) + + Res, Phase = appres(F,H,sign,chg,taux,c,mun,epsn,n) + + if plotIt: + + PlotAppRes(F, H, sign, chg, taux, c, mun, epsn, n, fenvelope=1000., PlotEnvelope=True) + + return Res, Phase + +if __name__ == '__main__': + run(3) + + + + + From 4cf5d49524d7d44e390ddd2dc8a12113dcc0b002 Mon Sep 17 00:00:00 2001 From: Thibaut Astic Date: Fri, 12 Feb 2016 15:48:48 -0800 Subject: [PATCH 09/18] Ignoring non functioning examples --- SimPEG/Examples/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index df65e16b..73665644 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -1,10 +1,10 @@ # Run this file to add imports. ##### AUTOIMPORTS ##### -import DC_PseudoSection_Simulation +#import DC_PseudoSection_Simulation import EM_FDEM_1D_Inversion import EM_FDEM_Analytic_MagDipoleWholespace -import EM_FDEM_SusEffects +#import EM_FDEM_SusEffects import EM_TDEM_1D_Inversion import FLOW_Richards_1D_Celia1990 import Forward_BasicDirectCurrent From 6286e48830ff06a82153f973a39faae59f1457e8 Mon Sep 17 00:00:00 2001 From: Thibaut Astic Date: Mon, 15 Feb 2016 13:13:44 -0800 Subject: [PATCH 10/18] Sphere Electrostatic example. code cleaned, commented and updated. --- .../Examples/sphereElectrostatic_example.py | 776 ++++++++++++++++++ 1 file changed, 776 insertions(+) create mode 100644 SimPEG/Examples/sphereElectrostatic_example.py diff --git a/SimPEG/Examples/sphereElectrostatic_example.py b/SimPEG/Examples/sphereElectrostatic_example.py new file mode 100644 index 00000000..4b96f740 --- /dev/null +++ b/SimPEG/Examples/sphereElectrostatic_example.py @@ -0,0 +1,776 @@ +from scipy.constants import epsilon_0 +import matplotlib.pyplot as plt +import matplotlib.colors as colors +import numpy as np +from ipywidgets import * +from SimPEG.Utils import ndgrid, mkvc + +''' +Authors: Thibaut Astic, Lindsey Heagy, Sanna Tyrvainen, Ronghua Peng +Date: December 2015 + +This code defines function to resolve analytically the electrostatic sphere problem. +We first define a problem configuration, with a conductive or resistive sphere in a +wholespace background. +We then calculate the potential, then the electric field, then the current density and +finally the charges accumulation. + +Several plotting functions are defined for data visualisation. + + +''' + +# Plot options +ftsize_title = 18 #font size for titles +ftsize_axis = 14 #font size for axis ticks +ftsize_label = 14 #font size for axis labels + +# Radius function, useful sigma ratio, and log scale converter +r = lambda x,y,z: np.sqrt(x**2.+y**2.+z**2.) +sigf = lambda sig0,sig1: (sig1-sig0)/(sig1+2.*sig0) + +#tools to convert log conductivity in conductivity +def conductivity_log_wrapper(log_sig0,log_sig1): + sig0 = 10.**log_sig0 + sig1 = 10.**log_sig1 + + return sig0,sig1 + +# Examples +#Plot the configuration. Label=False is used to generate a general case figure +def get_Setup(XYZ,sig0,sig1,R,E0,ax,label,colorsphere): + ''' + XYZ: ndgrid + sig0: conductivity of the background + sig1: conductivity of the sphere + R: radius of the sphere + E0: Amplitude of the uniform electrostatic field + ax: ax where to plot the configuration + label: True: plot real values, False: plot general case + colorsphere: color of the sphere, format [x,x,x] + ''' + + xplt = np.linspace(-R, R, num=100) + xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2]) + dx = xr[1]-xr[0] + top = np.sqrt(R**2-xplt**2) + bot = -np.sqrt(R**2-xplt**2) + + if R != 0: + ax.plot(xplt, top, xplt, bot, color=colorsphere,linewidth=1.5) + ax.fill_between(xplt,bot,top,color=colorsphere,alpha=0.5 ) + ax.arrow(0.,0.,np.sqrt(2.)*R/2.,np.sqrt(2.)*R/2.,head_width=0.,head_length=0.) + + if label: + ax.annotate(("$\sigma_1$=%3.3f mS/m")%(sig1*10.**(3.)), + xy=(0.,-R/2.), xycoords='data', + xytext=(0.,-R/2.), textcoords='data', + fontsize=14.) + ax.annotate(("$\sigma_0$= %3.3f mS/m")%(sig0*10.**(3.)), + xy=(0.,-1.5*R), xycoords='data', + xytext=(0.,-1.5*R), textcoords='data', + fontsize=14.) + ax.annotate(('$\mathbf{E_0} = %1i \mathbf{\hat{x}}$ V/m')%(E0), + xy=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), xycoords='data', + xytext=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), textcoords='data', + fontsize=14.) + ax.annotate(('$R$ = %1i m')%(R), + xy=(R/4.+(xr[1]-xr[0]),R/4.), xycoords='data', + xytext=(R/4.+(xr[1]-xr[0]),R/4.), textcoords='data', + fontsize=14.) + ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label) + ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label) + ax.tick_params(labelsize=ftsize_axis) + + else: + ax.set_xticklabels([]) + ax.set_yticklabels([]) + ax.text(-1.,-np.sqrt(R)/2.-10.,'$\sigma_1$',fontsize=14) + ax.text(-0.05,-R-10,'$\sigma_0$',fontsize=14) + ax.annotate(('$\mathbf{E_0} = E_0 \mathbf{\hat{x}}$ V/m'), + xy=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), xycoords='data', + xytext=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), textcoords='data', + fontsize=14.) + ax.annotate(('$R$'), + xy=(R/4.+(xr[1]-xr[0]),R/4.), xycoords='data', + xytext=(R/4.+(xr[1]-xr[0]),R/4.), textcoords='data', + fontsize=14.) + ax.set_xlabel('x',fontsize=12) + ax.set_ylabel('y',fontsize=12) + + else: + if label: + ax.annotate(("$\sigma_0$= %3.3f mS/m")%(sig0*10.**(3.)), + xy=(0.,-1.5*R), xycoords='data', + xytext=(0.,-1.5*R), textcoords='data', + fontsize=14.) + ax.annotate(('$\mathbf{E_0} = %1i \mathbf{\hat{x}}$ V/m')%(E0), + xy=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), xycoords='data', + xytext=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), textcoords='data', + fontsize=14.) + ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label) + ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label) + ax.tick_params(labelsize=ftsize_axis) + + else: + ax.set_xticklabels([]) + ax.set_yticklabels([]) + ax.text(-0.05,-10,'$\sigma_0$',fontsize=14) + ax.text(xr.min()+np.abs(xr.max()-xr.min())/20., 0, '$\mathbf{E_0} = E_0 \mathbf{\hat{x}}$ V/m', fontsize=14) + ax.set_xlabel('x',fontsize=12) + ax.set_ylabel('y',fontsize=12) + + + ax.set_xlim([xr.min(),xr.max()]) + ax.set_ylim([yr.min(),yr.max()]) + [ax.arrow(xr.min(),_,np.abs(xr.max()-xr.min())/20.,0.,head_width=5.,head_length=2.,color='k') for _ in np.linspace(yr.min(),yr.max(),num=10)] + ax.patch.set_facecolor([0.4,0.7,0.4]) + ax.patch.set_alpha(0.2) + + ax.set_aspect('equal') + + + + return ax + +def get_Conductivity(XYZ,sig0,sig1,R): + ''' + Define the conductivity for each point of the space + ''' + x,y,z = XYZ[:,0],XYZ[:,1],XYZ[:,2] + r_view=r(x,y,z) + + ind0= (r_view>R) + ind1= (r_view<=R) + + assert (ind0 + ind1).all(), 'Some indicies not included' + + Sigma = np.zeros_like(x) + + Sigma[ind0] = sig0 + Sigma[ind1] = sig1 + + return Sigma + + +def get_Potential(XYZ,sig0,sig1,R,E0): + + ''' + Function that returns the total, the primary and the secondary potentials, assumes an x-oriented inducing field and that the sphere is at the origin + :input: grid, outer sigma, inner sigma, radius of the sphere, strength of the electric field + ''' + + x,y,z = XYZ[:,0],XYZ[:,1],XYZ[:,2] + + sig_cur = sigf(sig0,sig1) + + r_cur = r(x,y,z) # current radius + + ind0 = (r_cur > R) + ind1 = (r_cur <= R) + + assert (ind0 + ind1).all(), 'Some indicies not included' + + Vt = np.zeros_like(x) + Vp = np.zeros_like(x) + Vs = np.zeros_like(x) + + Vt[ind0] = -E0*x[ind0]*(1.-sig_cur*R**3./r_cur[ind0]**3.) # total potential outside the sphere + Vt[ind1] = -E0*x[ind1]*3.*sig0/(sig1+2.*sig0) # inside the sphere + + + Vp = - E0*x # primary potential + + Vs = Vt - Vp # secondary potential + + return Vt,Vp,Vs + +#plot the primary potential on ax +def Plot_Primary_Potential(XYZ,sig0,sig1,R,E0,ax): + + Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0) + + xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2]) + + xcirc = xr[np.abs(xr) <= R] + + Pplot = ax.pcolor(xr,yr,Vp.reshape(xr.size,yr.size)) + ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k') + ax.set_title('Primary Potential',fontsize=ftsize_title) + cb = plt.colorbar(Pplot,ax=ax) + cb.set_label(label= 'Potential ($V$)',size=ftsize_label) + cb.ax.tick_params(labelsize=ftsize_axis) + ax.set_xlim([xr.min(),xr.max()]) + ax.set_ylim([yr.min(),yr.max()]) + ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label) + ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label) + ax.set_aspect('equal') + ax.tick_params(labelsize=ftsize_axis) + + return ax + +#plot the total potential on ax +def Plot_Total_Potential(XYZ,sig0,sig1,R,E0,ax): + + Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0) + + xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2]) + + xcirc = xr[np.abs(xr) <= R] + + + Pplot = ax.pcolor(xr,yr,Vt.reshape(xr.size,yr.size)) + ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k') + ax.set_title('Total Potential',fontsize=ftsize_title) + cb = plt.colorbar(Pplot,ax=ax) + cb.set_label(label= 'Potential ($V$)',size=ftsize_label) + cb.ax.tick_params(labelsize=ftsize_axis) + ax.set_xlim([xr.min(),xr.max()]) + ax.set_ylim([yr.min(),yr.max()]) + ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label) + ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label) + ax.set_aspect('equal') + ax.tick_params(labelsize=ftsize_axis) + + return ax + +#plot the secondary potential on ax +def Plot_Secondary_Potential(XYZ,sig0,sig1,R,E0,ax): + + Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0) + + xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2]) + + xcirc = xr[np.abs(xr) <= R] + + Pplot = ax.pcolor(xr,yr,Vs.reshape(xr.size,yr.size)) + ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k') + ax.set_title('Secondary Potential',fontsize=ftsize_title) + cb = plt.colorbar(Pplot,ax=ax) + cb.set_label(label= 'Potential ($V$)',size=ftsize_label) + cb.ax.tick_params(labelsize=ftsize_axis) + ax.set_xlim([xr.min(),xr.max()]) + ax.set_ylim([yr.min(),yr.max()]) + ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label) + ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label) + ax.set_aspect('equal') + ax.tick_params(labelsize=ftsize_axis) + + return ax + + +def get_ElectricField(XYZ,sig0,sig1,R,E0): + ''' + Function that returns the total, the primary and the secondary electric fields, + input: grid, outer sigma, inner sigma, radius of the sphere, strength of the electric field + ''' + + x,y,z= XYZ[:,0], XYZ[:,1], XYZ[:,2] + + r_cur=r(x,y,z) # current radius + + ind0= (r_cur>R) + ind1= (r_cur<=R) + + assert (ind0 + ind1).all(), 'Some indicies not included' + + Ep = np.zeros(shape=(len(x),3)) + Ep[:,0] = E0 + + Et = np.zeros(shape=(len(x),3)) + + Et[ind0,0] = E0 + E0*R**3./(r_cur[ind0]**5.)*sigf(sig0,sig1)*(2.*x[ind0]**2.-y[ind0]**2.-z[ind0]**2.); + Et[ind0,1] = E0*R**3./(r_cur[ind0]**5.)*3.*x[ind0]*y[ind0]*sigf(sig0,sig1); + Et[ind0,2] = E0*R**3./(r_cur[ind0]**5.)*3.*x[ind0]*z[ind0]*sigf(sig0,sig1); + + Et[ind1,0] = 3.*sig0/(sig1+2.*sig0)*E0; + Et[ind1,1] = 0.; + Et[ind1,2] = 0.; + + Es = Et - Ep + + return Et, Ep, Es + +#plot the total electric field on ax +def Plot_Total_ElectricField(XYZ,sig0,sig1,R,E0,ax): + + Et, Ep, Es = get_ElectricField(XYZ,sig0,sig1,R,E0) + + xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2]) + + xcirc = xr[np.abs(xr) <= R] + + EtXr = Et[:,0].reshape(xr.size, yr.size) + EtYr = Et[:,1].reshape(xr.size, yr.size) + EtAmp = np.sqrt(Et[:,0]**2+Et[:,1]**2 + Et[:,2]**2).reshape(xr.size, yr.size) + + ax.set_xlim([xr.min(),xr.max()]) + ax.set_ylim([yr.min(),yr.max()]) + ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label) + ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label) + ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k') + ax.tick_params(labelsize=ftsize_axis) + ax.set_aspect('equal') + + Eplot = ax.pcolor(xr,yr,EtAmp) + cb = plt.colorbar(Eplot,ax=ax) + cb.set_label(label= 'Amplitude ($V/m$)',size=ftsize_label) #weight='bold') + cb.ax.tick_params(labelsize=ftsize_axis) + ax.streamplot(xr,yr,EtXr,EtYr,color='gray',linewidth=2.,density=0.75)#angles='xy',scale_units='xy',scale=0.05) + ax.set_title('Total Field',fontsize=ftsize_title) + + + return ax + +#plot the secondary electric field on ax +def Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax): + + Et, Ep, Es = get_ElectricField(XYZ,sig0,sig1,R,E0) + + xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2]) + + xcirc = xr[np.abs(xr) <= R] + + EsXr = Es[:,0].reshape(xr.size, yr.size) + EsYr = Es[:,1].reshape(xr.size, yr.size) + EsAmp = np.sqrt(Es[:,0]**2+Es[:,1]**2+Es[:,2]**2).reshape(xr.size, yr.size) + + ax.set_xlim([xr.min(),xr.max()]) + ax.set_ylim([yr.min(),yr.max()]) + ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label) + ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label) + ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k') + ax.tick_params(labelsize=ftsize_axis) + ax.set_aspect('equal') + + Eplot = ax.pcolor(xr,yr,EsAmp) + cb = plt.colorbar(Eplot,ax=ax) + cb.set_label(label= 'Amplitude ($V/m$)',size=ftsize_label) #weight='bold') + cb.ax.tick_params(labelsize=ftsize_axis) + ax.streamplot(xr,yr,EsXr,EsYr,color='gray',linewidth=2.,density=0.75)#,angles='xy',scale_units='xy',scale=0.05) + ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k') + ax.set_title('Secondary Field',fontsize=ftsize_title) + + return ax + + +def get_Current(XYZ,sig0,sig1,R,Et,Ep,Es): + ''' + Function that returns the total, the primary and the secondary current densities, + :input: grid, outer sigma, inner sigma, radius of the sphere, total, the primary and the seconadry electric fields, + ''' + + x,y,z= XYZ[:,0], XYZ[:,1], XYZ[:,2] + + r_cur=r(x,y,z) + + ind0= (r_cur>R) + ind1= (r_cur<=R) + + assert (ind0 + ind1).all(), 'Some indicies not included' + + Jt = np.zeros(shape=(len(x),3)) + J0 = np.zeros(shape=(len(x),3)) + Js = np.zeros(shape=(len(x),3)) + + + Jp = sig0*Ep + + Jt[ind0,:] = sig0*Et[ind0,:] + Jt[ind1,:] = sig1*Et[ind1,:] + + Js[ind0,:] = sig0*(Et[ind0,:]-Ep[ind0,:]) + Js[ind1,:] = sig1*Et[ind1,:]-sig0*Ep[ind1,:] + + return Jt,Jp,Js + +#plot the total currents density on ax +def Plot_Total_Currents(XYZ,sig0,sig1,R,E0,ax): + + Et,Ep,Es = get_ElectricField(XYZ,sig0,sig1,R,E0) + Jt,Jp,Js = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es) + + xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2]) + xcirc = xr[np.abs(xr) <= R] + + JtXr = Jt[:,0].reshape(xr.size, yr.size) + JtYr = Jt[:,1].reshape(xr.size, yr.size) + JtAmp = np.sqrt(Jt[:,0]**2+Jt[:,1]**2+Jt[:,2]**2).reshape(xr.size, yr.size) + + ax.set_xlim([xr.min(),xr.max()]) + ax.set_ylim([yr.min(),yr.max()]) + ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k') + ax.set_ylabel('Y coordinate ($m$)',fontsize=ftsize_label) + ax.set_xlabel('X coordinate ($m$)',fontsize=ftsize_label) + ax.tick_params(labelsize=ftsize_axis) + ax.set_aspect('equal') + + Jplot = ax.pcolor(xr,yr,JtAmp.reshape(xr.size,yr.size)) + cb = plt.colorbar(Jplot,ax=ax) + cb.set_label(label= 'Current Density ($A/m^2$)',size=ftsize_label) #weight='bold') + cb.ax.tick_params(labelsize=ftsize_axis) + ax.streamplot(xr,yr,JtXr,JtYr,color='gray',linewidth=2.,density=0.75)#,angles='xy',scale_units='xy',scale=1) + ax.set_title('Total Current Density',fontsize=ftsize_title) + + return ax + + +#plot the secondary currents density on ax +def Plot_Secondary_Currents(XYZ,sig0,sig1,R,E0,ax): + + Et,Ep,Es = get_ElectricField(XYZ,sig0,sig1,R,E0) + Jt,Jp,Js = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es) + + xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2]) + xcirc = xr[np.abs(xr) <= R] + + JsXr = Js[:,0].reshape(xr.size, yr.size) + JsYr = Js[:,1].reshape(xr.size, yr.size) + JsAmp = np.sqrt(Js[:,1]**2+Js[:,0]**2+Jt[:,2]**2).reshape(xr.size,yr.size) + + ax.set_xlim([xr.min(),xr.max()]) + ax.set_ylim([yr.min(),yr.max()]) + ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k') + ax.set_ylabel('Y coordinate ($m$)',fontsize=ftsize_label) + ax.set_xlabel('X coordinate ($m$)',fontsize=ftsize_label) + ax.tick_params(labelsize=ftsize_axis) + ax.set_aspect('equal') + + Jplot = ax.pcolor(xr,yr,JsAmp.reshape(xr.size,yr.size)) + cb = plt.colorbar(Jplot,ax=ax) + cb.set_label(label= 'Current Density ($A/m^2$)',size=ftsize_label) #weight='bold') + cb.ax.tick_params(labelsize=ftsize_axis) + ax.streamplot(xr,yr,JsXr,JsYr,color='gray',linewidth=2.,density=0.75)#,angles='xy',scale_units='xy',scale=1) + ax.set_title('Secondary Current Density',fontsize=ftsize_title) + + return ax + + +def get_ChargesDensity(XYZ,sig0,sig1,R,Et,Ep): + ''' + Function that returns the charges accumulation at the background/sphere interface, + :input: grid, outer sigma, inner sigma, radius of the sphere, total and the primary electric fields, + ''' + + x,y,z= XYZ[:,0], XYZ[:,1], XYZ[:,2] + + dx = x[1]-x[0] + + r_cur=r(x,y,z) + + ind0 = (r_cur > R) + ind1 = (r_cur < R) + ind2 = ((r_cur < (R+dx/2)) & (r_cur > (R-dx/2)) ) + + assert (ind0 + ind1 + ind2).all(), 'Some indicies not included' + + rho = np.zeros_like(x) + + rho[ind0] = 0 + rho[ind1] = 0 + rho[ind2] = epsilon_0*3.*Ep[ind2,0]*sigf(sig0,sig1)*x[ind2]/(np.sqrt(x[ind2]**2.+y[ind2]**2.)) + + return rho + +#Plot charges density on ax +def Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax): + + xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2]) + xcirc = xr[np.abs(xr) <= R] + + Et, Ep, Es = get_ElectricField(XYZ,sig0,sig1,R,E0) + rho = get_ChargesDensity(XYZ,sig0,sig1,R,Et,Ep) + + ax.set_xlim([xr.min(),xr.max()]) + ax.set_ylim([yr.min(),yr.max()]) + ax.set_aspect('equal') + Cplot = ax.pcolor(xr,yr,rho.reshape(xr.size, yr.size)) + cb1 = plt.colorbar(Cplot,ax=ax) + cb1.set_label(label= 'Charge Density ($C/m^2$)',size=ftsize_label) #weight='bold') + cb1.ax.tick_params(labelsize=ftsize_axis) + ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k') + ax.set_ylabel('Y coordinate ($m$)',fontsize=ftsize_label) + ax.set_xlabel('X coordinate ($m$)',fontsize=ftsize_label) + ax.tick_params(labelsize=ftsize_axis) + ax.set_title('Charges Density', fontsize=ftsize_title) + + return ax + +def MN_Potential_total(sig0,sig1,R,E0,start,end,nbmp,mn): + + ''' + Function that return array of midpoints electrodes, electrodes positions, + potentials differences for total and secondary potentials fields, unormalized and + normalized to electrodes distances. + sig0: background conductivity + sig1: sphere conductivity + R: Sphere's radius + E0: uniform E field value + start: start point for the profile start.shape = (2,) + end: end point for the profile end.shape = (2,) + nbmp: number of dipoles + mn: Space between the M and N electrodes + ''' + + #D: total distance from start to end + D = np.sqrt((start[0]-end[0])**2.+(start[1]-end[1])**2.) + + #MP: dipoles'midpoint positions (x,y) + MP = np.zeros(shape=(nbmp,2)) + MP[:,0] = np.linspace(start[0],end[0],nbmp) + MP[:,1] = np.linspace(start[1],end[1],nbmp) + + #Dipoles'Electrodes positions around each midpoints + EL = np.zeros(shape=(2*nbmp,2)) + for n in range(0,len(EL),2): + EL[n,0] = MP[n/2,0] - ((end[0]-start[0])/D)*mn/2. + EL[n+1,0] = MP[n/2,0] + ((end[0]-start[0])/D)*mn/2. + EL[n,1] = MP[n/2,1] - ((end[1]-start[1])/D)*mn/2. + EL[n+1,1] = MP[n/2,1] + ((end[1]-start[1])/D)*mn/2. + + VtEL = np.zeros(2*nbmp) #Total Potential (Vt-) at each electrode (-EL) + VsEL = np.zeros(2*nbmp) #Secondary Potential (Vt-) at each electrode (-EL) + dVtMP = np.zeros(nbmp) #Diffence (d-) of Total Potential (Vt-) at each dipole (-MP) + dVtMPn = np.zeros(nbmp) #Diffence (d-) of Total Potential (Vt-) at each dipole (-MP) normalized for the mn spacing (n) + dVsMP = np.zeros(nbmp) #Diffence (d-) of Secondaty Potential (Vt-) at each dipole (-MP) + dVsMPn = np.zeros(nbmp) #Diffence (d-) of Secondary Potential (Vt-) at each dipole (-MP) normalized for the mn spacing (n) + dVpMP = np.zeros(nbmp) #Diffence (d-) of Primary Potential (Vt-) at each dipole (-MP) + dVpMPn = np.zeros(nbmp) #Diffence (d-) of Primary Potential (Vt-) at each dipole (-MP) normalized for the mn spacing (n) + + #Computing VtEL + for m in range(0,2*nbmp): + if (r(EL[m,0],EL[m,1],0) > R): + VtEL[m] = -E0*EL[m,0]*(1.-sigf(sig0,sig1)*R**3./r(EL[m,0],EL[m,1],0)**3.) + else: + VtEL[m] = -E0*EL[m,0]*3.*sig0/(sig1+2.*sig0) + + #Computing VsEL + VsEL = VtEL + E0*EL[:,0] + + #Computing dVtMP, dVsMP + for p in range(0,nbmp): + dVtMP[p] = VtEL[2*p]-VtEL[2*p+1] + dVtMPn[p] = dVtMP[p]/mn + dVsMP[p] = VsEL[2*p]-VsEL[2*p+1] + dVsMPn[p] = dVsMP[p]/mn + + return MP,EL,dVtMP,dVtMPn,dVsMP,dVsMPn + +#Compare the DC response of two configurations +def two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,nb_dipole,electrode_spacing,PlotOpt):#,linearcolor): + + #Define the mesh + xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2]) + + #Defining the Profile + start = np.array([xstart,ystart]) + end = np.array([xend,yend]) + + #Calculating the data from the defined survey line for Configuration 0 and 1 + MP0,EL0,VtdMP0,VtdMPn0,VsdMP0,VsdMPn0 = MN_Potential_total(sig0,sig1,R0,E0,start,end,nb_dipole,electrode_spacing) + MP1,EL1,VtdMP1,VtdMPn1,VsdMP1,VsdMPn1 = MN_Potential_total(sig0,sig2,R1,E0,start,end,nb_dipole,electrode_spacing) + + + # Initializing the figure + fig = plt.figure(figsize=(20,20)) + ax0 = plt.subplot2grid((20,12), (0, 0),colspan=6,rowspan=6) + ax1 = plt.subplot2grid((20,12), (0, 6),colspan=6,rowspan=6) + ax2 = plt.subplot2grid((20,12), (16, 2), colspan=9,rowspan=4) + ax3 = plt.subplot2grid((20,12), (8, 0),colspan=6,rowspan=6) + ax4 = plt.subplot2grid((20,12), (8, 6),colspan=6,rowspan=6) + + #Plotting the Configuration 0 + ax0 = get_Setup(XYZ,sig0,sig1,R0,E0,ax0,True,[0.6,0.1,0.1]) + + #Plotting the Configuration 1 + ax1 = get_Setup(XYZ,sig0,sig2,R1,E0,ax1,True,[0.1,0.1,0.6]) + + #Plotting the Data (Legends) + ax2.set_title('Potential Differences',fontsize=ftsize_title) + ax2.set_ylabel('Potential difference ($V$)',fontsize=ftsize_label) + ax2.set_xlabel('Distance from start point ($m$)',fontsize=ftsize_label) + ax2.tick_params(labelsize=ftsize_axis) + ax2.grid() + + if PlotOpt == 'Total': + ax3= Plot_Total_Potential(XYZ,sig0,sig1,R0,E0,ax3) + ax4= Plot_Total_Potential(XYZ,sig0,sig2,R1,E0,ax4) + + #Plot the Data (from Configuration 0) + gphy0 = ax2.plot(np.sqrt((MP0[0,0]-MP0[:,0])**2+(MP0[:,1]-MP0[0,1])**2),VtdMP0 + ,marker='o',color='blue',linewidth=3.,label ='Left Model Response' ) + + #Plot the Data (from Configuration 1) + gphy1 = ax2.plot(np.sqrt((MP1[0,0]-MP1[:,0])**2+(MP1[:,1]-MP1[0,1])**2),VtdMP1 + ,marker='o',color='red',linewidth=2.,label ='Right Model Response' ) + ax2.legend(('Left Model Response','Right Model Response'),loc=4) + + elif PlotOpt == 'Secondary': + #plot the secondary potentials + ax3= Plot_Secondary_Potential(XYZ,sig0,sig1,R0,E0,ax3) + ax4= Plot_Secondary_Potential(XYZ,sig0,sig2,R1,E0,ax4) + + #Plot the data(from configuration 0) + gphy0 = ax2.plot(np.sqrt((MP0[0,0]-MP0[:,0])**2+(MP0[:,1]-MP0[0,1])**2),VsdMP0,color='blue' + ,marker='o',linewidth=3.,label ='Left Model Response' ) + + + #Plot the Data (from Configuration 1) + gphy1 = ax2.plot(np.sqrt((MP1[0,0]-MP1[:,0])**2+(MP1[:,1]-MP1[0,1])**2),VsdMP1 + ,marker='o',color='red',linewidth=2.,label ='Right Model Response' ) + ax2.legend(('Left Model Response','Right Model Response'),loc=4 ) + + else: + print('What dont you get? Total or Secondary?') + + #Legends + ax3.plot(MP0[:,0],MP0[:,1],color='gray') + Dip_Midpoint0 = ax3.scatter(MP0[:,0],MP0[:,1],color='black') + Electrodes0 = ax3.scatter(EL0[:,0],EL0[:,1],color='red') + ax3.legend([Dip_Midpoint0,Electrodes0], ["Dipole Midpoint", "Electrodes"],scatterpoints=1) + + ax4.plot(MP1[:,0],MP1[:,1],color='gray') + Dip_Midpoint1 = ax4.scatter(MP1[:,0],MP1[:,1],color='black') + Electrodes1 = ax4.scatter(EL1[:,0],EL1[:,1],color='red') + ax4.legend([Dip_Midpoint1,Electrodes1], ["Dipole Midpoint", "Electrodes"],scatterpoints=1) + + return fig + +#Function to visualise and compare any two meaningful plots for the sphere in a uniform backgound with an unifom Electric Field +def interact_conductiveSphere(R,log_sig0,log_sig1,Figure1a,Figure1b,Figure2a,Figure2b): + + sig0,sig1 = conductivity_log_wrapper(log_sig0,log_sig1) + E0 = 1. # inducing field strength in V/m + n = 100 #level of discretisation + xr = np.linspace(-200., 200., n) # X-axis discretization + yr = xr.copy() # Y-axis discretization + zr = np.r_[0] # identical to saying `zr = np.array([0])` + XYZ = ndgrid(xr,yr,zr) # Space Definition + + fig, ax = plt.subplots(1,2,figsize=(18,6)) + + #Setup figure 1 with options Configuration, Total or Secondary, + #then Potential, ElectricField, Current Density or Charges Density + if Figure1a == 'Configuration': + ax[0] = get_Setup(XYZ,sig0,sig1,R,E0,ax[0],True,[0.1,0.1,0.6]) + + elif Figure1a == 'Total': + + if Figure1b == 'Potential': + ax[0] = Plot_Total_Potential(XYZ,sig0,sig1,R,E0,ax[0]) + + elif Figure1b == 'ElectricField': + ax[0] = Plot_Total_ElectricField(XYZ,sig0,sig1,R,E0,ax[0]) + + elif Figure1b == 'CurrentDensity': + ax[0] = Plot_Total_Currents(XYZ,sig0,sig1,R,E0,ax[0]) + + elif Figure1b == 'ChargesDensity': + ax[0] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[0]) + + elif Figure1a == 'Secondary': + + if Figure1b == 'Potential': + ax[0] = Plot_Secondary_Potential(XYZ,sig0,sig1,R,E0,ax[0]) + + elif Figure1b == 'ElectricField': + ax[0] = Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax[0]) + + elif Figure1b == 'CurrentDensity': + ax[0] = Plot_Secondary_Currents(XYZ,sig0,sig1,R,E0,ax[0]) + + elif Figure1b == 'ChargesDensity': + ax[0] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[0]) + + + if Figure1a== 'Configuration': + ax[1] = Plot_Primary_Potential(XYZ,sig0,sig1,R,E0,ax[1]) + print 'While figure1 is plotting Configuration, figure2 plots the primary field' + + elif Figure2a == 'Total': + if Figure2b == 'Potential': + ax[1] = Plot_Total_Potential(XYZ,sig0,sig1,R,E0,ax[1]) + + elif Figure2b == 'ElectricField': + ax[1] = Plot_Total_ElectricField(XYZ,sig0,sig1,R,E0,ax[1]) + + elif Figure2b == 'CurrentDensity': + ax[1]=Plot_Total_Currents(XYZ,sig0,sig1,R,E0,ax[1]) + + elif Figure2b == 'ChargesDensity': + ax[1] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[1]) + + + elif Figure2a == 'Secondary': + if Figure2b == 'Potential': + ax[1] = Plot_Secondary_Potential(XYZ,sig0,sig1,R,E0,ax[1]) + + elif Figure2b == 'ElectricField': + ax[1] = Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax[1]) + + elif Figure2b == 'CurrentDensity': + ax[1] = Plot_Secondary_Currents(XYZ,sig0,sig1,R,E0,ax[1]) + + elif Figure2b == 'ChargesDensity': + ax[1] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[1]) + + plt.tight_layout(True) + plt.show() + +#Interactive Visualisation of the responses of two configurations to a (pseudo) DC resistivity survey +def interactive_two_configurations_comparison(log_sig0,log_sig1,log_sig2,R0,R1,xstart,ystart,xend,yend,dipole_number,electrode_spacing,matching_spheres_example): + + sig0,sig1 = conductivity_log_wrapper(log_sig0,log_sig1) + sig2 = 10.**log_sig2 + E0 = 1. # inducing field strength in V/m + n = 100 #level of discretisation + xr = np.linspace(-200., 200., n) # X-axis discretization + yr = xr.copy() # Y-axis discretization + zr = np.r_[0] # identical to saying `zr = np.array([0])` + XYZ = ndgrid(xr,yr,zr) # Space Definition + PlotOpt = 'Total' + + if matching_spheres_example: + sig0 = 10.**(-3) + sig1 = 10.**(-2) + sig2 = 1.310344828 * 10**(-3) + R0 = 20. + R1 = 40. + + two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,dipole_number,electrode_spacing,PlotOpt) + + else: + two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,dipole_number,electrode_spacing,PlotOpt) + + plt.tight_layout(True) + plt.show() + + + +if __name__ == '__main__': + sig0 = -3. # conductivity of the wholespace + sig1 = -1. # conductivity of the sphere + sig0, sig1 = conductivity_log_wrapper(sig0,sig1) + R = 50. # radius of the sphere + E0 = 1. # inducing field strength + n = 100 #level of discretisation + xr = np.linspace(-2.*R, 2.*R, n) # X-axis discretization + yr = xr.copy() # Y-axis discretization + zr = np.r_[0] # identical to saying `zr = np.array([0])` + XYZ = ndgrid(xr,yr,zr) # Space Definition + + fig, ax = plt.subplots(2,5,figsize=(50,10)) + ax[0,0] = get_Setup(XYZ,sig0,sig1,R,E0,ax[0,0],True,[0.6,0.1,0.1]) + ax[1,0] = Plot_Primary_Potential(XYZ,sig0,sig1,R,E0,ax[1,0]) + ax[0,1] = Plot_Total_Potential(XYZ,sig0,sig1,R,E0,ax[0,1]) + ax[1,1] = Plot_Secondary_Potential(XYZ,sig0,sig1,R,E0,ax[1,1]) + ax[0,2] = Plot_Total_ElectricField(XYZ,sig0,sig1,R,E0,ax[0,2]) + ax[1,2] = Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax[1,2]) + ax[0,3] = Plot_Total_Currents(XYZ,sig0,sig1,R,E0,ax[0,3]) + ax[1,3] = Plot_Secondary_Currents(XYZ,sig0,sig1,R,E0,ax[1,3]) + ax[0,4] = Plot_Primary_Potential(XYZ,sig0,sig1,R,E0,ax[0,4]) + ax[1,4] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[1,4]) + + + plt.show() + From a9c9ce6bc8658317f3ad3e763b44cdc9db884f29 Mon Sep 17 00:00:00 2001 From: Thibaut Astic Date: Mon, 15 Feb 2016 19:50:04 -0800 Subject: [PATCH 11/18] remove ipywidget --- SimPEG/Examples/sphereElectrostatic_example.py | 1 - 1 file changed, 1 deletion(-) diff --git a/SimPEG/Examples/sphereElectrostatic_example.py b/SimPEG/Examples/sphereElectrostatic_example.py index 4b96f740..7ff1ead1 100644 --- a/SimPEG/Examples/sphereElectrostatic_example.py +++ b/SimPEG/Examples/sphereElectrostatic_example.py @@ -2,7 +2,6 @@ from scipy.constants import epsilon_0 import matplotlib.pyplot as plt import matplotlib.colors as colors import numpy as np -from ipywidgets import * from SimPEG.Utils import ndgrid, mkvc ''' From b765699d2fb4cc40bf26ba909f0bc093eb8307fc Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 25 Mar 2016 23:26:52 -0700 Subject: [PATCH 12/18] seperated out smallness and smoothness contributions --- SimPEG/Regularization.py | 198 +++++++++++++++++++-------------------- 1 file changed, 98 insertions(+), 100 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index c20ed973..a91da269 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -344,7 +344,6 @@ class BaseRegularization(object): def W(self): """Full regularization weighting matrix W.""" return sp.identity(self.regmesh.nC) - # self.regmesh._Pac.T * sp.identity(self.regmesh.nC) * self.regmesh._Pac # or do we want sp.identity(self.mesh.nC) or even just Utils.Identity() ? @Utils.timeIt def eval(self, m): @@ -375,11 +374,12 @@ class BaseRegularization(object): @Utils.timeIt def eval2Deriv(self, m, v=None): """ + Second derivative - :param numpy.array m: geophysical model - :param numpy.array v: vector to multiply - :rtype: scipy.sparse.csr_matrix or numpy.ndarray - :return: WtW or WtW*v + :param numpy.array m: geophysical model + :param numpy.array v: vector to multiply + :rtype: scipy.sparse.csr_matrix or numpy.ndarray + :return: WtW or WtW*v The regularization is: @@ -402,25 +402,48 @@ class BaseRegularization(object): class Tikhonov(BaseRegularization): """ + L2 Tikhonov regularization with both smallness and smoothness (first order + derivative) contributions. + + .. math:: + \phi_m(\mathbf{m}) = \\alpha_s \| W_s (\mathbf{m} - \mathbf{m_{ref}} ) \|^2 + + \\alpha_x \| W_x \\frac{\partial}{\partial x} (\mathbf{m} - \mathbf{m_{ref}} ) \|^2 + + \\alpha_y \| W_y \\frac{\partial}{\partial y} (\mathbf{m} - \mathbf{m_{ref}} ) \|^2 + + \\alpha_z \| W_z \\frac{\partial}{\partial z} (\mathbf{m} - \mathbf{m_{ref}} ) \|^2 + + Note if the key word argument `mrefInSmooth` is False, then mref is not + included in the smoothness contribution. + + :param Mesh mesh: SimPEG mesh + :param Maps mapping: regularization mapping, takes the model from model space to the thing you want to regularize + :param numpy.ndarray indActive: active cell indices for reducing the size of differential operators in the definition of a regularization mesh + :param bool mrefInSmooth: (default = False) put mref in the smoothness component? + :param float alpha_s: (default 1e-6) smallness weight + :param float alpha_x: (default 1) smoothness weight for first derivative in the x-direction + :param float alpha_y: (default 1) smoothness weight for first derivative in the y-direction + :param float alpha_z: (default 1) smoothness weight for first derivative in the z-direction + :param float alpha_xx: (default 1) smoothness weight for second derivative in the x-direction + :param float alpha_yy: (default 1) smoothness weight for second derivative in the y-direction + :param float alpha_zz: (default 1) smoothness weight for second derivative in the z-direction """ - mrefInSmooth = True #: SMOOTH and SMOOTH_MOD_DIF options - alpha_s = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Ws'], "Smallness weight") - alpha_x = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction") - alpha_y = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction") - alpha_z = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction") - alpha_xx = Utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'], "Weight for the second derivative in the x direction") - alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction") - alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction") + mrefInSmooth = False # put mref in the smoothness contribution + alpha_s = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Wsmall'], "Smallness weight") + alpha_x = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction") + alpha_y = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction") + alpha_z = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction") + alpha_xx = Utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'], "Weight for the second derivative in the x direction") + alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction") + alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction") def __init__(self, mesh, mapping=None, indActive = None, **kwargs): BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) @property - def Ws(self): - """Regularization matrix Ws""" - if getattr(self,'_Ws', None) is None: - self._Ws = Utils.sdiag((self.regmesh.vol*self.alpha_s)**0.5) - return self._Ws + def Wsmall(self): + """Regularization matrix Wsmall""" + if getattr(self,'_Wsmall', None) is None: + self._Wsmall = Utils.sdiag((self.regmesh.vol*self.alpha_s)**0.5) + return self._Wsmall @property def Wx(self): @@ -483,25 +506,44 @@ class Tikhonov(BaseRegularization): def W(self): """Full regularization matrix W""" if getattr(self, '_W', None) is None: - wlist = (self.Ws, self.Wsmooth) + wlist = (self.Wsmall, self.Wsmooth) self._W = sp.vstack(wlist) return self._W @Utils.timeIt - def eval(self, m): - if self.mrefInSmooth == True: - r1 = self.Wsmooth * ( self.mapping * (m) ) - r2 = self.Ws * ( self.mapping * (m - self.mref) ) - return 0.5*(r1.dot(r1)+r2.dot(r2)) - elif self.mrefInSmooth == False: - r = self.W * ( self.mapping * (m - self.mref) ) - return 0.5*r.dot(r) + def _evalSmall(self, m): + r = self.Wsmall * ( self.mapping * (m - self.mref) ) + return 0.5 * r.dot(r) + @Utils.timeIt + def _evalSmooth(self, m): + if self.mrefInSmooth == True: + r = self.Wsmooth * ( self.mapping * (m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wsmooth * ( self.mapping * (m) ) + return 0.5 * r.dot(r) + + @Utils.timeIt + def eval(self, m): + return self._evalSmall(m) + self._evalSmooth(m) + + @Utils.timeIt + def _evalSmallDeriv(self,m): + r = self.Wsmall * ( self.mapping * (m - self.mref) ) + return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) ) + + @Utils.timeIt + def _evalSmoothDeriv(self,m): + if self.mrefInSmooth == True: + r = self.Wsmooth * ( self.mapping * ( m - self.mref ) ) + return r.T * ( self.Wsmooth * self.mapping.deriv(m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wsmooth * ( self.mapping * m ) + return r.T * ( self.Wsmooth * self.mapping.deriv(m) ) @Utils.timeIt def evalDeriv(self, m): """ - The regularization is: .. math:: @@ -515,45 +557,33 @@ class Tikhonov(BaseRegularization): R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} """ - if self.mrefInSmooth == True: - mD1 = self.mapping.deriv(m) - mD2 = self.mapping.deriv(m - self.mref) - r1 = self.Wsmooth * ( self.mapping * (m)) - r2 = self.Ws * ( self.mapping * (m - self.mref) ) - out1 = mD1.T * ( self.Wsmooth.T * r1 ) - out2 = mD2.T * ( self.Ws.T * r2 ) - out = out1+out2 - elif self.mrefInSmooth == False: - mD = self.mapping.deriv(m - self.mref) - r = self.W * ( self.mapping * (m - self.mref) ) - out = mD.T * ( self.W.T * r ) - return out + return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) -class Simple(BaseRegularization): +class Simple(Tikhonov): """ - Only for tensor mesh + Simple regularization that does not include length scales in the derivatives. """ mrefInSmooth = True #: SMOOTH and SMOOTH_MOD_DIF options - alpha_s = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Ws'], "Smallness weight") + alpha_s = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "Smallness weight") alpha_x = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction") alpha_y = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction") alpha_z = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction") wght = 1. - + def __init__(self, mesh, mapping=None, indActive=None, **kwargs): BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) - + if isinstance(self.wght,float): self.wght = np.ones(self.regmesh.nC) * self.wght @property - def Ws(self): - """Regularization matrix Ws""" - if getattr(self,'_Ws', None) is None: - self._Ws = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5) - return self._Ws + def Wsmall(self): + """Regularization matrix Wsmall""" + if getattr(self,'_Wsmall', None) is None: + self._Wsmall = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5) + return self._Wsmall @property def Wx(self): @@ -592,54 +622,22 @@ class Simple(BaseRegularization): def W(self): """Full regularization matrix W""" if getattr(self, '_W', None) is None: - wlist = (self.Ws, self.Wsmooth) + wlist = (self.Wsmall, self.Wsmooth) self._W = sp.vstack(wlist) return self._W + @Utils.timeIt + def _evalSmall(self, m): + r = self.Wsmall * ( self.mapping * (m - self.mref) ) + return 0.5 * r.dot(r) @Utils.timeIt - def eval(self, m): + def _evalSmooth(self, m): if self.mrefInSmooth == True: - r1 = self.Wsmooth * ( self.mapping * (m) ) - r2 = self.Ws * ( self.mapping * (m - self.mref) ) - return 0.5*(r1.dot(r1)+r2.dot(r2)) + r = self.Wsmooth * ( self.mapping * (m - self.mref) ) elif self.mrefInSmooth == False: - r = self.W * ( self.mapping * (m - self.mref) ) - return 0.5*r.dot(r) - return phim - - - - @Utils.timeIt - def evalDeriv(self, m): - """ - - The regularization is: - - .. math:: - - R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})} - - So the derivative is straight forward: - - .. math:: - - R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} - - """ - if self.mrefInSmooth == True: - mD1 = self.mapping.deriv(m) - mD2 = self.mapping.deriv(m - self.mref) - r1 = self.Wsmooth * ( self.mapping * (m)) - r2 = self.Ws * ( self.mapping * (m - self.mref) ) - out1 = mD1.T * ( self.Wsmooth.T * r1 ) - out2 = mD2.T * ( self.Ws.T * r2 ) - out = out1+out2 - elif self.mrefInSmooth == False: - mD = self.mapping.deriv(m - self.mref) - r = self.W * ( self.mapping * (m - self.mref) ) - out = mD.T * ( self.W.T * r ) - return out + r = self.Wsmooth * ( self.mapping * m) + return 0.5 * r.dot(r) class Sparse(Simple): @@ -656,13 +654,13 @@ class Sparse(Simple): def __init__(self, mesh, mapping=None, indActive=None, **kwargs): Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) - + if isinstance(self.wght,float): self.wght = np.ones(self.regmesh.nC) * self.wght @property - def Ws(self): - """Regularization matrix Ws""" + def Wsmall(self): + """Regularization matrix Wsmall""" if getattr(self, 'curModel', None) is None: self.Rs = Utils.speye(self.regmesh.nC) @@ -671,14 +669,14 @@ class Sparse(Simple): self.rs = self.R(f_m , self.p) #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs)) self.Rs = Utils.sdiag( self.rs ) - + return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.wght)**0.5)*self.Rs @property def Wx(self): """Regularization matrix Wx""" - + if getattr(self, 'curModel', None) is None: self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0]) @@ -700,7 +698,7 @@ class Sparse(Simple): f_m = self.regmesh.cellDiffyStencil * self.curModel self.ry = self.R( f_m , self.qy) self.Ry = Utils.sdiag( self.ry ) - + return Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.Ry*self.regmesh.cellDiffyStencil @property @@ -733,10 +731,10 @@ class Sparse(Simple): def W(self): """Full regularization matrix W""" #if getattr(self, '_W', None) is None: - wlist = (self.Ws, self.Wsmooth) + wlist = (self.Wsmall, self.Wsmooth) #self._W = sp.vstack(wlist) return sp.vstack(wlist) - + def R(self, f_m , exponent): eta = (self.eps**(1-exponent/2.))**0.5 From 9c220ef37cde8c7faf1f67dfff622687c22e98dd Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 29 Mar 2016 23:05:04 -0700 Subject: [PATCH 13/18] default is mrefInSmooth = False --- SimPEG/Regularization.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index a91da269..e43670d2 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -565,7 +565,7 @@ class Simple(Tikhonov): Simple regularization that does not include length scales in the derivatives. """ - mrefInSmooth = True #: SMOOTH and SMOOTH_MOD_DIF options + mrefInSmooth = False #: SMOOTH and SMOOTH_MOD_DIF options alpha_s = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "Smallness weight") alpha_x = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction") alpha_y = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction") From 8f73b2e7be524a15411b4ab6765fae8b9fc5b81d Mon Sep 17 00:00:00 2001 From: D Fournier Date: Sun, 3 Apr 2016 17:18:39 -0700 Subject: [PATCH 14/18] Update directives Add IRLS example --- SimPEG/Directives.py | 101 +++++++++++++--------- SimPEG/Examples/Inversion_IRLS.py | 134 ++++++++++++++++++++---------- SimPEG/Regularization.py | 22 +++-- SimPEG/Survey.py | 2 +- 4 files changed, 160 insertions(+), 99 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index e5a63547..de7e0b6b 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -237,39 +237,41 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration): # Save the file as a npz np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred) -class SaveOutputDictEveryIteration(_SaveEveryIteration): - """SaveOutputDictEveryIteration - A directive that saves some relevant information from the inversion run to a numpy .npz dictionary file (see numpy.savez function for further info). - """ - - def initialize(self): - print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '%s-###.npz'"%self.fileName - - def endIter(self): - # Save the data. - ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) ) - phi_ms = 0.5*ms.dot(ms) - if self.reg.mrefInSmooth == True: - mref = self.reg.mref - else: - mref = 0 - mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) ) - phi_mx = 0.5 * mx.dot(mx) - if self.prob.mesh.dim==2: - my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) ) - phi_my = 0.5 * my.dot(my) - else: - phi_my = 'NaN' - if self.prob.mesh.dim==3 and 'CYL' not in self.prob.mesh._meshType: - mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) ) - phi_mz = 0.5 * mz.dot(mz) - else: - phi_mz = 'NaN' - - - # Save the file as a npz - np.savez('{:s}-{:03d}'.format(self.fileName,self.opt.iter), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred) - +#============================================================================== +# class SaveOutputDictEveryIteration(_SaveEveryIteration): +# """SaveOutputDictEveryIteration +# A directive that saves some relevant information from the inversion run to a numpy .npz dictionary file (see numpy.savez function for further info). +# """ +# +# def initialize(self): +# print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '%s-###.npz'"%self.fileName +# +# def endIter(self): +# # Save the data. +# ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) ) +# phi_ms = 0.5*ms.dot(ms) +# if self.reg.mrefInSmooth == True: +# mref = self.reg.mref +# else: +# mref = 0 +# mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) ) +# phi_mx = 0.5 * mx.dot(mx) +# if self.prob.mesh.dim==2: +# my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) ) +# phi_my = 0.5 * my.dot(my) +# else: +# phi_my = 'NaN' +# if self.prob.mesh.dim==3 and 'CYL' not in self.prob.mesh._meshType: +# mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) ) +# phi_mz = 0.5 * mz.dot(mz) +# else: +# phi_mz = 'NaN' +# +# +# # Save the file as a npz +# np.savez('{:s}-{:03d}'.format(self.fileName,self.opt.iter), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred) +# +#============================================================================== # class UpdateReferenceModel(Parameter): @@ -295,6 +297,8 @@ class update_IRLS(InversionDirective): # Scale the regularization for changes in norm if getattr(self, 'phi_m_last', None) is not None: + + self.reg.curModel = self.invProb.curModel self.reg.gamma = 1. phim_new = self.reg.eval(self.invProb.curModel) self.gamma = self.phi_m_last / phim_new @@ -320,22 +324,39 @@ class update_IRLS(InversionDirective): # Update the model used for the IRLS weights self.reg.curModel = self.invProb.curModel - - # Update the pre-conditioner - diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() * (self.reg.mapping * np.ones(self.reg.curModel.size))**2. - PC = Utils.sdiag(diagA**-1.) - self.opt.approxHinv = PC # Temporarely set gamma to 1. self.reg.gamma = 1. # Compute change in model objective function and update scaling phim_new = self.reg.eval(self.invProb.curModel) + self.reg.gamma = self.phi_m_last / phim_new - # TO DO: Re-scale beta if too much change in misfit - self.invProb.beta = self.invProb.beta * self.phi_d_last / self.invProb.phi_d + # TO DO: Check optimization class, data misfit not matching reality + #dpred = self.prob.fields(self.invProb.curModel) + #phid = self.invProb.dmisfit.eval(self.invProb.curModel) + #print self.survey.std[0] + #print phid + #print self.invProb.phi_d + #print self.invProb.phi_d_last + + self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d + +class update_lin_PreCond(InversionDirective): + + + def endIter(self): + # Cool the threshold parameter + + if getattr(self.opt, 'approxHinv', None) is not None: + # Update the pre-conditioner + diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() * (self.reg.mapping * np.ones(self.reg.curModel.size))**2. + PC = Utils.sdiag(diagA**-1.) + self.opt.approxHinv = PC + print 'Updated pre-cond' + #============================================================================== # import pylab as plt # plt.figure() diff --git a/SimPEG/Examples/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py index a7f5bd54..53240139 100644 --- a/SimPEG/Examples/Inversion_IRLS.py +++ b/SimPEG/Examples/Inversion_IRLS.py @@ -1,7 +1,7 @@ from SimPEG import * -def run(N=100, plotIt=True): +def run(N=200, plotIt=True): """ Inversion: Linear Problem ========================= @@ -9,39 +9,21 @@ def run(N=100, plotIt=True): Here we go over the basics of creating a linear problem and inversion. """ - - class LinearSurvey(Survey.BaseSurvey): - def projectFields(self, u): - return u - - class LinearProblem(Problem.BaseProblem): - - surveyPair = LinearSurvey - - def __init__(self, mesh, G, **kwargs): - Problem.BaseProblem.__init__(self, mesh, **kwargs) - self.G = G - - def fields(self, m, u=None): - return self.G.dot(m) - - def Jvec(self, m, v, u=None): - return self.G.dot(v) - - def Jtvec(self, m, v, u=None): - return self.G.T.dot(v) - - + + np.random.seed(1) + std_noise = 1e-2 + mesh = Mesh.TensorMesh([N]) + + m0 = np.ones(mesh.nC) * 1e-4 + nk = 10 + jk = np.linspace(1.,nk,nk) + p = -2. + q = 1. - nk = 20 - jk = np.linspace(1.,20.,nk) - p = -0.25 - q = 0.25 - - g = lambda k: np.exp(p*jk[k]*mesh.vectorCCx)*np.cos(2*np.pi*q*jk[k]*mesh.vectorCCx) + g = lambda k: np.exp(p*jk[k]*mesh.vectorCCx)*np.cos(np.pi*q*jk[k]*mesh.vectorCCx) G = np.empty((nk, mesh.nC)) @@ -52,25 +34,43 @@ def run(N=100, plotIt=True): mtrue[mesh.vectorCCx > 0.3] = 1. mtrue[mesh.vectorCCx > 0.45] = -0.5 mtrue[mesh.vectorCCx > 0.6] = 0 + - prob = LinearProblem(mesh, G) - survey = LinearSurvey() + prob = Problem.LinearProblem(mesh, G) + survey = Survey.LinearSurvey() survey.pair(prob) - survey.makeSyntheticData(mtrue, std=0.01) - - M = prob.mesh - - reg = Regularization.Tikhonov(mesh) + survey.dobs = prob.fields(mtrue) + std_noise * np.random.randn(nk) + #survey.makeSyntheticData(mtrue, std=std_noise) + + wd = np.ones(nk) * std_noise + + #print survey.std[0] + #M = prob.mesh + # Distance weighting + wr = np.sum(prob.G**2.,axis=0)**0.5 + wr = ( wr/np.max(wr) )**0 + + reg = Regularization.Simple(mesh) + reg.wght = wr + dmis = DataMisfit.l2_DataMisfit(survey) - opt = Optimization.ProjectedGNCG(maxIter=20,lower=-1.,upper=1., maxIterCG= 20, tolCG = 1e-3) + dmis.Wd = 1./wd + + opt = Optimization.ProjectedGNCG(maxIter=30,lower=-2.,upper=2., maxIterCG= 20, tolCG = 1e-4) invProb = InvProblem.BaseInvProblem(dmis, reg, opt) - beta = Directives.BetaSchedule() + invProb.curModel = m0 + + beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1) + target = Directives.TargetMisfit() + betaest = Directives.BetaEstimate_ByEig() - inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest]) - m0 = np.zeros_like(survey.mtrue) + inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target]) + mrec = inv.run(m0) - + + print "Final misfit:" + str(invProb.dmisfit.eval(mrec)) + if plotIt: import matplotlib.pyplot as plt @@ -78,12 +78,54 @@ def run(N=100, plotIt=True): for i in range(prob.G.shape[0]): axes[0].plot(prob.G[i,:]) axes[0].set_title('Columns of matrix G') - - axes[1].plot(M.vectorCCx, survey.mtrue, 'b-') - axes[1].plot(M.vectorCCx, mrec, 'r-') - axes[1].legend(('True Model', 'Recovered Model')) + + axes[1].plot(mesh.vectorCCx, mtrue, 'b-') + axes[1].plot(mesh.vectorCCx, mrec, 'r-') + #axes[1].legend(('True Model', 'Recovered Model')) + axes[1].set_ylim(-1.0,1.25) plt.show() + + # Switch regularization to sparse + phim = invProb.phi_m_last + phid = invProb.phi_d + + reg = Regularization.Sparse(mesh) + +#============================================================================== +# fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2)) +# dmdx = reg.mesh.cellDiffxStencil * mrec +# plt.plot(np.sort(dmdx)) +#============================================================================== + + #reg.recModel = mrec + reg.wght = np.ones(mesh.nC) + reg.mref = np.zeros(mesh.nC) + reg.eps_p = 2e-3 + reg.eps_q = 2e-3 + reg.norms = [0., 0., 2., 2.] + reg.wght = wr + + opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterCG= 200, tolCG = 1e-3) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.) + beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1) + #betaest = Directives.BetaEstimate_ByEig() + target = Directives.TargetMisfit() + IRLS =Directives.update_IRLS( phi_m_last = phim, phi_d_last = phid ) + + inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS]) + m0 = mrec + + # Run inversion + mrec = inv.run(m0) + + print "Final misfit:" + str(invProb.dmisfit.eval(mrec)) + + if plotIt: + axes[1].plot(mesh.vectorCCx, mrec, 'k-',lw = 2) + axes[1].legend(('True Model', 'Smooth l2-l2', + 'Sparse lp:' + str(reg.norms[0]) + ', lqx:' + str(reg.norms[1]) ), fontsize = 12) + return prob, survey, mesh, mrec if __name__ == '__main__': diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index ed75ada6..195e6203 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -643,13 +643,11 @@ class Simple(Tikhonov): class Sparse(Simple): # set default values - eps = 1e-1 + eps_p = 1e-1 + eps_q = 1e-1 curModel = None # use a model to compute the weights gamma = 1. - norms = [0., .2, 2., 2., 1.] - qx = 2. - qy = 2. - qz = 2. + norms = [0., 2., 2., 2.] wght = 1. def __init__(self, mesh, mapping=None, indActive=None, **kwargs): @@ -666,7 +664,7 @@ class Sparse(Simple): else: f_m = self.curModel - self.reg.mref - self.rs = self.R(f_m , self.norms[0]) + self.rs = self.R(f_m , self.eps_p, self.norms[0]) #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs)) self.Rs = Utils.sdiag( self.rs ) @@ -682,7 +680,7 @@ class Sparse(Simple): else: f_m = self.regmesh.cellDiffxStencil * self.curModel - self.rx = self.R( f_m , self.qx) + self.rx = self.R( f_m , self.eps_q, self.norms[1]) self.Rx = Utils.sdiag( self.rx ) return Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.Rx*self.regmesh.cellDiffxStencil @@ -696,7 +694,7 @@ class Sparse(Simple): else: f_m = self.regmesh.cellDiffyStencil * self.curModel - self.ry = self.R( f_m , self.qy) + self.ry = self.R( f_m , self.eps_q, self.norms[2]) self.Ry = Utils.sdiag( self.ry ) return Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.Ry*self.regmesh.cellDiffyStencil @@ -710,7 +708,7 @@ class Sparse(Simple): else: f_m = self.regmesh.cellDiffzStencil * self.curModel - self.rz = self.R( f_m , self.qz) + self.rz = self.R( f_m , self.eps_q, self.norms[3]) self.Rz = Utils.sdiag( self.rz ) return Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.Rz*self.regmesh.cellDiffzStencil @@ -735,9 +733,9 @@ class Sparse(Simple): #self._W = sp.vstack(wlist) return sp.vstack(wlist) - def R(self, f_m , exponent): + def R(self, f_m , eps, exponent): - eta = (self.eps**(1-exponent/2.))**0.5 - r = eta / (f_m**2.+self.eps**2.)**((1-exponent/2.)/2.) + eta = (eps**(1-exponent/2.))**0.5 + r = eta / (f_m**2.+ eps**2.)**((1-exponent/2.)/2.) return r diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 9f307c3f..eee2d538 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -372,7 +372,7 @@ class BaseSurvey(object): self.dtrue = self.dpred(m, u=u) noise = std*abs(self.dtrue)*np.random.randn(*self.dtrue.shape) self.dobs = self.dtrue+noise - self.std = self.dobs*0 + std + self.std = self.dobs*0. + std return self.dobs class LinearSurvey(BaseSurvey): From 2ee158e5d77f74603f58b98e054fdb4f343fa685 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Sun, 3 Apr 2016 17:26:04 -0700 Subject: [PATCH 15/18] Add distance weighting to example TO DO: Create example with and without distance weights --- SimPEG/Examples/Inversion_IRLS.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Examples/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py index 53240139..2cb75f31 100644 --- a/SimPEG/Examples/Inversion_IRLS.py +++ b/SimPEG/Examples/Inversion_IRLS.py @@ -48,7 +48,7 @@ def run(N=200, plotIt=True): #M = prob.mesh # Distance weighting wr = np.sum(prob.G**2.,axis=0)**0.5 - wr = ( wr/np.max(wr) )**0 + wr = ( wr/np.max(wr) ) reg = Regularization.Simple(mesh) reg.wght = wr From df620b42bd97fc8398806050b05881fbc995bbaa Mon Sep 17 00:00:00 2001 From: D Fournier Date: Tue, 5 Apr 2016 13:23:08 -0700 Subject: [PATCH 16/18] Fix plotting for Linear_IRLS example --- SimPEG/Examples/Inversion_IRLS.py | 34 +++++++++++++++---------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/SimPEG/Examples/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py index 2cb75f31..f5475680 100644 --- a/SimPEG/Examples/Inversion_IRLS.py +++ b/SimPEG/Examples/Inversion_IRLS.py @@ -68,23 +68,9 @@ def run(N=200, plotIt=True): mrec = inv.run(m0) - + ml2 = mrec print "Final misfit:" + str(invProb.dmisfit.eval(mrec)) - - if plotIt: - import matplotlib.pyplot as plt - - fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2)) - for i in range(prob.G.shape[0]): - axes[0].plot(prob.G[i,:]) - axes[0].set_title('Columns of matrix G') - - axes[1].plot(mesh.vectorCCx, mtrue, 'b-') - axes[1].plot(mesh.vectorCCx, mrec, 'r-') - #axes[1].legend(('True Model', 'Recovered Model')) - axes[1].set_ylim(-1.0,1.25) - plt.show() - + # Switch regularization to sparse phim = invProb.phi_m_last phid = invProb.phi_d @@ -105,7 +91,7 @@ def run(N=200, plotIt=True): reg.norms = [0., 0., 2., 2.] reg.wght = wr - opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterCG= 200, tolCG = 1e-3) + opt = Optimization.ProjectedGNCG(maxIter=5 ,lower=-2.,upper=2., maxIterCG= 100, tolCG = 1e-3) invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.) beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1) #betaest = Directives.BetaEstimate_ByEig() @@ -121,10 +107,24 @@ def run(N=200, plotIt=True): print "Final misfit:" + str(invProb.dmisfit.eval(mrec)) + if plotIt: + import matplotlib.pyplot as plt + + fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2)) + for i in range(prob.G.shape[0]): + axes[0].plot(prob.G[i,:]) + axes[0].set_title('Columns of matrix G') + + axes[1].plot(mesh.vectorCCx, mtrue, 'b-') + axes[1].plot(mesh.vectorCCx, ml2, 'r-') + #axes[1].legend(('True Model', 'Recovered Model')) + axes[1].set_ylim(-1.0,1.25) + axes[1].plot(mesh.vectorCCx, mrec, 'k-',lw = 2) axes[1].legend(('True Model', 'Smooth l2-l2', 'Sparse lp:' + str(reg.norms[0]) + ', lqx:' + str(reg.norms[1]) ), fontsize = 12) + plt.show() return prob, survey, mesh, mrec From ba977206eaa43d10b5df49c99631b4c78db9948c Mon Sep 17 00:00:00 2001 From: D Fournier Date: Wed, 6 Apr 2016 07:22:18 -0700 Subject: [PATCH 17/18] Create sensitivity re-weighting directive Adapt Map/polymap for actInd (topography) --- SimPEG/Directives.py | 31 ++++++++++++++++++++++++------ SimPEG/Maps.py | 45 ++++++++++++++++++++++++++++---------------- 2 files changed, 54 insertions(+), 22 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index de7e0b6b..ec77b7a1 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -356,10 +356,29 @@ class update_lin_PreCond(InversionDirective): PC = Utils.sdiag(diagA**-1.) self.opt.approxHinv = PC print 'Updated pre-cond' + +class update_Wj(InversionDirective): + """ + Create approx-sensitivity base weighting + """ + k = None -#============================================================================== -# import pylab as plt -# plt.figure() -# ax = plt.subplot(221) -# self.prob.mesh.plotSlice(self.invProb.curModel, ax = ax, normal = 'Z', ind=-5, clim = (0, 0.005)) -#============================================================================== + def endIter(self): + + if self.opt.iter == 2: + m = self.invProb.curModel + if self.k is None: + self.k = int(self.survey.nD/10) + + def JtJv(v): + + Jv = self.prob.Jvec(m, v) + + return self.prob.Jtvec(m,Jv) + + JtJdiag = Utils.diagEst(JtJv,6,k=100) + JtJdiag = JtJdiag + JtJdiag = JtJdiag / max(JtJdiag) + + + self.reg.wght = JtJdiag diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index 5f71d87c..f6b24f10 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -759,15 +759,27 @@ class PolyMap(IdentityMap): m = [\sigma_1, \sigma_2, c] - """ - def __init__(self, mesh, order, logSigma=True, normal='X'): + """ + def __init__(self, mesh, order, logSigma=True, normal='X', actInd = None): IdentityMap.__init__(self, mesh) self.logSigma = logSigma self.order = order - self.normal = normal - - slope = 1e4 - + self.normal = normal + self.actInd = actInd + + if getattr(self, 'actInd', None) is None: + self.actInd = range(self.mesh.nC) + self.nC = self.mesh.nC + + else: + self.nC = len(self.actInd) + + slope = 1e4 + + @property + def shape(self): + return (self.nC, self.nP) + @property def nP(self): if np.isscalar(self.order): @@ -785,8 +797,8 @@ class PolyMap(IdentityMap): sig1, sig2 = np.exp(sig1), np.exp(sig2) #2D if self.mesh.dim == 2: - X = self.mesh.gridCC[:,0] - Y = self.mesh.gridCC[:,1] + X = self.mesh.gridCC[self.actInd,0] + Y = self.mesh.gridCC[self.actInd,1] if self.normal =='X': f = polynomial.polyval(Y, c) - X elif self.normal =='Y': @@ -795,9 +807,9 @@ class PolyMap(IdentityMap): raise(Exception("Input for normal = X or Y or Z")) #3D elif self.mesh.dim == 3: - X = self.mesh.gridCC[:,0] - Y = self.mesh.gridCC[:,1] - Z = self.mesh.gridCC[:,2] + X = self.mesh.gridCC[self.actInd,0] + Y = self.mesh.gridCC[self.actInd,1] + Z = self.mesh.gridCC[self.actInd,2] if self.normal =='X': f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X elif self.normal =='Y': @@ -806,6 +818,7 @@ class PolyMap(IdentityMap): f = polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,self.order[1]+1))) - Z else: raise(Exception("Input for normal = X or Y or Z")) + else: raise(Exception("Only supports 2D")) @@ -819,8 +832,8 @@ class PolyMap(IdentityMap): sig1, sig2 = np.exp(sig1), np.exp(sig2) #2D if self.mesh.dim == 2: - X = self.mesh.gridCC[:,0] - Y = self.mesh.gridCC[:,1] + X = self.mesh.gridCC[self.actInd,0] + Y = self.mesh.gridCC[self.actInd,1] if self.normal =='X': f = polynomial.polyval(Y, c) - X @@ -832,9 +845,9 @@ class PolyMap(IdentityMap): raise(Exception("Input for normal = X or Y or Z")) #3D elif self.mesh.dim == 3: - X = self.mesh.gridCC[:,0] - Y = self.mesh.gridCC[:,1] - Z = self.mesh.gridCC[:,2] + X = self.mesh.gridCC[self.actInd,0] + Y = self.mesh.gridCC[self.actInd,1] + Z = self.mesh.gridCC[self.actInd,2] if self.normal =='X': f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X From 15d59a5b50511ae6396d7b74c455e68c1da44878 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Wed, 6 Apr 2016 09:00:13 -0700 Subject: [PATCH 18/18] Fix typos in Directives.update_Wj. overseer @lheagy --- SimPEG/Directives.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index ec77b7a1..9926055a 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -376,9 +376,7 @@ class update_Wj(InversionDirective): return self.prob.Jtvec(m,Jv) - JtJdiag = Utils.diagEst(JtJv,6,k=100) - JtJdiag = JtJdiag + JtJdiag = Utils.diagEst(JtJv,len(m),k=self.k) JtJdiag = JtJdiag / max(JtJdiag) - self.reg.wght = JtJdiag