diff --git a/SimPEG/DCIP/DCIPUtils.py b/SimPEG/DCIP/DCIPUtils.py index 97c30e89..8e1ca5e9 100644 --- a/SimPEG/DCIP/DCIPUtils.py +++ b/SimPEG/DCIP/DCIPUtils.py @@ -1,12 +1,16 @@ -from SimPEG import np +from SimPEG import np, Utils import BaseDC as DC import BaseDC as IP +import warnings def getActiveindfromTopo(mesh, topo): # def genActiveindfromTopo(mesh, topo): """ Get active indices from topography """ + warnings.warn( + "`getActiveindfromTopo` is deprecated and will be removed in future versions. Use `SimPEG.Utils.surface2ind_topo` instead", + FutureWarning) from scipy.interpolate import NearestNDInterpolator if mesh.dim==3: nCxy = mesh.nCx*mesh.nCy @@ -28,6 +32,9 @@ def gettopoCC(mesh, airind): """ Get topography from active indices of mesh. """ + warnings.warn( + "`gettopoCC` is deprecated and will be removed in future versions. Use `SimPEG.Utils.surface2ind_topo` instead", + FutureWarning) mesh2D = Mesh.TensorMesh([mesh.hx, mesh.hy], mesh.x0[:2]) zc = mesh.gridCC[:,2] AIRIND = airind.reshape((mesh.vnC[0]*mesh.vnC[1],mesh.vnC[2]), order='F') @@ -118,34 +125,27 @@ def readUBC_DC3Dobstopo(filename,mesh,topo,probType="CC"): def readUBC_DC2DModel(fileName): """ - Read UBC GIF 2DTensor model and generate 2D Tensor model in simpeg + Read UBC GIF 2DTensor model and generate 2D Tensor model in simpeg - Input: - :param fileName, path to the UBC GIF 2D model file - - Output: - :param SimPEG TensorMesh 2D object - :return - - Created on Thu Nov 12 13:14:10 2015 - - @author: dominiquef + :param string fileName: path to the UBC GIF 2D model file + :rtype: TensorMesh + :return: SimPEG TensorMesh 2D object """ from SimPEG import np, mkvc # Open fileand skip header... assume that we know the mesh already - obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!') + obsfile = np.genfromtxt(fileName, delimiter=' \n', dtype=np.str, comments='!') - dim = np.array(obsfile[0].split(),dtype=float) + dim = np.array(obsfile[0].split(), dtype=float) - temp = np.array(obsfile[1].split(),dtype=float) + temp = np.array(obsfile[1].split(), dtype=float) if len(temp) > 1: model = np.zeros(dim) for ii in range(len(obsfile)-1): - mm = np.array(obsfile[ii+1].split(),dtype=float) + mm = np.array(obsfile[ii+1].split(), dtype=float) model[:,ii] = mm model = model[:,::-1] @@ -153,10 +153,10 @@ def readUBC_DC2DModel(fileName): else: if len(obsfile[1:])==1: - mm = np.array(obsfile[1:].split(),dtype=float) + mm = np.array(obsfile[1:].split(), dtype=float) else: - mm = np.array(obsfile[1:],dtype=float) + mm = np.array(obsfile[1:], dtype=float) # Permute the second dimension to flip the order model = mm.reshape(dim[1],dim[0]) @@ -169,23 +169,19 @@ def readUBC_DC2DModel(fileName): return model -def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cblabel=True, axlabel = True, colorbar = True, contour = None): + +def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt', clim=None, cblabel=True, axlabel = True, colorbar = True, contour = None): """ - Read list of 2D tx-rx location and plot a speudo-section of apparent - resistivity. + Read list of 2D tx-rx location and plot a speudo-section of apparent + resistivity. - Assumes flat topo for now... + Assumes flat topo for now... - Input: - :param d2D, z0 - :switch stype -> Either 'pdp' (pole-dipole) | 'dpdp' (dipole-dipole) - :switch dtype=-> Either 'appr' (app. res) | 'appc' (app. con) | 'volt' (potential) - Output: - :figure scatter plot overlayed on image - - Edited Feb 17th, 2016 - - @author: dominiquef + :param SurveyDC DCsurvey: + :param string surveyType: Either 'pole-dipole' | 'dipole-dipole' + :param string unitType: Either 'appResistivity' | 'appConductivity' | 'volt' + :rtype: matplotlib.plt + :return: figure scatter plot overlayed on image """ from SimPEG import np @@ -218,39 +214,39 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cbl Cmid = (Tx[0][0] + Tx[1][0])/2 Pmid = (Rx[0][:,0] + Rx[1][:,0])/2 - # Change output for dtype - if dtype == 'volt': + # Change output for unitType + if unitType == 'volt': rho = np.hstack([rho,data]) else: # Compute pant leg of apparent rho - if stype == 'pdp': + if surveyType == 'pole-dipole': leg = data * 2*np.pi * MA * ( MA + MN ) / MN - elif stype == 'dpdp': + elif surveyType == 'dipole-dipole': leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA ) else: - print """dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """ + print """unitType must be 'pole-dipole' | 'dipole-dipole' """ break - if dtype == 'appc': + if unitType == 'appConductivity': leg = np.log10(abs(1./leg)) rho = np.hstack([rho,leg]) - elif dtype == 'appr': + elif unitType == 'appResistivity': leg = np.log10(abs(leg)) rho = np.hstack([rho,leg]) else: - print """dtype must be 'appr' | 'appc' | 'volt' """ + print """unitType must be 'appResistivity' | 'appConductivity' | 'volt' """ break midx = np.hstack([midx, ( Cmid + Pmid )/2 ]) @@ -259,7 +255,7 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cbl # Grid points grid_x, grid_z = np.mgrid[np.min(midx):np.max(midx), np.min(midz):np.max(midz)] grid_rho = griddata(np.c_[midx,midz], rho.T, (grid_x, grid_z), method='linear') - + # Scale the color scheme if clim == None: vmin, vmax = rho.min(), rho.max() @@ -268,36 +264,37 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cbl # Plot data grid_rho = np.ma.masked_where(np.isnan(grid_rho), grid_rho) + + + ph = plt.pcolormesh(grid_x[:,0],grid_z[0,:],grid_rho.T, vmin = vmin, vmax = vmax) plt.gca().tick_params(axis='both', which='major', labelsize=8) - + if contour is not None: plt.contour(grid_x,grid_z,grid_rho,levels = contour,colors = 'r', vmin = vmin, vmax = vmax) - + # Add scatter points axs.scatter(midx,midz,s=10,c=rho.T, vmin = vmin, vmax = vmax) - + if colorbar: - - if dtype == 'volt': + + if unitType == 'volt': cbar = plt.colorbar(ph, ax = axs, format="%4.1f",fraction=0.04,orientation="horizontal") - else: + else: cbar = plt.colorbar(ph, ax = axs, format="$10^{%.1f}$",fraction=0.04,orientation="horizontal") - - cmin,cmax = cbar.get_clim() - ticks = np.linspace(cmin,cmax,3) - cbar.set_ticks(ticks) - cbar.ax.tick_params(labelsize=10) - - if cblabel: - if dtype == 'appc': - cbar.set_label("App.Cond",size=12) - elif dtype == 'appr': - cbar.set_label("App.Res.",size=12) - elif dtype == 'volt': - cbar.set_label("Potential (V)",size=12) + cmin,cmax = cbar.get_clim() + ticks = np.linspace(cmin,cmax,3) + cbar.set_ticks(ticks) + cbar.ax.tick_params(labelsize=10) + + if unitType == 'appConductivity': + cbar.set_label("App.Cond",size=12) + elif unitType == 'appResistivity': + cbar.set_label("App.Res.",size=12) + elif unitType == 'volt': + cbar.set_label("Potential (V)",size=12) if not axlabel: @@ -310,27 +307,24 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cbl return ph -def gen_DCIPsurvey(endl, mesh, stype, a, b, n): +def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx): """ - Load in endpoints and survey specifications to generate Tx, Rx location - stations. + Load in endpoints and survey specifications to generate Tx, Rx location + stations. - Assumes flat topo for now... + Assumes flat topo for now... - Input: - :param endl -> input endpoints [x1, y1, z1, x2, y2, z2] - :object mesh -> SimPEG mesh object - :switch stype -> "dpdp" (dipole-dipole) | "pdp" (pole-dipole) | 'gradient' - : param a, n -> pole seperation, number of rx dipoles per tx + :param numpy.array endl: input endpoints [[x1, y1] , [x2, y2]] + :param Mesh mesh: SimPEG mesh object + :param string surveyType: 'dipole-dipole' | 'pole-dipole' | 'gradient' + :param float AM_sep: transmitter (A) - receiver (M) seperation + :param float b: receiver dipole seperation + :param float nrx: pole seperation, number of rx dipoles per tx - Output: - :param Tx, Rx -> List objects for each tx location - Lines: P1x, P1y, P1z, P2x, P2y, P2z + :rtype: DC.Survey, Src, Rx + :returns: DC survey, Source - Created on Wed December 9th, 2015 - - @author: dominiquef - !! Require clean up to deal with DCsurvey + !! Require clean up to deal with DCsurvey """ from SimPEG import np @@ -346,17 +340,17 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n): dl_x = ( endl[1,0] - endl[0,0] ) / dl_len dl_y = ( endl[1,1] - endl[0,1] ) / dl_len - nstn = np.floor( dl_len / a ) + nstn = np.floor( dl_len / AM_sep ) # Compute discrete pole location along line - stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*a - stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*a + stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*AM_sep + stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*AM_sep # Create line of P1 locations M = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]] # Create line of P2 locations - N = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]] + N = np.c_[stn_x+AM_sep*dl_x, stn_y+AM_sep*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]] ## Build list of Tx-Rx locations depending on survey type # Dipole-dipole: Moving tx with [a] spacing -> [AB a MN1 a MN2 ... a MNn] @@ -366,14 +360,14 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n): SrcList = [] - if stype != 'gradient': + if surveyType != 'gradient': for ii in range(0, int(nstn)-1): - if stype == 'dpdp': + if surveyType == 'dipole-dipole': tx = np.c_[M[ii,:],N[ii,:]] - elif stype == 'pdp': + elif surveyType == 'pole-dipole': tx = np.c_[M[ii,:],M[ii,:]] # Rx.append(np.c_[M[ii+1:indx,:],N[ii+1:indx,:]]) @@ -382,33 +376,33 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n): AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1]) # Number of receivers to fit - nstn = np.min([np.floor( (AB - b) / a ) , n]) + nstn = np.min([np.floor( (AB - MN_sep) / AM_sep ) , nrx]) # Check if there is enough space, else break the loop if nstn <= 0: continue # Compute discrete pole location along line - stn_x = N[ii,0] + dl_x*b + np.array(range(int(nstn)))*dl_x*a - stn_y = N[ii,1] + dl_y*b + np.array(range(int(nstn)))*dl_y*a + stn_x = N[ii,0] + dl_x*MN_sep + np.array(range(int(nstn)))*dl_x*AM_sep + stn_y = N[ii,1] + dl_y*MN_sep + np.array(range(int(nstn)))*dl_y*AM_sep # Create receiver poles # Create line of P1 locations P1 = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]] # Create line of P2 locations - P2 = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]] + P2 = np.c_[stn_x+AM_sep*dl_x, stn_y+AM_sep*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]] Rx.append(np.c_[P1,P2]) rxClass = DC.RxDipole(P1, P2) Tx.append(tx) - if stype == 'dpdp': + if surveyType == 'dipole-dipole': srcClass = DC.SrcDipole([rxClass], M[ii,:],N[ii,:]) - elif stype == 'pdp': + elif surveyType == 'pole-dipole': srcClass = DC.SrcDipole([rxClass], M[ii,:],M[ii,:]) SrcList.append(srcClass) - elif stype == 'gradient': + elif surveyType == 'gradient': # Gradient survey only requires Tx at end of line and creates a square # grid of receivers at in the middle at a pre-set minimum distance @@ -416,23 +410,23 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n): Tx.append(np.c_[M[0,:],N[-1,:]]) # Get the edge limit of survey area - min_x = endl[0,0] + dl_x * b - min_y = endl[0,1] + dl_y * b + min_x = endl[0,0] + dl_x * MN_sep + min_y = endl[0,1] + dl_y * MN_sep - max_x = endl[1,0] - dl_x * b - max_y = endl[1,1] - dl_y * b + max_x = endl[1,0] - dl_x * MN_sep + max_y = endl[1,1] - dl_y * MN_sep box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 ) box_w = box_l/2. - nstn = np.floor( box_l / a ) + nstn = np.floor( box_l / AM_sep ) # Compute discrete pole location along line - stn_x = min_x + np.array(range(int(nstn)))*dl_x*a - stn_y = min_y + np.array(range(int(nstn)))*dl_y*a + stn_x = min_x + np.array(range(int(nstn)))*dl_x*AM_sep + stn_y = min_y + np.array(range(int(nstn)))*dl_y*AM_sep # Define number of cross lines - nlin = int(np.floor( box_w / a )) + nlin = int(np.floor( box_w / AM_sep )) lind = range(-nlin,nlin+1) ngrad = nstn * len(lind) @@ -441,12 +435,12 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n): for ii in range( len(lind) ): # Move line in perpendicular direction by dipole spacing - lxx = stn_x - lind[ii]*a*dl_y - lyy = stn_y + lind[ii]*a*dl_x + lxx = stn_x - lind[ii]*AM_sep*dl_y + lyy = stn_y + lind[ii]*AM_sep*dl_x M = np.c_[ lxx, lyy , np.ones(nstn).T*mesh.vectorNz[-1]] - N = np.c_[ lxx+a*dl_x, lyy+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]] + N = np.c_[ lxx+AM_sep*dl_x, lyy+AM_sep*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]] rx[(ii*nstn):((ii+1)*nstn),:] = np.c_[M,N] @@ -455,44 +449,38 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n): srcClass = DC.SrcDipole([rxClass], M[0,:], N[-1,:]) SrcList.append(srcClass) else: - print """stype must be either 'pdp', 'dpdp' or 'gradient'. """ + print """surveyType must be either 'pole-dipole', 'dipole-dipole' or 'gradient'. """ survey = DC.SurveyDC(SrcList) return survey, Tx, Rx -def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0): + +def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0): """ Write UBC GIF DCIP 2D or 3D observation file - Input: - :string fileName -> including path where the file is written out - :DCsurvey DC survey class object - :string dtype -> either '2D' | '3D' - :string stype -> either 'SURFACE' | 'GENERAL' - - Output: - :param UBC2D-Data file - :return - - Last edit: February 16th, 2016 - - @author: dominiquef - + :param string fileName: including path where the file is written out + :param Survey DCsurvey: DC survey class object + :param string dim: either '2D' | '3D' + :param string surveyType: either 'SURFACE' | 'GENERAL' + :rtype: file + :return: UBC2D-Data file """ + from SimPEG import mkvc - assert (dtype=='2D') | (dtype=='3D'), "Data must be either '2D' | '3D'" - assert (stype=='SURFACE') | (stype=='GENERAL') | (stype=='SIMPLE'), "Data must be either 'SURFACE' | 'GENERAL' | 'SIMPLE'" - + assert (dim=='2D') | (dim=='3D'), "Data must be either '2D' | '3D'" + assert (surveyType=='SURFACE') | (surveyType=='GENERAL') | (surveyType=='SIMPLE'), "Data must be either 'SURFACE' | 'GENERAL' | 'SIMPLE'" + fid = open(fileName,'w') - - + fid.write('! ' + surveyType + ' FORMAT\n') + if iptype!=0: fid.write('IPTYPE=%i\n'%iptype) - + else: fid.write('! ' + stype + ' FORMAT\n') - + count = 0 for ii in range(DCsurvey.nSrc): @@ -506,33 +494,33 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0): M = rx[0] N = rx[1] - # Adapt source-receiver location for dtype and stype - if dtype=='2D': + # Adapt source-receiver location for dim and surveyType + if dim=='2D': - if stype == 'SIMPLE': + if surveyType == 'SIMPLE': #fid.writelines("%e " % ii for ii in mkvc(tx[0,:])) A = np.repeat(tx[0,0],M.shape[0],axis=0) B = np.repeat(tx[0,1],M.shape[0],axis=0) M = M[:,0] N = N[:,0] - + np.savetxt(fid, np.c_[A, B, M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n') else: - if stype == 'SURFACE': + if surveyType == 'SURFACE': fid.writelines("%f " % ii for ii in mkvc(tx[0,:])) M = M[:,0] N = N[:,0] - if stype == 'GENERAL': + if surveyType == 'GENERAL': # Flip sign for z-elevation to depth tx[2::2,:] = -tx[2::2,:] - + fid.writelines("%e " % ii for ii in mkvc(tx[::2,:])) M = M[:,0::2] N = N[:,0::2] @@ -540,31 +528,31 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0): # Flip sign for z-elevation to depth M[:,1::2] = -M[:,1::2] N[:,1::2] = -N[:,1::2] - + fid.write('%i\n'% nD) np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%f',delimiter=' ',newline='\n') - if dtype=='3D': + if dim=='3D': - if stype == 'SURFACE': + if surveyType == 'SURFACE': fid.writelines("%e " % ii for ii in mkvc(tx[0:2,:])) M = M[:,0:2] N = N[:,0:2] - if stype == 'GENERAL': + if surveyType == 'GENERAL': fid.writelines("%e " % ii for ii in mkvc(tx[0:3,:])) fid.write('%i\n'% nD) np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n') fid.write('\n') - + count += nD fid.close() -def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'): +def convertObs_DC3D_to_2D(DCsurvey, lineID, flag='local'): """ Read DC survey and projects the coordinate system according to the flag = 'Xloc' | 'Yloc' | 'local' (default) @@ -573,15 +561,9 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'): The Z value is preserved, but Y coordinates zeroed. - Input: - :param survey3D - - Output: - :figure survey2D - - Edited April 6th, 2016 - - @author: dominiquef + :param DC.Survey survey3D: 3D simpeg DC survey + :rtype: DC.Survey + :return: survey2D """ from SimPEG import np @@ -666,39 +648,34 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'): DCsurvey2D.std = np.asarray(DCsurvey.std) return DCsurvey2D - -def readUBC_DC3Dobs(fileName, dtype = 'DC'): + +def readUBC_DC3Dobs(fileName, rtype = 'DC'): """ Read UBC GIF IP 3D observation file and generate survey - Input: - :param fileName, path to the UBC GIF 3D obs file - - Output: - :param IPsurvey - :return - - @author: dominiquef + :param string fileName:, path to the UBC GIF 3D obs file + :rtype: Survey + :return: DCIPsurvey """ zflag = True # Flag for z value provided - + # Load file - if dtype == 'IP': + if rtype == 'IP': obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='IPTYPE') - - elif dtype == 'DC': + + elif rtype == 'DC': obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!') - + else: - print "dtype must be 'DC'(default) | 'IP'" - + print "rtype must be 'DC'(default) | 'IP'" + # Pre-allocate srcLists = [] Rx = [] d = [] wd = [] - + # Countdown for number of obs/tx count = 0 @@ -717,7 +694,7 @@ def readUBC_DC3Dobs(fileName, dtype = 'DC'): # Check if z value is provided, if False -> nan if len(temp)==5: tx = np.r_[temp[0:2],np.nan,temp[2:4],np.nan] - + zflag = False # Pass on the flag to the receiver loc else: @@ -729,12 +706,12 @@ def readUBC_DC3Dobs(fileName, dtype = 'DC'): temp = np.fromstring(obsfile[ii], dtype=float,sep=' ') # Get the string # Filter out negative IP -# if temp[-2] < 0: +# if temp[-2] < 0: # count = count -1 # print "Negative!" -# +# # else: - + # If the Z-location is provided, otherwise put nan if zflag: @@ -772,17 +749,9 @@ def readUBC_DC2Dobs(fileName): ------- NEEDS TO BE UPDATED ------ Read UBC GIF 2D observation file and generate arrays for tx-rx location - Input: - :param fileName, path to the UBC GIF 2D model file - - Output: - :param rx, tx - :return - - Created on Thu Nov 12 13:14:10 2015 - - @author: dominiquef - + :param string fileName: path to the UBC GIF 2D model file + :rtype: (DC.Src, DC.Rx, ??, ??) + :return: source_locs, rx_locs, ??, ?? """ from SimPEG import np @@ -822,11 +791,9 @@ def readUBC_DC2Dpre(fileName): Read UBC GIF DCIP 2D observation file and generate arrays for tx-rx location Input: - :param fileName, path to the UBC GIF 3D obs file - - Output: - DCsurvey - :return + :param string fileName: path to the UBC GIF 3D obs file + :rtype: DC.Survey + :return: DCsurvey Created on Mon March 9th, 2016 << Doug's 70th Birthday !! >> @@ -888,12 +855,9 @@ def readUBC_DC2DMesh(fileName): """ Read UBC GIF 2DTensor mesh and generate 2D Tensor mesh in simpeg - Input: - :param fileName, path to the UBC GIF mesh file - - Output: - :param SimPEG TensorMesh 2D object - :return + :param string fileName: path to the UBC GIF mesh file + :rtype: Mesh.TensorMesh + :return: SimPEG TensorMesh 2D object Created on Thu Nov 12 13:14:10 2015 @@ -959,12 +923,9 @@ def xy_2_lineID(DCsurvey): they were collected. May need to generalize for random point locations, but will be more expensive - Input: - :param DCdict Vectors of station location - - Output: - :param LineID Vector of integers - :return + :param numpy.array DCdict: Vectors of station location + :rtype: numpy.array + :return: LineID Vector of integers Created on Thu Feb 11, 2015 diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index e51c904e..93e51f62 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -146,10 +146,15 @@ class BetaSchedule(InversionDirective): class TargetMisfit(InversionDirective): + chifact = 1. + phi_d_star = None + @property def target(self): if getattr(self, '_target', None) is None: - self._target = self.survey.nD*0.5 + if self.phi_d_star is None: + self.phi_d_star = 0.5 * self.survey.nD + self._target = self.chifact * self.phi_d_star # the factor of 0.5 is because we do phid = 0.5*|| dpred - dobs||^2 return self._target @target.setter def target(self, val): @@ -309,24 +314,24 @@ class Update_lin_PreCond(InversionDirective): Create a Jacobi preconditioner for the linear problem """ onlyOnStart=False - + def initialize(self): - + if getattr(self.opt, 'approxHinv', None) is 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.) + PC = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.) self.opt.approxHinv = PC - + def endIter(self): # Cool the threshold parameter if self.onlyOnStart==True: return - + 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.) + PC = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.) self.opt.approxHinv = PC diff --git a/SimPEG/EM/Analytics/DC.py b/SimPEG/EM/Analytics/DC.py new file mode 100644 index 00000000..f67843d6 --- /dev/null +++ b/SimPEG/EM/Analytics/DC.py @@ -0,0 +1,118 @@ +import numpy as np +from scipy.constants import mu_0, pi +from scipy import special + +def DCAnalyticHalf(txloc, rxlocs, sigma, earth_type="wholespace"): + """ + Analytic solution for electric potential from a postive pole + + :param array txloc: a xyz location of A (+) electrode (np.r_[xa, ya, za]) + :param list rxlocs: xyz locations of M (+) and N (-) electrodes [M, N] + + e.g. + rxlocs = [M, N] + M: xyz locations of M (+) electrode (np.c_[xmlocs, ymlocs, zmlocs]) + N: xyz locations of N (-) electrode (np.c_[xnlocs, ynlocs, znlocs]) + + :param float or complex sigma: values of conductivity + :param string earth_type: values of conductivity ("wholsespace" or "halfspace") + + """ + M = rxlocs[0] + N = rxlocs[1] + + rM = np.sqrt( (M[:,0]-txloc[0])**2 + (M[:,1]-txloc[1])**2 + (M[:,2]-txloc[1])**2 ) + rN = np.sqrt( (N[:,0]-txloc[0])**2 + (N[:,1]-txloc[1])**2 + (N[:,2]-txloc[1])**2 ) + + phiM = 1./(4*np.pi*rM*sigma) + phiN = 1./(4*np.pi*rN*sigma) + phi = phiM - phiN + + if earth_type == "halfspace": + phi *= 2 + + return phi + +deg2rad = lambda deg: deg/180.*np.pi +rad2deg = lambda rad: rad*180./np.pi + +def DCAnalyticSphere(txloc, rxloc, xc, radius, sigma, sigma1, \ + field_type = "secondary", order=12, halfspace=False): +# def DCSpherePointCurrent(txloc, rxloc, xc, radius, rho, rho1, \ +# field_type = "secondary", order=12): + """ + + Parameters: + + :param array txloc: A (+) current electrode location (x,y,z) + :param array xc: x center of depressed sphere + :param array rxloc: M(+) electrode locations / (Nx3 array, # of electrodes) + + :param float radius: radius (float): radius of the sphere (m) + :param float rho: resistivity of the background (ohm-m) + :param float rho1: resistivity of the sphere + :param string field_type: : "secondary", "total", "primary" + (default="secondary") + "secondary": secondary potential only due to sphere + "primary": primary potential from the point source + "total": "secondary"+"primary" + :param float order: maximum order of Legendre polynomial (default=12) + + Written by Seogi Kang (skang@eos.ubc.ca) + Ph.D. Candidate of University of British Columbia, Canada + + """ + + Pleg = [] + # Compute Legendre Polynomial + for i in range(order): + Pleg.append(special.legendre(i, monic=0)) + + + rho = 1./sigma + rho1 = 1./sigma1 + + # Center of the sphere should be aligned in txloc in y-direction + yc = txloc[1] + xyz = np.c_[rxloc[:,0]-xc, rxloc[:,1]-yc, rxloc[:,2]] + r = np.sqrt( (xyz**2).sum(axis=1) ) + + x0 = abs(txloc[0]-xc) + + costheta = xyz[:,0]/r * (txloc[0]-xc)/x0 + phi = np.zeros_like(r) + R = (r**2+x0**2.-2.*r*x0*costheta)**0.5 + # primary potential in a whole space + prim = rho*1./(4*np.pi*R) + + if field_type =="primary": + return prim + + sphind = r < radius + out = np.zeros_like(r) + for n in range(order): + An, Bn = AnBnfun(n, radius, x0, rho, rho1) + dumout = An*r[~sphind]**(-n-1.)*Pleg[n](costheta[~sphind]) + out[~sphind] += dumout + dumin = Bn*r[sphind]**(n)*Pleg[n](costheta[sphind]) + out[sphind] += dumin + + out[~sphind] += prim[~sphind] + + if halfspace: + scale = 2 + else: + scale = 1 + + if field_type == "secondary": + return scale*(out-prim) + elif field_type == "total": + return scale*out + +def AnBnfun(n, radius, x0, rho, rho1, I=1.): + const = I*rho/(4*np.pi) + bunmo = n*rho + (n+1)*rho1 + An = const * radius**(2*n+1) / x0 ** (n+1.) * n * \ + (rho1-rho) / bunmo + Bn = const * 1. / x0 ** (n+1.) * (2*n+1) * (rho1) / bunmo + return An, Bn diff --git a/SimPEG/EM/Analytics/__init__.py b/SimPEG/EM/Analytics/__init__.py index 5b7a8851..9df2aef7 100644 --- a/SimPEG/EM/Analytics/__init__.py +++ b/SimPEG/EM/Analytics/__init__.py @@ -1,3 +1,4 @@ from TDEM import hzAnalyticDipoleT from FDEM import hzAnalyticDipoleF from FDEMcasing import * +from DC import DCAnalyticHalf, DCAnalyticSphere diff --git a/SimPEG/EM/Base.py b/SimPEG/EM/Base.py index a16cdb91..2a2ce363 100644 --- a/SimPEG/EM/Base.py +++ b/SimPEG/EM/Base.py @@ -1,6 +1,7 @@ from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solver as SimpegSolver from scipy.constants import mu_0 + class EMPropMap(Maps.PropMap): """ Property Map for EM Problems. The electrical conductivity (\\(\\sigma\\)) is the default inversion property, and the default value of the magnetic permeability is that of free space (\\(\\mu = 4\\pi\\times 10^{-7} \\) H/m) @@ -61,6 +62,15 @@ class BaseEMProblem(Problem.BaseProblem): self._Me = self.mesh.getEdgeInnerProduct() return self._Me + @property + def MeI(self): + """ + Edge inner product matrix + """ + if getattr(self, '_MeI', None) is None: + self._MeI = self.mesh.getEdgeInnerProduct(invMat=True) + return self._MeI + @property def Mf(self): """ @@ -70,6 +80,20 @@ class BaseEMProblem(Problem.BaseProblem): self._Mf = self.mesh.getFaceInnerProduct() return self._Mf + @property + def MfI(self): + """ + Face inner product matrix + """ + if getattr(self, '_MfI', None) is None: + self._MfI = self.mesh.getFaceInnerProduct(invMat=True) + return self._MfI + + @property + def Vol(self): + if getattr(self, '_Vol', None) is None: + self._Vol = Utils.sdiag(self.mesh.vol) + return self._Vol # ----- Magnetic Permeability ----- # @property @@ -127,7 +151,6 @@ class BaseEMProblem(Problem.BaseProblem): """ return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) * self.curModel.sigmaDeriv - @property def MeSigmaI(self): """ @@ -146,10 +169,7 @@ class BaseEMProblem(Problem.BaseProblem): dMeSigmaI_dI = -self.MeSigmaI**2 dMe_dsig = self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) - dsig_dm = self.curModel.sigmaDeriv - return dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm)) - # return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma, invMat=True)(u) - + return dMeSigmaI_dI * ( dMe_dsig * self.curModel.sigmaDeriv ) @property def MfRho(self): @@ -165,8 +185,7 @@ class BaseEMProblem(Problem.BaseProblem): """ Derivative of :code:`MfRho` with respect to the model. """ - return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv) - # self.curModel.rhoDeriv + return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * self.curModel.rhoDeriv @property def MfRhoI(self): @@ -183,7 +202,10 @@ class BaseEMProblem(Problem.BaseProblem): """ Derivative of :code:`MfRhoI` with respect to the model. """ - return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv + + dMfRhoI_dI = -self.MfRhoI**2 + dMf_drho = self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) + return dMfRhoI_dI * ( dMf_drho * self.curModel.rhoDeriv ) class BaseEMSurvey(Survey.BaseSurvey): @@ -192,7 +214,7 @@ class BaseEMSurvey(Survey.BaseSurvey): self.srcList = srcList Survey.BaseSurvey.__init__(self, **kwargs) - def eval(self, u): + def eval(self, f): """ Project fields to receiver locations :param Fields u: fields object @@ -202,8 +224,8 @@ class BaseEMSurvey(Survey.BaseSurvey): data = Survey.Data(self) for src in self.srcList: for rx in src.rxList: - data[src, rx] = rx.eval(src, self.mesh, u) + data[src, rx] = rx.eval(src, self.mesh, f) return data - def evalDeriv(self, u): + def evalDeriv(self, f): raise Exception('Use Receivers to project fields deriv.') diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index e2193973..253b9984 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -160,9 +160,9 @@ class Fields(SimPEG.Problem.Fields): return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint) return np.array(self._jDeriv_u(src, du_dm_v, adjoint) + self._jDeriv_m(src, v, adjoint), dtype = complex) -class Fields_e(Fields): +class Fields3D_e(Fields): """ - Fields object for Problem_e. + Fields object for Problem3D_e. :param Mesh mesh: mesh :param Survey survey: survey @@ -181,7 +181,7 @@ class Fields_e(Fields): } def __init__(self, mesh, survey, **kwargs): - Fields.__init__(self,mesh,survey,**kwargs) + Fields.__init__(self, mesh, survey, **kwargs) def startup(self): self.prob = self.survey.prob @@ -426,9 +426,9 @@ class Fields_e(Fields): -class Fields_b(Fields): +class Fields3D_b(Fields): """ - Fields object for Problem_b. + Fields object for Problem3D_b. :param Mesh mesh: mesh :param Survey survey: survey @@ -693,9 +693,9 @@ class Fields_b(Fields): return Zero() -class Fields_j(Fields): +class Fields3D_j(Fields): """ - Fields object for Problem_j. + Fields object for Problem3D_j. :param Mesh mesh: mesh :param Survey survey: survey @@ -988,9 +988,9 @@ class Fields_j(Fields): return 1./(1j * omega(src.freq)) * VI * (self._aveE2CCV * ( s_mDeriv(v) - self._edgeCurl.T * ( self._MfRhoDeriv(jSolution) * v ) ) ) -class Fields_h(Fields): +class Fields3D_h(Fields): """ - Fields object for Problem_h. + Fields object for Problem3D_h. :param Mesh mesh: mesh :param Survey survey: survey diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/ProblemFDEM.py similarity index 95% rename from SimPEG/EM/FDEM/FDEM.py rename to SimPEG/EM/FDEM/ProblemFDEM.py index caca7602..9eb2e915 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/ProblemFDEM.py @@ -1,7 +1,7 @@ from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver from scipy.constants import mu_0 from SurveyFDEM import Survey as SurveyFDEM -from FieldsFDEM import Fields, Fields_e, Fields_b, Fields_h, Fields_j +from FieldsFDEM import Fields, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j from SimPEG.EM.Base import BaseEMProblem from SimPEG.EM.Utils import omega @@ -17,8 +17,8 @@ class BaseFDEMProblem(BaseEMProblem): \mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\ {\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}} - if using the E-B formulation (:code:`Problem_e` - or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity. + if using the E-B formulation (:code:`Problem3D_e` + or :code:`Problem3D_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity. If we write Maxwell's equations in terms of \\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\) @@ -28,7 +28,7 @@ class BaseFDEMProblem(BaseEMProblem): \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{s_m} \\\\ \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e} - if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity. + if using the H-J formulation (:code:`Problem3D_j` or :code:`Problem3D_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity. The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\) """ @@ -87,7 +87,7 @@ class BaseFDEMProblem(BaseEMProblem): du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v ) for rx in src.rxList: - df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_dmFun = getattr(f, '_{0}Deriv'.format(rx.projField), None) df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False) Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v) Ainv.clean() @@ -125,7 +125,7 @@ class BaseFDEMProblem(BaseEMProblem): for rx in src.rxList: PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m - df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_duTFun = getattr(f, '_{0}Deriv'.format(rx.projField), None) df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True) ATinvdf_duT = ATinv * df_duT @@ -137,10 +137,9 @@ class BaseFDEMProblem(BaseEMProblem): df_dmT = df_dmT + du_dmT # TODO: this should be taken care of by the reciever? - real_or_imag = rx.projComp - if real_or_imag is 'real': + if rx.component is 'real': Jtv += np.array(df_dmT, dtype=complex).real - elif real_or_imag is 'imag': + elif rx.component is 'imag': Jtv += - np.array(df_dmT, dtype=complex).real else: raise Exception('Must be real or imag') @@ -167,6 +166,7 @@ class BaseFDEMProblem(BaseEMProblem): for i, src in enumerate(Srcs): smi, sei = src.eval(self) + #Why are you adding? s_m[:,i] = s_m[:,i] + smi s_e[:,i] = s_e[:,i] + sei @@ -177,7 +177,7 @@ class BaseFDEMProblem(BaseEMProblem): ################################ E-B Formulation ######################################### ########################################################################################## -class Problem_e(BaseFDEMProblem): +class Problem3D_e(BaseFDEMProblem): """ By eliminating the magnetic flux density using @@ -199,7 +199,7 @@ class Problem_e(BaseFDEMProblem): _solutionType = 'eSolution' _formulation = 'EB' - fieldsPair = Fields_e + fieldsPair = Fields3D_e def __init__(self, mesh, **kwargs): BaseFDEMProblem.__init__(self, mesh, **kwargs) @@ -288,7 +288,7 @@ class Problem_e(BaseFDEMProblem): return C.T * (MfMui * s_mDeriv(v)) -1j * omega(freq) * s_eDeriv(v) -class Problem_b(BaseFDEMProblem): +class Problem3D_b(BaseFDEMProblem): """ We eliminate :math:`\mathbf{e}` using @@ -310,7 +310,7 @@ class Problem_b(BaseFDEMProblem): _solutionType = 'bSolution' _formulation = 'EB' - fieldsPair = Fields_b + fieldsPair = Fields3D_b def __init__(self, mesh, **kwargs): BaseFDEMProblem.__init__(self, mesh, **kwargs) @@ -436,7 +436,7 @@ class Problem_b(BaseFDEMProblem): ########################################################################################## -class Problem_j(BaseFDEMProblem): +class Problem3D_j(BaseFDEMProblem): """ We eliminate \\\(\\\mathbf{h}\\\) using @@ -458,7 +458,7 @@ class Problem_j(BaseFDEMProblem): _solutionType = 'jSolution' _formulation = 'HJ' - fieldsPair = Fields_j + fieldsPair = Fields3D_j def __init__(self, mesh, **kwargs): BaseFDEMProblem.__init__(self, mesh, **kwargs) @@ -577,7 +577,7 @@ class Problem_j(BaseFDEMProblem): -class Problem_h(BaseFDEMProblem): +class Problem3D_h(BaseFDEMProblem): """ We eliminate \\\(\\\mathbf{j}\\\) using @@ -596,7 +596,7 @@ class Problem_h(BaseFDEMProblem): _solutionType = 'hSolution' _formulation = 'HJ' - fieldsPair = Fields_h + fieldsPair = Fields3D_h def __init__(self, mesh, **kwargs): BaseFDEMProblem.__init__(self, mesh, **kwargs) diff --git a/SimPEG/EM/FDEM/RxFDEM.py b/SimPEG/EM/FDEM/RxFDEM.py new file mode 100644 index 00000000..53d6c722 --- /dev/null +++ b/SimPEG/EM/FDEM/RxFDEM.py @@ -0,0 +1,126 @@ +import SimPEG +from SimPEG import sp + +class BaseRx(SimPEG.Survey.BaseRx): + """ + Frequency domain receiver base class + + :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`) + :param string orientation: receiver orientation 'x', 'y' or 'z' + :param string component: real or imaginary component 'real' or 'imag' + """ + + def __init__(self, locs, orientation=None, component=None): + assert(orientation in ['x','y','z']), "Orientation %s not known. Orientation must be in 'x', 'y', 'z'. Arbitrary orientations have not yet been implemented."%orientation + assert(component in ['real', 'imag']), "'component' must be 'real' or 'imag', not %s"%component + + self.projComp = orientation + self.component = component + + SimPEG.Survey.BaseRx.__init__(self, locs, rxType=None) #TODO: remove rxType from baseRx + + def projGLoc(self, u): + """Grid Location projection (e.g. Ex Fy ...)""" + return u._GLoc(self.projField) + self.projComp + + def eval(self, src, mesh, f): + """ + Project fields to recievers to get data. + + :param Source src: FDEM source + :param Mesh mesh: mesh used + :param Fields f: fields object + :rtype: numpy.ndarray + :return: fields projected to recievers + """ + + P = self.getP(mesh, self.projGLoc(f)) + f_part_complex = f[src, self.projField] + f_part = getattr(f_part_complex, self.component) # get the real or imag component + + return P*f_part + + def evalDeriv(self, src, mesh, f, v, adjoint=False): + """ + Derivative of projected fields with respect to the inversion model times a vector. + + :param Source src: FDEM source + :param Mesh mesh: mesh used + :param Fields f: fields object + :param numpy.ndarray v: vector to multiply + :rtype: numpy.ndarray + :return: fields projected to recievers + """ + + P = self.getP(mesh, self.projGLoc(f)) + + if not adjoint: + Pv_complex = P * v + Pv = getattr(Pv_complex, self.component) + elif adjoint: + Pv_real = P.T * v + + if self.component == 'imag': + Pv = 1j*Pv_real + elif self.component == 'real': + Pv = Pv_real.astype(complex) + else: + raise NotImplementedError('must be real or imag') + + return Pv + + +class Point_e(BaseRx): + """ + Electric field FDEM receiver + + :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`) + :param string orientation: receiver orientation 'x', 'y' or 'z' + :param string component: real or imaginary component 'real' or 'imag' + """ + + def __init__(self, locs, orientation=None, component=None): + self.projField = 'e' + super(Point_e, self).__init__(locs, orientation, component) + + +class Point_b(BaseRx): + """ + Magnetic flux FDEM receiver + + :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`) + :param string orientation: receiver orientation 'x', 'y' or 'z' + :param string component: real or imaginary component 'real' or 'imag' + """ + + def __init__(self, locs, orientation=None, component=None): + self.projField = 'b' + super(Point_b, self).__init__(locs, orientation, component) + + +class Point_h(BaseRx): + """ + Magnetic field FDEM receiver + + :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`) + :param string orientation: receiver orientation 'x', 'y' or 'z' + :param string component: real or imaginary component 'real' or 'imag' + """ + + def __init__(self, locs, orientation=None, component=None): + self.projField = 'h' + super(Point_h, self).__init__(locs, orientation, component) + + +class Point_j(BaseRx): + """ + Current density FDEM receiver + + :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`) + :param string orientation: receiver orientation 'x', 'y' or 'z' + :param string component: real or imaginary component 'real' or 'imag' + """ + + def __init__(self, locs, orientation=None, component=None): + self.projField = 'j' + super(Point_j, self).__init__(locs, orientation, component) diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py index 87967dd5..e6c8f005 100644 --- a/SimPEG/EM/FDEM/SrcFDEM.py +++ b/SimPEG/EM/FDEM/SrcFDEM.py @@ -9,12 +9,17 @@ class BaseSrc(Survey.BaseSrc): """ freq = None - # rxPair = RxFDEM - integrate = True + integrate = False + _ePrimary = None + _bPrimary = None + _hPrimary = None + _jPrimary = None + + def __init__(self, rxList, **kwargs): + Survey.BaseSrc.__init__(self, rxList, **kwargs) def eval(self, prob): """ - Evaluate the source terms. - :math:`s_m` : magnetic source term - :math:`s_e` : electric source term @@ -51,7 +56,9 @@ class BaseSrc(Survey.BaseSrc): :rtype: numpy.ndarray :return: primary magnetic flux density """ - return Zero() + if self._bPrimary is None: + return Zero() + return self._bPrimary def hPrimary(self, prob): """ @@ -61,7 +68,9 @@ class BaseSrc(Survey.BaseSrc): :rtype: numpy.ndarray :return: primary magnetic field """ - return Zero() + if self._hPrimary is None: + return Zero() + return self._hPrimary def ePrimary(self, prob): """ @@ -71,7 +80,9 @@ class BaseSrc(Survey.BaseSrc): :rtype: numpy.ndarray :return: primary electric field """ - return Zero() + if self._ePrimary is None: + return Zero() + return self._ePrimary def jPrimary(self, prob): """ @@ -81,7 +92,9 @@ class BaseSrc(Survey.BaseSrc): :rtype: numpy.ndarray :return: primary current density """ - return Zero() + if self._jPrimary is None: + return Zero() + return self._jPrimary def s_m(self, prob): """ @@ -136,15 +149,14 @@ class RawVec_e(BaseSrc): :param list rxList: receiver list :param float freq: frequency :param numpy.array s_e: electric source term - :param bool integrate: Integrate the source term (multiply by Me) [True] + :param bool integrate: Integrate the source term (multiply by Me) [False] """ - def __init__(self, rxList, freq, s_e, integrate=True): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): + def __init__(self, rxList, freq, s_e, **kwargs): self._s_e = np.array(s_e, dtype=complex) self.freq = float(freq) - self.integrate = integrate - BaseSrc.__init__(self, rxList) + BaseSrc.__init__(self, rxList, **kwargs) def s_e(self, prob): """ @@ -166,15 +178,14 @@ class RawVec_m(BaseSrc): :param float freq: frequency :param rxList: receiver list :param numpy.array s_m: magnetic source term - :param bool integrate: Integrate the source term (multiply by Me) [True] + :param bool integrate: Integrate the source term (multiply by Me) [False] """ - def __init__(self, rxList, freq, s_m, integrate=True): #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()): + def __init__(self, rxList, freq, s_m, **kwargs): #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()): self._s_m = np.array(s_m, dtype=complex) self.freq = float(freq) - self.integrate = integrate - BaseSrc.__init__(self, rxList) + BaseSrc.__init__(self, rxList, **kwargs) def s_m(self, prob): """ @@ -197,14 +208,13 @@ class RawVec(BaseSrc): :param float freq: frequency :param numpy.array s_m: magnetic source term :param numpy.array s_e: electric source term - :param bool integrate: Integrate the source term (multiply by Me) [True] + :param bool integrate: Integrate the source term (multiply by Me) [False] """ - def __init__(self, rxList, freq, s_m, s_e, integrate=True): + def __init__(self, rxList, freq, s_m, s_e, **kwargs): self._s_m = np.array(s_m, dtype=complex) self._s_e = np.array(s_e, dtype=complex) self.freq = float(freq) - self.integrate = integrate - BaseSrc.__init__(self, rxList) + BaseSrc.__init__(self, rxList, **kwargs) def s_m(self, prob): """ @@ -278,14 +288,13 @@ class MagDipole(BaseSrc): :param float mu: background magnetic permeability """ - def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu=mu_0): + def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu=mu_0, **kwargs): self.freq = float(freq) self.loc = loc self.orientation = orientation assert orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..." self.moment = moment self.mu = mu - self.integrate = False BaseSrc.__init__(self, rxList) def bPrimary(self, prob): @@ -543,7 +552,7 @@ class CircularLoop(BaseSrc): if not prob.mesh.isSymmetric: # TODO ? raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') - a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu) + a = MagneticLoopVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu) else: srcfct = MagneticDipoleVectorPotential diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index 1552a12c..46ae2523 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -4,126 +4,9 @@ from SimPEG.EM.Base import BaseEMSurvey from scipy.constants import mu_0 from SimPEG.Utils import Zero, Identity import SrcFDEM as Src +import RxFDEM as Rx from SimPEG import sp - -#################################################### -# Receivers -#################################################### - -class Rx(SimPEG.Survey.BaseRx): - """ - Frequency domain receivers - - :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`) - :param string rxType: reciever type from knownRxTypes - """ - - knownRxTypes = { - 'exr':['e', 'x', 'real'], - 'eyr':['e', 'y', 'real'], - 'ezr':['e', 'z', 'real'], - 'exi':['e', 'x', 'imag'], - 'eyi':['e', 'y', 'imag'], - 'ezi':['e', 'z', 'imag'], - - 'bxr':['b', 'x', 'real'], - 'byr':['b', 'y', 'real'], - 'bzr':['b', 'z', 'real'], - 'bxi':['b', 'x', 'imag'], - 'byi':['b', 'y', 'imag'], - 'bzi':['b', 'z', 'imag'], - - 'jxr':['j', 'x', 'real'], - 'jyr':['j', 'y', 'real'], - 'jzr':['j', 'z', 'real'], - 'jxi':['j', 'x', 'imag'], - 'jyi':['j', 'y', 'imag'], - 'jzi':['j', 'z', 'imag'], - - 'hxr':['h', 'x', 'real'], - 'hyr':['h', 'y', 'real'], - 'hzr':['h', 'z', 'real'], - 'hxi':['h', 'x', 'imag'], - 'hyi':['h', 'y', 'imag'], - 'hzi':['h', 'z', 'imag'], - } - radius = None - - def __init__(self, locs, rxType): - SimPEG.Survey.BaseRx.__init__(self, locs, rxType) - - @property - def projField(self): - """Field Type projection (e.g. e b ...)""" - return self.knownRxTypes[self.rxType][0] - - @property - def projComp(self): - """Component projection (real/imag)""" - return self.knownRxTypes[self.rxType][2] - - def projGLoc(self, u): - """Grid Location projection (e.g. Ex Fy ...)""" - return u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1] - - def eval(self, src, mesh, f): - """ - Project fields to recievers to get data. - - :param Source src: FDEM source - :param Mesh mesh: mesh used - :param Fields f: fields object - :rtype: numpy.ndarray - :return: fields projected to recievers - """ - # projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0]) - # projGLoc += self.knownRxTypes[self.rxType][1] - - P = self.getP(mesh, self.projGLoc(f)) - f_part_complex = f[src, self.projField] - # get the real or imag component - real_or_imag = self.projComp - f_part = getattr(f_part_complex, real_or_imag) - - return P*f_part - - def evalDeriv(self, src, mesh, f, v, adjoint=False): - """ - Derivative of projected fields with respect to the inversion model times a vector. - - :param Source src: FDEM source - :param Mesh mesh: mesh used - :param Fields f: fields object - :param numpy.ndarray v: vector to multiply - :rtype: numpy.ndarray - :return: fields projected to recievers - """ - - P = self.getP(mesh, self.projGLoc(f)) - - if not adjoint: - Pv_complex = P * v - real_or_imag = self.projComp - Pv = getattr(Pv_complex, real_or_imag) - elif adjoint: - Pv_real = P.T * v - - real_or_imag = self.projComp - if real_or_imag == 'imag': - Pv = 1j*Pv_real - elif real_or_imag == 'real': - Pv = Pv_real.astype(complex) - else: - raise NotImplementedError('must be real or imag') - - return Pv - - -#################################################### -# Survey -#################################################### - class Survey(BaseEMSurvey): """ Frequency domain electromagnetic survey @@ -132,7 +15,7 @@ class Survey(BaseEMSurvey): """ srcPair = Src.BaseSrc - rxPair = Rx + rxPair = Rx.BaseRx def __init__(self, srcList, **kwargs): # Sort these by frequency diff --git a/SimPEG/EM/FDEM/__init__.py b/SimPEG/EM/FDEM/__init__.py index 978972f5..c4ff6451 100644 --- a/SimPEG/EM/FDEM/__init__.py +++ b/SimPEG/EM/FDEM/__init__.py @@ -1,3 +1,5 @@ -from SurveyFDEM import Rx, Src, Survey -from FDEM import BaseFDEMProblem, Problem_e, Problem_b, Problem_j, Problem_h -from FieldsFDEM import * \ No newline at end of file +from SurveyFDEM import Survey +import SrcFDEM as Src +import RxFDEM as Rx +from ProblemFDEM import Problem3D_e, Problem3D_b, Problem3D_j, Problem3D_h +from FieldsFDEM import Fields3D_e, Fields3D_b, Fields3D_j, Fields3D_h diff --git a/SimPEG/EM/Static/DC/BoundaryUtils.py b/SimPEG/EM/Static/DC/BoundaryUtils.py new file mode 100644 index 00000000..3967eb46 --- /dev/null +++ b/SimPEG/EM/Static/DC/BoundaryUtils.py @@ -0,0 +1,160 @@ +import numpy as np + +def getxBCyBC_CC(mesh, alpha, beta, gamma): +# def getxBCyBC(mesh, alpha, beta, gamma): + """ + This is a subfunction generating mixed-boundary condition: + + .. math:: + + \nabla \cdot \vec{j} = -\nabla \cdot \vec{j}_s = q + + \rho \vec{j} = -\nabla \phi \phi + + \alpha \phi + \beta \frac{\partial \phi}{\partial r} = \gamma \ at \ r = \partial \Omega + + xBC = f_1(\alpha, \beta, \gamma) + yBC = f(\alpha, \beta, \gamma) + + Computes xBC and yBC for cell-centered discretizations + """ + if mesh.dim == 1: #1D + if (len(alpha) != 2 or len(beta) != 2 or len(gamma) != 2): + raise Exception("Lenght of list, alpha should be 2") + fCCxm,fCCxp = mesh.cellBoundaryInd + nBC = fCCxm.sum()+fCCxp.sum() + h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + + alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] + alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] + + # h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + h_xm, h_xp = mesh.hx[0], mesh.hx[-1] + + a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) + b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) + a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp) + b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp) + + xBC_xm = 0.5*a_xm + xBC_xp = 0.5*a_xp/b_xp + yBC_xm = 0.5*(1.-b_xm) + yBC_xp = 0.5*(1.-1./b_xp) + + xBC = np.r_[xBC_xm, xBC_xp] + yBC = np.r_[yBC_xm, yBC_xp] + + elif mesh.dim == 2: #2D + if (len(alpha) != 4 or len(beta) != 4 or len(gamma) != 4): + raise Exception("Lenght of list, alpha should be 4") + + fxm,fxp,fym,fyp = mesh.faceBoundaryInd + nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum() + + alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] + alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] + alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2] + alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3] + + # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] + # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + + h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp) + h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp) + + a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) + b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) + a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp) + b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp) + + a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym) + b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym) + a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp) + b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp) + + xBC_xm = 0.5*a_xm + xBC_xp = 0.5*a_xp/b_xp + yBC_xm = 0.5*(1.-b_xm) + yBC_xp = 0.5*(1.-1./b_xp) + xBC_ym = 0.5*a_ym + xBC_yp = 0.5*a_yp/b_yp + yBC_ym = 0.5*(1.-b_ym) + yBC_yp = 0.5*(1.-1./b_yp) + + sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]]) + sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]]) + + xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx] + xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy] + yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx] + yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy] + + xBC = np.r_[xBC_x, xBC_y] + yBC = np.r_[yBC_x, yBC_y] + + elif mesh.dim == 3: #3D + if (len(alpha) != 6 or len(beta) != 6 or len(gamma) != 6): + raise Exception("Lenght of list, alpha should be 6") + # fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd + fxm,fxp,fym,fyp,fzm,fzp = mesh.faceBoundaryInd + nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum() + + alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] + alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] + alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2] + alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3] + alpha_zm, beta_zm, gamma_zm = alpha[4], beta[4], gamma[4] + alpha_zp, beta_zp, gamma_zp = alpha[5], beta[5], gamma[5] + + # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] + # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + # h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2] + + h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp) + h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp) + h_zm, h_zp = mesh.hz[0]*np.ones_like(alpha_zm), mesh.hz[-1]*np.ones_like(alpha_zp) + + a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) + b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) + a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp) + b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp) + + a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym) + b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym) + a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp) + b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp) + + a_zm = gamma_zm/(0.5*alpha_zm-beta_zm/h_zm) + b_zm = (0.5*alpha_zm+beta_zm/h_zm)/(0.5*alpha_zm-beta_zm/h_zm) + a_zp = gamma_zp/(0.5*alpha_zp-beta_zp/h_zp) + b_zp = (0.5*alpha_zp+beta_zp/h_zp)/(0.5*alpha_zp-beta_zp/h_zp) + + xBC_xm = 0.5*a_xm + xBC_xp = 0.5*a_xp/b_xp + yBC_xm = 0.5*(1.-b_xm) + yBC_xp = 0.5*(1.-1./b_xp) + xBC_ym = 0.5*a_ym + xBC_yp = 0.5*a_yp/b_yp + yBC_ym = 0.5*(1.-b_ym) + yBC_yp = 0.5*(1.-1./b_yp) + xBC_zm = 0.5*a_zm + xBC_zp = 0.5*a_zp/b_zp + yBC_zm = 0.5*(1.-b_zm) + yBC_zp = 0.5*(1.-1./b_zp) + + sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]]) + sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]]) + sortindsfz = np.argsort(np.r_[np.arange(mesh.nFz)[fzm], np.arange(mesh.nFz)[fzp]]) + + xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx] + xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy] + xBC_z = np.r_[xBC_zm, xBC_zp][sortindsfz] + + yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx] + yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy] + yBC_z = np.r_[yBC_zm, yBC_zp][sortindsfz] + + xBC = np.r_[xBC_x, xBC_y, xBC_z] + yBC = np.r_[yBC_x, yBC_y, yBC_z] + + return xBC, yBC diff --git a/SimPEG/EM/Static/DC/FieldsDC.py b/SimPEG/EM/Static/DC/FieldsDC.py new file mode 100644 index 00000000..6fcea083 --- /dev/null +++ b/SimPEG/EM/Static/DC/FieldsDC.py @@ -0,0 +1,148 @@ +import SimPEG +from SimPEG.Utils import Identity, Zero +import numpy as np +from scipy.constants import epsilon_0 + +class Fields(SimPEG.Problem.Fields): + knownFields = {} + dtype = float + + def _phiDeriv(self, src, du_dm_v, v, adjoint=False): + if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None: + raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0]) + + if adjoint: + return self._phiDeriv_u(src, v, adjoint=adjoint), self._phiDeriv_m(src, v, adjoint=adjoint) + + return np.array(self._phiDeriv_u(src, du_dm_v, adjoint) + self._phiDeriv_m(src, v, adjoint), dtype = float) + + def _eDeriv(self, src, du_dm_v, v, adjoint=False): + if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None: + raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0]) + + if adjoint: + return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint) + return np.array(self._eDeriv_u(src, du_dm_v, adjoint) + self._eDeriv_m(src, v, adjoint), dtype = float) + + def _jDeriv(self, src, du_dm_v, v, adjoint=False): + if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None: + raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0]) + + if adjoint: + return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint) + return np.array(self._jDeriv_u(src, du_dm_v, adjoint) + self._jDeriv_m(src, v, adjoint), dtype = float) + + +class Fields_CC(Fields): + knownFields = {'phiSolution':'CC'} + aliasFields = { + 'phi': ['phiSolution','CC','_phi'], + 'j' : ['phiSolution','F','_j'], + 'e' : ['phiSolution','F','_e'], + 'charge' : ['phiSolution','CC','_charge'], + } + # primary - secondary + # CC variables + + def __init__(self, mesh, survey, **kwargs): + Fields.__init__(self, mesh, survey, **kwargs) + mesh.setCellGradBC("neumann") + cellGrad = mesh.cellGrad + def startup(self): + self.prob = self.survey.prob + + def _GLoc(self, fieldType): + if fieldType == 'phi': + return 'CC' + elif fieldType == 'e' or fieldType == 'j': + return 'F' + else: + raise Exception('Field type must be phi, e, j') + + def _phi(self, phiSolution, srcList): + return phiSolution + + def _phiDeriv_u(self, src, v, adjoint = False): + return Identity()*v + + def _phiDeriv_m(self, src, v, adjoint = False): + return Zero() + + def _j(self, phiSolution, srcList): + """ + .. math:: + \mathbf{j} = \mathbf{M}^{f \ -1}_{\rho} \mathbf{G} \phi + """ + return self.prob.MfRhoI*self.prob.Grad*phiSolution + + def _e(self, phiSolution, srcList): + """ + In HJ formulation e is not well-defined!! + .. math:: + \vec{e} = -\nabla \phi + """ + return -self.mesh.cellGrad*phiSolution + + def _charge(self, phiSolution, srcList): + """ + .. math:: + \int \nabla \codt \vec{e} = \int \frac{\rho_v }{\epsillon_0} + """ + return epsilon_0*self.prob.Vol*(self.mesh.faceDiv*self._e(phiSolution, srcList)) + +class Fields_N(Fields): + knownFields = {'phiSolution':'N'} + aliasFields = { + 'phi': ['phiSolution','N','_phi'], + 'j' : ['phiSolution','E','_j'], + 'e' : ['phiSolution','E','_e'], + 'charge' : ['phiSolution','N','_charge'], + } + # primary - secondary + # N variables + + def __init__(self, mesh, survey, **kwargs): + Fields.__init__(self, mesh, survey, **kwargs) + + def startup(self): + self.prob = self.survey.prob + + def _GLoc(self, fieldType): + if fieldType == 'phi': + return 'N' + elif fieldType == 'e' or fieldType == 'j': + return 'E' + else: + raise Exception('Field type must be phi, e, j') + + def _phi(self, phiSolution, srcList): + return phiSolution + + def _phiDeriv_u(self, src, v, adjoint = False): + return Identity()*v + + def _phiDeriv_m(self, src, v, adjoint = False): + return Zero() + + def _j(self, phiSolution, srcList): + """ + In EB formulation j is not well-defined!! + .. math:: + \mathbf{j} = - \mathbf{M}^{e}_{\sigma} \mathbf{G} \phi + """ + return self.prob.MeSigma * self._e(phiSolution, srcList) + + def _e(self, phiSolution, srcList): + """ + In HJ formulation e is not well-defined!! + .. math:: + \vec{e} = -\nabla \phi + """ + return -self.mesh.nodalGrad * phiSolution + + def _charge(self, phiSolution, srcList): + """ + .. math:: + \int \nabla \codt \vec{e} = \int \frac{\rho_v }{\epsillon_0} + """ + return - epsilon_0*(self.mesh.nodalGrad.T*self.mesh.getEdgeInnerProduct()*self._e(phiSolution, srcList)) diff --git a/SimPEG/EM/Static/DC/FieldsDC_2D.py b/SimPEG/EM/Static/DC/FieldsDC_2D.py new file mode 100644 index 00000000..da1fbf97 --- /dev/null +++ b/SimPEG/EM/Static/DC/FieldsDC_2D.py @@ -0,0 +1,146 @@ +import SimPEG +from SimPEG.Utils import Identity, Zero +import numpy as np + +class Fields_ky(SimPEG.Problem.TimeFields): + + """ + + Fancy Field Storage for a 2.5D code. + + u[:,'phi', kyInd] = phi + print u[src0,'phi'] + + Only one field type is stored for + each problem, the rest are computed. The fields obejct acts like an array and is indexed by + .. code-block:: python + f = problem.fields(m) + e = f[srcList,'e'] + j = f[srcList,'j'] + + If accessing all sources for a given field, use the :code:`:` + .. code-block:: python + f = problem.fields(m) + phi = f[:,'phi'] + e = f[:,'e'] + b = f[:,'b'] + The array returned will be size (nE or nF, nSrcs :math:`\\times` nFrequencies) + """ + + knownFields = {} + dtype = float + + def _phiDeriv(self,kyInd, src, du_dm_v, v, adjoint=False): + if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None: + raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0]) + + if adjoint: + return self._phiDeriv_u(kyInd, src, v, adjoint=adjoint), self._phiDeriv_m(kyInd, src, v, adjoint=adjoint) + + return np.array(self._phiDeriv_u(kyInd, src, du_dm_v, adjoint) + self._phiDeriv_m(kyInd, src, v, adjoint), dtype = float) + + def _eDeriv(self,kyInd, src, du_dm_v, v, adjoint=False): + if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None: + raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0]) + + if adjoint: + return self._eDeriv_u(kyInd, src, v, adjoint), self._eDeriv_m(kyInd, src, v, adjoint) + return np.array(self._eDeriv_u(kyInd, src, du_dm_v, adjoint) + self._eDeriv_m(kyInd, src, v, adjoint), dtype = float) + + def _jDeriv(self,kyInd, src, du_dm_v, v, adjoint=False): + if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None: + raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0]) + + if adjoint: + return self._jDeriv_u(kyInd, src, v, adjoint), self._jDeriv_m(kyInd, src, v, adjoint) + return np.array(self._jDeriv_u(kyInd, src, du_dm_v, adjoint) + self._jDeriv_m(kyInd, src, v, adjoint), dtype = float) + + + # def _eDeriv(self, tInd, src, dun_dm_v, v, adjoint=False): + # if adjoint is True: + # return self._eDeriv_u(tInd, src, v, adjoint), self._eDeriv_m(tInd, src, v, adjoint) + # return self._eDeriv_u(tInd, src, dun_dm_v) + self._eDeriv_m(tInd, src, v) + + # def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False): + # if adjoint is True: + # return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint) + # return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v) + + +class Fields_ky_CC(Fields_ky): + knownFields = {'phiSolution':'CC'} + aliasFields = { + 'phi': ['phiSolution','CC','_phi'], + 'j' : ['phiSolution','F','_j'], + 'e' : ['phiSolution','F','_e'], + } + # primary - secondary + # CC variables + + def __init__(self, mesh, survey, **kwargs): + Fields_ky.__init__(self, mesh, survey, **kwargs) + + def startup(self): + self.prob = self.survey.prob + + def _GLoc(self, fieldType): + if fieldType == 'phi': + return 'CC' + elif fieldType == 'e' or fieldType == 'j': + return 'F' + else: + raise Exception('Field type must be phi, e, j') + + def _phi(self, phiSolution, src, kyInd): + return phiSolution + + def _phiDeriv_u(self, kyInd, src, v, adjoint = False): + return Identity()*v + + def _phiDeriv_m(self, kyInd, src, v, adjoint = False): + return Zero() + + def _j(self, phiSolution, srcList): + raise NotImplementedError + + def _e(self, phiSolution, srcList): + raise NotImplementedError + +class Fields_ky_N(Fields_ky): + knownFields = {'phiSolution':'N'} + aliasFields = { + 'phi': ['phiSolution','N','_phi'], + 'j' : ['phiSolution','E','_j'], + 'e' : ['phiSolution','E','_e'], + } + # primary - secondary + # CC variables + + def __init__(self, mesh, survey, **kwargs): + Fields_ky.__init__(self, mesh, survey, **kwargs) + + def startup(self): + self.prob = self.survey.prob + + def _GLoc(self, fieldType): + if fieldType == 'phi': + return 'N' + elif fieldType == 'e' or fieldType == 'j': + return 'E' + else: + raise Exception('Field type must be phi, e, j') + + def _phi(self, phiSolution, src, kyInd): + return phiSolution + + def _phiDeriv_u(self, kyInd, src, v, adjoint = False): + return Identity()*v + + def _phiDeriv_m(self, kyInd, src, v, adjoint = False): + return Zero() + + def _j(self, phiSolution, srcList): + raise NotImplementedError + + def _e(self, phiSolution, srcList): + raise NotImplementedError diff --git a/SimPEG/EM/Static/DC/ProblemDC.py b/SimPEG/EM/Static/DC/ProblemDC.py new file mode 100644 index 00000000..81653e83 --- /dev/null +++ b/SimPEG/EM/Static/DC/ProblemDC.py @@ -0,0 +1,296 @@ +from SimPEG import Problem, Utils +from SimPEG.EM.Base import BaseEMProblem +from SurveyDC import Survey +from FieldsDC import Fields, Fields_CC, Fields_N +from SimPEG.Utils import sdiag +import numpy as np +from SimPEG.Utils import Zero +from BoundaryUtils import getxBCyBC_CC + +class BaseDCProblem(BaseEMProblem): + + surveyPair = Survey + fieldsPair = Fields + Ainv = None + + def fields(self, m): + self.curModel = m + + if not self.Ainv == None: + self.Ainv.clean() + + f = self.fieldsPair(self.mesh, self.survey) + A = self.getA() + self.Ainv = self.Solver(A, **self.solverOpts) + RHS = self.getRHS() + u = self.Ainv * RHS + Srcs = self.survey.srcList + f[Srcs, self._solutionType] = u + return f + + def Jvec(self, m, v, f=None): + + if f is None: + f = self.fields(m) + + self.curModel = m + + Jv = self.dataPair(self.survey) #same size as the data + + A = self.getA() + + for src in self.survey.srcList: + u_src = f[src, self._solutionType] # solution vector + dA_dm_v = self.getADeriv(u_src, v) + dRHS_dm_v = self.getRHSDeriv(src, v) + du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v ) + + for rx in src.rxList: + df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False) + Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v) + return Utils.mkvc(Jv) + + def Jtvec(self, m, v, f=None): + if f is None: + f = self.fields(m) + + self.curModel = m + + # Ensure v is a data object. + if not isinstance(v, self.dataPair): + v = self.dataPair(self.survey, v) + + Jtv = np.zeros(m.size) + AT = self.getA() + + + for src in self.survey.srcList: + u_src = f[src, self._solutionType] + for rx in src.rxList: + PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m + df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True) + + ATinvdf_duT = self.Ainv * df_duT + + dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True) + dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True) + du_dmT = -dA_dmT + dRHS_dmT + Jtv += (df_dmT + du_dmT).astype(float) + + return Utils.mkvc(Jtv) + + def getSourceTerm(self): + """ + takes concept of source and turns it into a matrix + """ + """ + Evaluates the sources, and puts them in matrix form + + :rtype: (numpy.ndarray, numpy.ndarray) + :return: q (nC or nN, nSrc) + """ + + Srcs = self.survey.srcList + + if self._formulation is 'EB': + n = self.mesh.nN + # return NotImplementedError + + elif self._formulation is 'HJ': + n = self.mesh.nC + + q = np.zeros((n, len(Srcs))) + + for i, src in enumerate(Srcs): + q[:,i] = src.eval(self) + return q + +class Problem3D_CC(BaseDCProblem): + + _solutionType = 'phiSolution' + _formulation = 'HJ' # CC potentials means J is on faces + fieldsPair = Fields_CC + + def __init__(self, mesh, **kwargs): + BaseDCProblem.__init__(self, mesh, **kwargs) + self.setBC() + + def getA(self): + """ + + Make the A matrix for the cell centered DC resistivity problem + + A = D MfRhoI G + + """ + + D = self.Div + G = self.Grad + MfRhoI = self.MfRhoI + A = D * MfRhoI * G + + # I think we should deprecate this for DC problem. + # if self._makeASymmetric is True: + # return V.T * A + return A + + def getADeriv(self, u, v, adjoint= False): + + D = self.Div + G = self.Grad + MfRhoIDeriv = self.MfRhoIDeriv + + if adjoint: + return(MfRhoIDeriv( G * u ).T) * ( D.T * v) + + return D * (MfRhoIDeriv( G * u ) * v) + + def getRHS(self): + """ + RHS for the DC problem + + q + """ + + RHS = self.getSourceTerm() + + return RHS + + def getRHSDeriv(self, src, v, adjoint=False): + """ + Derivative of the right hand side with respect to the model + """ + # TODO: add qDeriv for RHS depending on m + # qDeriv = src.evalDeriv(self, adjoint=adjoint) + # return qDeriv + return Zero() + + def setBC(self): + if self.mesh.dim==3: + fxm,fxp,fym,fyp,fzm,fzp = self.mesh.faceBoundaryInd + gBFxm = self.mesh.gridFx[fxm,:] + gBFxp = self.mesh.gridFx[fxp,:] + gBFym = self.mesh.gridFy[fym,:] + gBFyp = self.mesh.gridFy[fyp,:] + gBFzm = self.mesh.gridFz[fzm,:] + gBFzp = self.mesh.gridFz[fzp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + temp_zm, temp_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2]) + + alpha_xm, alpha_xp = temp_xm*0., temp_xp*0. + alpha_ym, alpha_yp = temp_ym*0., temp_yp*0. + alpha_zm, alpha_zp = temp_zm*0., temp_zp*0. + + beta_xm, beta_xp = temp_xm, temp_xp + beta_ym, beta_yp = temp_ym, temp_yp + beta_zm, beta_zp = temp_zm, temp_zp + + gamma_xm, gamma_xp = temp_xm*0., temp_xp*0. + gamma_ym, gamma_yp = temp_ym*0., temp_yp*0. + gamma_zm, gamma_zp = temp_zm*0., temp_zp*0. + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp] + + elif self.mesh.dim==2: + + fxm,fxp,fym,fyp = self.mesh.faceBoundaryInd + gBFxm = self.mesh.gridFx[fxm,:] + gBFxp = self.mesh.gridFx[fxp,:] + gBFym = self.mesh.gridFy[fym,:] + gBFyp = self.mesh.gridFy[fyp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + + alpha_xm, alpha_xp = temp_xm*0., temp_xp*0. + alpha_ym, alpha_yp = temp_ym*0., temp_yp*0. + + beta_xm, beta_xp = temp_xm, temp_xp + beta_ym, beta_yp = temp_ym, temp_yp + + gamma_xm, gamma_xp = temp_xm*0., temp_xp*0. + gamma_ym, gamma_yp = temp_ym*0., temp_yp*0. + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp] + + x_BC, y_BC = getxBCyBC_CC(self.mesh, alpha, beta, gamma) + V = self.Vol + self.Div = V * self.mesh.faceDiv + P_BC, B = self.mesh.getBCProjWF_simple() + M = B*self.mesh.aveCC2F + self.Grad = self.Div.T - P_BC*Utils.sdiag(y_BC)*M + + +class Problem3D_N(BaseDCProblem): + + _solutionType = 'phiSolution' + _formulation = 'EB' # N potentials means B is on faces + fieldsPair = Fields_N + + def __init__(self, mesh, **kwargs): + BaseDCProblem.__init__(self, mesh, **kwargs) + + def getA(self): + """ + + Make the A matrix for the cell centered DC resistivity problem + + A = G.T MeSigma G + + """ + + MeSigma = self.MeSigma + Grad = self.mesh.nodalGrad + A = Grad.T * MeSigma * Grad + + # Handling Null space of A + A[0,0] = A[0,0] + 1. + + return A + + def getADeriv(self, u, v, adjoint=False): + """ + + Product of the derivative of our system matrix with respect to the model and a vector + + """ + MeSigma = self.MeSigma + Grad = self.mesh.nodalGrad + if not adjoint: + return Grad.T*(self.MeSigmaDeriv(Grad*u)*v) + elif adjoint: + return self.MeSigmaDeriv(Grad*u).T * (Grad*v) + + + def getRHS(self): + """ + RHS for the DC problem + + q + """ + + RHS = self.getSourceTerm() + return RHS + + def getRHSDeriv(self, src, v, adjoint=False): + """ + Derivative of the right hand side with respect to the model + """ + # TODO: add qDeriv for RHS depending on m + # qDeriv = src.evalDeriv(self, adjoint=adjoint) + # return qDeriv + return Zero() + + + + diff --git a/SimPEG/EM/Static/DC/ProblemDC_2D.py b/SimPEG/EM/Static/DC/ProblemDC_2D.py new file mode 100644 index 00000000..8e0561f9 --- /dev/null +++ b/SimPEG/EM/Static/DC/ProblemDC_2D.py @@ -0,0 +1,349 @@ +from SimPEG import Problem, Utils +from SimPEG.EM.Base import BaseEMProblem +from SurveyDC import Survey, Survey_ky +from FieldsDC_2D import Fields_ky, Fields_ky_CC, Fields_ky_N +from SimPEG.Utils import sdiag +import numpy as np +from SimPEG.Utils import Zero +from BoundaryUtils import getxBCyBC_CC + +class BaseDCProblem_2D(BaseEMProblem): + + surveyPair = Survey_ky + fieldsPair = Fields_ky + nky = 15 + kys = np.logspace(-4, 1, nky) + Ainv = [None for i in range(nky)] + nT = nky # Only for using TimeFields + + def fields(self, m): + self.curModel = m + + if not self.Ainv[0] == None: + for i in range(self.nky): + self.Ainv[i].clean() + + f = self.fieldsPair(self.mesh, self.survey) + Srcs = self.survey.srcList + for iky in range(self.nky): + ky = self.kys[iky] + A = self.getA(ky) + self.Ainv[iky] = self.Solver(A, **self.solverOpts) + RHS = self.getRHS(ky) + u = self.Ainv[iky] * RHS + f[Srcs, self._solutionType, iky] = u + return f + + def Jvec(self, m, v, f=None): + + if f is None: + f = self.fields(m) + + self.curModel = m + + Jv = self.dataPair(self.survey) #same size as the data + Jv0 = self.dataPair(self.survey) + + # Assume y=0. + # This needs some thoughts to implement in general when src is dipole + dky = np.diff(self.kys) + dky = np.r_[dky[0], dky] + y = 0. + + #TODO: this loop is pretty slow .. (Parellize) + for iky in range(self.nky): + ky = self.kys[iky] + A = self.getA(ky) + for src in self.survey.srcList: + u_src = f[src, self._solutionType, iky] # solution vector + dA_dm_v = self.getADeriv(ky, u_src, v) + dRHS_dm_v = self.getRHSDeriv(ky, src, v) + du_dm_v = self.Ainv[iky] * ( - dA_dm_v + dRHS_dm_v ) + for rx in src.rxList: + df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_dm_v = df_dmFun(iky, src, du_dm_v, v, adjoint=False) + # Trapezoidal intergration + Jv1_temp = 1./np.pi*rx.evalDeriv(ky, src, self.mesh, f, df_dm_v) + if iky==0: + #First assigment + Jv[src, rx] = Jv1_temp*dky[iky]*np.cos(ky*y) + else: + Jv[src, rx] += Jv1_temp*dky[iky] /2.*np.cos(ky*y) + Jv[src, rx] += Jv0[src, rx]*dky[iky]/2.*np.cos(ky*y) + Jv0[src, rx] = Jv1_temp.copy() + return Utils.mkvc(Jv) + + def Jtvec(self, m, v, f=None): + if f is None: + f = self.fields(m) + + self.curModel = m + + # Ensure v is a data object. + if not isinstance(v, self.dataPair): + v = self.dataPair(self.survey, v) + + Jtv = np.zeros(m.size, dtype=float) + + # Assume y=0. + # This needs some thoughts to implement in general when src is dipole + dky = np.diff(self.kys) + dky = np.r_[dky[0], dky] + y = 0. + + for src in self.survey.srcList: + for rx in src.rxList: + Jtv_temp1 = np.zeros(m.size, dtype=float) + Jtv_temp0 = np.zeros(m.size, dtype=float) + #TODO: this loop is pretty slow .. (Parellize) + for iky in range(self.nky): + u_src = f[src, self._solutionType, iky] + ky = self.kys[iky] + AT = self.getA(ky) + PTv = rx.evalDeriv(ky, src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m + df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_duT, df_dmT = df_duTFun(iky, src, None, PTv, adjoint=True) + + ATinvdf_duT = self.Ainv[iky] * df_duT + + dA_dmT = self.getADeriv(ky, u_src, ATinvdf_duT, adjoint=True) + dRHS_dmT = self.getRHSDeriv(ky, src, ATinvdf_duT, adjoint=True) + du_dmT = -dA_dmT + dRHS_dmT + Jtv_temp1 = 1./np.pi*(df_dmT + du_dmT).astype(float) + # Trapezoidal intergration + if iky==0: + #First assigment + Jtv += Jtv_temp1*dky[iky]*np.cos(ky*y) + else: + Jtv += Jtv_temp1*dky[iky]/2.*np.cos(ky*y) + Jtv += Jtv_temp0*dky[iky]/2.*np.cos(ky*y) + Jtv_temp0 = Jtv_temp1.copy() + return Utils.mkvc(Jtv) + + def getSourceTerm(self, ky): + """ + takes concept of source and turns it into a matrix + """ + """ + Evaluates the sources, and puts them in matrix form + + :rtype: (numpy.ndarray, numpy.ndarray) + :return: q (nC or nN, nSrc) + """ + + Srcs = self.survey.srcList + + if self._formulation is 'EB': + n = self.mesh.nN + # return NotImplementedError + + elif self._formulation is 'HJ': + n = self.mesh.nC + + q = np.zeros((n, len(Srcs))) + + for i, src in enumerate(Srcs): + q[:,i] = src.eval(self) + return q + +class Problem2D_CC(BaseDCProblem_2D): + + _solutionType = 'phiSolution' + _formulation = 'HJ' # CC potentials means J is on faces + fieldsPair = Fields_ky_CC + + def __init__(self, mesh, **kwargs): + BaseDCProblem_2D.__init__(self, mesh, **kwargs) + self.setBC() + + def getA(self, ky): + """ + + Make the A matrix for the cell centered DC resistivity problem + + A = D MfRhoI G + + """ + + D = self.Div + G = self.Grad + vol = self.mesh.vol + MfRhoI = self.MfRhoI + # Get resistivity rho + rho = self.curModel.rho + A = D * MfRhoI * G + Utils.sdiag(ky**2*vol/rho) + return A + + def getADeriv(self, ky, u, v, adjoint= False): + + D = self.Div + G = self.Grad + vol = self.mesh.vol + MfRhoIDeriv = self.MfRhoIDeriv + rho = self.curModel.rho + if adjoint: + return(MfRhoIDeriv( G * u ).T) * ( D.T * v) + ky**2*Utils.sdiag(u.flatten()*vol*(-1./rho**2))*v + return D * ((MfRhoIDeriv( G * u )) * v) + ky**2*Utils.sdiag(u.flatten()*vol*(-1./rho**2))*v + + def getRHS(self, ky): + """ + RHS for the DC problem + + q + """ + + RHS = self.getSourceTerm(ky) + return RHS + + def getRHSDeriv(self, ky, src, v, adjoint=False): + """ + Derivative of the right hand side with respect to the model + """ + # TODO: add qDeriv for RHS depending on m + # qDeriv = src.evalDeriv(self, ky, adjoint=adjoint) + # return qDeriv + return Zero() + + def setBC(self): + if self.mesh.dim==3: + fxm,fxp,fym,fyp,fzm,fzp = self.mesh.faceBoundaryInd + gBFxm = self.mesh.gridFx[fxm,:] + gBFxp = self.mesh.gridFx[fxp,:] + gBFym = self.mesh.gridFy[fym,:] + gBFyp = self.mesh.gridFy[fyp,:] + gBFzm = self.mesh.gridFz[fzm,:] + gBFzp = self.mesh.gridFz[fzp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + temp_zm, temp_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2]) + + alpha_xm, alpha_xp = temp_xm*0., temp_xp*0. + alpha_ym, alpha_yp = temp_ym*0., temp_yp*0. + alpha_zm, alpha_zp = temp_zm*0., temp_zp*0. + + beta_xm, beta_xp = temp_xm, temp_xp + beta_ym, beta_yp = temp_ym, temp_yp + beta_zm, beta_zp = temp_zm, temp_zp + + gamma_xm, gamma_xp = temp_xm*0., temp_xp*0. + gamma_ym, gamma_yp = temp_ym*0., temp_yp*0. + gamma_zm, gamma_zp = temp_zm*0., temp_zp*0. + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp] + + elif self.mesh.dim==2: + + fxm,fxp,fym,fyp = self.mesh.faceBoundaryInd + gBFxm = self.mesh.gridFx[fxm,:] + gBFxp = self.mesh.gridFx[fxp,:] + gBFym = self.mesh.gridFy[fym,:] + gBFyp = self.mesh.gridFy[fyp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + + alpha_xm, alpha_xp = temp_xm*0., temp_xp*0. + alpha_ym, alpha_yp = temp_ym*0., temp_yp*0. + + beta_xm, beta_xp = temp_xm, temp_xp + beta_ym, beta_yp = temp_ym, temp_yp + + gamma_xm, gamma_xp = temp_xm*0., temp_xp*0. + gamma_ym, gamma_yp = temp_ym*0., temp_yp*0. + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp] + + x_BC, y_BC = getxBCyBC_CC(self.mesh, alpha, beta, gamma) + V = self.Vol + self.Div = V * self.mesh.faceDiv + P_BC, B = self.mesh.getBCProjWF_simple() + M = B*self.mesh.aveCC2F + self.Grad = self.Div.T - P_BC*Utils.sdiag(y_BC)*M + +class Problem2D_N(BaseDCProblem_2D): + + _solutionType = 'phiSolution' + _formulation = 'EB' # CC potentials means J is on faces + fieldsPair = Fields_ky_N + + def __init__(self, mesh, **kwargs): + BaseDCProblem_2D.__init__(self, mesh, **kwargs) + # self.setBC() + + @property + def MnSigma(self): + """ + Node inner product matrix for \\(\\sigma\\). Used in the E-B formulation + """ + # TODO: only works isotropic sigma + sigma = self.curModel.sigma + vol = self.mesh.vol + MnSigma = Utils.sdiag(self.mesh.aveN2CC.T*(Utils.sdiag(vol)*sigma)) + + return MnSigma + + def MnSigmaDeriv(self, u): + """ + Derivative of MnSigma with respect to the model + """ + sigma = self.curModel.sigma + sigmaderiv = self.curModel.sigmaDeriv + vol = self.mesh.vol + return Utils.sdiag(u)*self.mesh.aveN2CC.T*Utils.sdiag(vol) * self.curModel.sigmaDeriv + + def getA(self, ky): + """ + + Make the A matrix for the cell centered DC resistivity problem + + A = D MfRhoI G + + """ + + MeSigma = self.MeSigma + MnSigma = self.MnSigma + Grad = self.mesh.nodalGrad + # Get conductivity sigma + sigma = self.curModel.sigma + A = Grad.T * MeSigma * Grad + ky**2*MnSigma + + # Handling Null space of A + A[0,0] = A[0,0] + 1. + return A + + def getADeriv(self, ky, u, v, adjoint= False): + + MeSigma = self.MeSigma + Grad = self.mesh.nodalGrad + sigma = self.curModel.sigma + vol = self.mesh.vol + + if adjoint: + return self.MeSigmaDeriv(Grad*u).T * (Grad*v) + ky**2*self.MnSigmaDeriv(u).T*v + return Grad.T*(self.MeSigmaDeriv(Grad*u)*v) + ky**2*self.MnSigmaDeriv(u)*v + + def getRHS(self, ky): + """ + RHS for the DC problem + + q + """ + + RHS = self.getSourceTerm(ky) + return RHS + + def getRHSDeriv(self, ky, src, v, adjoint=False): + """ + Derivative of the right hand side with respect to the model + """ + # TODO: add qDeriv for RHS depending on m + # qDeriv = src.evalDeriv(self, ky, adjoint=adjoint) + # return qDeriv + return Zero() diff --git a/SimPEG/EM/Static/DC/RxDC.py b/SimPEG/EM/Static/DC/RxDC.py new file mode 100644 index 00000000..d2ec6098 --- /dev/null +++ b/SimPEG/EM/Static/DC/RxDC.py @@ -0,0 +1,129 @@ +import SimPEG +import numpy as np +from SimPEG.Utils import Zero, closestPoints + +class BaseRx(SimPEG.Survey.BaseRx): + locs = None + rxType = None + + knownRxTypes = { + 'phi':['phi',None], + 'ex':['e','x'], + 'ey':['e','y'], + 'ez':['e','z'], + 'jx':['j','x'], + 'jy':['j','y'], + 'jz':['j','z'], + } + + def __init__(self, locs, rxType, **kwargs): + SimPEG.Survey.BaseRx.__init__(self, locs, rxType, **kwargs) + + + @property + def projField(self): + """Field Type projection (e.g. e b ...)""" + return self.knownRxTypes[self.rxType][0] + + def projGLoc(self, f): + """Grid Location projection (e.g. Ex Fy ...)""" + comp = self.knownRxTypes[self.rxType][1] + if comp is not None: + return f._GLoc(self.rxType) + comp + return f._GLoc(self.rxType) + + def eval(self, src, mesh, f): + P = self.getP(mesh, self.projGLoc(f)) + return P*f[src, self.projField] + + def evalDeriv(self, src, mesh, f, v, adjoint=False): + P = self.getP(mesh, self.projGLoc(f)) + if not adjoint: + return P*v + elif adjoint: + return P.T*v + +# DC.Rx.Dipole(locs) +class Dipole(BaseRx): + + def __init__(self, locsM, locsN, rxType = 'phi', **kwargs): + assert locsM.shape == locsN.shape, 'locsM and locsN need to be the same size' + locs = [locsM, locsN] + # We may not need this ... + BaseRx.__init__(self, locs, rxType) + + @property + def nD(self): + """Number of data in the receiver.""" + return self.locs[0].shape[0] + + # Not sure why ... + # return int(self.locs[0].size / 2) + + + def getP(self, mesh, Gloc): + if mesh in self._Ps: + return self._Ps[mesh] + + P0 = mesh.getInterpolationMat(self.locs[0], Gloc) + P1 = mesh.getInterpolationMat(self.locs[1], Gloc) + P = P0 - P1 + + if self.storeProjections: + self._Ps[mesh] = P + + return P + + +class Dipole_ky(BaseRx): + + def __init__(self, locsM, locsN, rxType = 'phi', **kwargs): + assert locsM.shape == locsN.shape, 'locsM and locsN need to be the same size' + locs = [locsM, locsN] + # We may not need this ... + BaseRx.__init__(self, locs, rxType) + + @property + def nD(self): + """Number of data in the receiver.""" + return self.locs[0].shape[0] + + # Not sure why ... + # return int(self.locs[0].size / 2) + + def getP(self, mesh, Gloc): + if mesh in self._Ps: + return self._Ps[mesh] + + P0 = mesh.getInterpolationMat(self.locs[0], Gloc) + P1 = mesh.getInterpolationMat(self.locs[1], Gloc) + P = P0 - P1 + if self.storeProjections: + self._Ps[mesh] = P + return P + + def eval(self, kys, src, mesh, f): + P = self.getP(mesh, self.projGLoc(f)) + Pf = P*f[src, self.projField,:] + return self.IntTrapezoidal(kys, Pf, y=0.) + + def evalDeriv(self, ky, src, mesh, f, v, adjoint=False): + P = self.getP(mesh, self.projGLoc(f)) + if not adjoint: + return P*v + elif adjoint: + return P.T*v + + def IntTrapezoidal(self, kys, Pf, y=0.): + phi = np.zeros(Pf.shape[0]) + nky = kys.size + dky = np.diff(kys) + dky = np.r_[dky[0], dky] + phi0 = 1./np.pi*Pf[:,0] + for iky in range(nky): + phi1 = 1./np.pi*Pf[:,iky] + phi += phi1*dky[iky]/2.*np.cos(kys[iky]*y) + phi += phi0*dky[iky]/2.*np.cos(kys[iky]*y) + phi0 = phi1.copy() + return phi + diff --git a/SimPEG/EM/Static/DC/SrcDC.py b/SimPEG/EM/Static/DC/SrcDC.py new file mode 100644 index 00000000..ae52f376 --- /dev/null +++ b/SimPEG/EM/Static/DC/SrcDC.py @@ -0,0 +1,86 @@ +import SimPEG +# from SimPEG.EM.Base import BaseEMSurvey +from SimPEG.Utils import Zero, closestPoints, mkvc +import numpy as np + +class BaseSrc(SimPEG.Survey.BaseSrc): + + current = 1.0 + loc = None + + def __init__(self, rxList, **kwargs): + SimPEG.Survey.BaseSrc.__init__(self, rxList, **kwargs) + + def eval(self, prob): + raise NotImplementedError + + def evalDeriv(self, prob): + return Zero() + + +class Dipole(BaseSrc): + + def __init__(self, rxList, locA, locB, **kwargs): + assert locA.shape == locB.shape, 'Shape of locA and locB should be the same' + self.loc = [locA, locB] + BaseSrc.__init__(self, rxList, **kwargs) + + def eval(self, prob): + if prob._formulation == 'HJ': + inds = closestPoints(prob.mesh, self.loc, gridLoc='CC') + q = np.zeros(prob.mesh.nC) + q[inds] = self.current * np.r_[1., -1.] + elif prob._formulation == 'EB': + qa = prob.mesh.getInterpolationMat(self.loc[0], locType='N').todense() + qb = -prob.mesh.getInterpolationMat(self.loc[1], locType='N').todense() + q = self.current * mkvc(qa+qb) + return q + +class Pole(BaseSrc): + + def __init__(self, rxList, loc, **kwargs): + BaseSrc.__init__(self, rxList, loc=loc, **kwargs) + + def eval(self, prob): + if prob._formulation == 'HJ': + inds = closestPoints(prob.mesh, self.loc) + q = np.zeros(prob.mesh.nC) + q[inds] = self.current * np.r_[1.] + elif prob._formulation == 'EB': + q = prob.mesh.getInterpolationMat(self.loc, locType='N').todense() + q = self.current * mkvc(q) + return q + + +# class Dipole_ky(BaseSrc): + +# def __init__(self, rxList, locA, locB, **kwargs): +# assert locA.shape == locB.shape, 'Shape of locA and locB should be the same' +# self.loc = [locA[[0,2]], locB[[0,2]]] +# BaseSrc.__init__(self, rxList, **kwargs) + +# def eval(self, prob): +# if prob._formulation == 'HJ': +# inds = closestPoints(prob.mesh, self.loc, gridLoc='CC') +# q = np.zeros(prob.mesh.nC) +# q[inds] = self.current * np.r_[1., -1.] +# elif prob._formulation == 'EB': +# qa = prob.mesh.getInterpolationMat(self.loc[0], locType='N').todense() +# qb = -prob.mesh.getInterpolationMat(self.loc[1], locType='N').todense() +# q = self.current * mkvc(qa+qb) +# return q + +# class Pole_ky(BaseSrc): + +# def __init__(self, rxList, loc, **kwargs): +# BaseSrc.__init__(self, rxList, loc=loc, **kwargs) + +# def eval(self, prob): +# if prob._formulation == 'HJ': +# inds = closestPoints(prob.mesh, self.loc[[0,2]]) +# q = np.zeros(prob.mesh.nC) +# q[inds] = self.current * np.r_[1.] +# elif prob._formulation == 'EB': +# q = prob.mesh.getInterpolationMat(self.loc[[0,2]], locType='N').todense() +# q = self.current * mkvc(q) +# return q diff --git a/SimPEG/EM/Static/DC/SurveyDC.py b/SimPEG/EM/Static/DC/SurveyDC.py new file mode 100644 index 00000000..d9c493a2 --- /dev/null +++ b/SimPEG/EM/Static/DC/SurveyDC.py @@ -0,0 +1,38 @@ +import SimPEG +from SimPEG.EM.Base import BaseEMSurvey +from SimPEG import sp, Survey +from SimPEG.Utils import Zero, Identity +from RxDC import BaseRx +from SrcDC import BaseSrc + +class Survey(BaseEMSurvey): + rxPair = BaseRx + srcPair = BaseSrc + + def __init__(self, srcList, **kwargs): + self.srcList = srcList + BaseEMSurvey.__init__(self, srcList, **kwargs) + +class Survey_ky(BaseEMSurvey): + rxPair = BaseRx + srcPair = BaseSrc + + def __init__(self, srcList, **kwargs): + self.srcList = srcList + BaseEMSurvey.__init__(self, srcList, **kwargs) + + def eval(self, f): + """ + Project fields to receiver locations + :param Fields u: fields object + :rtype: numpy.ndarray + :return: data + """ + data = SimPEG.Survey.Data(self) + kys = self.prob.kys + for src in self.srcList: + for rx in src.rxList: + data[src, rx] = rx.eval(kys, src, self.mesh, f) + return data + + diff --git a/SimPEG/EM/Static/DC/Utils.py b/SimPEG/EM/Static/DC/Utils.py new file mode 100644 index 00000000..626ae368 --- /dev/null +++ b/SimPEG/EM/Static/DC/Utils.py @@ -0,0 +1,38 @@ +import numpy as np + +def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False): + + import SimPEG.EM.Static.DC as DC + + elocs = np.arange(0,aSpacing*nElecs,aSpacing) + elocs -= (nElecs*aSpacing - aSpacing)/2 + space = 1 + WENNER = np.zeros((0,),dtype=int) + for ii in range(nElecs): + for jj in range(nElecs): + test = np.r_[jj,jj+space,jj+space*2,jj+space*3] + if np.any(test >= nElecs): + break + WENNER = np.r_[WENNER, test] + space += 1 + WENNER = WENNER.reshape((-1,4)) + + + if plotIt: + for i, s in enumerate('rbkg'): + plt.plot(elocs[WENNER[:,i]],s+'.') + plt.show() + + # Create sources and receivers + i = 0 + if in2D: + getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0] + else: + getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0, 0] + srcList = [] + for i in range(WENNER.shape[0]): + rx = DC.Rx.Dipole(getLoc(i,1).reshape([1,-1]),getLoc(i,2).reshape([1,-1])) + src = DC.Src.Dipole([rx], getLoc(i,0),getLoc(i,3)) + srcList += [src] + + return srcList diff --git a/SimPEG/EM/Static/DC/__init__.py b/SimPEG/EM/Static/DC/__init__.py new file mode 100644 index 00000000..1080e391 --- /dev/null +++ b/SimPEG/EM/Static/DC/__init__.py @@ -0,0 +1,8 @@ +from ProblemDC import Problem3D_CC, Problem3D_N +from ProblemDC_2D import Problem2D_CC, Problem2D_N +from SurveyDC import Survey, Survey_ky +import SrcDC as Src #Pole +import RxDC as Rx +from FieldsDC import Fields_CC +from BoundaryUtils import getxBCyBC_CC +import Utils diff --git a/SimPEG/EM/Static/IP/ProblemIP.py b/SimPEG/EM/Static/IP/ProblemIP.py new file mode 100644 index 00000000..77767452 --- /dev/null +++ b/SimPEG/EM/Static/IP/ProblemIP.py @@ -0,0 +1,372 @@ +from SimPEG import Problem, Utils, Maps, Mesh +from SimPEG.EM.Base import BaseEMProblem +from SimPEG.EM.Static.DC.FieldsDC import Fields, Fields_CC, Fields_N +from SimPEG.Utils import sdiag +import numpy as np +from SimPEG.Utils import Zero +from SimPEG.EM.Static.DC import getxBCyBC_CC +from SurveyIP import Survey + +class IPPropMap(Maps.PropMap): + """ + Property Map for IP Problems. The electrical chargeability, + (\\(\\eta\\)) is the default inversion property + """ + eta = Maps.Property("Electrical Chargeability", defaultInvProp = True) + +class BaseIPProblem(BaseEMProblem): + + surveyPair = Survey + fieldsPair = Fields + PropMap = IPPropMap + Ainv = None + sigma = None + rho = None + f = None + Ainv = None + + def fields(self, m): + self.curModel = m + if self.f is None: + self.f = self.fieldsPair(self.mesh, self.survey) + if self.Ainv == None: + A = self.getA() + self.Ainv = self.Solver(A, **self.solverOpts) + RHS = self.getRHS() + u = self.Ainv * RHS + Srcs = self.survey.srcList + self.f[Srcs, self._solutionType] = u + return self.f + + def Jvec(self, m, v, f=None): + + if f is None: + f = self.fields(m) + + self.curModel = m + + Jv = self.dataPair(self.survey) #same size as the data + + A = self.getA() + + for src in self.survey.srcList: + u_src = f[src, self._solutionType] # solution vector + dA_dm_v = self.getADeriv(u_src, v) + dRHS_dm_v = self.getRHSDeriv(src, v) + du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v ) + + for rx in src.rxList: + df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False) + Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v) + # Conductivity (d u / d log sigma) + if self._formulation is 'EB': + return -Utils.mkvc(Jv) + # Conductivity (d u / d log rho) + if self._formulation is 'HJ': + return Utils.mkvc(Jv) + + def Jtvec(self, m, v, f=None): + if f is None: + f = self.fields(m) + + self.curModel = m + + # Ensure v is a data object. + if not isinstance(v, self.dataPair): + v = self.dataPair(self.survey, v) + + Jtv = np.zeros(m.size) + AT = self.getA() + + for src in self.survey.srcList: + u_src = f[src, self._solutionType] + for rx in src.rxList: + PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m + df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True) + ATinvdf_duT = self.Ainv * df_duT + dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True) + dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True) + du_dmT = -dA_dmT + dRHS_dmT + Jtv += (df_dmT + du_dmT).astype(float) + # Conductivity ((d u / d log sigma).T) + if self._formulation is 'EB': + return -Utils.mkvc(Jtv) + # Conductivity ((d u / d log rho).T) + if self._formulation is 'HJ': + return Utils.mkvc(Jtv) + + def getSourceTerm(self): + """ + takes concept of source and turns it into a matrix + """ + """ + Evaluates the sources, and puts them in matrix form + + :rtype: (numpy.ndarray, numpy.ndarray) + :return: q (nC or nN, nSrc) + """ + + Srcs = self.survey.srcList + + if self._formulation is 'EB': + n = self.mesh.nN + # return NotImplementedError + + elif self._formulation is 'HJ': + n = self.mesh.nC + + q = np.zeros((n, len(Srcs))) + + for i, src in enumerate(Srcs): + q[:,i] = src.eval(self) + return q + + @property + def deleteTheseOnModelUpdate(self): + toDelete = [] + return toDelete + + # assume log rho or log cond + @property + def MeSigma(self): + """ + Edge inner product matrix for \\(\\sigma\\). Used in the E-B formulation + """ + if getattr(self, '_MeSigma', None) is None: + self._MeSigma = self.mesh.getEdgeInnerProduct(self.sigma) + return self._MeSigma + + @property + def MfRhoI(self): + """ + Inverse of :code:`MfRho` + """ + if getattr(self, '_MfRhoI', None) is None: + self._MfRhoI = self.mesh.getFaceInnerProduct(self.rho, invMat=True) + return self._MfRhoI + + def MfRhoIDeriv(self,u): + """ + Derivative of :code:`MfRhoI` with respect to the model. + """ + + dMfRhoI_dI = -self.MfRhoI**2 + dMf_drho = self.mesh.getFaceInnerProductDeriv(self.rho)(u) + drho_dlogrho = Utils.sdiag(self.rho)*self.curModel.etaDeriv + return dMfRhoI_dI * ( dMf_drho * ( drho_dlogrho)) + + # TODO: This should take a vector + def MeSigmaDeriv(self, u): + """ + Derivative of MeSigma with respect to the model + """ + dsigma_dlogsigma = Utils.sdiag(self.sigma)*self.curModel.etaDeriv + return self.mesh.getEdgeInnerProductDeriv(self.sigma)(u) * dsigma_dlogsigma + +class Problem3D_CC(BaseIPProblem): + + _solutionType = 'phiSolution' + _formulation = 'HJ' # CC potentials means J is on faces + fieldsPair = Fields_CC + + def __init__(self, mesh, **kwargs): + BaseIPProblem.__init__(self, mesh, **kwargs) + self.setBC() + + def getA(self): + """ + + Make the A matrix for the cell centered DC resistivity problem + + A = D MfRhoI G + + """ + + D = self.Div + G = self.Grad + MfRhoI = self.MfRhoI + A = D * MfRhoI * G + + # I think we should deprecate this for DC problem. + # if self._makeASymmetric is True: + # return V.T * A + return A + + def getADeriv(self, u, v, adjoint= False): + + D = self.Div + G = self.Grad + MfRhoIDeriv = self.MfRhoIDeriv + + if adjoint: + # if self._makeASymmetric is True: + # v = V * v + return(MfRhoIDeriv( G * u ).T) * ( D.T * v) + + # I think we should deprecate this for DC problem. + # if self._makeASymmetric is True: + # return V.T * ( D * ( MfRhoIDeriv( D.T * ( V * u ) ) * v ) ) + return D * (MfRhoIDeriv( G * u ) * v) + + def getRHS(self): + """ + RHS for the DC problem + + q + """ + + RHS = self.getSourceTerm() + + # I think we should deprecate this for DC problem. + # if self._makeASymmetric is True: + # return self.Vol.T * RHS + + return RHS + + def getRHSDeriv(self, src, v, adjoint=False): + """ + Derivative of the right hand side with respect to the model + """ + # TODO: add qDeriv for RHS depending on m + # qDeriv = src.evalDeriv(self, adjoint=adjoint) + # return qDeriv + return Zero() + + def setBC(self): + if self.mesh.dim==3: + fxm,fxp,fym,fyp,fzm,fzp = self.mesh.faceBoundaryInd + gBFxm = self.mesh.gridFx[fxm,:] + gBFxp = self.mesh.gridFx[fxp,:] + gBFym = self.mesh.gridFy[fym,:] + gBFyp = self.mesh.gridFy[fyp,:] + gBFzm = self.mesh.gridFz[fzm,:] + gBFzp = self.mesh.gridFz[fzp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + temp_zm, temp_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2]) + + alpha_xm, alpha_xp = temp_xm*0., temp_xp*0. + alpha_ym, alpha_yp = temp_ym*0., temp_yp*0. + alpha_zm, alpha_zp = temp_zm*0., temp_zp*0. + + beta_xm, beta_xp = temp_xm, temp_xp + beta_ym, beta_yp = temp_ym, temp_yp + beta_zm, beta_zp = temp_zm, temp_zp + + gamma_xm, gamma_xp = temp_xm*0., temp_xp*0. + gamma_ym, gamma_yp = temp_ym*0., temp_yp*0. + gamma_zm, gamma_zp = temp_zm*0., temp_zp*0. + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp] + + elif self.mesh.dim==2: + + fxm,fxp,fym,fyp = self.mesh.faceBoundaryInd + gBFxm = self.mesh.gridFx[fxm,:] + gBFxp = self.mesh.gridFx[fxp,:] + gBFym = self.mesh.gridFy[fym,:] + gBFyp = self.mesh.gridFy[fyp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + + alpha_xm, alpha_xp = temp_xm*0., temp_xp*0. + alpha_ym, alpha_yp = temp_ym*0., temp_yp*0. + + beta_xm, beta_xp = temp_xm, temp_xp + beta_ym, beta_yp = temp_ym, temp_yp + + gamma_xm, gamma_xp = temp_xm*0., temp_xp*0. + gamma_ym, gamma_yp = temp_ym*0., temp_yp*0. + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp] + + x_BC, y_BC = getxBCyBC_CC(self.mesh, alpha, beta, gamma) + V = self.Vol + self.Div = V * self.mesh.faceDiv + P_BC, B = self.mesh.getBCProjWF_simple() + M = B*self.mesh.aveCC2F + self.Grad = self.Div.T - P_BC*Utils.sdiag(y_BC)*M + + +class Problem3D_N(BaseIPProblem): + + _solutionType = 'phiSolution' + _formulation = 'EB' # N potentials means B is on faces + fieldsPair = Fields_N + + def __init__(self, mesh, **kwargs): + BaseIPProblem.__init__(self, mesh, **kwargs) + + def getA(self): + """ + + Make the A matrix for the cell centered DC resistivity problem + + A = G.T MeSigma G + + """ + + MeSigma = self.MeSigma + Grad = self.mesh.nodalGrad + A = Grad.T * MeSigma * Grad + + # Handling Null space of A + A[0,0] = A[0,0] + 1. + + return A + + def getADeriv(self, u, v, adjoint=False): + """ + + Product of the derivative of our system matrix with respect to the model and a vector + + """ + MeSigma = self.MeSigma + Grad = self.mesh.nodalGrad + if not adjoint: + return Grad.T*(self.MeSigmaDeriv(Grad*u)*v) + elif adjoint: + return self.MeSigmaDeriv(Grad*u).T * (Grad*v) + + + def getRHS(self): + """ + RHS for the DC problem + + q + """ + + RHS = self.getSourceTerm() + return RHS + + def getRHSDeriv(self, src, v, adjoint=False): + """ + Derivative of the right hand side with respect to the model + """ + # TODO: add qDeriv for RHS depending on m + # qDeriv = src.evalDeriv(self, adjoint=adjoint) + # return qDeriv + return Zero() + +if __name__ == '__main__': + + + cs = 12.5 + hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)] + hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)] + hz = [(cs,7, -1.3),(cs,20)] + mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN") + sigma = np.ones(mesh.nC) + prob = BaseIPProblem(mesh, sigma=sigma) + + diff --git a/SimPEG/EM/Static/IP/SurveyIP.py b/SimPEG/EM/Static/IP/SurveyIP.py new file mode 100644 index 00000000..c980d927 --- /dev/null +++ b/SimPEG/EM/Static/IP/SurveyIP.py @@ -0,0 +1,23 @@ +import SimPEG +from SimPEG.EM.Base import BaseEMSurvey +from SimPEG import sp, Survey +from SimPEG.Utils import Zero, Identity +from SimPEG.EM.Static.DC.SrcDC import BaseSrc +from SimPEG.EM.Static.DC.RxDC import BaseRx + +class Survey(BaseEMSurvey): + rxPair = BaseRx + srcPair = BaseSrc + + def __init__(self, srcList, **kwargs): + self.srcList = srcList + BaseEMSurvey.__init__(self, srcList, **kwargs) + + def dpred(self, m, f=None): + """ + Predicted data. + + .. math:: + d_\\text{pred} = Pf(m) + """ + return self.prob.Jvec(m, m, f=f) diff --git a/SimPEG/EM/Static/IP/__init__.py b/SimPEG/EM/Static/IP/__init__.py new file mode 100644 index 00000000..663117d3 --- /dev/null +++ b/SimPEG/EM/Static/IP/__init__.py @@ -0,0 +1,2 @@ +from ProblemIP import Problem3D_CC, Problem3D_N +from SurveyIP import Survey diff --git a/SimPEG/EM/Static/SIP/ProblemSIP.py b/SimPEG/EM/Static/SIP/ProblemSIP.py new file mode 100644 index 00000000..97be624d --- /dev/null +++ b/SimPEG/EM/Static/SIP/ProblemSIP.py @@ -0,0 +1,445 @@ +from SimPEG import Problem, Utils, Maps, Mesh +from SimPEG.EM.Base import BaseEMProblem +from SimPEG.EM.Static.DC.FieldsDC import Fields, Fields_CC, Fields_N +from SimPEG.Utils import sdiag +import numpy as np +from SimPEG.Utils import Zero +from SimPEG.EM.Static.DC import getxBCyBC_CC +from SurveySIP import Survey, Data + +class ColeColePropMap(Maps.PropMap): + """ + Property Map for EM Problems. The electrical conductivity (\\(\\sigma\\)) is the default inversion property, and the default value of the magnetic permeability is that of free space (\\(\\mu = 4\\pi\\times 10^{-7} \\) H/m) + """ + + eta = Maps.Property("Electrical Conductivity", defaultInvProp=True) + tau = Maps.Property("Electrical Conductivity", defaultVal=0.1, propertyLink=('taui', Maps.ReciprocalMap)) + taui = Maps.Property("Electrical Conductivity", defaultVal=1., propertyLink=('tau', Maps.ReciprocalMap)) + c = Maps.Property("Electrical Conductivity", defaultVal=1.) + + +class BaseSIPProblem(BaseEMProblem): + + surveyPair = Survey + fieldsPair = Fields + dataPair = Data + PropMap = ColeColePropMap + Ainv = None + sigma = None + rho = None + f = None + Ainv = None + + def DebyeTime(self, t): + peta = self.curModel.eta*np.exp(-self.curModel.taui*t) + return peta + + def EtaDeriv(self, t, v, adjoint=False): + v = np.array(v, dtype=float) + if adjoint: + return self.curModel.etaDeriv.T * (np.exp(-self.curModel.taui*t)*v) + else: + return np.exp(-self.curModel.taui*t) * (self.curModel.etaDeriv*v) + + + def TauiDeriv(self, t, v, adjoint=False): + v = np.array(v, dtype=float) + if adjoint: + return -self.curModel.tauiDeriv.T * (self.curModel.eta*t*np.exp(-self.curModel.taui*t)*v) + else: + return -self.curModel.eta*t*np.exp(-self.curModel.taui*t) * (self.curModel.tauiDeriv*v) + + def fields(self, m): + self.curModel = m + if self.f is None: + self.f = self.fieldsPair(self.mesh, self.survey) + if self.Ainv == None: + A = self.getA() + self.Ainv = self.Solver(A, **self.solverOpts) + RHS = self.getRHS() + u = self.Ainv * RHS + Srcs = self.survey.srcList + self.f[Srcs, self._solutionType] = u + return self.f + + def forward(self, m, f=None): + + if f is None: + f = self.fields(m) + + self.curModel = m + Jv = self.dataPair(self.survey) #same size as the data + # A = self.getA() + JvAll = [] + for tind in range(len(self.survey.times)): + #Pseudo-chareability + t = self.survey.times[tind] + v = self.DebyeTime(t) + for src in self.survey.srcList: + u_src = f[src, self._solutionType] # solution vector + dA_dm_v = self.getADeriv(u_src, v) + dRHS_dm_v = self.getRHSDeriv(src, v) + du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v ) + for rx in src.rxList: + timeindex = rx.getTimeP(self.survey.times) + if timeindex[tind]: + df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False) + Jv[src, rx, t] = rx.evalDeriv(src, self.mesh, f, df_dm_v) + + # Conductivity (d u / d log sigma) + if self._formulation is 'EB': + return -Utils.mkvc(Jv) + # Resistivity (d u / d log rho) + if self._formulation is 'HJ': + return Utils.mkvc(Jv) + + def Jvec(self, m, v, f=None): + + if f is None: + f = self.fields(m) + + self.curModel = m + Jv = self.dataPair(self.survey) #same size as the data + # A = self.getA() + JvAll = [] + #Assume only eta and tau (eta first then tau) + # v = [2*Mx1] + v = v.reshape((int(v.size/2), 2), order='F') + + for tind in range(len(self.survey.times)): + t = self.survey.times[tind] + v0 = self.EtaDeriv(t, v[:,0]) + v1 = self.TauiDeriv(t, v[:,1]) + for src in self.survey.srcList: + u_src = f[src, self._solutionType] # solution vector + dA_dm_v0 = self.getADeriv(u_src, v0) + dRHS_dm_v0 = self.getRHSDeriv(src, v0) + du_dm_v0 = self.Ainv * ( - dA_dm_v0 + dRHS_dm_v0 ) + dA_dm_v1 = self.getADeriv(u_src, v1) + dRHS_dm_v1 = self.getRHSDeriv(src, v1) + du_dm_v1 = self.Ainv * ( - dA_dm_v1 + dRHS_dm_v1 ) + for rx in src.rxList: + timeindex = rx.getTimeP(self.survey.times) + if timeindex[tind]: + df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_dm_v0 = df_dmFun(src, du_dm_v0, v0, adjoint=False) + df_dm_v1 = df_dmFun(src, du_dm_v1, v1, adjoint=False) + Jv[src, rx, t] = rx.evalDeriv(src, self.mesh, f, df_dm_v0) + Jv[src, rx, t] += rx.evalDeriv(src, self.mesh, f, df_dm_v1) + # Conductivity (d u / d log sigma) + if self._formulation is 'EB': + return -Jv.tovec() + # Resistivity (d u / d log rho) + if self._formulation is 'HJ': + return Jv.tovec() + + def Jtvec(self, m, v, f=None): + if f is None: + f = self.fields(m) + + self.curModel = m + + # Ensure v is a data object. + if not isinstance(v, self.dataPair): + v = self.dataPair(self.survey, v) + + Jtv= np.zeros(m.size) + for tind in range(len(self.survey.times)): + t = self.survey.times[tind] + for src in self.survey.srcList: + u_src = f[src, self._solutionType] + for rx in src.rxList: + timeindex = rx.getTimeP(self.survey.times) + if timeindex[tind]: + PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx, t], adjoint=True) # wrt f, need possibility wrt m + df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True) + ATinvdf_duT = self.Ainv * df_duT + dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True) + dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True) + du_dmT = -dA_dmT + dRHS_dmT + Jtv += np.r_[self.EtaDeriv(self.survey.times[tind], du_dmT, adjoint=True), self.TauiDeriv(self.survey.times[tind], du_dmT, adjoint=True)] + + # Conductivity ((d u / d log sigma).T) + if self._formulation is 'EB': + return -Jtv + # Conductivity ((d u / d log rho).T) + if self._formulation is 'HJ': + return Jtv + + def getSourceTerm(self): + """ + takes concept of source and turns it into a matrix + """ + """ + Evaluates the sources, and puts them in matrix form + + :rtype: (numpy.ndarray, numpy.ndarray) + :return: q (nC or nN, nSrc) + """ + + Srcs = self.survey.srcList + + if self._formulation is 'EB': + n = self.mesh.nN + # return NotImplementedError + + elif self._formulation is 'HJ': + n = self.mesh.nC + + q = np.zeros((n, len(Srcs))) + + for i, src in enumerate(Srcs): + q[:,i] = src.eval(self) + return q + + @property + def deleteTheseOnModelUpdate(self): + toDelete = [] + return toDelete + + # assume log rho or log cond + @property + def MeSigma(self): + """ + Edge inner product matrix for \\(\\sigma\\). Used in the E-B formulation + """ + if getattr(self, '_MeSigma', None) is None: + self._MeSigma = self.mesh.getEdgeInnerProduct(self.sigma) + return self._MeSigma + + @property + def MfRhoI(self): + """ + Inverse of :code:`MfRho` + """ + if getattr(self, '_MfRhoI', None) is None: + self._MfRhoI = self.mesh.getFaceInnerProduct(self.rho, invMat=True) + return self._MfRhoI + + def MfRhoIDeriv(self,u): + """ + Derivative of :code:`MfRhoI` with respect to the model. + """ + + dMfRhoI_dI = -self.MfRhoI**2 + dMf_drho = self.mesh.getFaceInnerProductDeriv(self.rho)(u) + drho_dlogrho = Utils.sdiag(self.rho) + return dMfRhoI_dI * ( dMf_drho * ( drho_dlogrho)) + + # TODO: This should take a vector + def MeSigmaDeriv(self, u): + """ + Derivative of MeSigma with respect to the model + """ + dsigma_dlogsigma = Utils.sdiag(self.sigma) + return self.mesh.getEdgeInnerProductDeriv(self.sigma)(u) * dsigma_dlogsigma + +class Problem3D_CC(BaseSIPProblem): + + _solutionType = 'phiSolution' + _formulation = 'HJ' # CC potentials means J is on faces + fieldsPair = Fields_CC + + def __init__(self, mesh, **kwargs): + BaseSIPProblem.__init__(self, mesh, **kwargs) + self.setBC() + + def getA(self): + """ + + Make the A matrix for the cell centered DC resistivity problem + + A = D MfRhoI G + + """ + + D = self.Div + G = self.Grad + # TODO: this won't work for full anisotropy + MfRhoI = self.MfRhoI + A = D * MfRhoI * G + + # I think we should deprecate this for DC problem. + # if self._makeASymmetric is True: + # return V.T * A + return A + + def getADeriv(self, u, v, adjoint= False): + + D = self.Div + G = self.Grad + MfRhoIDeriv = self.MfRhoIDeriv + + if adjoint: + # if self._makeASymmetric is True: + # v = V * v + return(MfRhoIDeriv( G * u ).T) * ( D.T * v) + + # I think we should deprecate this for DC problem. + # if self._makeASymmetric is True: + # return V.T * ( D * ( MfRhoIDeriv( D.T * ( V * u ) ) * v ) ) + return D * (MfRhoIDeriv( G * u ) * v) + + def getRHS(self): + """ + RHS for the DC problem + + q + """ + + RHS = self.getSourceTerm() + + # I think we should deprecate this for DC problem. + # if self._makeASymmetric is True: + # return self.Vol.T * RHS + + return RHS + + def getRHSDeriv(self, src, v, adjoint=False): + """ + Derivative of the right hand side with respect to the model + """ + # TODO: add qDeriv for RHS depending on m + # qDeriv = src.evalDeriv(self, adjoint=adjoint) + # return qDeriv + return Zero() + + def setBC(self): + if self.mesh.dim==3: + fxm,fxp,fym,fyp,fzm,fzp = self.mesh.faceBoundaryInd + gBFxm = self.mesh.gridFx[fxm,:] + gBFxp = self.mesh.gridFx[fxp,:] + gBFym = self.mesh.gridFy[fym,:] + gBFyp = self.mesh.gridFy[fyp,:] + gBFzm = self.mesh.gridFz[fzm,:] + gBFzp = self.mesh.gridFz[fzp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + temp_zm, temp_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2]) + + alpha_xm, alpha_xp = temp_xm*0., temp_xp*0. + alpha_ym, alpha_yp = temp_ym*0., temp_yp*0. + alpha_zm, alpha_zp = temp_zm*0., temp_zp*0. + + beta_xm, beta_xp = temp_xm, temp_xp + beta_ym, beta_yp = temp_ym, temp_yp + beta_zm, beta_zp = temp_zm, temp_zp + + gamma_xm, gamma_xp = temp_xm*0., temp_xp*0. + gamma_ym, gamma_yp = temp_ym*0., temp_yp*0. + gamma_zm, gamma_zp = temp_zm*0., temp_zp*0. + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp] + + elif self.mesh.dim==2: + + fxm,fxp,fym,fyp = self.mesh.faceBoundaryInd + gBFxm = self.mesh.gridFx[fxm,:] + gBFxp = self.mesh.gridFx[fxp,:] + gBFym = self.mesh.gridFy[fym,:] + gBFyp = self.mesh.gridFy[fyp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + + alpha_xm, alpha_xp = temp_xm*0., temp_xp*0. + alpha_ym, alpha_yp = temp_ym*0., temp_yp*0. + + beta_xm, beta_xp = temp_xm, temp_xp + beta_ym, beta_yp = temp_ym, temp_yp + + gamma_xm, gamma_xp = temp_xm*0., temp_xp*0. + gamma_ym, gamma_yp = temp_ym*0., temp_yp*0. + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp] + + x_BC, y_BC = getxBCyBC_CC(self.mesh, alpha, beta, gamma) + V = self.Vol + self.Div = V * self.mesh.faceDiv + P_BC, B = self.mesh.getBCProjWF_simple() + M = B*self.mesh.aveCC2F + self.Grad = self.Div.T - P_BC*Utils.sdiag(y_BC)*M + + +class Problem3D_N(BaseSIPProblem): + + _solutionType = 'phiSolution' + _formulation = 'EB' # N potentials means B is on faces + fieldsPair = Fields_N + + def __init__(self, mesh, **kwargs): + BaseSIPProblem.__init__(self, mesh, **kwargs) + + def getA(self): + """ + + Make the A matrix for the cell centered DC resistivity problem + + A = G.T MeSigma G + + """ + + # TODO: this won't work for full anisotropy + MeSigma = self.MeSigma + Grad = self.mesh.nodalGrad + A = Grad.T * MeSigma * Grad + + # Handling Null space of A + A[0,0] = A[0,0] + 1. + + return A + + def getADeriv(self, u, v, adjoint=False): + """ + + Product of the derivative of our system matrix with respect to the model and a vector + + """ + MeSigma = self.MeSigma + Grad = self.mesh.nodalGrad + if not adjoint: + return Grad.T*(self.MeSigmaDeriv(Grad*u)*v) + elif adjoint: + return self.MeSigmaDeriv(Grad*u).T * (Grad*v) + + + def getRHS(self): + """ + RHS for the DC problem + + q + """ + + RHS = self.getSourceTerm() + return RHS + + def getRHSDeriv(self, src, v, adjoint=False): + """ + Derivative of the right hand side with respect to the model + """ + # TODO: add qDeriv for RHS depending on m + # qDeriv = src.evalDeriv(self, adjoint=adjoint) + # return qDeriv + return Zero() + +if __name__ == '__main__': + + + cs = 12.5 + hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)] + hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)] + hz = [(cs,7, -1.3),(cs,20)] + mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN") + sigma = np.ones(mesh.nC) + prob = BaseSIPProblem(mesh, sigma=sigma) + + diff --git a/SimPEG/EM/Static/SIP/Regularization.py b/SimPEG/EM/Static/SIP/Regularization.py new file mode 100644 index 00000000..3c20e6fe --- /dev/null +++ b/SimPEG/EM/Static/SIP/Regularization.py @@ -0,0 +1,204 @@ +from SimPEG import Utils, Maps, Mesh, sp, np +from SimPEG.Regularization import BaseRegularization, Simple + +class MultiRegularization(Simple): + """ + **MultiRegularization Class** + + This is used to regularize the model space + having multiple models [m1, m2, m3, ...] :: + + reg = Regularization(mesh) + + """ + nModels = None # Number of models + ratios = None # Ratio for different models + crossgrad = False # Use cross gradient or not + betacross = 1. + wx = [] + wy = [] + wz = [] + + def __init__(self, mesh, mapping=None, indActive=None, **kwargs): + BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) + if self.nModels == None: + raise Exception("Put nModels as a initial input!") + if self.ratios == None: + self.ratios = [1. for imodel in range(self.nModels)] + + @property + def Wsmall(self): + """Regularization matrix Wsmall""" + if getattr(self,'_Wsmall', None) is None: + vecs = [] + for imodel in range(self.nModels): + vecs.append((self.regmesh.vol*self.alpha_s*self.wght*self.ratios[imodel])**0.5) + self._Wsmall = Utils.sdiag(np.hstack(vecs)) + return self._Wsmall + + @property + def Wx(self): + """Regularization matrix Wx""" + if getattr(self, '_Wx', None) is None: + mats = [] + for imodel in range(self.nModels): + self.wx.append(Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x*self.ratios[imodel]*(self.regmesh.aveCC2Fx*self.wght))**0.5)) + mats.append(self.wx[imodel]*self.regmesh.cellDiffxStencil) + self._Wx = sp.block_diag(mats) + return self._Wx + + @property + def Wy(self): + """Regularization matrix Wy""" + if getattr(self, '_Wy', None) is None: + mats = [] + for imodel in range(self.nModels): + self.wy.append(Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol*self.alpha_y*self.ratios[imodel]*(self.regmesh.aveCC2Fy*self.wght))**0.5)) + mats.append(self.wy[imodel]*self.regmesh.cellDiffyStencil) + self._Wy = sp.block_diag(mats) + return self._Wy + + @property + def Wz(self): + """Regularization matrix Wz""" + if getattr(self, '_Wz', None) is None: + mats = [] + for imodel in range(self.nModels): + self.wz.append(Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z*self.ratios[imodel]*(self.regmesh.aveCC2Fz*self.wght))**0.5)) + mats.append(self.wz[imodel]*self.regmesh.cellDiffzStencil) + self._Wz = sp.block_diag(mats) + return self._Wz + + @property + def Wsmooth(self): + """Full smoothness regularization matrix W""" + if getattr(self, '_Wsmooth', None) is None: + wlist = (self.Wx,) + if self.regmesh.dim > 1: + wlist += (self.Wy,) + if self.regmesh.dim > 2: + wlist += (self.Wz,) + self._Wsmooth = sp.vstack(wlist) + return self._Wsmooth + + @property + def W(self): + """Full regularization matrix W""" + if getattr(self, '_W', None) is None: + wlist = (self.Wsmall, self.Wsmooth) + self._W = sp.vstack(wlist) + return self._W + + + @Utils.timeIt + def eval(self, m): + return self._evalSmall(m) + self._evalSmooth(m) + + @Utils.timeIt + 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) + + def cross(a,b): + ax, ay, az = a[0], a[1], a[2] + bx, by, bz = b[0], b[1], b[2] + cx = ay*bz - az*by + cy = az*bx - ax*bz + cz = ax*by - ay*bx + return [cx, cy, cz] + + # TODO: Implement Cross Gradients.. + @Utils.timeIt + def _evalCross(self, m): + if self.crossgrad == False: + return 0. + elif self.crossgrad == True: + M = (self.mapping * m).reshape((self.regmesh.nC, self.nModels), order="F") + + ax = self.regmesh.aveFx2CC*self.regmesh.wx[0]*M[:,0] + ay = self.regmesh.aveFy2CC*self.regmesh.wy[0]*M[:,0] + az = self.regmesh.aveFz2CC*self.regmesh.wz[0]*M[:,0] + bx = self.regmesh.aveFx2CC*self.regmesh.wx[1]*M[:,1] + by = self.regmesh.aveFy2CC*self.regmesh.wy[1]*M[:,1] + bz = self.regmesh.aveFz2CC*self.regmesh.wz[1]*M[:,1] + #ab + out_ab = cross([ax, ay, az], [bx, by, bz]) + r = np.r_[out_ab[0], out_ab[1], out_ab[2]]*np.sqrt(self.betacross) + + if self.nModels == 3: + cx = self.regmesh.aveFx2CC*self.regmesh.wx[1]*M[:,1] + cy = self.regmesh.aveFy2CC*self.regmesh.wy[1]*M[:,1] + cz = self.regmesh.aveFz2CC*self.regmesh.wz[1]*M[:,1] + #ac + out_ac = cross([ax, ay, az], [cx, cy, cz]) + #bc + out_bc = cross([bx, by, bz], [cx, cy, cz]) + r = np.r_[r, np.hstack(out_ac)*np.sqrt(self.betacross), np.hstack(out_bc)*np.sqrt(self.betacross)] + + return 0.5 * r.dot(r) + + @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})} + + """ + deriv = self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + if self.crossgrad==True: + deriv += self._evalCrossDeriv(m) + return deriv + + @Utils.timeIt + def _evalCrossDeriv(self,m): + r = self.Wsmall * ( self.mapping * (m - self.mref) ) + return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) ) + + @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 + + The regularization is: + + .. math:: + + R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})} + + So the second derivative is straight forward: + + .. math:: + + R(m) = \mathbf{W^\\top W} + + """ + mD = self.mapping.deriv(m - self.mref) + if v is None: + return mD.T * self.W.T * self.W * mD + + return mD.T * ( self.W.T * ( self.W * ( mD * v) ) ) + + + diff --git a/SimPEG/EM/Static/SIP/RxSIP.py b/SimPEG/EM/Static/SIP/RxSIP.py new file mode 100644 index 00000000..b8144d0f --- /dev/null +++ b/SimPEG/EM/Static/SIP/RxSIP.py @@ -0,0 +1,88 @@ +import SimPEG +import numpy as np +from SimPEG.Utils import Zero, closestPoints + +class BaseRx(SimPEG.Survey.BaseTimeRx): + locs = None + rxType = None + + knownRxTypes = { + 'phi':['phi',None], + 'ex':['e','x'], + 'ey':['e','y'], + 'ez':['e','z'], + 'jx':['j','x'], + 'jy':['j','y'], + 'jz':['j','z'], + } + + def __init__(self, locs, times, rxType, **kwargs): + SimPEG.Survey.BaseTimeRx.__init__(self, locs, times, rxType, **kwargs) + + @property + def projField(self): + """Field Type projection (e.g. e b ...)""" + return self.knownRxTypes[self.rxType][0] + + def projGLoc(self, f): + """Grid Location projection (e.g. Ex Fy ...)""" + comp = self.knownRxTypes[self.rxType][1] + if comp is not None: + return f._GLoc(self.rxType) + comp + return f._GLoc(self.rxType) + + def getTimeP(self, timesall): + """ + Returns the time projection matrix. + + .. note:: + + This is not stored in memory, but is created on demand. + """ + time_inds = np.in1d(timesall, self.times) + return time_inds + + def evalDeriv(self, src, mesh, f, v, adjoint=False): + P = self.getP(mesh, self.projGLoc(f)) + if not adjoint: + return P*v + elif adjoint: + return P.T*v + + +# DC.Rx.Dipole(locs) +class Dipole(BaseRx): + + def __init__(self, locsM, locsN, times, rxType = 'phi', **kwargs): + assert locsM.shape == locsN.shape, 'locsM and locsN need to be the same size' + locs = [locsM, locsN] + # We may not need this ... + BaseRx.__init__(self, locs, times, rxType) + + @property + def nD(self): + """Number of data in the receiver.""" + # return self.locs[0].shape[0] * len(self.times) + return self.locs[0].shape[0] + + @property + def nRx(self): + """Number of data in the receiver.""" + return self.locs[0].shape[0] + + # Not sure why ... + # return int(self.locs[0].size / 2) + + + def getP(self, mesh, Gloc): + if mesh in self._Ps: + return self._Ps[mesh] + + P0 = mesh.getInterpolationMat(self.locs[0], Gloc) + P1 = mesh.getInterpolationMat(self.locs[1], Gloc) + P = P0 - P1 + + if self.storeProjections: + self._Ps[mesh] = P + + return P diff --git a/SimPEG/EM/Static/SIP/SrcSIP.py b/SimPEG/EM/Static/SIP/SrcSIP.py new file mode 100644 index 00000000..b1f0a452 --- /dev/null +++ b/SimPEG/EM/Static/SIP/SrcSIP.py @@ -0,0 +1,64 @@ +import SimPEG +# from SimPEG.EM.Base import BaseEMSurvey +from SimPEG.Utils import Zero, closestPoints, mkvc +import numpy as np + +class BaseSrc(SimPEG.Survey.BaseSrc): + + current = 1.0 + loc = None + + def __init__(self, rxList, **kwargs): + SimPEG.Survey.BaseSrc.__init__(self, rxList, **kwargs) + + def eval(self, prob): + raise NotImplementedError + + def evalDeriv(self, prob): + return Zero() + + @property + def nD(self): + """Number of data""" + return self.vnD.sum() + + @property + def vnD(self): + """Vector number of data""" + return np.array([rx.nD*len(rx.times) for rx in self.rxList]) + + + +class Dipole(BaseSrc): + + def __init__(self, rxList, locA, locB, **kwargs): + assert locA.shape == locB.shape, 'Shape of locA and locB should be the same' + self.loc = [locA, locB] + BaseSrc.__init__(self, rxList, **kwargs) + + def eval(self, prob): + if prob._formulation == 'HJ': + inds = closestPoints(prob.mesh, self.loc, gridLoc='CC') + q = np.zeros(prob.mesh.nC) + q[inds] = self.current * np.r_[1., -1.] + elif prob._formulation == 'EB': + qa = prob.mesh.getInterpolationMat(self.loc[0], locType='N').todense() + qb = -prob.mesh.getInterpolationMat(self.loc[1], locType='N').todense() + q = self.current * mkvc(qa+qb) + return q + +class Pole(BaseSrc): + + def __init__(self, rxList, loc, **kwargs): + BaseSrc.__init__(self, rxList, loc=loc, **kwargs) + + def eval(self, prob): + if prob._formulation == 'HJ': + inds = closestPoints(prob.mesh, self.loc) + q = np.zeros(prob.mesh.nC) + q[inds] = self.current * np.r_[1.] + elif prob._formulation == 'EB': + q = prob.mesh.getInterpolationMat(self.loc, locType='N').todense() + q = self.current * mkvc(q) + return q + diff --git a/SimPEG/EM/Static/SIP/SurveySIP.py b/SimPEG/EM/Static/SIP/SurveySIP.py new file mode 100644 index 00000000..9362a35c --- /dev/null +++ b/SimPEG/EM/Static/SIP/SurveySIP.py @@ -0,0 +1,102 @@ +import SimPEG +from SimPEG.EM.Base import BaseEMSurvey +from SimPEG import np, sp, Survey, Utils +from SimPEG.Utils import Zero, Identity +from SimPEG.EM.Static.SIP.SrcSIP import BaseSrc +from SimPEG.EM.Static.SIP.RxSIP import BaseRx +import uuid + + +class Survey(BaseEMSurvey): + rxPair = BaseRx + srcPair = BaseSrc + times = None + + def __init__(self, srcList, **kwargs): + self.srcList = srcList + BaseEMSurvey.__init__(self, srcList, **kwargs) + self.getUniqueTimes() + + def getUniqueTimes(self): + time_rx = [] + for src in self.srcList: + for rx in src.rxList: + time_rx.append(rx.times) + self.times = np.unique(np.hstack(time_rx)) + + def dpred(self, m, f=None): + """ + Predicted data. + + .. math:: + d_\\text{pred} = Pf(m) + """ + return self.prob.forward(m, f=f) + + +class Data(SimPEG.Survey.Data): + """Fancy data storage by Src and Rx""" + + def __init__(self, survey, v=None): + self.uid = str(uuid.uuid4()) + self.survey = survey + self._dataDict = {} + for src in self.survey.srcList: + self._dataDict[src] = {} + for rx in src.rxList: + self._dataDict[src][rx] = {} + + if v is not None: + self.fromvec(v) + + def _ensureCorrectKey(self, key): + if type(key) is tuple: + if len(key) is not 3: + raise KeyError('Key must be [Src, Rx, tInd]') + if key[0] not in self.survey.srcList: + raise KeyError('Src Key must be a source in the survey.') + if key[1] not in key[0].rxList: + raise KeyError('Rx Key must be a receiver for the source.') + return key + elif isinstance(key, self.survey.srcPair): + if key not in self.survey.srcList: + raise KeyError('Key must be a source in the survey.') + return key, None, None + else: + raise KeyError('Key must be [Src] or [Src,Rx] or [Src, Rx, tInd]') + + def __setitem__(self, key, value): + src, rx, t = self._ensureCorrectKey(key) + assert rx is not None, 'set data using [Src, Rx]' + assert isinstance(value, np.ndarray), 'value must by ndarray' + assert value.size == rx.nD, "value must have the same number of data as the source." + self._dataDict[src][rx][t] = Utils.mkvc(value) + + def __getitem__(self, key): + src, rx, t = self._ensureCorrectKey(key) + if rx is not None: + if rx not in self._dataDict[src]: + raise Exception('Data for receiver has not yet been set.') + return self._dataDict[src][rx][t] + + return np.concatenate([self[src,rx, t] for rx in src.rxList]) + + def tovec(self): + val = [] + for src in self.survey.srcList: + for rx in src.rxList: + for t in rx.times: + val.append(self[src, rx, t]) + return np.concatenate(val) + + + def fromvec(self, v): + v = Utils.mkvc(v) + assert v.size == self.survey.nD, 'v must have the correct number of data.' + indBot, indTop = 0, 0 + for src in self.survey.srcList: + for rx in src.rxList: + for t in rx.times: + indTop += rx.nRx + self[src, rx, t] = v[indBot:indTop] + indBot += rx.nRx diff --git a/SimPEG/EM/Static/SIP/__init__.py b/SimPEG/EM/Static/SIP/__init__.py new file mode 100644 index 00000000..1de46fcf --- /dev/null +++ b/SimPEG/EM/Static/SIP/__init__.py @@ -0,0 +1,5 @@ +from ProblemSIP import Problem3D_CC, Problem3D_N +from SurveySIP import Survey, Data +import SrcSIP as Src #Pole +import RxSIP as Rx +from Regularization import MultiRegularization diff --git a/SimPEG/EM/Static/Utils/StaticUtils.py b/SimPEG/EM/Static/Utils/StaticUtils.py new file mode 100644 index 00000000..43181c8a --- /dev/null +++ b/SimPEG/EM/Static/Utils/StaticUtils.py @@ -0,0 +1,317 @@ +from SimPEG import np +from SimPEG.EM.Static import DC, IP + +def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None): + """ + Read list of 2D tx-rx location and plot a speudo-section of apparent + resistivity. + + Assumes flat topo for now... + + Input: + :param d2D, z0 + :switch stype -> Either 'pdp' (pole-dipole) | 'dpdp' (dipole-dipole) + :switch dtype=-> Either 'appr' (app. res) | 'appc' (app. con) | 'volt' (potential) + Output: + :figure scatter plot overlayed on image + + Edited Feb 17th, 2016 + + @author: dominiquef + + """ + from SimPEG import np + from scipy.interpolate import griddata + import pylab as plt + + # Set depth to 0 for now + z0 = 0. + + # Pre-allocate + midx = [] + midz = [] + rho = [] + LEG = [] + count = 0 # Counter for data + for ii in range(DCsurvey.nSrc): + + Tx = DCsurvey.srcList[ii].loc + Rx = DCsurvey.srcList[ii].rxList[0].locs + + nD = DCsurvey.srcList[ii].rxList[0].nD + + data = DCsurvey.dobs[count:count+nD] + count += nD + + # Get distances between each poles A-B-M-N + if stype == 'pdp': + MA = np.abs(Tx[0] - Rx[0][:,0]) + NA = np.abs(Tx[0] - Rx[1][:,0]) + MN = np.abs(Rx[1][:,0] - Rx[0][:,0]) + + # Create mid-point location + Cmid = Tx[0] + Pmid = (Rx[0][:,0] + Rx[1][:,0])/2 + if DCsurvey.mesh.dim == 2: + zsrc = Tx[1] + elif DCsurvey.mesh.dim ==3: + zsrc = Tx[2] + + elif stype == 'dpdp': + MA = np.abs(Tx[0][0] - Rx[0][:,0]) + MB = np.abs(Tx[1][0] - Rx[0][:,0]) + NA = np.abs(Tx[0][0] - Rx[1][:,0]) + NB = np.abs(Tx[1][0] - Rx[1][:,0]) + + # Create mid-point location + Cmid = (Tx[0][0] + Tx[1][0])/2 + Pmid = (Rx[0][:,0] + Rx[1][:,0])/2 + if DCsurvey.mesh.dim == 2: + zsrc = (Tx[0][1] + Tx[1][1])/2 + elif DCsurvey.mesh.dim ==3: + zsrc = (Tx[0][2] + Tx[1][2])/2 + + # Change output for dtype + if dtype == 'volt': + + rho = np.hstack([rho,data]) + + else: + + # Compute pant leg of apparent rho + if stype == 'pdp': + + leg = data * 2*np.pi * MA * ( MA + MN ) / MN + + elif stype == 'dpdp': + + leg = data * 2*np.pi / ( 1/MA - 1/MB + 1/NB - 1/NA ) + LEG.append(1./(2*np.pi) *( 1/MA - 1/MB + 1/NB - 1/NA )) + else: + print """dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """ + break + + + if dtype == 'appc': + + leg = np.log10(abs(1./leg)) + rho = np.hstack([rho,leg]) + + elif dtype == 'appr': + + leg = np.log10(abs(leg)) + rho = np.hstack([rho,leg]) + + else: + print """dtype must be 'appr' | 'appc' | 'volt' """ + break + + + midx = np.hstack([midx, ( Cmid + Pmid )/2 ]) + if DCsurvey.mesh.dim==3: + midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + zsrc ]) + elif DCsurvey.mesh.dim==2: + midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + zsrc ]) + ax = axs + + # Grid points + grid_x, grid_z = np.mgrid[np.min(midx):np.max(midx), np.min(midz):np.max(midz)] + grid_rho = griddata(np.c_[midx,midz], rho.T, (grid_x, grid_z), method='linear') + + if clim == None: + vmin, vmax = rho.min(), rho.max() + else: + vmin, vmax = clim[0], clim[1] + + grid_rho = np.ma.masked_where(np.isnan(grid_rho), grid_rho) + ph = plt.pcolormesh(grid_x[:,0],grid_z[0,:],grid_rho.T, clim=(vmin, vmax), vmin=vmin, vmax=vmax) + cbar = plt.colorbar(format="$10^{%.1f}$",fraction=0.04,orientation="horizontal") + + cmin,cmax = cbar.get_clim() + ticks = np.linspace(cmin,cmax,3) + cbar.set_ticks(ticks) + cbar.ax.tick_params(labelsize=10) + + if dtype == 'appc': + cbar.set_label("App.Cond",size=12) + elif dtype == 'appr': + cbar.set_label("App.Res.",size=12) + elif dtype == 'volt': + cbar.set_label("Potential (V)",size=12) + + # Plot apparent resistivity + ax.scatter(midx,midz,s=10,c=rho.T, vmin =vmin, vmax = vmax, clim=(vmin, vmax)) + + #ax.set_xticklabels([]) + #ax.set_yticklabels([]) + + plt.gca().set_aspect('equal', adjustable='box') + + + + return ph, LEG + +def gen_DCIPsurvey(endl, mesh, stype, a, b, n): + """ + Load in endpoints and survey specifications to generate Tx, Rx location + stations. + + Assumes flat topo for now... + + Input: + :param endl -> input endpoints [x1, y1, z1, x2, y2, z2] + :object mesh -> SimPEG mesh object + :switch stype -> "dpdp" (dipole-dipole) | "pdp" (pole-dipole) | 'gradient' + : param a, n -> pole seperation, number of rx dipoles per tx + + Output: + :param Tx, Rx -> List objects for each tx location + Lines: P1x, P1y, P1z, P2x, P2y, P2z + + Created on Wed December 9th, 2015 + + @author: dominiquef + !! Require clean up to deal with DCsurvey + """ + + from SimPEG import np + + def xy_2_r(x1,x2,y1,y2): + r = np.sqrt( np.sum((x2 - x1)**2 + (y2 - y1)**2) ) + return r + + ## Evenly distribute electrodes and put on surface + # Mesure survey length and direction + dl_len = xy_2_r(endl[0,0],endl[1,0],endl[0,1],endl[1,1]) + + dl_x = ( endl[1,0] - endl[0,0] ) / dl_len + dl_y = ( endl[1,1] - endl[0,1] ) / dl_len + + nstn = np.floor( dl_len / a ) + + # Compute discrete pole location along line + stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*a + stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*a + + if mesh.dim==2: + ztop = mesh.vectorNy[-1] + # Create line of P1 locations + M = np.c_[stn_x, np.ones(nstn).T*ztop] + # Create line of P2 locations + N = np.c_[stn_x+a*dl_x, np.ones(nstn).T*ztop] + + elif mesh.dim==3: + ztop = mesh.vectorNz[-1] + # Create line of P1 locations + M = np.c_[stn_x, stn_y, np.ones(nstn).T*ztop] + # Create line of P2 locations + N = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*ztop] + + + ## Build list of Tx-Rx locations depending on survey type + # Dipole-dipole: Moving tx with [a] spacing -> [AB a MN1 a MN2 ... a MNn] + # Pole-dipole: Moving pole on one end -> [A a MN1 a MN2 ... MNn a B] + SrcList = [] + + + if stype != 'gradient': + + for ii in range(0, int(nstn)-1): + + + if stype == 'dpdp': + tx = np.c_[M[ii,:],N[ii,:]] + elif stype == 'pdp': + tx = np.c_[M[ii,:],M[ii,:]] + + # Rx.append(np.c_[M[ii+1:indx,:],N[ii+1:indx,:]]) + + # Current elctrode seperation + AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1]) + + # Number of receivers to fit + nstn = np.min([np.floor( (AB - b) / a ) , n]) + + # Check if there is enough space, else break the loop + if nstn <= 0: + continue + + # Compute discrete pole location along line + stn_x = N[ii,0] + dl_x*b + np.array(range(int(nstn)))*dl_x*a + stn_y = N[ii,1] + dl_y*b + np.array(range(int(nstn)))*dl_y*a + + # Create receiver poles + + if mesh.dim==3: + # Create line of P1 locations + P1 = np.c_[stn_x, stn_y, np.ones(nstn).T*ztop] + # Create line of P2 locations + P2 = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*ztop] + rxClass = DC.Rx.Dipole(P1, P2) + + elif mesh.dim==2: + # Create line of P1 locations + P1 = np.c_[stn_x, np.ones(nstn).T*ztop] + # Create line of P2 locations + P2 = np.c_[stn_x+a*dl_x, np.ones(nstn).T*ztop] + rxClass = DC.Rx.Dipole_ky(P1, P2) + + if stype == 'dpdp': + srcClass = DC.Src.Dipole([rxClass], M[ii,:],N[ii,:]) + elif stype == 'pdp': + srcClass = DC.Src.Pole([rxClass], M[ii,:]) + SrcList.append(srcClass) + + elif stype == 'gradient': + + # Gradient survey only requires Tx at end of line and creates a square + # grid of receivers at in the middle at a pre-set minimum distance + + # Get the edge limit of survey area + min_x = endl[0,0] + dl_x * b + min_y = endl[0,1] + dl_y * b + + max_x = endl[1,0] - dl_x * b + max_y = endl[1,1] - dl_y * b + + box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 ) + box_w = box_l/2. + + nstn = np.floor( box_l / a ) + + # Compute discrete pole location along line + stn_x = min_x + np.array(range(int(nstn)))*dl_x*a + stn_y = min_y + np.array(range(int(nstn)))*dl_y*a + + # Define number of cross lines + nlin = int(np.floor( box_w / a )) + lind = range(-nlin,nlin+1) + + ngrad = nstn * len(lind) + + rx = np.zeros([ngrad,6]) + for ii in range( len(lind) ): + + # Move line in perpendicular direction by dipole spacing + lxx = stn_x - lind[ii]*a*dl_y + lyy = stn_y + lind[ii]*a*dl_x + + + M = np.c_[ lxx, lyy , np.ones(nstn).T*ztop] + N = np.c_[ lxx+a*dl_x, lyy+a*dl_y, np.ones(nstn).T*ztop] + rx[(ii*nstn):((ii+1)*nstn),:] = np.c_[M,N] + + if mesh.dim==3: + rxClass = DC.Rx.Dipole(rx[:,:3], rx[:,3:]) + elif mesh.dim==2: + M = M[:,[0,2]] + N = N[:,[0,2]] + rxClass = DC.Rx.Dipole_ky(rx[:,[0,2]], rx[:,[3,5]]) + srcClass = DC.Src.Dipole([rxClass], M[0,:], N[-1,:]) + SrcList.append(srcClass) + else: + print """stype must be either 'pdp', 'dpdp' or 'gradient'. """ + + + return SrcList + diff --git a/SimPEG/EM/Static/Utils/__init__.py b/SimPEG/EM/Static/Utils/__init__.py new file mode 100644 index 00000000..adc317ae --- /dev/null +++ b/SimPEG/EM/Static/Utils/__init__.py @@ -0,0 +1 @@ +from StaticUtils import * diff --git a/SimPEG/EM/Static/__init__.py b/SimPEG/EM/Static/__init__.py new file mode 100644 index 00000000..f6e0e757 --- /dev/null +++ b/SimPEG/EM/Static/__init__.py @@ -0,0 +1,3 @@ +import DC +import IP +import SIP diff --git a/SimPEG/EM/TDEM/SurveyTDEM.py b/SimPEG/EM/TDEM/SurveyTDEM.py index 4d04f0ae..9c9ed4b9 100644 --- a/SimPEG/EM/TDEM/SurveyTDEM.py +++ b/SimPEG/EM/TDEM/SurveyTDEM.py @@ -87,7 +87,7 @@ class SrcTDEM_VMD_MVP(SrcTDEM): def getInitialFields(self, mesh): """Vertical magnetic dipole, magnetic vector potential""" if self.waveformType == "STEPOFF": - print ">> Step waveform: Non-zero initial condition" + print ">> Step waveform: Non-zero initial condition" if mesh._meshType is 'CYL': if mesh.isSymmetric: MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey') @@ -96,8 +96,8 @@ class SrcTDEM_VMD_MVP(SrcTDEM): elif mesh._meshType is 'TENSOR': MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez']) else: - raise Exception('Unknown mesh for VMD') - return {"b": mesh.edgeCurl*MVP} + raise Exception('Unknown mesh for VMD') + return {"b": mesh.edgeCurl*MVP} elif self.waveformType == "GENERAL": print ">> General waveform: Zero initial condition" return {"b": np.zeros(mesh.nF)} @@ -113,7 +113,7 @@ class SrcTDEM_VMD_MVP(SrcTDEM): elif mesh._meshType is 'TENSOR': MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez']) else: - raise Exception('Unknown mesh for VMD') + raise Exception('Unknown mesh for VMD') return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP @@ -122,7 +122,7 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM): self.loc = loc self.radius = radius self.waveformType = waveformType - SrcTDEM.__init__(self,rxList) + SrcTDEM.__init__(self,rxList) def getInitialFields(self, mesh): """Circular Loop, magnetic vector potential""" @@ -153,7 +153,7 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM): elif mesh._meshType is 'TENSOR': MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius) else: - raise Exception('Unknown mesh for CircularLoop') + raise Exception('Unknown mesh for CircularLoop') return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP diff --git a/SimPEG/EM/Utils/testingUtils.py b/SimPEG/EM/Utils/testingUtils.py index 569f8e6d..22c925b6 100644 --- a/SimPEG/EM/Utils/testingUtils.py +++ b/SimPEG/EM/Utils/testingUtils.py @@ -20,56 +20,61 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False): mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C']) if useMu is True: - mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))] + mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))] else: mapping = Maps.ExpMap(mesh) x = np.array([np.linspace(-5.*cs,-2.*cs,3),np.linspace(5.*cs,2.*cs,3)]) + cs/4. #don't sample right by the source, slightly off alignment from either staggered grid XYZ = Utils.ndgrid(x,x,np.linspace(-2.*cs,2.*cs,5)) - Rx0 = EM.FDEM.Rx(XYZ, comp) + Rx0 = getattr(EM.FDEM.Rx, 'Point_' + comp[0]) + if comp[2] == 'r': + real_or_imag = 'real' + elif comp[2] == 'i': + real_or_imag = 'imag' + rx0 = Rx0(XYZ, comp[1], 'imag') Src = [] for SrcType in SrcList: if SrcType is 'MagDipole': - Src.append(EM.FDEM.Src.MagDipole([Rx0], freq=freq, loc=np.r_[0.,0.,0.])) + Src.append(EM.FDEM.Src.MagDipole([rx0], freq=freq, loc=np.r_[0.,0.,0.])) elif SrcType is 'MagDipole_Bfield': - Src.append(EM.FDEM.Src.MagDipole_Bfield([Rx0], freq=freq, loc=np.r_[0.,0.,0.])) + Src.append(EM.FDEM.Src.MagDipole_Bfield([rx0], freq=freq, loc=np.r_[0.,0.,0.])) elif SrcType is 'CircularLoop': - Src.append(EM.FDEM.Src.CircularLoop([Rx0], freq=freq, loc=np.r_[0.,0.,0.])) + Src.append(EM.FDEM.Src.CircularLoop([rx0], freq=freq, loc=np.r_[0.,0.,0.])) elif SrcType is 'RawVec': if fdemType is 'e' or fdemType is 'b': S_m = np.zeros(mesh.nF) S_e = np.zeros(mesh.nE) S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3 S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3 - Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e)) + Src.append(EM.FDEM.Src.RawVec([rx0], freq, S_m, mesh.getEdgeInnerProduct()*S_e)) elif fdemType is 'h' or fdemType is 'j': S_m = np.zeros(mesh.nE) S_e = np.zeros(mesh.nF) S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3 S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3 - Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e)) + Src.append(EM.FDEM.Src.RawVec([rx0], freq, mesh.getEdgeInnerProduct()*S_m, S_e)) if verbose: print ' Fetching %s problem' % (fdemType) if fdemType == 'e': survey = EM.FDEM.Survey(Src) - prb = EM.FDEM.Problem_e(mesh, mapping=mapping) + prb = EM.FDEM.Problem3D_e(mesh, mapping=mapping) elif fdemType == 'b': survey = EM.FDEM.Survey(Src) - prb = EM.FDEM.Problem_b(mesh, mapping=mapping) + prb = EM.FDEM.Problem3D_b(mesh, mapping=mapping) elif fdemType == 'j': survey = EM.FDEM.Survey(Src) - prb = EM.FDEM.Problem_j(mesh, mapping=mapping) + prb = EM.FDEM.Problem3D_j(mesh, mapping=mapping) elif fdemType == 'h': survey = EM.FDEM.Survey(Src) - prb = EM.FDEM.Problem_h(mesh, mapping=mapping) + prb = EM.FDEM.Problem3D_h(mesh, mapping=mapping) else: raise NotImplementedError() @@ -90,7 +95,7 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose) mesh = prb1.mesh print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp) - + logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY) mu = np.ones(mesh.nC)*MU diff --git a/SimPEG/EM/__init__.py b/SimPEG/EM/__init__.py index 565f63a8..e42afbc5 100644 --- a/SimPEG/EM/__init__.py +++ b/SimPEG/EM/__init__.py @@ -1,5 +1,6 @@ import TDEM import FDEM +import Static import Base import Analytics import Utils diff --git a/SimPEG/Examples/DC_Forward_PseudoSection.py b/SimPEG/Examples/DC_Forward_PseudoSection.py index 53467826..240cbb33 100644 --- a/SimPEG/Examples/DC_Forward_PseudoSection.py +++ b/SimPEG/Examples/DC_Forward_PseudoSection.py @@ -2,7 +2,7 @@ from SimPEG import Mesh, Utils, np, sp import SimPEG.DCIP as DC import time -def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', plotIt=True): +def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', unitType='appConductivity', plotIt=True): """ DC Forward Simulation ===================== @@ -15,14 +15,14 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p loc = Location of spheres [[x1,y1,z1],[x2,y2,z2]] radi = Radius of spheres [r1,r2] param = Conductivity of background and two spheres [m0,m1,m2] - stype = survey type "pdp" (pole dipole) or "dpdp" (dipole dipole) - dtype = Data type "appr" (app res) | "appc" (app cond) | "volt" (potential) + surveyType = survey type 'pole-dipole' or 'dipole-dipole' + unitType = Data type "appResistivity" | "appConductivity" | "volt" Created by @fourndo """ - assert stype in ['pdp', 'dpdp'], "Source type (stype) must be pdp or dpdp (pole dipole or dipole dipole)" - assert dtype in ['appr', 'appc', 'volt'], "Data type (dtype) must be appr (app res) or appc (app cond) or volt (potential)" + assert surveyType in ['pole-dipole', 'dipole-dipole'], "Source type (surveyType) must be pdp or dpdp (pole dipole or dipole dipole)" + assert unitType in ['appResistivity', 'appConductivity', 'volt'], "Unit type (unitType) must be appResistivity or appConductivity or volt (potential)" if loc is None: loc = np.c_[[-50.,0.,-50.],[50.,0.,-50.]] @@ -73,8 +73,8 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p 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]) - survey, Tx, Rx = DC.gen_DCIPsurvey(locs, mesh, stype, param[0], param[1], param[2]) + # [Tx, Rx] = DC.gen_DCIPsurvey(locs, mesh, surveyType, param[0], param[1], param[2]) + survey, Tx, Rx = DC.gen_DCIPsurvey(locs, mesh, surveyType, param[0], param[1], param[2]) # Define some global geometry dl_len = np.sqrt( np.sum((locs[0,:] - locs[1,:])**2) ) @@ -118,8 +118,8 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p rxloc_N = np.asarray(Rx[ii][:,3:]) - # For usual cases "dpdp" or "gradient" - if stype == 'pdp': + # For usual cases 'dipole-dipole' or "gradient" + if surveyType == 'pole-dipole': # Create an "inifinity" pole tx = np.squeeze(Tx[ii][:,0:1]) tinf = tx + np.array([dl_x,dl_y,0])*dl_len*2 @@ -157,12 +157,12 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p fig = plt.figure(figsize=(7,7)) ax = plt.subplot(2,1,1, aspect='equal') # Plot the location of the spheres for reference - circle1=plt.Circle((loc[0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3) - circle2=plt.Circle((loc[0,1],loc[2,1]),radi[1],color='k',fill=False, lw=3) + circle1=plt.Circle((loc[0,0], loc[2,0]), radi[0], color='w', fill=False, lw=3) + circle2=plt.Circle((loc[0,1], loc[2,1]), radi[1], color='k', fill=False, lw=3) ax.add_artist(circle1) ax.add_artist(circle2) - dat = mesh.plotSlice(np.log10(model), ax =ax, normal = 'Y', + dat = mesh.plotSlice(np.log10(model), ax = ax, normal = 'Y', ind = indy,grid=True, clim = np.log10([sig.min(),sig.max()])) ax.set_title('3-D model') @@ -188,15 +188,13 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p ax2 = plt.subplot(2,1,2, aspect='equal') # Plot the location of the spheres for reference - circle1=plt.Circle((loc[0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3) - circle2=plt.Circle((loc[0,1],loc[2,1]),radi[1],color='k',fill=False, lw=3) + circle1=plt.Circle((loc[0,0], loc[2,0]), radi[0], color='w', fill=False, lw=3) + circle2=plt.Circle((loc[0,1], loc[2,1]), radi[1], color='k', fill=False, lw=3) ax2.add_artist(circle1) ax2.add_artist(circle2) # Add the speudo section - dat = DC.plot_pseudoSection(survey2D,ax2,stype=stype, dtype = dtype) - - # plt.scatter(Tx2d[0][:],Tx[0][2,:],s=40,c='g', marker='v') + dat = DC.plot_pseudoSection(survey2D, ax2, surveyType=surveyType, unitType=unitType) # 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') ax2.set_title('Apparent Conductivity data') diff --git a/SimPEG/Examples/EM_FDEM_1D_Inversion.py b/SimPEG/Examples/EM_FDEM_1D_Inversion.py index 29f51ed4..4275903d 100644 --- a/SimPEG/Examples/EM_FDEM_1D_Inversion.py +++ b/SimPEG/Examples/EM_FDEM_1D_Inversion.py @@ -42,8 +42,8 @@ def run(plotIt=True): ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5) - rxOffset=10. - bzi = EM.FDEM.Rx(np.array([[rxOffset, 0., 1e-3]]), 'bzi') + rxOffset=10. + bzi = EM.FDEM.Rx.Point_b(np.array([[rxOffset, 0., 1e-3]]), orientation='z', component='imag') freqs = np.logspace(1,3,10) srcLoc = np.array([0., 0., 10.]) @@ -51,7 +51,7 @@ def run(plotIt=True): srcList = [EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z') for freq in freqs] survey = EM.FDEM.Survey(srcList) - prb = EM.FDEM.Problem_b(mesh, mapping=mapping) + prb = EM.FDEM.Problem3D_b(mesh, mapping=mapping) try: from pymatsolver import MumpsSolver diff --git a/SimPEG/Examples/EM_Schenkel_Morrison_Casing.py b/SimPEG/Examples/EM_Schenkel_Morrison_Casing.py index 76af4a3d..930b452a 100644 --- a/SimPEG/Examples/EM_Schenkel_Morrison_Casing.py +++ b/SimPEG/Examples/EM_Schenkel_Morrison_Casing.py @@ -215,7 +215,7 @@ def run(plotIt=True): # ------------ Problem and Survey --------------- survey = FDEM.Survey(sg_p + dg_p) mapping = [('sigma', Maps.IdentityMap(mesh))] - problem = FDEM.Problem_h(mesh, mapping=mapping) + problem = FDEM.Problem3D_h(mesh, mapping=mapping) problem.pair(survey) # ------------- Solve --------------------------- diff --git a/SimPEG/Examples/Utils_surface2ind_topo.py b/SimPEG/Examples/Utils_surface2ind_topo.py new file mode 100644 index 00000000..4ed1cdaf --- /dev/null +++ b/SimPEG/Examples/Utils_surface2ind_topo.py @@ -0,0 +1,41 @@ +from SimPEG import * +from SimPEG.Utils import surface2ind_topo + + +def run(plotIt=False, nx = 5, ny = 5): + """ + Here we show how to use :code:`Utils.surface2ind_topo` to identify cells below + a topographic surface. + + """ + + mesh = Mesh.TensorMesh([nx,ny], x0='CC') # 2D mesh + xtopo = np.linspace(mesh.gridN[:,0].min(), mesh.gridN[:,0].max()) + topo = 0.4*np.sin(xtopo*5) # define a topographic surface + + Topo = np.hstack([Utils.mkvc(xtopo,2),Utils.mkvc(topo,2)]) #make it an array + + indcc = surface2ind_topo(mesh, Topo,'CC') + + if plotIt: + from matplotlib.pylab import plt + from scipy.interpolate import interp1d + fig, ax = plt.subplots(1,1,figsize=(6,6)) + mesh.plotGrid(ax=ax, nodes=True, centers=True) + ax.plot(xtopo,topo,'k',linewidth=1) + # ax.plot(mesh.vectorNx, interp1d(xtopo,topo)(mesh.vectorNx),'--k',linewidth=3) + ax.plot(mesh.vectorCCx, interp1d(xtopo,topo)(mesh.vectorCCx),'--k',linewidth=3) + + + aveN2CC = Utils.sdiag(mesh.aveN2CC.T.sum(1))*mesh.aveN2CC.T + a = aveN2CC * indcc + a[a > 0] = 1. + a[a < 0.25] = np.nan + a = a.reshape(mesh.vnN, order='F') + masked_array = np.ma.array(a, mask=np.isnan(a)) + ax.pcolor(mesh.vectorNx,mesh.vectorNy,masked_array.T, cmap = plt.cm.gray,alpha=0.2) + plt.show() + + +if __name__ == '__main__': + run(plotIt=True) diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index c592e454..3c771f77 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -20,8 +20,9 @@ import Mesh_QuadTree_HangingNodes import Mesh_Tensor_Creation import MT_1D_ForwardAndInversion import MT_3D_Foward +import Utils_surface2ind_topo -__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_ForwardDC", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward"] +__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_ForwardDC", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "Utils_surface2ind_topo"] ##### AUTOIMPORTS ##### diff --git a/SimPEG/MT/BaseMT.py b/SimPEG/MT/BaseMT.py index c201dfb0..579f590f 100644 --- a/SimPEG/MT/BaseMT.py +++ b/SimPEG/MT/BaseMT.py @@ -1,5 +1,5 @@ from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np -from SimPEG.EM.FDEM.FDEM import BaseFDEMProblem +from SimPEG.EM.FDEM.ProblemFDEM import BaseFDEMProblem from SurveyMT import Survey, Data from FieldsMT import BaseMTFields diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index f17d2314..3e97499f 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -533,83 +533,6 @@ class ActiveCells(InjectActiveCells): FutureWarning) InjectActiveCells.__init__(self, mesh, indActive, valInactive, nC) -class InjectActiveCellsTopo(IdentityMap): - """ - Active model parameters. Extend for cells on topography to air cell (only works for tensor mesh) - - """ - - indActive = None #: Active Cells - valInactive = None #: Values of inactive Cells - nC = None #: Number of cells in the full model - - def __init__(self, mesh, indActive, nC=None): - self.mesh = mesh - - self.nC = nC or mesh.nC - - if indActive.dtype is not bool: - z = np.zeros(self.nC,dtype=bool) - z[indActive] = True - indActive = z - self.indActive = indActive - - self.indInactive = np.logical_not(indActive) - inds = np.nonzero(self.indActive)[0] - self.P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))), shape=(self.nC, self.nP)) - - @property - def shape(self): - return (self.nC, self.nP) - - @property - def nP(self): - """Number of parameters in the model.""" - return self.indActive.sum() - - def _transform(self, m): - val_temp = np.zeros(self.mesh.nC) - val_temp[self.indActive] = m - valInactive = np.zeros(self.mesh.nC) - #1D - if self.mesh.dim == 1: - z_temp = self.mesh.gridCC - val_temp[~self.indActive] = val_temp[np.argmax(z_temp[self.indActive])] - #2D - elif self.mesh.dim == 2: - act_temp = self.indActive.reshape((self.mesh.nCx, self.mesh.nCy), order = 'F') - val_temp = val_temp.reshape((self.mesh.nCx, self.mesh.nCy), order = 'F') - y_temp = self.mesh.gridCC[:,1].reshape((self.mesh.nCx, self.mesh.nCy), order = 'F') - for i in range(self.mesh.nCx): - act_tempx = act_temp[i,:] == 1 - val_temp[i,~act_tempx] = val_temp[i,np.argmax(y_temp[i,act_tempx])] - valInactive[~self.indActive] = Utils.mkvc(val_temp)[~self.indActive] - #3D - elif self.mesh.dim == 3: - act_temp = self.indActive.reshape((self.mesh.nCx*self.mesh.nCy, self.mesh.nCz), order = 'F') - val_temp = val_temp.reshape((self.mesh.nCx*self.mesh.nCy, self.mesh.nCz), order = 'F') - z_temp = self.mesh.gridCC[:,2].reshape((self.mesh.nCx*self.mesh.nCy, self.mesh.nCz), order = 'F') - for i in range(self.mesh.nCx*self.mesh.nCy): - act_tempxy = act_temp[i,:] == 1 - val_temp[i,~act_tempxy] = val_temp[i,np.argmax(z_temp[i,act_tempxy])] - valInactive[~self.indActive] = Utils.mkvc(val_temp)[~self.indActive] - - self.valInactive = valInactive - - return self.P*m + self.valInactive - - def inverse(self, D): - return self.P.T*D - - def deriv(self, m): - return self.P - -class ActiveCellsTopo(InjectActiveCellsTopo): - def __init__(self, mesh, indActive, valInactive, nC=None): - warnings.warn( - "`ActiveCellsTopo` is deprecated and will be removed in future versions. Use `InjectActiveCellsTopo` instead", - FutureWarning) - InjectActiveCellsTopo.__init__(self, mesh, indActive, valInactive, nC) class Weighting(IdentityMap): """ diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index ecdf36ac..f0b6751e 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -330,7 +330,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView): raise NotImplementedError('wrapping in the averaging is not yet implemented') return self._aveF2CCV - def getInterpolationMatCartMesh(self, Mrect, locType='CC'): + def getInterpolationMatCartMesh(self, Mrect, locType='CC', locTypeTo=None): """ Takes a cartesian mesh and returns a projection to translate onto the cartesian grid. """ @@ -338,19 +338,22 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView): assert self.isSymmetric, "Currently we have not taken into account other projections for more complicated CylMeshes" + if locTypeTo is None: + locTypeTo = locType + if locType == 'F': # do this three times for each component - X = self.getInterpolationMatCartMesh(Mrect, locType='Fx') - Y = self.getInterpolationMatCartMesh(Mrect, locType='Fy') - Z = self.getInterpolationMatCartMesh(Mrect, locType='Fz') + X = self.getInterpolationMatCartMesh(Mrect, locType='Fx', locTypeTo=locTypeTo+'x') + Y = self.getInterpolationMatCartMesh(Mrect, locType='Fy', locTypeTo=locTypeTo+'y') + Z = self.getInterpolationMatCartMesh(Mrect, locType='Fz', locTypeTo=locTypeTo+'z') return sp.vstack((X,Y,Z)) if locType == 'E': - X = self.getInterpolationMatCartMesh(Mrect, locType='Ex') - Y = self.getInterpolationMatCartMesh(Mrect, locType='Ey') - Z = spzeros(Mrect.nEz, self.nE) + X = self.getInterpolationMatCartMesh(Mrect, locType='Ex', locTypeTo=locTypeTo+'x') + Y = self.getInterpolationMatCartMesh(Mrect, locType='Ey', locTypeTo=locTypeTo+'y') + Z = spzeros(getattr(Mrect, 'n' + locTypeTo + 'z'), self.nE) return sp.vstack((X,Y,Z)) - grid = getattr(Mrect, 'grid' + locType) + grid = getattr(Mrect, 'grid' + locTypeTo) # This is unit circle stuff, 0 to 2*pi, starting at x-axis, rotating counter clockwise in an x-y slice theta = - np.arctan2(grid[:,0] - self.cartesianOrigin[0], grid[:,1] - self.cartesianOrigin[1]) + np.pi/2 theta[theta < 0] += np.pi*2.0 @@ -366,7 +369,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView): 'Ex': Mrect.tangents[:Mrect.nEx,:], 'Ey': Mrect.tangents[Mrect.nEx:(Mrect.nEx+Mrect.nEy),:], 'Ez': Mrect.tangents[-Mrect.nEz:,:], - }[locType] + }[locTypeTo] if 'F' in locType: normals = np.c_[np.cos(theta), np.sin(theta), np.zeros(theta.size)] proj = ( normals * dotMe ).sum(axis=1) diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index d0363001..b4a4b1ba 100644 --- a/SimPEG/Mesh/DiffOperators.py +++ b/SimPEG/Mesh/DiffOperators.py @@ -584,7 +584,67 @@ class DiffOperators(object): return Pbc, Pin, Pout + def getBCProjWF_simple(self, discretization='CC'): + """ + The weak form boundary condition projection matrices + when mixed boundary condition is used + + + """ + + if discretization is not 'CC': + raise NotImplementedError('Boundary conditions only implemented for CC discretization.') + + def projBC(n): + ij = ([0,n], [0,1]) + vals = [0,0] + vals[0] = 1 + vals[1] = 1 + return sp.csr_matrix((vals, ij), shape=(n+1,2)) + + def projDirichlet(n, bc): + bc = checkBC(bc) + ij = ([0,n], [0,1]) + vals = [0,0] + if(bc[0] == 'dirichlet'): + vals[0] = -1 + if(bc[1] == 'dirichlet'): + vals[1] = 1 + return sp.csr_matrix((vals, ij), shape=(n+1,2)) + + BC = [['dirichlet','dirichlet'],['dirichlet','dirichlet'],['dirichlet','dirichlet']] + n = self.vnC + indF = self.faceBoundaryInd + if(self.dim == 1): + Pbc = projDirichlet(n[0], BC[0]) + B = projBC(n[0]) + indF = indF[0] | indF[1] + Pbc = Pbc*sdiag(self.area[indF]) + + elif(self.dim == 2): + Pbc1 = sp.kron(speye(n[1]), projDirichlet(n[0], BC[0])) + Pbc2 = sp.kron(projDirichlet(n[1], BC[1]), speye(n[0])) + Pbc = sp.block_diag((Pbc1, Pbc2), format="csr") + B1 = sp.kron(speye(n[1]), projBC(n[0])) + B2 = sp.kron(projBC(n[1]), speye(n[0])) + B = sp.block_diag((B1, B2), format="csr") + indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3])] + Pbc = Pbc*sdiag(self.area[indF]) + + elif(self.dim == 3): + Pbc1 = kron3(speye(n[2]), speye(n[1]), projDirichlet(n[0], BC[0])) + Pbc2 = kron3(speye(n[2]), projDirichlet(n[1], BC[1]), speye(n[0])) + Pbc3 = kron3(projDirichlet(n[2], BC[2]), speye(n[1]), speye(n[0])) + Pbc = sp.block_diag((Pbc1, Pbc2, Pbc3), format="csr") + B1 = kron3(speye(n[2]), speye(n[1]), projBC(n[0])) + B2 = kron3(speye(n[2]), projBC(n[1]), speye(n[0])) + B3 = kron3(projBC(n[2]), speye(n[1]), speye(n[0])) + B = sp.block_diag((B1, B2, B3), format="csr") + indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3]), (indF[4] | indF[5])] + Pbc = Pbc*sdiag(self.area[indF]) + + return Pbc, B.T # --------------- Averaging --------------------- @property diff --git a/SimPEG/Mesh/View.py b/SimPEG/Mesh/View.py index 72225040..32e97539 100644 --- a/SimPEG/Mesh/View.py +++ b/SimPEG/Mesh/View.py @@ -218,7 +218,7 @@ class TensorView(object): return out viewOpts = ['real','imag','abs','vec'] normalOpts = ['X', 'Y', 'Z'] - vTypeOpts = ['CC', 'CCv','F','E','Fx','Fy','Fz','E','Ex','Ey','Ez'] + vTypeOpts = ['CC', 'CCv','N','F','E','Fx','Fy','Fz','E','Ex','Ey','Ez'] # Some user error checking assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts) diff --git a/SimPEG/PropMaps.py b/SimPEG/PropMaps.py index 527a6f7e..e4f65973 100644 --- a/SimPEG/PropMaps.py +++ b/SimPEG/PropMaps.py @@ -74,7 +74,7 @@ class Property(object): if linkedMap is None: return None linkMap = linkMapClass(None) * linkedMap - m = getattr(self, '%s'%linkName) + m = getattr(self, '%sModel'%linkName) return linkMap.deriv( m ) m = getattr(self, '%sModel'%prop.name) @@ -239,7 +239,7 @@ class PropMap(object): setattr(self, '%sMap'%name, mapping) setattr(self, '%sIndex'%name, slices.get(name, slice(nP, nP + mapping.nP))) nP += mapping.nP - self.nP = nP + self.nP = nP @property def defaultInvProp(self): diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index ed1039ec..fc101a61 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -311,6 +311,9 @@ class BaseRegularization(object): tmp = indActive indActive = np.zeros(mesh.nC, dtype=bool) indActive[tmp] = True + if indActive is not None and mapping is None: + mapping = Maps.IdentityMap(nP=indActive.nonzero()[0].size) + self.regmesh = RegularizationMesh(mesh,indActive) self.mapping = mapping or self.mapPair(mesh) self.mapping._assertMatchesPair(self.mapPair) diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 18c1994f..c7597fee 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -7,3 +7,4 @@ from CounterUtils import * import ModelBuilder import SolverUtils from coordutils import * +from modelutils import * diff --git a/SimPEG/Utils/modelutils.py b/SimPEG/Utils/modelutils.py new file mode 100644 index 00000000..dad92fae --- /dev/null +++ b/SimPEG/Utils/modelutils.py @@ -0,0 +1,63 @@ +from matutils import mkvc, ndgrid +import numpy as np + +def surface2ind_topo(mesh, topo, gridLoc='CC'): +# def genActiveindfromTopo(mesh, topo): + """ + Get active indices from topography + """ + + + if mesh.dim == 3: + from scipy.interpolate import NearestNDInterpolator + Ftopo = NearestNDInterpolator(topo[:,:2], topo[:,2]) + + if gridLoc == 'CC': + XY = ndgrid(mesh.vectorCCx, mesh.vectorCCy) + Zcc = mesh.gridCC[:,2].reshape((np.prod(mesh.vnC[:2]), mesh.nCz), order='F') + + gridTopo = Ftopo(XY) + actind = [gridTopo[ixy] <= Zcc[ixy,:] for ixy in range(np.prod(mesh.vnC[0]))] + actind = np.hstack(actind) + + elif gridLoc == 'N': + + XY = ndgrid(mesh.vectorNx, mesh.vectorNy) + gridTopo = Ftopo(XY).reshape(mesh.vnN[:2], order='F') + + if mesh._meshType not in ['TENSOR', 'CYL', 'BASETENSOR']: + raise NotImplementedError('Nodal surface2ind_topo not implemented for %s mesh'%mesh._meshType) + + Nz = mesh.vectorNz[1:] # TODO: this will only work for tensor meshes + actind = np.array([False]*mesh.nC).reshape(mesh.vnC, order='F') + + for ii in range(mesh.nCx): + for jj in range(mesh.nCy): + actind[ii,jj,:] = [np.all(gridTopo[ii:ii+2, jj:jj+2] >= Nz[kk]) for kk in range(len(Nz)) ] + + elif mesh.dim == 2: + from scipy.interpolate import interp1d + Ftopo = interp1d(topo[:,0], topo[:,1]) + + if gridLoc == 'CC': + gridTopo = Ftopo(mesh.gridCC[:,0]) + actind = mesh.gridCC[:,1] <= gridTopo + + elif gridLoc == 'N': + + gridTopo = Ftopo(mesh.vectorNx) + if mesh._meshType not in ['TENSOR', 'CYL', 'BASETENSOR']: + raise NotImplementedError('Nodal surface2ind_topo not implemented for %s mesh'%mesh._meshType) + + Ny = mesh.vectorNy[1:] # TODO: this will only work for tensor meshes + actind = np.array([False]*mesh.nC).reshape(mesh.vnC, order='F') + + for ii in range(mesh.nCx): + actind[ii,:] = [np.all(gridTopo[ii:ii+2] > Ny[kk]) for kk in range(len(Ny)) ] + + else: + raise NotImplementedError('surface2ind_topo not implemented for 1D mesh') + + return mkvc(actind) + + diff --git a/docs/examples/DC_Forward_PseudoSection.rst b/docs/examples/DC_Forward_PseudoSection.rst index 80cf0307..4231e944 100644 --- a/docs/examples/DC_Forward_PseudoSection.rst +++ b/docs/examples/DC_Forward_PseudoSection.rst @@ -20,9 +20,9 @@ INPUT: loc = Location of spheres [[x1,y1,z1],[x2,y2,z2]] radi = Radius of spheres [r1,r2] param = Conductivity of background and two spheres [m0,m1,m2] -stype = survey type "pdp" (pole dipole) or "dpdp" (dipole dipole) -dtype = Data type "appr" (app res) | "appc" (app cond) | "volt" (potential) -Created by @fourndo on Mon Feb 01 19:28:06 2016 +surveyType = survey type 'pole-dipole' or 'dipole-dipole' +unitType = Data type "appResistivity" | "appConductivity" | "volt" +Created by @fourndo diff --git a/docs/examples/Utils_surface2ind_topo.rst b/docs/examples/Utils_surface2ind_topo.rst new file mode 100644 index 00000000..04b27d5d --- /dev/null +++ b/docs/examples/Utils_surface2ind_topo.rst @@ -0,0 +1,24 @@ +.. _examples_Utils_surface2ind_topo: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +Here we show how to use :code:`Utils.surface2ind_topo` to identify cells below +a topographic surface. + + + +.. plot:: + + from SimPEG import Examples + Examples.Utils_surface2ind_topo.run() + +.. literalinclude:: ../../SimPEG/Examples/Utils_surface2ind_topo.py + :language: python + :linenos: diff --git a/setup.py b/setup.py index 3383c7f7..5308db47 100644 --- a/setup.py +++ b/setup.py @@ -5,16 +5,17 @@ SimPEG is a python package for simulation and gradient based parameter estimation in the context of geophysical applications. """ -import numpy as np - import os import sys import subprocess from distutils.core import setup +from distutils.command.build_ext import build_ext from setuptools import find_packages from distutils.extension import Extension + + CLASSIFIERS = [ 'Development Status :: 4 - Beta', 'Intended Audience :: Developers', @@ -51,11 +52,16 @@ if args.count("build_ext") > 0 and args.count("--inplace") == 0: try: from Cython.Build import cythonize from Cython.Distutils import build_ext - cythonKwargs = dict(cmdclass={'build_ext': build_ext}) USE_CYTHON = True except Exception, e: USE_CYTHON = False - cythonKwargs = dict() + +class NumpyBuild(build_ext): + def finalize_options(self): + build_ext.finalize_options(self) + __builtins__.__NUMPY_SETUP__ = False + import numpy + self.include_dirs.append(numpy.get_include()) ext = '.pyx' if USE_CYTHON else '.c' @@ -94,8 +100,8 @@ setup( classifiers=CLASSIFIERS, platforms = ["Windows", "Linux", "Solaris", "Mac OS-X", "Unix"], use_2to3 = False, - include_dirs=[np.get_include()], + cmdclass={'build_ext':NumpyBuild}, + setup_requires=['numpy'], ext_modules = extensions, scripts=scripts, - **cythonKwargs ) diff --git a/tests/base/test_PropMaps.py b/tests/base/test_PropMaps.py index ef22aaad..e012949f 100644 --- a/tests/base/test_PropMaps.py +++ b/tests/base/test_PropMaps.py @@ -1,6 +1,7 @@ import unittest from SimPEG import * from scipy.constants import mu_0 +from SimPEG import Tests class MyPropMap(Maps.PropMap): @@ -187,6 +188,34 @@ class TestPropMaps(unittest.TestCase): MyReciprocalPropMap([('sigma', iMap), ('mu', iMap)]) # This should be fine + def test_linked_derivs_sigma(self): + mesh = Mesh.TensorMesh([4,5], x0='CC') + + mapping = Maps.ExpMap(mesh) + propmap = MyReciprocalPropMap([('rho', mapping)]) + + x0 = np.random.rand(mesh.nC) + m = propmap(x0) + + # test Sigma + testme = lambda v: [1./(m.rhoMap*v), m.sigmaDeriv] + print 'Testing Rho from Sigma' + Tests.checkDerivative(testme, x0, dx=0.01*x0, num=5, plotIt=False) + + def test_linked_derivs_rho(self): + mesh = Mesh.TensorMesh([4,5], x0='CC') + + mapping = Maps.ExpMap(mesh) + propmap = MyReciprocalPropMap([('sigma', mapping)]) + + x0 = np.random.rand(mesh.nC) + m = propmap(x0) + + # test Sigma + testme = lambda v: [1./(m.sigmaMap*v), m.rhoDeriv] + print 'Testing Rho from Sigma' + Tests.checkDerivative(testme, x0, dx=0.01*x0, num=5, plotIt=False) + if __name__ == '__main__': unittest.main() diff --git a/tests/base/test_regularization.py b/tests/base/test_regularization.py index 5aa21ad5..6dddfb74 100644 --- a/tests/base/test_regularization.py +++ b/tests/base/test_regularization.py @@ -65,10 +65,8 @@ class RegularizationTests(unittest.TestCase): elif mesh.dim == 3: indActive = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5) - mapping = Maps.IdentityMap(nP=indActive.nonzero()[0].size) - for indAct in [indActive, indActive.nonzero()[0]]: # test both bool and integers - reg = r(mesh, mapping=mapping, indActive=indAct) + reg = r(mesh, indActive=indAct) m = np.random.rand(mesh.nC)[indAct] reg.mref = np.ones_like(m)*np.mean(m) diff --git a/tests/em/fdem/forward/test_FDEM_analytics.py b/tests/em/fdem/forward/test_FDEM_analytics.py index 9786e7c8..6f283666 100644 --- a/tests/em/fdem/forward/test_FDEM_analytics.py +++ b/tests/em/fdem/forward/test_FDEM_analytics.py @@ -28,12 +28,12 @@ class FDEM_analyticTests(unittest.TestCase): x = np.linspace(-10,10,5) XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0]) - rxList = EM.FDEM.Rx(XYZ, 'exi') + rxList = EM.FDEM.Rx.Point_e(XYZ, orientation='x', component='imag') Src0 = EM.FDEM.Src.MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq) survey = EM.FDEM.Survey([Src0]) - prb = EM.FDEM.Problem_b(mesh, mapping=mapping) + prb = EM.FDEM.Problem3D_b(mesh, mapping=mapping) prb.pair(survey) try: @@ -125,8 +125,8 @@ class FDEM_analyticTests(unittest.TestCase): mapping = [('sigma', Maps.IdentityMap(mesh)),('mu', Maps.IdentityMap(mesh))] - prbe = EM.FDEM.Problem_h(mesh, mapping=mapping) - prbm = EM.FDEM.Problem_e(mesh, mapping=mapping) + prbe = EM.FDEM.Problem3D_h(mesh, mapping=mapping) + prbm = EM.FDEM.Problem3D_e(mesh, mapping=mapping) prbe.pair(surveye) # pair problem and survey prbm.pair(surveym) diff --git a/tests/em/fdem/forward/test_FDEM_forwardHB.py b/tests/em/fdem/forward/test_FDEM_forwardHB.py index 545a5014..8fce615b 100644 --- a/tests/em/fdem/forward/test_FDEM_forwardHB.py +++ b/tests/em/fdem/forward/test_FDEM_forwardHB.py @@ -12,7 +12,7 @@ testBH = True verbose = False TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.) -#TODO: choose better testing parameters to lower this +#TODO: choose better testing parameters to lower this SrcList = ['RawVec', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop'] @@ -125,4 +125,4 @@ class FDEM_CrossCheck(unittest.TestCase): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzi', verbose=verbose, TOL=TOLEJHB)) if __name__ == '__main__': - unittest.main() \ No newline at end of file + unittest.main() diff --git a/tests/em/static/__init__.py b/tests/em/static/__init__.py new file mode 100644 index 00000000..420388ef --- /dev/null +++ b/tests/em/static/__init__.py @@ -0,0 +1,12 @@ +import os +import glob +import unittest + +if __name__ == '__main__': + test_file_strings = glob.glob('test_*.py') + module_strings = [str[0:len(str)-3] for str in test_file_strings] + suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str + in module_strings] + testSuite = unittest.TestSuite(suites) + + unittest.TextTestRunner(verbosity=2).run(testSuite) diff --git a/tests/em/static/test_DC_2D_analytic.py b/tests/em/static/test_DC_2D_analytic.py new file mode 100644 index 00000000..699bcbce --- /dev/null +++ b/tests/em/static/test_DC_2D_analytic.py @@ -0,0 +1,69 @@ +import unittest +from SimPEG import Mesh, Utils, EM, Maps, np +import SimPEG.EM.Static.DC as DC + +class DCProblemAnalyticTests(unittest.TestCase): + + def setUp(self): + + cs = 12.5 + hx = [(cs,7, -1.3),(cs,61),(cs,7, 1.3)] + hy = [(cs,7, -1.3),(cs,20)] + mesh = Mesh.TensorMesh([hx, hy],x0="CN") + sighalf = 1e-2 + sigma = np.ones(mesh.nC)*sighalf + x = np.linspace(-135, 250., 20) + M = Utils.ndgrid(x-12.5, np.r_[0.]) + N = Utils.ndgrid(x+12.5, np.r_[0.]) + A0loc = np.r_[-150, 0.] + A1loc = np.r_[-130, 0.] + rxloc = [np.c_[M, np.zeros(20)], np.c_[N, np.zeros(20)]] + data_anal = EM.Analytics.DCAnalyticHalf(np.r_[A0loc, 0.], rxloc, sighalf, earth_type="halfspace") + + rx = DC.Rx.Dipole_ky(M, N) + src0 = DC.Src.Pole([rx], A0loc) + survey = DC.Survey_ky([src0]) + + self.survey = survey + self.mesh = mesh + self.sigma = sigma + self.data_anal = data_anal + + try: + from pymatsolver import MumpsSolver + self.Solver = MumpsSolver + except ImportError, e: + self.Solver = SolverLU + + def test_Problem3D_N(self): + + problem = DC.Problem2D_N(self.mesh) + problem.Solver = self.Solver + problem.pair(self.survey) + data = self.survey.dpred(self.sigma) + err= np.linalg.norm((data-self.data_anal)/self.data_anal)**2 / self.data_anal.size + if err < 0.05: + passed = True + print ">> DC analytic test for Problem3D_N is passed" + else: + passed = False + print ">> DC analytic test for Problem3D_N is failed" + self.assertTrue(passed) + + def test_Problem3D_CC(self): + problem = DC.Problem2D_CC(self.mesh) + problem.Solver = self.Solver + problem.pair(self.survey) + data = self.survey.dpred(self.sigma) + err= np.linalg.norm((data-self.data_anal)/self.data_anal)**2 / self.data_anal.size + if err < 0.05: + passed = True + print ">> DC analytic test for Problem3D_CC is passed" + else: + passed = False + print ">> DC analytic test for Problem3D_CC is failed" + self.assertTrue(passed) + +if __name__ == '__main__': + unittest.main() + diff --git a/tests/em/static/test_DC_2D_jvecjtvecadj.py b/tests/em/static/test_DC_2D_jvecjtvecadj.py new file mode 100644 index 00000000..b96529cf --- /dev/null +++ b/tests/em/static/test_DC_2D_jvecjtvecadj.py @@ -0,0 +1,127 @@ +import unittest +from SimPEG import * +import SimPEG.EM.Static.DC as DC + + +class DCProblem_2DTestsCC(unittest.TestCase): + + def setUp(self): + + cs = 12.5 + hx = [(cs,7, -1.3),(cs,61),(cs,7, 1.3)] + hy = [(cs,7, -1.3),(cs,20)] + mesh = Mesh.TensorMesh([hx, hy],x0="CN") + x = np.linspace(-135, 250., 20) + M = Utils.ndgrid(x-12.5, np.r_[0.]) + N = Utils.ndgrid(x+12.5, np.r_[0.]) + A0loc = np.r_[-150, 0.] + A1loc = np.r_[-130, 0.] + rxloc = [np.c_[M, np.zeros(20)], np.c_[N, np.zeros(20)]] + rx = DC.Rx.Dipole_ky(M, N) + src0 = DC.Src.Pole([rx], A0loc) + src1 = DC.Src.Pole([rx], A1loc) + survey = DC.Survey_ky([src0, src1]) + problem = DC.Problem2D_CC(mesh, mapping=[('rho', Maps.IdentityMap(mesh))]) + problem.pair(survey) + + mSynth = np.ones(mesh.nC)*1. + survey.makeSyntheticData(mSynth) + + # Now set up the problem to do some minimization + dmis = DataMisfit.l2_DataMisfit(survey) + reg = Regularization.Tikhonov(mesh) + opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e0) + inv = Inversion.BaseInversion(invProb) + + self.inv = inv + self.reg = reg + self.p = problem + self.mesh = mesh + self.m0 = mSynth + self.survey = survey + self.dmis = dmis + + def test_misfit(self): + derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + + def test_adjoint(self): + # Adjoint Test + u = np.random.rand(self.mesh.nC*self.survey.nSrc) + v = np.random.rand(self.mesh.nC) + w = np.random.rand(self.survey.dobs.shape[0]) + wtJv = w.dot(self.p.Jvec(self.m0, v)) + vtJtw = v.dot(self.p.Jtvec(self.m0, w)) + passed = np.abs(wtJv - vtJtw) < 1e-10 + print 'Adjoint Test', np.abs(wtJv - vtJtw), passed + self.assertTrue(passed) + + def test_dataObj(self): + derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + +class DCProblemTestsN(unittest.TestCase): + + def setUp(self): + + cs = 12.5 + hx = [(cs,7, -1.3),(cs,61),(cs,7, 1.3)] + hy = [(cs,7, -1.3),(cs,20)] + mesh = Mesh.TensorMesh([hx, hy],x0="CN") + x = np.linspace(-135, 250., 20) + M = Utils.ndgrid(x-12.5, np.r_[0.]) + N = Utils.ndgrid(x+12.5, np.r_[0.]) + A0loc = np.r_[-150, 0.] + A1loc = np.r_[-130, 0.] + rxloc = [np.c_[M, np.zeros(20)], np.c_[N, np.zeros(20)]] + rx = DC.Rx.Dipole_ky(M, N) + src0 = DC.Src.Pole([rx], A0loc) + src1 = DC.Src.Pole([rx], A1loc) + survey = DC.Survey_ky([src0, src1]) + problem = DC.Problem2D_N(mesh, mapping=[('rho', Maps.IdentityMap(mesh))]) + problem.pair(survey) + + mSynth = np.ones(mesh.nC)*1. + survey.makeSyntheticData(mSynth) + + # Now set up the problem to do some minimization + dmis = DataMisfit.l2_DataMisfit(survey) + reg = Regularization.Tikhonov(mesh) + opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e0) + inv = Inversion.BaseInversion(invProb) + + self.inv = inv + self.reg = reg + self.p = problem + self.mesh = mesh + self.m0 = mSynth + self.survey = survey + self.dmis = dmis + + def test_misfit(self): + derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + + def test_adjoint(self): + # Adjoint Test + u = np.random.rand(self.mesh.nC*self.survey.nSrc) + v = np.random.rand(self.mesh.nC) + w = np.random.rand(self.survey.dobs.shape[0]) + wtJv = w.dot(self.p.Jvec(self.m0, v)) + vtJtw = v.dot(self.p.Jtvec(self.m0, w)) + passed = np.abs(wtJv - vtJtw) < 1e-8 + print 'Adjoint Test', np.abs(wtJv - vtJtw), passed + self.assertTrue(passed) + + def test_dataObj(self): + derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + +if __name__ == '__main__': + unittest.main() diff --git a/tests/em/static/test_DC_analytic.py b/tests/em/static/test_DC_analytic.py new file mode 100644 index 00000000..b8ebfc81 --- /dev/null +++ b/tests/em/static/test_DC_analytic.py @@ -0,0 +1,71 @@ +import unittest +from SimPEG import Mesh, Utils, EM, Maps, np +import SimPEG.EM.Static.DC as DC + +class DCProblemAnalyticTests(unittest.TestCase): + + def setUp(self): + + cs = 25. + hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)] + hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)] + hz = [(cs,7, -1.3),(cs,20)] + mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN") + sigma = np.ones(mesh.nC)*1e-2 + + x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)] + y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)] + Aloc = np.r_[-200., 0., 0.] + Bloc = np.r_[200., 0., 0.] + M = Utils.ndgrid(x-25.,y, np.r_[0.]) + N = Utils.ndgrid(x+25.,y, np.r_[0.]) + phiA = EM.Analytics.DCAnalyticHalf(Aloc, [M,N], 1e-2, earth_type="halfspace") + phiB = EM.Analytics.DCAnalyticHalf(Bloc, [M,N], 1e-2, earth_type="halfspace") + data_anal = phiA-phiB + + rx = DC.Rx.Dipole(M, N) + src = DC.Src.Dipole([rx], Aloc, Bloc) + survey = DC.Survey([src]) + + self.survey = survey + self.mesh = mesh + self.sigma = sigma + self.data_anal = data_anal + + try: + from pymatsolver import MumpsSolver + self.Solver = MumpsSolver + except ImportError, e: + self.Solver = SolverLU + + def test_Problem3D_N(self): + problem = DC.Problem3D_N(self.mesh) + problem.Solver = self.Solver + problem.pair(self.survey) + data = self.survey.dpred(self.sigma) + err= np.linalg.norm(data-self.data_anal)/np.linalg.norm(self.data_anal) + if err < 0.2: + passed = True + print ">> DC analytic test for Problem3D_N is passed" + else: + passed = False + print ">> DC analytic test for Problem3D_N is failed" + self.assertTrue(passed) + + def test_Problem3D_CC(self): + problem = DC.Problem3D_CC(self.mesh) + problem.Solver = self.Solver + problem.pair(self.survey) + data = self.survey.dpred(self.sigma) + err= np.linalg.norm(data-self.data_anal)/np.linalg.norm(self.data_anal) + if err < 0.2: + passed = True + print ">> DC analytic test for Problem3D_CC is passed" + else: + passed = False + print ">> DC analytic test for Problem3D_CC is failed" + self.assertTrue(passed) + +if __name__ == '__main__': + unittest.main() + diff --git a/tests/em/static/test_DC_jvecjtvecadj.py b/tests/em/static/test_DC_jvecjtvecadj.py new file mode 100644 index 00000000..227d4614 --- /dev/null +++ b/tests/em/static/test_DC_jvecjtvecadj.py @@ -0,0 +1,127 @@ +import unittest +from SimPEG import * +import SimPEG.EM.Static.DC as DC + + +class DCProblemTestsCC(unittest.TestCase): + + def setUp(self): + + aSpacing=2.5 + nElecs=5 + + surveySize = nElecs*aSpacing - aSpacing + cs = surveySize/nElecs/4 + + mesh = Mesh.TensorMesh([ + [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)], + [(cs,3, -1.3),(cs,3,1.3)], + # [(cs,5, -1.3),(cs,10)] + ],'CN') + + srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True) + survey = DC.Survey(srcList) + problem = DC.Problem3D_CC(mesh, mapping=[('rho', Maps.IdentityMap(mesh))]) + problem.pair(survey) + + mSynth = np.ones(mesh.nC) + survey.makeSyntheticData(mSynth) + + # Now set up the problem to do some minimization + dmis = DataMisfit.l2_DataMisfit(survey) + reg = Regularization.Tikhonov(mesh) + opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4) + inv = Inversion.BaseInversion(invProb) + + self.inv = inv + self.reg = reg + self.p = problem + self.mesh = mesh + self.m0 = mSynth + self.survey = survey + self.dmis = dmis + + def test_misfit(self): + derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + + def test_adjoint(self): + # Adjoint Test + u = np.random.rand(self.mesh.nC*self.survey.nSrc) + v = np.random.rand(self.mesh.nC) + w = np.random.rand(self.survey.dobs.shape[0]) + wtJv = w.dot(self.p.Jvec(self.m0, v)) + vtJtw = v.dot(self.p.Jtvec(self.m0, w)) + passed = np.abs(wtJv - vtJtw) < 1e-10 + print 'Adjoint Test', np.abs(wtJv - vtJtw), passed + self.assertTrue(passed) + + def test_dataObj(self): + derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + +class DCProblemTestsN(unittest.TestCase): + + def setUp(self): + + aSpacing=2.5 + nElecs=10 + + surveySize = nElecs*aSpacing - aSpacing + cs = surveySize/nElecs/4 + + mesh = Mesh.TensorMesh([ + [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)], + [(cs,3, -1.3),(cs,3,1.3)], + # [(cs,5, -1.3),(cs,10)] + ],'CN') + + srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True) + survey = DC.Survey(srcList) + problem = DC.Problem3D_N(mesh, mapping=[('rho', Maps.IdentityMap(mesh))]) + problem.pair(survey) + + mSynth = np.ones(mesh.nC) + survey.makeSyntheticData(mSynth) + + # Now set up the problem to do some minimization + dmis = DataMisfit.l2_DataMisfit(survey) + reg = Regularization.Tikhonov(mesh) + opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4) + inv = Inversion.BaseInversion(invProb) + + self.inv = inv + self.reg = reg + self.p = problem + self.mesh = mesh + self.m0 = mSynth + self.survey = survey + self.dmis = dmis + + def test_misfit(self): + derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) + self.assertTrue(passed) + + def test_adjoint(self): + # Adjoint Test + u = np.random.rand(self.mesh.nC*self.survey.nSrc) + v = np.random.rand(self.mesh.nC) + w = np.random.rand(self.survey.dobs.shape[0]) + wtJv = w.dot(self.p.Jvec(self.m0, v)) + vtJtw = v.dot(self.p.Jtvec(self.m0, w)) + passed = np.abs(wtJv - vtJtw) < 1e-8 + print 'Adjoint Test', np.abs(wtJv - vtJtw), passed + self.assertTrue(passed) + + def test_dataObj(self): + derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False) + self.assertTrue(passed) + +if __name__ == '__main__': + unittest.main() diff --git a/tests/em/static/test_IP_fwd.py b/tests/em/static/test_IP_fwd.py new file mode 100644 index 00000000..98bb7e90 --- /dev/null +++ b/tests/em/static/test_IP_fwd.py @@ -0,0 +1,96 @@ +import unittest +from SimPEG import Mesh, Utils, EM, Maps, np +import SimPEG.EM.Static.DC as DC +import SimPEG.EM.Static.IP as IP + +class IPProblemAnalyticTests(unittest.TestCase): + + def setUp(self): + + cs = 12.5 + npad=2 + hx = [(cs,npad, -1.3),(cs,21),(cs,npad, 1.3)] + hy = [(cs,npad, -1.3),(cs,21),(cs,npad, 1.3)] + hz = [(cs,npad, -1.3),(cs,20)] + mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN") + + x = mesh.vectorCCx[(mesh.vectorCCx>-80.)&(mesh.vectorCCx<80.)] + y = mesh.vectorCCx[(mesh.vectorCCy>-80.)&(mesh.vectorCCy<80.)] + Aloc = np.r_[-100., 0., 0.] + Bloc = np.r_[100., 0., 0.] + M = Utils.ndgrid(x-12.5,y, np.r_[0.]) + N = Utils.ndgrid(x+12.5,y, np.r_[0.]) + radius = 50. + xc = np.r_[0., 0., -100] + blkind = Utils.ModelBuilder.getIndicesSphere(xc, radius, mesh.gridCC) + sigmaInf = np.ones(mesh.nC)*1e-2 + eta = np.zeros(mesh.nC) + eta[blkind] = 0.1 + sigma0 = sigmaInf*(1.-eta) + + rx = DC.Rx.Dipole(M, N) + src = DC.Src.Dipole([rx], Aloc, Bloc) + surveyDC = DC.Survey([src]) + + self.surveyDC = surveyDC + self.mesh = mesh + self.sigmaInf = sigmaInf + self.sigma0 = sigma0 + self.src = src + self.eta = eta + + try: + from pymatsolver import MumpsSolver + self.Solver = MumpsSolver + except ImportError, e: + self.Solver = SolverLU + + def test_Problem3D_N(self): + + problemDC = DC.Problem3D_N(self.mesh) + problemDC.Solver = self.Solver + problemDC.pair(self.surveyDC) + data0 = self.surveyDC.dpred(self.sigma0) + finf = problemDC.fields(self.sigmaInf) + datainf = self.surveyDC.dpred(self.sigmaInf, f=finf) + problemIP = IP.Problem3D_N(self.mesh, sigma=self.sigmaInf, Ainv=problemDC.Ainv, f=finf) + problemIP.Solver = self.Solver + surveyIP = IP.Survey([self.src]) + problemIP.pair(surveyIP) + data_full = data0 - datainf + data = surveyIP.dpred(self.eta) + err= np.linalg.norm((data-data_full)/data_full)**2 / data_full.size + if err < 0.05: + passed = True + print ">> IP forward test for Problem3D_N is passed" + else: + passed = False + print ">> IP forward test for Problem3D_N is failed" + self.assertTrue(passed) + + def test_Problem3D_CC(self): + + problemDC = DC.Problem3D_CC(self.mesh) + problemDC.Solver = self.Solver + problemDC.pair(self.surveyDC) + data0 = self.surveyDC.dpred(self.sigma0) + finf = problemDC.fields(self.sigmaInf) + datainf = self.surveyDC.dpred(self.sigmaInf, f=finf) + problemIP = IP.Problem3D_CC(self.mesh, rho=1./self.sigmaInf, Ainv=problemDC.Ainv, f=finf) + problemIP.Solver = self.Solver + surveyIP = IP.Survey([self.src]) + problemIP.pair(surveyIP) + data_full = data0 - datainf + data = surveyIP.dpred(self.eta) + err= np.linalg.norm((data-data_full)/data_full)**2 / data_full.size + if err < 0.05: + passed = True + print ">> IP forward test for Problem3D_CC is passed" + else: + passed = False + print ">> IP forward test for Problem3D_CC is failed" + self.assertTrue(passed) + +if __name__ == '__main__': + unittest.main() + diff --git a/tests/em/static/test_IP_jvecjtvecadj.py b/tests/em/static/test_IP_jvecjtvecadj.py new file mode 100644 index 00000000..7a455784 --- /dev/null +++ b/tests/em/static/test_IP_jvecjtvecadj.py @@ -0,0 +1,126 @@ +import unittest +from SimPEG import * +import SimPEG.EM.Static.DC as DC +import SimPEG.EM.Static.IP as IP + + +class IPProblemTestsCC(unittest.TestCase): + + def setUp(self): + + aSpacing=2.5 + nElecs=5 + + surveySize = nElecs*aSpacing - aSpacing + cs = surveySize/nElecs/4 + + mesh = Mesh.TensorMesh([ + [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)], + [(cs,3, -1.3),(cs,3,1.3)], + # [(cs,5, -1.3),(cs,10)] + ],'CN') + + srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True) + survey = IP.Survey(srcList) + sigma = np.ones(mesh.nC) + problem = IP.Problem3D_CC(mesh, rho=1./sigma) + problem.pair(survey) + mSynth = np.ones(mesh.nC)*0.1 + survey.makeSyntheticData(mSynth) + # Now set up the problem to do some minimization + dmis = DataMisfit.l2_DataMisfit(survey) + reg = Regularization.Tikhonov(mesh) + opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4) + inv = Inversion.BaseInversion(invProb) + + self.inv = inv + self.reg = reg + self.p = problem + self.mesh = mesh + self.m0 = mSynth + self.survey = survey + self.dmis = dmis + + def test_misfit(self): + derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + + def test_adjoint(self): + # Adjoint Test + u = np.random.rand(self.mesh.nC*self.survey.nSrc) + v = np.random.rand(self.mesh.nC) + w = np.random.rand(self.survey.dobs.shape[0]) + wtJv = w.dot(self.p.Jvec(self.m0, v)) + vtJtw = v.dot(self.p.Jtvec(self.m0, w)) + passed = np.abs(wtJv - vtJtw) < 1e-10 + print 'Adjoint Test', np.abs(wtJv - vtJtw), passed + self.assertTrue(passed) + + def test_dataObj(self): + derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + +class IPProblemTestsN(unittest.TestCase): + + def setUp(self): + + aSpacing=2.5 + nElecs=5 + + surveySize = nElecs*aSpacing - aSpacing + cs = surveySize/nElecs/4 + + mesh = Mesh.TensorMesh([ + [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)], + [(cs,3, -1.3),(cs,3,1.3)], + # [(cs,5, -1.3),(cs,10)] + ],'CN') + + srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True) + survey = IP.Survey(srcList) + sigma = np.ones(mesh.nC) + problem = IP.Problem3D_N(mesh, sigma=sigma) + problem.pair(survey) + mSynth = np.ones(mesh.nC)*0.1 + survey.makeSyntheticData(mSynth) + # Now set up the problem to do some minimization + dmis = DataMisfit.l2_DataMisfit(survey) + reg = Regularization.Tikhonov(mesh) + opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4) + inv = Inversion.BaseInversion(invProb) + + self.inv = inv + self.reg = reg + self.p = problem + self.mesh = mesh + self.m0 = mSynth + self.survey = survey + self.dmis = dmis + + def test_misfit(self): + derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + + def test_adjoint(self): + # Adjoint Test + u = np.random.rand(self.mesh.nC*self.survey.nSrc) + v = np.random.rand(self.mesh.nC) + w = np.random.rand(self.survey.dobs.shape[0]) + wtJv = w.dot(self.p.Jvec(self.m0, v)) + vtJtw = v.dot(self.p.Jtvec(self.m0, w)) + passed = np.abs(wtJv - vtJtw) < 1e-8 + print 'Adjoint Test', np.abs(wtJv - vtJtw), passed + self.assertTrue(passed) + + def test_dataObj(self): + derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + +if __name__ == '__main__': + unittest.main() diff --git a/tests/em/static/test_SIP_jvecjtvecadj.py b/tests/em/static/test_SIP_jvecjtvecadj.py new file mode 100644 index 00000000..7cdb6def --- /dev/null +++ b/tests/em/static/test_SIP_jvecjtvecadj.py @@ -0,0 +1,232 @@ +import unittest +from SimPEG import * +import SimPEG +from SimPEG import Mesh, Utils, EM, Maps, np, Survey +from SimPEG.EM.Static import SIP, DC, IP +from pymatsolver import MumpsSolver + + +class IPProblemTestsCC(unittest.TestCase): + + def setUp(self): + + cs = 25. + hx = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)] + hy = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)] + hz = [(cs,0, -1.3),(cs,20)] + mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN") + blkind0 = Utils.ModelBuilder.getIndicesSphere(np.r_[-100., -100., -200.], 75., mesh.gridCC) + blkind1 = Utils.ModelBuilder.getIndicesSphere(np.r_[100., 100., -200.], 75., mesh.gridCC) + sigma = np.ones(mesh.nC)*1e-2 + eta = np.zeros(mesh.nC) + tau = np.ones_like(sigma)*1. + eta[blkind0] = 0.1 + eta[blkind1] = 0.1 + tau[blkind0] = 0.1 + tau[blkind1] = 0.01 + + x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)] + y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)] + Aloc = np.r_[-200., 0., 0.] + Bloc = np.r_[200., 0., 0.] + M = Utils.ndgrid(x-25.,y, np.r_[0.]) + N = Utils.ndgrid(x+25.,y, np.r_[0.]) + + times = np.arange(10)*1e-3 + 1e-3 + rx = SIP.Rx.Dipole(M, N, times) + src = SIP.Src.Dipole([rx], Aloc, Bloc) + survey = SIP.Survey([src]) + colemap = [("eta", Maps.IdentityMap(mesh)), ("taui", Maps.IdentityMap(mesh))] + problem = SIP.Problem3D_CC(mesh, rho=1./sigma, mapping=colemap) + problem.Solver = MumpsSolver + problem.pair(survey) + mSynth = np.r_[eta, 1./tau] + survey.makeSyntheticData(mSynth) + # Now set up the problem to do some minimization + dmis = DataMisfit.l2_DataMisfit(survey) + reg = Regularization.Tikhonov(mesh) + opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4) + inv = Inversion.BaseInversion(invProb) + + self.inv = inv + self.reg = reg + self.p = problem + self.mesh = mesh + self.m0 = mSynth + self.survey = survey + self.dmis = dmis + + def test_misfit(self): + derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + + def test_adjoint(self): + # Adjoint Test + u = np.random.rand(self.mesh.nC*self.survey.nSrc) + v = np.random.rand(self.mesh.nC*2) + w = np.random.rand(self.survey.dobs.shape[0]) + wtJv = w.dot(self.p.Jvec(self.m0, v)) + vtJtw = v.dot(self.p.Jtvec(self.m0, w)) + passed = np.abs(wtJv - vtJtw) < 1e-10 + print 'Adjoint Test', np.abs(wtJv - vtJtw), passed + self.assertTrue(passed) + + def test_dataObj(self): + derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + +class IPProblemTestsN(unittest.TestCase): + + def setUp(self): + + cs = 25. + hx = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)] + hy = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)] + hz = [(cs,0, -1.3),(cs,20)] + mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN") + blkind0 = Utils.ModelBuilder.getIndicesSphere(np.r_[-100., -100., -200.], 75., mesh.gridCC) + blkind1 = Utils.ModelBuilder.getIndicesSphere(np.r_[100., 100., -200.], 75., mesh.gridCC) + sigma = np.ones(mesh.nC)*1e-2 + eta = np.zeros(mesh.nC) + tau = np.ones_like(sigma)*1. + eta[blkind0] = 0.1 + eta[blkind1] = 0.1 + tau[blkind0] = 0.1 + tau[blkind1] = 0.01 + + x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)] + y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)] + Aloc = np.r_[-200., 0., 0.] + Bloc = np.r_[200., 0., 0.] + M = Utils.ndgrid(x-25.,y, np.r_[0.]) + N = Utils.ndgrid(x+25.,y, np.r_[0.]) + + times = np.arange(10)*1e-3 + 1e-3 + rx = SIP.Rx.Dipole(M, N, times) + src = SIP.Src.Dipole([rx], Aloc, Bloc) + survey = SIP.Survey([src]) + colemap = [("eta", Maps.IdentityMap(mesh)), ("taui", Maps.IdentityMap(mesh))] + problem = SIP.Problem3D_N(mesh, sigma=sigma, mapping=colemap) + problem.Solver = MumpsSolver + problem.pair(survey) + mSynth = np.r_[eta, 1./tau] + survey.makeSyntheticData(mSynth) + # Now set up the problem to do some minimization + dmis = DataMisfit.l2_DataMisfit(survey) + reg = Regularization.Tikhonov(mesh) + opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4) + inv = Inversion.BaseInversion(invProb) + + self.inv = inv + self.reg = reg + self.p = problem + self.mesh = mesh + self.m0 = mSynth + self.survey = survey + self.dmis = dmis + + def test_misfit(self): + derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + + def test_adjoint(self): + # Adjoint Test + u = np.random.rand(self.mesh.nC*self.survey.nSrc) + v = np.random.rand(self.mesh.nC*2) + w = np.random.rand(self.survey.dobs.shape[0]) + wtJv = w.dot(self.p.Jvec(self.m0, v)) + vtJtw = v.dot(self.p.Jtvec(self.m0, w)) + passed = np.abs(wtJv - vtJtw) < 1e-8 + print 'Adjoint Test', np.abs(wtJv - vtJtw), passed + self.assertTrue(passed) + + def test_dataObj(self): + derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + +class IPProblemTestsN_air(unittest.TestCase): + + def setUp(self): + + cs = 25. + hx = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)] + hy = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)] + hz = [(cs,0, -1.3),(cs,20),(cs,0, 1.3)] + mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCC") + blkind0 = Utils.ModelBuilder.getIndicesSphere(np.r_[-100., -100., -200.], 75., mesh.gridCC) + blkind1 = Utils.ModelBuilder.getIndicesSphere(np.r_[100., 100., -200.], 75., mesh.gridCC) + sigma = np.ones(mesh.nC)*1e-2 + airind = mesh.gridCC[:,2]>0. + sigma[airind] = 1e-8 + eta = np.zeros(mesh.nC) + tau = np.ones_like(sigma)*1. + eta[blkind0] = 0.1 + eta[blkind1] = 0.1 + tau[blkind0] = 0.1 + tau[blkind1] = 0.01 + + actmapeta = Maps.InjectActiveCells(mesh, ~airind, 0.) + actmaptau = Maps.InjectActiveCells(mesh, ~airind, 1.) + + x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)] + y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)] + Aloc = np.r_[-200., 0., 0.] + Bloc = np.r_[200., 0., 0.] + M = Utils.ndgrid(x-25.,y, np.r_[0.]) + N = Utils.ndgrid(x+25.,y, np.r_[0.]) + + times = np.arange(10)*1e-3 + 1e-3 + rx = SIP.Rx.Dipole(M, N, times) + src = SIP.Src.Dipole([rx], Aloc, Bloc) + survey = SIP.Survey([src]) + colemap = [("eta", Maps.IdentityMap(mesh)*actmapeta), ("taui", Maps.IdentityMap(mesh)*actmaptau)] + problem = SIP.Problem3D_N(mesh, sigma=sigma, mapping=colemap) + problem.Solver = MumpsSolver + problem.pair(survey) + mSynth = np.r_[eta[~airind], 1./tau[~airind]] + survey.makeSyntheticData(mSynth) + # Now set up the problem to do some minimization + dmis = DataMisfit.l2_DataMisfit(survey) + regmap = Maps.IdentityMap(nP=int(mSynth[~airind].size*2)) + reg = SIP.MultiRegularization(mesh, mapping=regmap, nModels=2, indActive=~airind) + opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4) + inv = Inversion.BaseInversion(invProb) + + self.inv = inv + self.reg = reg + self.p = problem + self.mesh = mesh + self.m0 = mSynth + self.survey = survey + self.dmis = dmis + + def test_misfit(self): + derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + + def test_adjoint(self): + # Adjoint Test + u = np.random.rand(self.mesh.nC*self.survey.nSrc) + v = np.random.rand(self.mesh.nC) + w = np.random.rand(self.survey.dobs.shape[0]) + wtJv = w.dot(self.p.Jvec(self.m0, v)) + vtJtw = v.dot(self.p.Jtvec(self.m0, w)) + passed = np.abs(wtJv - vtJtw) < 1e-8 + print 'Adjoint Test', np.abs(wtJv - vtJtw), passed + self.assertTrue(passed) + + def test_dataObj(self): + derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + +if __name__ == '__main__': + unittest.main() diff --git a/tests/mesh/test_Mixed_boundaryPoisson.py b/tests/mesh/test_Mixed_boundaryPoisson.py new file mode 100644 index 00000000..3aa1dbcd --- /dev/null +++ b/tests/mesh/test_Mixed_boundaryPoisson.py @@ -0,0 +1,411 @@ +import numpy as np +import scipy.sparse as sp +import unittest +import matplotlib.pyplot as plt +from SimPEG import * + +MESHTYPES = ['uniformTensorMesh'] + +def getxBCyBC_CC(mesh, alpha, beta, gamma): +# def getxBCyBC(mesh, alpha, beta, gamma): + """ + This is a subfunction generating mixed-boundary condition: + + .. math:: + + \nabla \cdot \vec{j} = -\nabla \cdot \vec{j}_s = q + + \rho \vec{j} = -\nabla \phi \phi + + \alpha \phi + \beta \frac{\partial \phi}{\partial r} = \gamma \ at \ r = \partial \Omega + + xBC = f_1(\alpha, \beta, \gamma) + yBC = f(\alpha, \beta, \gamma) + + Computes xBC and yBC for cell-centered discretizations + """ + if mesh.dim == 1: #1D + if (len(alpha) != 2 or len(beta) != 2 or len(gamma) != 2): + raise Exception("Lenght of list, alpha should be 2") + fCCxm,fCCxp = mesh.cellBoundaryInd + nBC = fCCxm.sum()+fCCxp.sum() + h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + + alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] + alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] + + # h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp] + h_xm, h_xp = mesh.hx[0], mesh.hx[-1] + + a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) + b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) + a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp) + b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp) + + xBC_xm = 0.5*a_xm + xBC_xp = 0.5*a_xp/b_xp + yBC_xm = 0.5*(1.-b_xm) + yBC_xp = 0.5*(1.-1./b_xp) + + xBC = np.r_[xBC_xm, xBC_xp] + yBC = np.r_[yBC_xm, yBC_xp] + + elif mesh.dim == 2: #2D + if (len(alpha) != 4 or len(beta) != 4 or len(gamma) != 4): + raise Exception("Lenght of list, alpha should be 4") + + fxm,fxp,fym,fyp = mesh.faceBoundaryInd + nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum() + + alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] + alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] + alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2] + alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3] + + # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] + # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + + h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp) + h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp) + + a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) + b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) + a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp) + b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp) + + a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym) + b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym) + a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp) + b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp) + + xBC_xm = 0.5*a_xm + xBC_xp = 0.5*a_xp/b_xp + yBC_xm = 0.5*(1.-b_xm) + yBC_xp = 0.5*(1.-1./b_xp) + xBC_ym = 0.5*a_ym + xBC_yp = 0.5*a_yp/b_yp + yBC_ym = 0.5*(1.-b_ym) + yBC_yp = 0.5*(1.-1./b_yp) + + sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]]) + sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]]) + + xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx] + xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy] + yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx] + yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy] + + xBC = np.r_[xBC_x, xBC_y] + yBC = np.r_[yBC_x, yBC_y] + + elif mesh.dim == 3: #3D + if (len(alpha) != 6 or len(beta) != 6 or len(gamma) != 6): + raise Exception("Lenght of list, alpha should be 6") + # fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd + fxm,fxp,fym,fyp,fzm,fzp = mesh.faceBoundaryInd + nBC = fxm.sum()+fxp.sum()+fxm.sum()+fxp.sum() + + alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0] + alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1] + alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2] + alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3] + alpha_zm, beta_zm, gamma_zm = alpha[4], beta[4], gamma[4] + alpha_zp, beta_zp, gamma_zp = alpha[5], beta[5], gamma[5] + + # h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] + # h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] + # h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2] + + h_xm, h_xp = mesh.hx[0]*np.ones_like(alpha_xm), mesh.hx[-1]*np.ones_like(alpha_xp) + h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp) + h_zm, h_zp = mesh.hz[0]*np.ones_like(alpha_zm), mesh.hz[-1]*np.ones_like(alpha_zp) + + a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) + b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm) + a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp) + b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp) + + a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym) + b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym) + a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp) + b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp) + + a_zm = gamma_zm/(0.5*alpha_zm-beta_zm/h_zm) + b_zm = (0.5*alpha_zm+beta_zm/h_zm)/(0.5*alpha_zm-beta_zm/h_zm) + a_zp = gamma_zp/(0.5*alpha_zp-beta_zp/h_zp) + b_zp = (0.5*alpha_zp+beta_zp/h_zp)/(0.5*alpha_zp-beta_zp/h_zp) + + xBC_xm = 0.5*a_xm + xBC_xp = 0.5*a_xp/b_xp + yBC_xm = 0.5*(1.-b_xm) + yBC_xp = 0.5*(1.-1./b_xp) + xBC_ym = 0.5*a_ym + xBC_yp = 0.5*a_yp/b_yp + yBC_ym = 0.5*(1.-b_ym) + yBC_yp = 0.5*(1.-1./b_yp) + xBC_zm = 0.5*a_zm + xBC_zp = 0.5*a_zp/b_zp + yBC_zm = 0.5*(1.-b_zm) + yBC_zp = 0.5*(1.-1./b_zp) + + sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]]) + sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]]) + sortindsfz = np.argsort(np.r_[np.arange(mesh.nFz)[fzm], np.arange(mesh.nFz)[fzp]]) + + xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx] + xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy] + xBC_z = np.r_[xBC_zm, xBC_zp][sortindsfz] + + yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx] + yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy] + yBC_z = np.r_[yBC_zm, yBC_zp][sortindsfz] + + xBC = np.r_[xBC_x, xBC_y, xBC_z] + yBC = np.r_[yBC_x, yBC_y, yBC_z] + + return xBC, yBC + +class Test1D_InhomogeneousMixed(Tests.OrderTest): + name = "1D - Mixed" + meshTypes = MESHTYPES + meshDimension = 1 + expectedOrders = 2 + meshSizes = [4, 8, 16, 32] + + def getError(self): + #Test function + phi_fun = lambda x: np.cos(np.pi*x) + j_fun = lambda x: np.pi*np.sin(np.pi*x) + phi_deriv = lambda x: -j_fun(x) + q_fun = lambda x: (np.pi**2)*np.cos(np.pi*x) + + xc_ana = phi_fun(self.M.gridCC) + q_ana = q_fun(self.M.gridCC) + j_ana = j_fun(self.M.gridFx) + + # Get boundary locations + vecN = self.M.vectorNx + vecC = self.M.vectorCCx + + # Setup Mixed B.C (alpha, beta, gamma) + alpha_xm, alpha_xp = 1., 1. + beta_xm, beta_xp = 1., 1. + alpha = np.r_[alpha_xm, alpha_xp] + beta = np.r_[beta_xm, beta_xp] + vecN = self.M.vectorNx + vecC = self.M.vectorCCx + phi_bc = phi_fun(vecN[[0,-1]]) + phi_deriv_bc = phi_deriv(vecN[[0,-1]]) + gamma = alpha*phi_bc + beta*phi_deriv_bc + x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma) + + + sigma = np.ones(self.M.nC) + Mfrho = self.M.getFaceInnerProduct(1./sigma) + MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True) + V = Utils.sdiag(self.M.vol) + Div = V*self.M.faceDiv + P_BC, B = self.M.getBCProjWF_simple() + q = q_fun(self.M.gridCC) + M = B*self.M.aveCC2F + G = Div.T - P_BC*Utils.sdiag(y_BC)*M + # Mrhoj = D.T V phi + P_BC*Utils.sdiag(y_BC)*M phi - P_BC*x_BC + rhs = V*q + Div*MfrhoI*P_BC*x_BC + A = Div*MfrhoI*G + + if self.myTest == 'xc': + #TODO: fix the null space + Ainv = Solver(A) + xc = Ainv*rhs + err = np.linalg.norm((xc-xc_ana), np.inf) + else: + NotImplementedError + return err + + + def test_order(self): + print "==== Testing Mixed boudary conduction for CC-problem ====" + self.name = "1D" + self.myTest = 'xc' + self.orderTest() + +class Test2D_InhomogeneousMixed(Tests.OrderTest): + name = "2D - Mixed" + meshTypes = MESHTYPES + meshDimension = 2 + expectedOrders = 2 + meshSizes = [4, 8, 16, 32] + + def getError(self): + #Test function + phi_fun = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1]) + j_funX = lambda x: +np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1]) + j_funY = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1]) + phideriv_funX = lambda x: -j_funX(x) + phideriv_funY = lambda x: -j_funY(x) + q_fun = lambda x: +2*(np.pi**2)*phi_fun(x) + + xc_ana = phi_fun(self.M.gridCC) + q_ana = q_fun(self.M.gridCC) + jX_ana = j_funX(self.M.gridFx) + jY_ana = j_funY(self.M.gridFy) + j_ana = np.r_[jX_ana,jY_ana] + + # Get boundary locations + fxm,fxp,fym,fyp = self.M.faceBoundaryInd + gBFxm = self.M.gridFx[fxm,:] + gBFxp = self.M.gridFx[fxp,:] + gBFym = self.M.gridFy[fym,:] + gBFyp = self.M.gridFy[fyp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + alpha_xm, alpha_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + beta_xm, beta_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + beta_ym, beta_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + + phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp) + phi_bc_ym, phi_bc_yp = phi_fun(gBFym), phi_fun(gBFyp) + + phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp) + phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp) + + gamma_fun = lambda alpha, beta, phi, phi_deriv: alpha*phi + beta*phi_deriv + gamma_xm = gamma_fun(alpha_xm, beta_xm, phi_bc_xm, phiderivX_bc_xm) + gamma_xp = gamma_fun(alpha_xp, beta_xp, phi_bc_xp, phiderivX_bc_xp) + gamma_ym = gamma_fun(alpha_ym, beta_ym, phi_bc_ym, phiderivY_bc_ym) + gamma_yp = gamma_fun(alpha_yp, beta_yp, phi_bc_yp, phiderivY_bc_yp) + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp] + + x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma) + + + sigma = np.ones(self.M.nC) + Mfrho = self.M.getFaceInnerProduct(1./sigma) + MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True) + V = Utils.sdiag(self.M.vol) + Div = V*self.M.faceDiv + P_BC, B = self.M.getBCProjWF_simple() + q = q_fun(self.M.gridCC) + M = B*self.M.aveCC2F + G = Div.T - P_BC*Utils.sdiag(y_BC)*M + rhs = V*q + Div*MfrhoI*P_BC*x_BC + A = Div*MfrhoI*G + + if self.myTest == 'xc': + Ainv = Solver(A) + xc = Ainv*rhs + err = np.linalg.norm((xc-xc_ana), np.inf) + else: + NotImplementedError + return err + + + def test_order(self): + print "==== Testing Mixed boudary conduction for CC-problem ====" + self.name = "2D" + self.myTest = 'xc' + self.orderTest() + +class Test3D_InhomogeneousMixed(Tests.OrderTest): + name = "3D - Mixed" + meshTypes = MESHTYPES + meshDimension = 3 + expectedOrders = 2 + meshSizes = [4, 8, 16] + + def getError(self): + #Test function + phi_fun = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.cos(np.pi*x[:,2]) + j_funX = lambda x: +np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.cos(np.pi*x[:,2]) + j_funY = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])*np.cos(np.pi*x[:,2]) + j_funZ = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.sin(np.pi*x[:,2]) + + phideriv_funX = lambda x: -j_funX(x) + phideriv_funY = lambda x: -j_funY(x) + phideriv_funZ = lambda x: -j_funZ(x) + + q_fun = lambda x: 3*(np.pi**2)*phi_fun(x) + + xc_ana = phi_fun(self.M.gridCC) + q_ana = q_fun(self.M.gridCC) + jX_ana = j_funX(self.M.gridFx) + jY_ana = j_funY(self.M.gridFy) + j_ana = np.r_[jX_ana,jY_ana,jY_ana] + + # Get boundary locations + fxm,fxp,fym,fyp,fzm,fzp = self.M.faceBoundaryInd + gBFxm = self.M.gridFx[fxm,:] + gBFxp = self.M.gridFx[fxp,:] + gBFym = self.M.gridFy[fym,:] + gBFyp = self.M.gridFy[fyp,:] + gBFzm = self.M.gridFz[fzm,:] + gBFzp = self.M.gridFz[fzp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + alpha_xm, alpha_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + beta_xm, beta_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + beta_ym, beta_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + alpha_zm, alpha_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2]) + beta_zm, beta_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2]) + + + phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp) + phi_bc_ym, phi_bc_yp = phi_fun(gBFym), phi_fun(gBFyp) + phi_bc_zm, phi_bc_zp = phi_fun(gBFzm), phi_fun(gBFzp) + + phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp) + phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp) + phiderivY_bc_zm, phiderivY_bc_zp = phideriv_funZ(gBFzm), phideriv_funZ(gBFzp) + + gamma_fun = lambda alpha, beta, phi, phi_deriv: alpha*phi + beta*phi_deriv + gamma_xm = gamma_fun(alpha_xm, beta_xm, phi_bc_xm, phiderivX_bc_xm) + gamma_xp = gamma_fun(alpha_xp, beta_xp, phi_bc_xp, phiderivX_bc_xp) + gamma_ym = gamma_fun(alpha_ym, beta_ym, phi_bc_ym, phiderivY_bc_ym) + gamma_yp = gamma_fun(alpha_yp, beta_yp, phi_bc_yp, phiderivY_bc_yp) + gamma_zm = gamma_fun(alpha_zm, beta_zm, phi_bc_zm, phiderivY_bc_zm) + gamma_zp = gamma_fun(alpha_zp, beta_zp, phi_bc_zp, phiderivY_bc_zp) + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp] + + x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma) + + + sigma = np.ones(self.M.nC) + Mfrho = self.M.getFaceInnerProduct(1./sigma) + MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True) + V = Utils.sdiag(self.M.vol) + Div = V*self.M.faceDiv + P_BC, B = self.M.getBCProjWF_simple() + q = q_fun(self.M.gridCC) + M = B*self.M.aveCC2F + G = Div.T - P_BC*Utils.sdiag(y_BC)*M + rhs = V*q + Div*MfrhoI*P_BC*x_BC + A = Div*MfrhoI*G + + if self.myTest == 'xc': + #TODO: fix the null space + Ainv = Solver(A) + xc = Ainv*rhs + err = np.linalg.norm((xc-xc_ana), np.inf) + else: + NotImplementedError + return err + + + def test_order(self): + print "==== Testing Mixed boudary conduction for CC-problem ====" + self.name = "3D" + self.myTest = 'xc' + self.orderTest() + + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/mesh/test_cylMesh.py b/tests/mesh/test_cylMesh.py index 5a963398..c7cd44ae 100644 --- a/tests/mesh/test_cylMesh.py +++ b/tests/mesh/test_cylMesh.py @@ -146,6 +146,20 @@ class TestCyl2DMesh(unittest.TestCase): assert np.abs(Pr*(Pc2r*mc) - Pc*mc).max() < 1e-3 + def test_getInterpMatCartMesh_Cells2Nodes(self): + + Mr = Mesh.TensorMesh([100,100,2], x0='CC0') + Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0]) + + mc = np.arange(Mc.nC) + xr = np.linspace(0,0.4,50) + xc = np.linspace(0,0.4,50) + 0.2 + Pr = Mr.getInterpolationMat(np.c_[xr,np.ones(50)*-0.2,np.ones(50)*0.5],'N') + Pc = Mc.getInterpolationMat(np.c_[xc,np.zeros(50),np.ones(50)*0.5],'CC') + Pc2r = Mc.getInterpolationMatCartMesh(Mr, 'CC', locTypeTo='N') + + assert np.abs(Pr*(Pc2r*mc) - Pc*mc).max() < 1e-3 + def test_getInterpMatCartMesh_Faces(self): Mr = Mesh.TensorMesh([100,100,2], x0='CC0') @@ -177,6 +191,37 @@ class TestCyl2DMesh(unittest.TestCase): assert np.abs(mag[dist > 0.1].min() - 1) < TOL + def test_getInterpMatCartMesh_Faces2Edges(self): + + Mr = Mesh.TensorMesh([100,100,2], x0='CC0') + Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0]) + + Pf2e = Mc.getInterpolationMatCartMesh(Mr, 'F', locTypeTo='E') + mf = np.ones(Mc.nF) + + ecart = Pf2e * mf + + excc = Mr.aveEx2CC*Mr.r(ecart, 'E', 'Ex') + eycc = Mr.aveEy2CC*Mr.r(ecart, 'E', 'Ey') + ezcc = Mr.r(ecart, 'E', 'Ez') + + indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5]) + indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5]) + + TOL = 1e-2 + assert np.abs(float(excc[indX]) - 1) < TOL + assert np.abs(float(excc[indY]) - 0) < TOL + assert np.abs(float(eycc[indX]) - 0) < TOL + assert np.abs(float(eycc[indY]) - 1) < TOL + assert np.abs((ezcc - 1).sum()) < TOL + + mag = (excc**2 + eycc**2)**0.5 + dist = ((Mr.gridCC[:,0] + 0.2)**2 + (Mr.gridCC[:,1] + 0.2)**2)**0.5 + + assert np.abs(mag[dist > 0.1].max() - 1) < TOL + assert np.abs(mag[dist > 0.1].min() - 1) < TOL + + def test_getInterpMatCartMesh_Edges(self): Mr = Mesh.TensorMesh([100,100,2], x0='CC0') @@ -185,11 +230,42 @@ class TestCyl2DMesh(unittest.TestCase): Pe = Mc.getInterpolationMatCartMesh(Mr, 'E') me = np.ones(Mc.nE) - erect = Pe * me + ecart = Pe * me - excc = Mr.aveEx2CC*Mr.r(erect, 'E', 'Ex') - eycc = Mr.aveEy2CC*Mr.r(erect, 'E', 'Ey') - ezcc = Mr.r(erect, 'E', 'Ez') + excc = Mr.aveEx2CC*Mr.r(ecart, 'E', 'Ex') + eycc = Mr.aveEy2CC*Mr.r(ecart, 'E', 'Ey') + ezcc = Mr.aveEz2CC*Mr.r(ecart, 'E', 'Ez') + + indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5]) + indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5]) + + TOL = 1e-2 + assert np.abs(float(excc[indX]) - 0) < TOL + assert np.abs(float(excc[indY]) + 1) < TOL + assert np.abs(float(eycc[indX]) - 1) < TOL + assert np.abs(float(eycc[indY]) - 0) < TOL + assert np.abs(ezcc.sum()) < TOL + + mag = (excc**2 + eycc**2)**0.5 + dist = ((Mr.gridCC[:,0] + 0.2)**2 + (Mr.gridCC[:,1] + 0.2)**2)**0.5 + + assert np.abs(mag[dist > 0.1].max() - 1) < TOL + assert np.abs(mag[dist > 0.1].min() - 1) < TOL + + + def test_getInterpMatCartMesh_Edges2Faces(self): + + Mr = Mesh.TensorMesh([100,100,2], x0='CC0') + Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0]) + + Pe2f = Mc.getInterpolationMatCartMesh(Mr, 'E', locTypeTo='F') + me = np.ones(Mc.nE) + + frect = Pe2f * me + + excc = Mr.aveFx2CC*Mr.r(frect, 'F', 'Fx') + eycc = Mr.aveFy2CC*Mr.r(frect, 'F', 'Fy') + ezcc = Mr.r(frect, 'F', 'Fz') indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5]) indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])