diff --git a/simpegMT/Analytics/MT1Danalytic.py b/simpegMT/Analytics/MT1Danalytic.py deleted file mode 100644 index 66db5df5..00000000 --- a/simpegMT/Analytics/MT1Danalytic.py +++ /dev/null @@ -1,101 +0,0 @@ -# Analytic solution of EM fields due to a plane wave - -import numpy as np, SimPEG as simpeg - -def getEHfields(m1d,sigma,freq,zd): - '''Analytic solution for MT 1D layered earth. Returns E and H fields. - - :param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information. - :param numpy.array, vector sigma: Physical property of conductivity corresponding with the mesh. - :param float, freq: Frequency to calculate data at. - :param numpy array, vector zd: location to calculate EH fields at - - Assumes a halfspace with the same conductive as the last cell below. - - ''' - # Note add an error check for the mesh and sigma are the same size. - # Need make the solution e^-iwt dependent - - # Constants: Assume constant - mu = 4*np.pi*1e-7*np.ones((m1d.nC+1)) - eps = 8.85*1e-12*np.ones((m1d.nC+1)) - # Angular freq - w = 2*np.pi*freq - # Add the halfspace value to the property - sig = np.concatenate((np.array([sigma[0]]),sigma)) - # Calculate the wave number - k = np.sqrt(eps*mu*w**2-1j*mu*sig*w) - - # Initiate the propagation matrix, in the order down up. - UDp = np.zeros((2,m1d.nC+1),dtype=complex) - UDp[1,0] = 1 # Set the wave amplitude as 1 into the half-space at the bottom of the mesh - # Loop over all the layers, starting at the bottom layer - for lnr, h in enumerate(m1d.hx): # lnr-number of layer, h-thickness of the layer - # Calculate - yp1 = k[lnr]/(w*mu[lnr]) # Admittance of the layer below the current layer - zp = (w*mu[lnr+1])/k[lnr+1] # Impedance in the current layer - # Build the propagation matrix - - # Convert fields to down/up going components in layer below current layer - Pj1 = np.array([[1,1],[yp1,-yp1]]) - # Convert fields to down/up going components in current layer - Pjinv = 1./2*np.array([[1,zp],[1,-zp]]) - # Propagate down and up components through the current layer - elamh = np.array([[np.exp(-1j*k[lnr+1]*h),0],[0,np.exp(1j*k[lnr+1]*h)]]) - - # The down and up component in current layer. - UDp[:,lnr+1] = elamh.dot(Pjinv.dot(Pj1)).dot(UDp[:,lnr]) - - # Calculate the fields - Ed = np.zeros((zd.size,),dtype=complex) - Eu = np.zeros((zd.size,),dtype=complex) - Hd = np.zeros((zd.size,),dtype=complex) - Hu = np.zeros((zd.size,),dtype=complex) - - # Loop over the layers and calculate the fields - for ki,mui,epsi,dlow,dup,Up,Dp in zip(k[1::],mu[1::],eps[1::],m1d.vectorNx[:-1],m1d.vectorNx[1::],UDp[0,1::],UDp[1,1::]): - dind = np.logical_and(dup >= zd, zd > dlow) - Ed[dind] = Dp*np.exp(-1j*ki*(dup-zd[dind])) - Eu[dind] = Up*np.exp(1j*ki*(dup-zd[dind])) - Hd[dind] = (ki/(w*mui))*Dp*np.exp(-1j*ki*(dup-zd[dind])) - Hu[dind] = -(ki/(w*mui))*Up*np.exp(1j*ki*(dup-zd[dind])) - - # Return return the fields - return Ed, Eu, Hd, Hu - -def getImpedance(m1d,sigma,freq): - """Analytic solution for MT 1D layered earth. Returns the impedance at the surface. - - :param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information. - :param numpy.array, vector sigma: Physical property corresponding with the mesh. - :param numpy.array, vector freq: Frequencies to calculate data at. - - - """ - - # Define constants - mu0 = 4*np.pi*1e-7 - eps0 = 8.85e-12 - - # Initiate the impedances - Z1d = np.empty(len(freq) , dtype='complex') - h = m1d.hx #vectorNx[:-1] - # Start the process - for nrFr, fr in enumerate(freq): - om = 2*np.pi*fr - Zall = np.empty(len(h)+1,dtype='complex') - # Calculate the impedance for the bottom layer - Zall[0] = (mu0*om)/np.sqrt(mu0*eps0*(om)**2 - 1j*mu0*sigma[0]*om) - - for nr,hi in enumerate(h): - # Calculate the wave number - # print nr,sigma[nr] - k = np.sqrt(mu0*eps0*om**2 - 1j*mu0*sigma[nr]*om) - Z = (mu0*om)/k - - Zall[nr+1] = Z *((Zall[nr] + Z*np.tanh(1j*k*hi))/(Z + Zall[nr]*np.tanh(1j*k*hi))) - - #pdb.set_trace() - Z1d[nrFr] = Zall[-1] - - return Z1d \ No newline at end of file diff --git a/simpegMT/Analytics/MT1Dsolutions.py b/simpegMT/Analytics/MT1Dsolutions.py deleted file mode 100644 index a2900240..00000000 --- a/simpegMT/Analytics/MT1Dsolutions.py +++ /dev/null @@ -1,43 +0,0 @@ -import numpy as np, SimPEG as simpeg -from MT1Danalytic import getEHfields -from scipy.constants import mu_0 - -def get1DEfields(m1d,sigma,freq,sourceAmp=1.0): - """Function to get 1D electrical fields""" - - # Get the gradient - G = m1d.nodalGrad - # Mass matrices - # Magnetic permeability - Mmu = simpeg.Utils.sdiag(m1d.vol*(1.0/mu_0)) - # Conductivity - Msig = m1d.getFaceInnerProduct(sigma) - # Set up the solution matrix - A = G.T*Mmu*G - 1j*2.*np.pi*freq*Msig - # Define the inner part of the solution matrix - Aii = A[1:-1,1:-1] - # Define the outer part of the solution matrix - Aio = A[1:-1,[0,-1]] - - # Set the boundary conditions - Ed_low, Eu_low, Hd_low, Hu_low = getEHfields(m1d,sigma,freq,np.array([m1d.vectorNx[0]])) - Etot_low = Ed_low + Eu_low - ## Note: need to use conjugate of the analytic solution. It is derived with e^iwt - bc = np.r_[Etot_low.conj(),sourceAmp] - # The right hand side - rhs = -Aio*bc - # Solve the system - Aii_inv = simpeg.Solver(Aii) - eii = Aii_inv*rhs - # Assign the boundary conditions - e = np.r_[bc[0],eii,bc[1]] - # Return the electrical fields - return e - - -if __name__ == '__main__': - - hz = [(100.,18)] - M = simpeg.Mesh.TensorMesh([hz],'C') - sig = np.zeros(M.nC) + 1e-8 - sig[M.vectorCCx<=0] = sigHalf \ No newline at end of file diff --git a/simpegMT/Analytics/__init__.py b/simpegMT/Analytics/__init__.py deleted file mode 100644 index 92f785d8..00000000 --- a/simpegMT/Analytics/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from MT1Dsolutions import * # Add the names of the functions diff --git a/simpegMT/Utils/MT1Danalytic.py b/simpegMT/Utils/MT1Danalytic.py index f5605125..02406fc4 100644 --- a/simpegMT/Utils/MT1Danalytic.py +++ b/simpegMT/Utils/MT1Danalytic.py @@ -3,98 +3,99 @@ import numpy as np, SimPEG as simpeg def getEHfields(m1d,sigma,freq,zd): - '''Analytic solution for MT 1D layered earth. Returns E and H fields. + '''Analytic solution for MT 1D layered earth. Returns E and H fields. - :param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information. - :param numpy.array, vector sigma: Physical property of conductivity corresponding with the mesh. - :param float, freq: Frequency to calculate data at. - :param numpy array, vector zd: location to calculate EH fields at + :param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information. + :param numpy.array, vector sigma: Physical property of conductivity corresponding with the mesh. + :param float, freq: Frequency to calculate data at. + :param numpy array, vector zd: location to calculate EH fields at - Assumes a halfspace with the same conductive as the last cell below. + Assumes a halfspace with the same conductive as the last cell below. - ''' - # Note add an error check for the mesh. + ''' + # Note add an error check for the mesh and sigma are the same size. + # Need make the solution e^-iwt dependent - # Constants: Assume constant - mu = 4*np.pi*1e-7*np.ones((m1d.nC+1)) - eps = 8.85*1e-12*np.ones((m1d.nC+1)) - # Angular freq - w = 2*np.pi*freq - # Add the halfspace value to the property - sig = np.concatenate((np.array([sigma[0]]),sigma)) - # Calculate the wave number - k = np.sqrt(eps*mu*w**2-1j*mu*sig*w) + # Constants: Assume constant + mu = 4*np.pi*1e-7*np.ones((m1d.nC+1)) + eps = 8.85*1e-12*np.ones((m1d.nC+1)) + # Angular freq + w = 2*np.pi*freq + # Add the halfspace value to the property + sig = np.concatenate((np.array([sigma[0]]),sigma)) + # Calculate the wave number + k = np.sqrt(eps*mu*w**2-1j*mu*sig*w) - # Initiate the propagation matrix, in the order down up. - UDp = np.zeros((2,m1d.nC+1),dtype=complex) - UDp[1,0] = 1 # Set the wave amplitude as 1 into the half-space at the bottom of the mesh - # Loop over all the layers, starting at the bottom layer - for lnr, h in enumerate(m1d.hx): # lnr-number of layer, h-thickness of the layer - # Calculate - yp1 = k[lnr]/(w*mu[lnr]) # Admittance of the layer below the current layer - zp = (w*mu[lnr+1])/k[lnr+1] # Impedance in the current layer - # Build the propagation matrix - - # Convert fields to down/up going components in layer below current layer - Pj1 = np.array([[1,1],[yp1,-yp1]]) - # Convert fields to down/up going components in current layer - Pjinv = 1./2*np.array([[1,zp],[1,-zp]]) - # Propagate down and up components through the current layer - elamh = np.array([[np.exp(-1j*k[lnr+1]*h),0],[0,np.exp(1j*k[lnr+1]*h)]]) + # Initiate the propagation matrix, in the order down up. + UDp = np.zeros((2,m1d.nC+1),dtype=complex) + UDp[1,0] = 1 # Set the wave amplitude as 1 into the half-space at the bottom of the mesh + # Loop over all the layers, starting at the bottom layer + for lnr, h in enumerate(m1d.hx): # lnr-number of layer, h-thickness of the layer + # Calculate + yp1 = k[lnr]/(w*mu[lnr]) # Admittance of the layer below the current layer + zp = (w*mu[lnr+1])/k[lnr+1] # Impedance in the current layer + # Build the propagation matrix - # The down and up component in current layer. - UDp[:,lnr+1] = elamh.dot(Pjinv.dot(Pj1)).dot(UDp[:,lnr]) + # Convert fields to down/up going components in layer below current layer + Pj1 = np.array([[1,1],[yp1,-yp1]]) + # Convert fields to down/up going components in current layer + Pjinv = 1./2*np.array([[1,zp],[1,-zp]]) + # Propagate down and up components through the current layer + elamh = np.array([[np.exp(-1j*k[lnr+1]*h),0],[0,np.exp(1j*k[lnr+1]*h)]]) - # Calculate the fields - Ed = np.zeros((zd.size,),dtype=complex) - Eu = np.zeros((zd.size,),dtype=complex) - Hd = np.zeros((zd.size,),dtype=complex) - Hu = np.zeros((zd.size,),dtype=complex) + # The down and up component in current layer. + UDp[:,lnr+1] = elamh.dot(Pjinv.dot(Pj1)).dot(UDp[:,lnr]) - # Loop over the layers and calculate the fields - for ki,mui,epsi,dlow,dup,Up,Dp in zip(k[1::],mu[1::],eps[1::],m1d.vectorNx[:-1],m1d.vectorNx[1::],UDp[0,1::],UDp[1,1::]): - dind = np.logical_and(dup >= zd, zd > dlow) - Ed[dind] = Dp*np.exp(-1j*ki*(dup-zd[dind])) - Eu[dind] = Up*np.exp(1j*ki*(dup-zd[dind])) - Hd[dind] = (ki/(w*mui))*Dp*np.exp(-1j*ki*(dup-zd[dind])) - Hu[dind] = -(ki/(w*mui))*Up*np.exp(1j*ki*(dup-zd[dind])) + # Calculate the fields + Ed = np.zeros((zd.size,),dtype=complex) + Eu = np.zeros((zd.size,),dtype=complex) + Hd = np.zeros((zd.size,),dtype=complex) + Hu = np.zeros((zd.size,),dtype=complex) - # Return return the fields - return Ed, Eu, Hd, Hu + # Loop over the layers and calculate the fields + for ki,mui,epsi,dlow,dup,Up,Dp in zip(k[1::],mu[1::],eps[1::],m1d.vectorNx[:-1],m1d.vectorNx[1::],UDp[0,1::],UDp[1,1::]): + dind = np.logical_and(dup >= zd, zd > dlow) + Ed[dind] = Dp*np.exp(-1j*ki*(dup-zd[dind])) + Eu[dind] = Up*np.exp(1j*ki*(dup-zd[dind])) + Hd[dind] = (ki/(w*mui))*Dp*np.exp(-1j*ki*(dup-zd[dind])) + Hu[dind] = -(ki/(w*mui))*Up*np.exp(1j*ki*(dup-zd[dind])) + + # Return return the fields + return Ed, Eu, Hd, Hu def getImpedance(m1d,sigma,freq): - """Analytic solution for MT 1D layered earth. Returns the impedance at the surface. + """Analytic solution for MT 1D layered earth. Returns the impedance at the surface. - :param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information. - :param numpy.array, vector sigma: Physical property corresponding with the mesh. - :param numpy.array, vector freq: Frequencies to calculate data at. + :param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information. + :param numpy.array, vector sigma: Physical property corresponding with the mesh. + :param numpy.array, vector freq: Frequencies to calculate data at. - """ + """ - # Define constants - mu0 = 4*np.pi*1e-7 - eps0 = 8.85e-12 + # Define constants + mu0 = 4*np.pi*1e-7 + eps0 = 8.85e-12 - # Initiate the impedances - Z1d = np.empty(len(freq) , dtype='complex') - h = m1d.hx #vectorNx[:-1] - # Start the process - for nrFr, fr in enumerate(freq): - om = 2*np.pi*fr - Zall = np.empty(len(h)+1,dtype='complex') - # Calculate the impedance for the bottom layer - Zall[0] = (mu0*om)/np.sqrt(mu0*eps0*(om)**2 - 1j*mu0*sigma[0]*om) + # Initiate the impedances + Z1d = np.empty(len(freq) , dtype='complex') + h = m1d.hx #vectorNx[:-1] + # Start the process + for nrFr, fr in enumerate(freq): + om = 2*np.pi*fr + Zall = np.empty(len(h)+1,dtype='complex') + # Calculate the impedance for the bottom layer + Zall[0] = (mu0*om)/np.sqrt(mu0*eps0*(om)**2 - 1j*mu0*sigma[0]*om) - for nr,hi in enumerate(h): - # Calculate the wave number - # print nr,sigma[nr] - k = np.sqrt(mu0*eps0*om**2 - 1j*mu0*sigma[nr]*om) - Z = (mu0*om)/k - - Zall[nr+1] = Z *((Zall[nr] + Z*np.tanh(1j*k*hi))/(Z + Zall[nr]*np.tanh(1j*k*hi))) + for nr,hi in enumerate(h): + # Calculate the wave number + # print nr,sigma[nr] + k = np.sqrt(mu0*eps0*om**2 - 1j*mu0*sigma[nr]*om) + Z = (mu0*om)/k - #pdb.set_trace() - Z1d[nrFr] = Zall[-1] + Zall[nr+1] = Z *((Zall[nr] + Z*np.tanh(1j*k*hi))/(Z + Zall[nr]*np.tanh(1j*k*hi))) - return Z1d \ No newline at end of file + #pdb.set_trace() + Z1d[nrFr] = Zall[-1] + + return Z1d diff --git a/simpegMT/Utils/MT1Dsolutions.py b/simpegMT/Utils/MT1Dsolutions.py index 3b1149b8..fd3347d5 100644 --- a/simpegMT/Utils/MT1Dsolutions.py +++ b/simpegMT/Utils/MT1Dsolutions.py @@ -3,40 +3,41 @@ from MT1Danalytic import getEHfields from scipy.constants import mu_0 def get1DEfields(m1d,sigma,freq,sourceAmp=1.0): - """Function to get 1D electrical fields""" - - # Get the gradient - G = m1d.nodalGrad - # Mass matrices - # Magnetic permeability - Mmu = simpeg.Utils.sdiag(m1d.vol*(1.0/mu_0)) - # Conductivity - Msig = m1d.getFaceInnerProduct(sigma) - # Set up the solution matrix - A = G.T*Mmu*G - 1j*2.*np.pi*freq*Msig - # Define the inner part of the solution matrix - Aii = A[1:-1,1:-1] - # Define the outer part of the solution matrix - Aio = A[1:-1,[0,-1]] + """Function to get 1D electrical fields""" - # Set the boundary conditions - Ed_low, Eu_low, Hd_low, Hu_low = getEHfields(m1d,sigma,freq,np.array([m1d.vectorNx[0]])) - Etot_low = Ed_low + Eu_low - bc = np.r_[Etot_low,sourceAmp] - # The right hand side - rhs = -Aio*bc - # Solve the system - Aii_inv = simpeg.Solver(Aii) - eii = Aii_inv*rhs - # Assign the boundary conditions - e = np.r_[bc[0],eii,bc[1]] - # Return the electrical fields - return e + # Get the gradient + G = m1d.nodalGrad + # Mass matrices + # Magnetic permeability + Mmu = simpeg.Utils.sdiag(m1d.vol*(1.0/mu_0)) + # Conductivity + Msig = m1d.getFaceInnerProduct(sigma) + # Set up the solution matrix + A = G.T*Mmu*G - 1j*2.*np.pi*freq*Msig + # Define the inner part of the solution matrix + Aii = A[1:-1,1:-1] + # Define the outer part of the solution matrix + Aio = A[1:-1,[0,-1]] + + # Set the boundary conditions + Ed_low, Eu_low, Hd_low, Hu_low = getEHfields(m1d,sigma,freq,np.array([m1d.vectorNx[0]])) + Etot_low = Ed_low + Eu_low + ## Note: need to use conjugate of the analytic solution. It is derived with e^iwt + bc = np.r_[Etot_low.conj(),sourceAmp] + # The right hand side + rhs = -Aio*bc + # Solve the system + Aii_inv = simpeg.Solver(Aii) + eii = Aii_inv*rhs + # Assign the boundary conditions + e = np.r_[bc[0],eii,bc[1]] + # Return the electrical fields + return e if __name__ == '__main__': - hz = [(100.,18)] - M = simpeg.Mesh.TensorMesh([hz],'C') - sig = np.zeros(M.nC) + 1e-8 - sig[M.vectorCCx<=0] = sigHalf \ No newline at end of file + hz = [(100.,18)] + M = simpeg.Mesh.TensorMesh([hz],'C') + sig = np.zeros(M.nC) + 1e-8 + sig[M.vectorCCx<=0] = sigHalf diff --git a/simpegMT/Utils/__init__.py b/simpegMT/Utils/__init__.py index 4414cd54..73005269 100644 --- a/simpegMT/Utils/__init__.py +++ b/simpegMT/Utils/__init__.py @@ -1,2 +1,2 @@ from MT1Dsolutions import * # Add the names of the functions -from MT1Danalytic import * +from MT1Danalytic import *