wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-06-27 18:25:42 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'em/dev' into em/ref/tdem

Browse Source 
# Conflicts:
#	SimPEG/EM/TDEM/SurveyTDEM.py
This commit is contained in:
Lindsey Heagy
2016-06-02 07:41:03 -07:00
parent 6f8f065643 382b31bd12
commit 9c7c89b8ba
65 changed files with 5600 additions and 818 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/DCIP/DCIPUtils.py
+157 -196
View File
@@ -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
SimPEG/Directives.py
+125 -59
View File
@@ -144,12 +144,18 @@ class BetaSchedule(InversionDirective):
if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
self.invProb.beta /= self.coolingFactor
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):
@@ -237,12 +243,6 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
# Save the file as a npz
np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
# class UpdateReferenceModel(Parameter):
# mref0 = None
# def nextIter(self):
# mref = getattr(self, 'm_prev', None)
# if mref is None:
# if self.debug: print 'UpdateReferenceModel is using mref0'
@@ -253,76 +253,158 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
class Update_IRLS(InversionDirective):
eps_min = None
eps_p = None
eps_q = None
norms = [2.,2.,2.,2.]
factor = None
gamma = None
phi_m_last = None
phi_d_last = None
f_old = None
f_min_change = 1e-2
beta_tol = 5e-2
# Solving parameter for IRLS (mode:2)
IRLSiter = 0
minGNiter = 5
maxIRLSiter = 10
iterStart = 0
# Beta schedule
coolingFactor = 2.
coolingRate = 1
mode = 1
@property
def target(self):
if getattr(self, '_target', None) is None:
self._target = self.survey.nD*0.5
return self._target
@target.setter
def target(self, val):
self._target = val
def initialize(self):
# Scale the regularization for changes in norm
if getattr(self, 'phi_m_last', None) is not None:
self.reg.curModel = self.invProb.curModel
self.reg.gamma = 1.
phim_new = self.reg.eval(self.invProb.curModel)
self.gamma = self.phi_m_last / phim_new
self.reg.curModel = self.invProb.curModel
self.reg.gamma = self.gamma
if getattr(self, 'phi_d_last', None) is None:
self.phi_d_last = self.invProb.phi_d
if self.mode == 1:
self.reg.norms = [2., 2., 2., 2.]
def endIter(self):
# Cool the threshold parameter if required
if getattr(self, 'factor', None) is not None:
eps = self.reg.eps / self.factor
if getattr(self, 'eps_min', None) is not None:
self.reg.eps = np.max([self.eps_min,eps])
# After reaching target misfit with l2-norm, switch to IRLS (mode:2)
if self.invProb.phi_d < self.target and self.mode == 1:
print "Convergence with smooth l2-norm regularization: Start IRLS steps..."
self.mode = 2
print self.eps_p, self.eps_q, self.norms
self.reg.eps_p = self.eps_p
self.reg.eps_q = self.eps_q
self.reg.norms = self.norms
self.coolingFactor = 1.
self.coolingRate = 1
self.iterStart = self.opt.iter
self.phi_d_last = self.invProb.phi_d
self.phi_m_last = self.invProb.phi_m_last
self.reg.l2model = self.invProb.curModel
self.reg.curModel = self.invProb.curModel
if getattr(self, 'f_old', None) is None:
self.f_old = self.reg.eval(self.invProb.curModel)#self.invProb.evalFunction(self.invProb.curModel, return_g=False, return_H=False)
# Beta Schedule
if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0:
if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
self.invProb.beta /= self.coolingFactor
# Only update after GN iterations
if (self.opt.iter-self.iterStart) % self.minGNiter == 0 and self.mode==2:
self.IRLSiter += 1
phim_new = self.reg.eval(self.invProb.curModel)
self.f_change = np.abs(self.f_old - phim_new) / self.f_old
print "Regularization decrease: %6.3e" % (self.f_change)
# Check for maximum number of IRLS cycles
if self.IRLSiter == self.maxIRLSiter:
print "Reach maximum number of IRLS cycles: %i" % self.maxIRLSiter
self.opt.stopNextIteration = True
return
# Check if the function has changed enough
if self.f_change < self.f_min_change and self.IRLSiter > 1:
print "Minimum decrease in regularization. End of IRLS"
self.opt.stopNextIteration = True
return
else:
self.reg.eps = eps
self.f_old = phim_new
# Get phi_m at the end of current iteration
self.phi_m_last = self.invProb.phi_m_last
# Cool the threshold parameter if required
if getattr(self, 'factor', None) is not None:
eps = self.reg.eps / self.factor
# Update the model used for the IRLS weights
self.reg.curModel = self.invProb.curModel
if getattr(self, 'eps_min', None) is not None:
self.reg.eps = np.max([self.eps_min,eps])
else:
self.reg.eps = eps
# Temporarely set gamma to 1. to get raw phi_m
self.reg.gamma = 1.
# Get phi_m at the end of current iteration
self.phi_m_last = self.invProb.phi_m_last
# Compute new model objective function value
phim_new = self.reg.eval(self.invProb.curModel)
# Reset the regularization matrices so that it is
# recalculated for current model
self.reg._Wsmall = None
self.reg._Wx = None
self.reg._Wy = None
self.reg._Wz = None
# Update gamma to scale the regularization between IRLS iterations
self.reg.gamma = self.phi_m_last / phim_new
# Update the model used for the IRLS weights
self.reg.curModel = self.invProb.curModel
# Set the weighting matrix to None so that it is recomputed next time
# it is called in the inversion
self.reg._W = None
# Temporarely set gamma to 1. to get raw phi_m
self.reg.gamma = 1.
# Compute new model objective function value
phim_new = self.reg.eval(self.invProb.curModel)
# Update gamma to scale the regularization between IRLS iterations
self.reg.gamma = self.phi_m_last / phim_new
# Reset the regularization matrices again for new gamma
self.reg._Wsmall = None
self.reg._Wx = None
self.reg._Wy = None
self.reg._Wz = None
# Check if misfit is within the tolerance, otherwise scale beta
val = self.invProb.phi_d / (self.survey.nD*0.5)
if np.abs(1.-val) > self.beta_tol:
self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d
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((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.
@@ -355,19 +437,3 @@ class Update_Wj(InversionDirective):
JtJdiag = JtJdiag / max(JtJdiag)
self.reg.wght = JtJdiag
class Scale_Beta(InversionDirective):
"""
Instead of a linear cooling schedule, beta is allowed to change based
on the ratio between the target misfit and the current data misfit. The
update is done only if the misfit is outside some threshold bounds.
"""
tol = 0.05
def endIter(self):
# Check if misfit is within the tolerance, otherwise adjust beta
val = self.invProb.phi_d / (self.survey.nD*0.5)
if np.abs(1.-val) > self.tol:
self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d
SimPEG/EM/Analytics/DC.py
+118
View File
@@ -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
SimPEG/EM/Analytics/__init__.py
+1
View File
@@ -1,3 +1,4 @@
from TDEM import hzAnalyticDipoleT
from FDEM import hzAnalyticDipoleF
from FDEMcasing import *
from DC import DCAnalyticHalf, DCAnalyticSphere
SimPEG/EM/Base.py
+15 -11
View File
@@ -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)
@@ -88,6 +89,11 @@ class BaseEMProblem(Problem.BaseProblem):
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
@@ -145,7 +151,6 @@ class BaseEMProblem(Problem.BaseProblem):
"""
return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) * self.curModel.sigmaDeriv
@property
def MeSigmaI(self):
"""
@@ -164,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):
@@ -183,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):
@@ -201,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):
@@ -210,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
@@ -220,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.')
SimPEG/EM/FDEM/FieldsFDEM.py
+1 -1
View File
@@ -181,7 +181,7 @@ class Fields3D_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
SimPEG/EM/FDEM/FDEM.py → SimPEG/EM/FDEM/ProblemFDEM.py
+5 -4
View File
@@ -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,9 +137,9 @@ class BaseFDEMProblem(BaseEMProblem):
df_dmT = df_dmT + du_dmT
# TODO: this should be taken care of by the reciever?
if rx.real_or_imag is 'real':
if rx.component is 'real':
Jtv += np.array(df_dmT, dtype=complex).real
elif rx.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')
@@ -166,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
SimPEG/EM/FDEM/RxFDEM.py
+24 -24
View File
@@ -7,15 +7,15 @@ class BaseRx(SimPEG.Survey.BaseRx):
:param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`)
:param string orientation: receiver orientation 'x', 'y' or 'z'
:param string real_or_imag: real or imaginary component 'real' or 'imag'
:param string component: real or imaginary component 'real' or 'imag'
"""
def __init__(self, locs, orientation=None, real_or_imag=None):
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(real_or_imag in ['real', 'imag']), "'real_or_imag' must be 'real' or 'imag', not %s"%real_or_imag
assert(component in ['real', 'imag']), "'component' must be 'real' or 'imag', not %s"%component
self.projComp = orientation
self.real_or_imag = real_or_imag
self.component = component
SimPEG.Survey.BaseRx.__init__(self, locs, rxType=None) #TODO: remove rxType from baseRx
@@ -36,7 +36,7 @@ class BaseRx(SimPEG.Survey.BaseRx):
P = self.getP(mesh, self.projGLoc(f))
f_part_complex = f[src, self.projField]
f_part = getattr(f_part_complex, self.real_or_imag) # get the real or imag component
f_part = getattr(f_part_complex, self.component) # get the real or imag component
return P*f_part
@@ -56,13 +56,13 @@ class BaseRx(SimPEG.Survey.BaseRx):
if not adjoint:
Pv_complex = P * v
Pv = getattr(Pv_complex, self.real_or_imag)
Pv = getattr(Pv_complex, self.component)
elif adjoint:
Pv_real = P.T * v
if self.real_or_imag == 'imag':
if self.component == 'imag':
Pv = 1j*Pv_real
elif self.real_or_imag == 'real':
elif self.component == 'real':
Pv = Pv_real.astype(complex)
else:
raise NotImplementedError('must be real or imag')
@@ -70,57 +70,57 @@ class BaseRx(SimPEG.Survey.BaseRx):
return Pv
class eField(BaseRx):
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 real_or_imag: real or imaginary component 'real' or 'imag'
:param string component: real or imaginary component 'real' or 'imag'
"""
def __init__(self, locs, orientation=None, real_or_imag=None):
def __init__(self, locs, orientation=None, component=None):
self.projField = 'e'
BaseRx.__init__(self, locs, orientation, real_or_imag)
super(Point_e, self).__init__(locs, orientation, component)
class bField(BaseRx):
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 real_or_imag: real or imaginary component 'real' or 'imag'
:param string component: real or imaginary component 'real' or 'imag'
"""
def __init__(self, locs, orientation=None, real_or_imag=None):
def __init__(self, locs, orientation=None, component=None):
self.projField = 'b'
BaseRx.__init__(self, locs, orientation, real_or_imag)
super(Point_b, self).__init__(locs, orientation, component)
class hField(BaseRx):
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 real_or_imag: real or imaginary component 'real' or 'imag'
:param string component: real or imaginary component 'real' or 'imag'
"""
def __init__(self, locs, orientation=None, real_or_imag=None):
def __init__(self, locs, orientation=None, component=None):
self.projField = 'h'
BaseRx.__init__(self, locs, orientation, real_or_imag)
super(Point_h, self).__init__(locs, orientation, component)
class jField(BaseRx):
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 real_or_imag: real or imaginary component 'real' or 'imag'
:param string component: real or imaginary component 'real' or 'imag'
"""
def __init__(self, locs, orientation=None, real_or_imag=None):
def __init__(self, locs, orientation=None, component=None):
self.projField = 'j'
BaseRx.__init__(self, locs, orientation, real_or_imag)
super(Point_j, self).__init__(locs, orientation, component)
SimPEG/EM/FDEM/SrcFDEM.py
+20 -9
View File
@@ -10,13 +10,16 @@ class BaseSrc(Survey.BaseSrc):
freq = None
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
@@ -53,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):
"""
@@ -63,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):
"""
@@ -73,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):
"""
@@ -83,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):
"""
@@ -141,11 +152,11 @@ class RawVec_e(BaseSrc):
:param bool integrate: Integrate the source term (multiply by Me) [False]
"""
def __init__(self, rxList, freq, s_e):
def __init__(self, rxList, freq, s_e, **kwargs):
self._s_e = np.array(s_e, dtype=complex)
self.freq = float(freq)
BaseSrc.__init__(self, rxList)
BaseSrc.__init__(self, rxList, **kwargs)
def s_e(self, prob):
"""
@@ -170,11 +181,11 @@ class RawVec_m(BaseSrc):
: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)
BaseSrc.__init__(self, rxList)
BaseSrc.__init__(self, rxList, **kwargs)
def s_m(self, prob):
"""
SimPEG/EM/FDEM/__init__.py
+1 -1
View File
@@ -1,5 +1,5 @@
from SurveyFDEM import Survey
import SrcFDEM as Src
import RxFDEM as Rx
from FDEM import Problem3D_e, Problem3D_b, Problem3D_j, Problem3D_h
from ProblemFDEM import Problem3D_e, Problem3D_b, Problem3D_j, Problem3D_h
from FieldsFDEM import Fields3D_e, Fields3D_b, Fields3D_j, Fields3D_h
SimPEG/EM/Static/DC/BoundaryUtils.py
+160
View File
@@ -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
SimPEG/EM/Static/DC/FieldsDC.py
+148
View File
@@ -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))
SimPEG/EM/Static/DC/FieldsDC_2D.py
+146
View File
@@ -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
SimPEG/EM/Static/DC/ProblemDC.py
+296
View File
@@ -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()
SimPEG/EM/Static/DC/ProblemDC_2D.py
+349
View File
@@ -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()
SimPEG/EM/Static/DC/RxDC.py
+129
View File
@@ -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
SimPEG/EM/Static/DC/SrcDC.py
+86
View File
@@ -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
SimPEG/EM/Static/DC/SurveyDC.py
+38
View File
@@ -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
SimPEG/EM/Static/DC/Utils.py
+38
View File
@@ -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
SimPEG/EM/Static/DC/__init__.py
+8
View File
@@ -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
SimPEG/EM/Static/IP/ProblemIP.py
+372
View File
@@ -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)
SimPEG/EM/Static/IP/SurveyIP.py
+23
View File
@@ -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)
SimPEG/EM/Static/IP/__init__.py
+2
View File
@@ -0,0 +1,2 @@
from ProblemIP import Problem3D_CC, Problem3D_N
from SurveyIP import Survey
SimPEG/EM/Static/SIP/ProblemSIP.py
+445
View File
@@ -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)
SimPEG/EM/Static/SIP/Regularization.py
+204
View File
@@ -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) ) )
SimPEG/EM/Static/SIP/RxSIP.py
+88
View File
@@ -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
SimPEG/EM/Static/SIP/SrcSIP.py
+64
View File
@@ -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
SimPEG/EM/Static/SIP/SurveySIP.py
+102
View File
@@ -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
SimPEG/EM/Static/SIP/__init__.py
+5
View File
@@ -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
SimPEG/EM/Static/Utils/StaticUtils.py
+317
View File
@@ -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
SimPEG/EM/Static/Utils/__init__.py
+1
View File
@@ -0,0 +1 @@
from StaticUtils import *
SimPEG/EM/Static/__init__.py
+3
View File
@@ -0,0 +1,3 @@
import DC
import IP
import SIP
SimPEG/EM/Utils/testingUtils.py
+1 -1
View File
@@ -26,7 +26,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
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 = getattr(EM.FDEM.Rx, comp[0] + 'Field')
Rx0 = getattr(EM.FDEM.Rx, 'Point_' + comp[0])
if comp[2] == 'r':
real_or_imag = 'real'
elif comp[2] == 'i':
SimPEG/EM/__init__.py
+1
View File
@@ -1,5 +1,6 @@
import TDEM
import FDEM
import Static
import Base
import Analytics
import Utils
SimPEG/Examples/DC_Forward_PseudoSection.py
+15 -17
View File
@@ -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')
SimPEG/Examples/EM_FDEM_1D_Inversion.py
+1 -1
View File
@@ -43,7 +43,7 @@ def run(plotIt=True):
rxOffset=10.
bzi = EM.FDEM.Rx.bField(np.array([[rxOffset, 0., 1e-3]]), orientation='z', real_or_imag='imag')
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.])
SimPEG/Examples/Inversion_IRLS.py
+36 -44
View File
@@ -1,7 +1,7 @@
from SimPEG import *
def run(N=200, plotIt=True):
def run(N=100, plotIt=True):
"""
Inversion: Linear Problem
=========================
@@ -18,6 +18,8 @@ def run(N=200, plotIt=True):
mesh = Mesh.TensorMesh([N])
m0 = np.ones(mesh.nC) * 1e-4
mref = np.zeros(mesh.nC)
nk = 10
jk = np.linspace(1.,nk,nk)
p = -2.
@@ -50,57 +52,47 @@ def run(N=200, plotIt=True):
wr = np.sum(prob.G**2.,axis=0)**0.5
wr = ( wr/np.max(wr) )
reg = Regularization.Simple(mesh)
reg.wght = wr
# reg = Regularization.Simple(mesh)
# reg.mref = mref
# reg.cell_weights = wr
#
dmis = DataMisfit.l2_DataMisfit(survey)
dmis.Wd = 1./wd
opt = Optimization.ProjectedGNCG(maxIter=30,lower=-2.,upper=2., maxIterCG= 20, tolCG = 1e-4)
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
invProb.curModel = m0
beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
target = Directives.TargetMisfit()
#
# opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, tolCG = 1e-4)
# invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
# invProb.curModel = m0
#
# beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
# target = Directives.TargetMisfit()
#
betaest = Directives.BetaEstimate_ByEig()
inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target])
mrec = inv.run(m0)
ml2 = mrec
print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
# Switch regularization to sparse
phim = invProb.phi_m_last
phid = invProb.phi_d
# inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target])
#
#
# mrec = inv.run(m0)
# ml2 = mrec
# print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
#
# # Switch regularization to sparse
# phim = invProb.phi_m_last
# phid = invProb.phi_d
reg = Regularization.Sparse(mesh)
reg.mref = mref
reg.cell_weights = wr
#==============================================================================
# fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
# dmdx = reg.mesh.cellDiffxStencil * mrec
# plt.plot(np.sort(dmdx))
#==============================================================================
#reg.recModel = mrec
reg.wght = np.ones(mesh.nC)
reg.mref = np.zeros(mesh.nC)
reg.eps_p = 5e-2
reg.eps_q = 1e-2
reg.norms = [0., 0., 2., 2.]
reg.wght = wr
eps_p = 5e-2
eps_q = 5e-2
norms = [0., 0., 2., 2.]
opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 20, tolCG = 1e-3)
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.)
beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1)
#betaest = Directives.BetaEstimate_ByEig()
target = Directives.TargetMisfit()
IRLS =Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid )
opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3)
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
update_Jacobi = Directives.Update_lin_PreCond()
IRLS = Directives.Update_IRLS( norms=norms, eps_p=eps_p, eps_q=eps_q)
inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS])
m0 = mrec
inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest,update_Jacobi])
# Run inversion
mrec = inv.run(m0)
@@ -117,7 +109,7 @@ def run(N=200, plotIt=True):
axes[0].set_title('Columns of matrix G')
axes[1].plot(mesh.vectorCCx, mtrue, 'b-')
axes[1].plot(mesh.vectorCCx, ml2, 'r-')
axes[1].plot(mesh.vectorCCx, reg.l2model, 'r-')
#axes[1].legend(('True Model', 'Recovered Model'))
axes[1].set_ylim(-1.0,1.25)
SimPEG/Examples/Forward_BasicDirectCurrent.py → SimPEG/Examples/Mesh_Basic_ForwardDC.py
+12 -31
View File
@@ -1,22 +1,25 @@
from SimPEG import Mesh, Utils, np, SolverLU
## 2D DC forward modeling example with Tensor and Curvilinear Meshes
def run(plotIt=True):
"""
Mesh: Basic Forward 2D DC Resistivity
=====================================
2D DC forward modeling example with Tensor and Curvilinear Meshes
"""
# Step1: Generate Tensor and Curvilinear Mesh
sz = [40,40]
# Tensor Mesh
tM = Mesh.TensorMesh(sz)
# Curvilinear Mesh
rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
# Step2: Direct Current (DC) operator
def DCfun(mesh, pts):
D = mesh.faceDiv
G = D.T
sigma = 1e-2*np.ones(mesh.nC)
Msigi = mesh.getFaceInnerProduct(1./sigma)
MsigI = Utils.sdInv(Msigi)
A = D*MsigI*G
MsigI = mesh.getFaceInnerProduct(sigma, invProp=True, invMat=True)
A = -D*MsigI*D.T
A[-1,-1] /= mesh.vol[-1] # Remove null space
rhs = np.zeros(mesh.nC)
txind = Utils.meshutils.closestPoints(mesh, pts)
@@ -37,39 +40,17 @@ def run(plotIt=True):
if not plotIt: return
import matplotlib.pyplot as plt
import matplotlib
from matplotlib.mlab import griddata
#Step4: Making Figure
fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
label = ["(a)", "(b)"]
opts = {}
vmin, vmax = phitM.min(), phitM.max()
dat = tM.plotImage(phitM, ax=axes[0], clim=(vmin, vmax), grid=True)
#TODO: At the moment Curvilinear Mesh do not have plotimage
Xi = tM.gridCC[:,0].reshape(sz[0], sz[1], order='F')
Yi = tM.gridCC[:,1].reshape(sz[0], sz[1], order='F')
PHIrM = griddata(rM.gridCC[:,0], rM.gridCC[:,1], phirM, Xi, Yi, interp='linear')
axes[1].contourf(Xi, Yi, PHIrM, 100, vmin=vmin, vmax=vmax)
dat = rM.plotImage(phirM, ax=axes[1], clim=(vmin, vmax), grid=True)
cb = plt.colorbar(dat[0], ax=axes[0]); cb.set_label("Voltage (V)")
cb = plt.colorbar(dat[0], ax=axes[1]); cb.set_label("Voltage (V)")
tM.plotGrid(ax=axes[0], **opts)
axes[0].set_title('TensorMesh')
rM.plotGrid(ax=axes[1], **opts)
axes[1].set_title('CurvilinearMesh')
for i in range(2):
axes[i].set_xlim(0.025, 0.975)
axes[i].set_ylim(0.025, 0.975)
axes[i].text(0., 1.0, label[i], fontsize=20)
if i==0:
axes[i].set_ylabel("y")
else:
axes[i].set_ylabel(" ")
axes[i].set_xlabel("x")
plt.show()
SimPEG/Examples/Utils_surface2ind_topo.py
+41
View File
@@ -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)
SimPEG/Examples/__init__.py
+3 -2
View File
@@ -8,9 +8,9 @@ import EM_FDEM_Analytic_MagDipoleWholespace
import EM_Schenkel_Morrison_Casing
import EM_TDEM_1D_Inversion
import FLOW_Richards_1D_Celia1990
import Forward_BasicDirectCurrent
import Inversion_IRLS
import Inversion_Linear
import Mesh_Basic_ForwardDC
import Mesh_Basic_PlotImage
import Mesh_Basic_Types
import Mesh_Operators_CahnHilliard
@@ -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", "Forward_BasicDirectCurrent", "Inversion_IRLS", "Inversion_Linear", "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 #####
SimPEG/MT/BaseMT.py
+1 -1
View File
@@ -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
SimPEG/Maps.py
-77
View File
@@ -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):
"""
SimPEG/Mesh/CurvilinearMesh.py
+2 -97
View File
@@ -2,6 +2,7 @@ from SimPEG import Utils, np
from BaseMesh import BaseRectangularMesh
from DiffOperators import DiffOperators
from InnerProducts import InnerProducts
from View import CurvView
# Some helper functions.
length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5
@@ -10,7 +11,7 @@ normalize2D = lambda x: x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2))
normalize3D = lambda x: x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2))
class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvView):
"""
CurvilinearMesh is a mesh class that deals with curvilinear meshes.
@@ -330,102 +331,6 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
#############################################
# Plotting Functions #
#############################################
def plotGrid(self, ax=None, nodes=False, faces=False, centers=False, edges=False, lines=True, showIt=False):
"""Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.
.. plot::
:include-source:
from SimPEG import Mesh, Utils
X, Y = Utils.exampleLrmGrid([3,3],'rotate')
M = Mesh.CurvilinearMesh([X, Y])
M.plotGrid(showIt=True)
"""
import matplotlib.pyplot as plt
import matplotlib
from mpl_toolkits.mplot3d import Axes3D
mkvc = Utils.mkvc
axOpts = {'projection':'3d'} if self.dim == 3 else {}
if ax is None: ax = plt.subplot(111, **axOpts)
NN = self.r(self.gridN, 'N', 'N', 'M')
if self.dim == 2:
if lines:
X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()
X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()
X = np.r_[X1, X2]
Y = np.r_[Y1, Y2]
ax.plot(X, Y, 'b-')
if centers:
ax.plot(self.gridCC[:,0],self.gridCC[:,1],'ro')
# Nx = self.r(self.normals, 'F', 'Fx', 'V')
# Ny = self.r(self.normals, 'F', 'Fy', 'V')
# Tx = self.r(self.tangents, 'E', 'Ex', 'V')
# Ty = self.r(self.tangents, 'E', 'Ey', 'V')
# ax.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
# nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
# nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
# ax.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
# ax.plot(nX, nY, 'r-')
# nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
# nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
# #ax.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
# ax.plot(nX, nY, 'g-')
# tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
# tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
# ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
# ax.plot(tX, tY, 'r-')
# nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
# nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
# #ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
# ax.plot(nX, nY, 'g-')
elif self.dim == 3:
X1 = np.c_[mkvc(NN[0][:-1, :, :]), mkvc(NN[0][1:, :, :]), mkvc(NN[0][:-1, :, :])*np.nan].flatten()
Y1 = np.c_[mkvc(NN[1][:-1, :, :]), mkvc(NN[1][1:, :, :]), mkvc(NN[1][:-1, :, :])*np.nan].flatten()
Z1 = np.c_[mkvc(NN[2][:-1, :, :]), mkvc(NN[2][1:, :, :]), mkvc(NN[2][:-1, :, :])*np.nan].flatten()
X2 = np.c_[mkvc(NN[0][:, :-1, :]), mkvc(NN[0][:, 1:, :]), mkvc(NN[0][:, :-1, :])*np.nan].flatten()
Y2 = np.c_[mkvc(NN[1][:, :-1, :]), mkvc(NN[1][:, 1:, :]), mkvc(NN[1][:, :-1, :])*np.nan].flatten()
Z2 = np.c_[mkvc(NN[2][:, :-1, :]), mkvc(NN[2][:, 1:, :]), mkvc(NN[2][:, :-1, :])*np.nan].flatten()
X3 = np.c_[mkvc(NN[0][:, :, :-1]), mkvc(NN[0][:, :, 1:]), mkvc(NN[0][:, :, :-1])*np.nan].flatten()
Y3 = np.c_[mkvc(NN[1][:, :, :-1]), mkvc(NN[1][:, :, 1:]), mkvc(NN[1][:, :, :-1])*np.nan].flatten()
Z3 = np.c_[mkvc(NN[2][:, :, :-1]), mkvc(NN[2][:, :, 1:]), mkvc(NN[2][:, :, :-1])*np.nan].flatten()
X = np.r_[X1, X2, X3]
Y = np.r_[Y1, Y2, Y3]
Z = np.r_[Z1, Z2, Z3]
ax.plot(X, Y, 'b', zs=Z)
ax.set_zlabel('x3')
ax.grid(True)
ax.set_xlabel('x1')
ax.set_ylabel('x2')
if showIt: plt.show()
if __name__ == '__main__':
nc = 5
h1 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])
SimPEG/Mesh/DiffOperators.py
+60
View File
@@ -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
SimPEG/Mesh/View.py
+79 -41
View File
@@ -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)
@@ -552,7 +552,8 @@ class CurvView(object):
def __init__(self):
pass
def plotGrid(self, length=0.05, showIt=False):
def plotGrid(self, ax=None, nodes=False, faces=False, centers=False, edges=False, lines=True, showIt=False):
"""Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.
@@ -560,60 +561,63 @@ class CurvView(object):
:include-source:
from SimPEG import Mesh, Utils
X, Y = Utils.exampleCurvGird([3,3],'rotate')
X, Y = Utils.exampleLrmGrid([3,3],'rotate')
M = Mesh.CurvilinearMesh([X, Y])
M.plotGrid(showIt=True)
"""
import matplotlib.pyplot as plt
import matplotlib
from mpl_toolkits.mplot3d import Axes3D
axOpts = {'projection':'3d'} if self.dim == 3 else {}
if ax is None: ax = plt.subplot(111, **axOpts)
NN = self.r(self.gridN, 'N', 'N', 'M')
if self.dim == 2:
fig = plt.figure(2)
fig.clf()
ax = plt.subplot(111)
X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()
X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()
if lines:
X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()
X = np.r_[X1, X2]
Y = np.r_[Y1, Y2]
X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()
plt.plot(X, Y)
X = np.r_[X1, X2]
Y = np.r_[Y1, Y2]
plt.hold(True)
Nx = self.r(self.normals, 'F', 'Fx', 'V')
Ny = self.r(self.normals, 'F', 'Fy', 'V')
Tx = self.r(self.tangents, 'E', 'Ex', 'V')
Ty = self.r(self.tangents, 'E', 'Ey', 'V')
ax.plot(X, Y, 'b-')
if centers:
ax.plot(self.gridCC[:,0],self.gridCC[:,1],'ro')
plt.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
# Nx = self.r(self.normals, 'F', 'Fx', 'V')
# Ny = self.r(self.normals, 'F', 'Fy', 'V')
# Tx = self.r(self.tangents, 'E', 'Ex', 'V')
# Ty = self.r(self.tangents, 'E', 'Ey', 'V')
nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
plt.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
plt.plot(nX, nY, 'r-')
# ax.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
#plt.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
plt.plot(nX, nY, 'g-')
# nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
# nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
# ax.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
# ax.plot(nX, nY, 'r-')
tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
plt.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
plt.plot(tX, tY, 'r-')
# nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
# nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
# #ax.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
# ax.plot(nX, nY, 'g-')
nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
#plt.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
plt.plot(nX, nY, 'g-')
plt.axis('equal')
# tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
# tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
# ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
# ax.plot(tX, tY, 'r-')
# nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
# nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
# #ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
# ax.plot(nX, nY, 'g-')
elif self.dim == 3:
fig = plt.figure(3)
fig.clf()
ax = fig.add_subplot(111, projection='3d')
X1 = np.c_[mkvc(NN[0][:-1, :, :]), mkvc(NN[0][1:, :, :]), mkvc(NN[0][:-1, :, :])*np.nan].flatten()
Y1 = np.c_[mkvc(NN[1][:-1, :, :]), mkvc(NN[1][1:, :, :]), mkvc(NN[1][:-1, :, :])*np.nan].flatten()
Z1 = np.c_[mkvc(NN[2][:-1, :, :]), mkvc(NN[2][1:, :, :]), mkvc(NN[2][:-1, :, :])*np.nan].flatten()
@@ -630,16 +634,50 @@ class CurvView(object):
Y = np.r_[Y1, Y2, Y3]
Z = np.r_[Z1, Z2, Z3]
plt.plot(X, Y, 'b', zs=Z)
ax.plot(X, Y, 'b', zs=Z)
ax.set_zlabel('x3')
ax.grid(True)
ax.hold(False)
ax.set_xlabel('x1')
ax.set_ylabel('x2')
if showIt: plt.show()
def plotImage(self, I, ax=None, showIt=False, grid=False, clim=None):
if self.dim == 3: raise NotImplementedError('This is not yet done!')
import matplotlib.pyplot as plt
import matplotlib
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.colors as colors
import matplotlib.cm as cmx
if ax is None: ax = plt.subplot(111)
jet = cm = plt.get_cmap('jet')
cNorm = colors.Normalize(
vmin=I.min() if clim is None else clim[0],
vmax=I.max() if clim is None else clim[1])
scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet)
# ax.set_xlim((self.x0[0], self.h[0].sum()))
# ax.set_ylim((self.x0[1], self.h[1].sum()))
Nx = self.r(self.gridN[:,0],'N','N','M')
Ny = self.r(self.gridN[:,1],'N','N','M')
cell = self.r(I,'CC','CC','M')
for ii in range(self.nCx):
for jj in range(self.nCy):
I = [ii,ii+1,ii+1,ii]
J = [jj,jj,jj+1,jj+1]
ax.add_patch(plt.Polygon(np.c_[Nx[I,J],Ny[I,J]], facecolor=scalarMap.to_rgba(cell[ii,jj]), edgecolor='k' if grid else 'none'))
scalarMap._A = [] # http://stackoverflow.com/questions/8342549/matplotlib-add-colorbar-to-a-sequence-of-line-plots
ax.set_xlabel('x')
ax.set_ylabel('y')
if showIt: plt.show()
return [scalarMap]
if __name__ == '__main__':
from SimPEG import *
SimPEG/Optimization.py
+1 -1
View File
@@ -1008,4 +1008,4 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0))
delx[indx] = 0.
return delx
return delx
SimPEG/PropMaps.py
+2 -2
View File
@@ -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):
SimPEG/Regularization.py
+442 -184
View File
@@ -1,4 +1,6 @@
import Utils, Maps, Mesh, numpy as np, scipy.sparse as sp
import Utils, Maps, Mesh
import numpy as np
import scipy.sparse as sp
class RegularizationMesh(object):
"""
@@ -403,7 +405,238 @@ class BaseRegularization(object):
return mD.T * ( self.W.T * ( self.W * ( mD * v) ) )
class Tikhonov(BaseRegularization):
class Simple(BaseRegularization):
"""
Simple regularization that does not include length scales in the derivatives.
"""
mrefInSmooth = False #: include mref in the smoothness?
alpha_s = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "Smallness weight")
alpha_x = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
alpha_y = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
alpha_z = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
cell_weights = 1.
def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
if isinstance(self.cell_weights,float):
self.cell_weights = np.ones(self.regmesh.nC) * self.cell_weights
@property
def Wsmall(self):
"""Regularization matrix Wsmall"""
if getattr(self,'_Wsmall', None) is None:
self._Wsmall = Utils.sdiag((self.alpha_s*self.cell_weights)**0.5)
return self._Wsmall
@property
def Wx(self):
"""Regularization matrix Wx"""
if getattr(self, '_Wx', None) is None:
self._Wx = Utils.sdiag((self.alpha_x * (self.regmesh.aveCC2Fx*self.cell_weights))**0.5)*self.regmesh.cellDiffxStencil
return self._Wx
@property
def Wy(self):
"""Regularization matrix Wy"""
if getattr(self, '_Wy', None) is None:
self._Wy = Utils.sdiag((self.alpha_y * (self.regmesh.aveCC2Fy*self.cell_weights))**0.5)*self.regmesh.cellDiffyStencil
return self._Wy
@property
def Wz(self):
"""Regularization matrix Wz"""
if getattr(self, '_Wz', None) is None:
self._Wz = Utils.sdiag((self.alpha_z * (self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.regmesh.cellDiffzStencil
return self._Wz
# @property
# def Wsmooth(self):
# """Full smoothness regularization matrix W"""
# print 'wtf why are we using Wsmooth'
# raise NotImplementedError
# 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"""
# print 'wtf why are we using W'
# if getattr(self, '_W', None) is None:
# wlist = (self.Wsmall, self.Wx)
# if self.regmesh.dim > 1:
# wlist += (self.Wy,)
# if self.regmesh.dim > 2:
# wlist += (self.Wz,)
# self._W = sp.vstack(wlist)
# return self._W
@Utils.timeIt
def _evalSmall(self, m):
r = self.Wsmall * ( self.mapping * (m - self.mref) )
return 0.5 * r.dot(r)
@Utils.timeIt
def _evalSmallDeriv(self, m):
r = self.Wsmall * ( self.mapping * (m - self.mref) )
return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) )
@Utils.timeIt
def _evalSmall2Deriv(self, m, v = None):
rDeriv = self.Wsmall * ( self.mapping.deriv(m - self.mref) )
if v is not None:
return rDeriv.T * (rDeriv * v)
return rDeriv.T * rDeriv
@Utils.timeIt
def _evalSmoothx(self, m):
if self.mrefInSmooth == True:
r = self.Wx * ( self.mapping * (m - self.mref) )
elif self.mrefInSmooth == False:
r = self.Wx * ( self.mapping * (m) )
return 0.5 * r.dot(r)
@Utils.timeIt
def _evalSmoothy(self, m):
if self.mrefInSmooth == True:
r = self.Wy * ( self.mapping * (m - self.mref) )
elif self.mrefInSmooth == False:
r = self.Wy * ( self.mapping * (m) )
return 0.5 * r.dot(r)
@Utils.timeIt
def _evalSmoothz(self, m):
if self.mrefInSmooth == True:
r = self.Wz * ( self.mapping * (m - self.mref) )
elif self.mrefInSmooth == False:
r = self.Wz * ( self.mapping * (m) )
return 0.5 * r.dot(r)
@Utils.timeIt
def _evalSmooth(self, m):
phiSmooth = self._evalSmoothx(m)
if self.regmesh.dim > 1:
phiSmooth += self._evalSmoothy(m)
if self.regmesh.dim > 2:
phiSmooth += self._evalSmoothz(m)
return phiSmooth
@Utils.timeIt
def _evalSmoothxDeriv(self, m):
if self.mrefInSmooth == True:
r = self.Wx * ( self.mapping * ( m - self.mref ) )
return r.T * ( self.Wx * self.mapping.deriv(m - self.mref) )
elif self.mrefInSmooth == False:
r = self.Wx * ( self.mapping * m )
return r.T * ( self.Wx * self.mapping.deriv(m) )
@Utils.timeIt
def _evalSmoothx2Deriv(self, m, v=None):
if self.mrefInSmooth == True:
rDeriv = self.Wx * ( self.mapping.deriv( m - self.mref ) )
elif self.mrefInSmooth == False:
rDeriv = self.Wx * ( self.mapping.deriv(m) )
if v is not None:
return rDeriv.T * ( rDeriv * v )
return rDeriv.T * rDeriv
@Utils.timeIt
def _evalSmoothyDeriv(self, m):
if self.mrefInSmooth == True:
r = self.Wy * ( self.mapping * ( m - self.mref ) )
return r.T * ( self.Wy * self.mapping.deriv(m - self.mref) )
elif self.mrefInSmooth == False:
r = self.Wy * ( self.mapping * m )
return r.T * ( self.Wy * self.mapping.deriv(m) )
@Utils.timeIt
def _evalSmoothy2Deriv(self, m, v=None):
if self.mrefInSmooth == True:
rDeriv = self.Wy * ( self.mapping.deriv( m - self.mref ) )
elif self.mrefInSmooth == False:
rDeriv = self.Wy * ( self.mapping.deriv(m) )
if v is not None:
return rDeriv.T * ( rDeriv * v )
return rDeriv.T * rDeriv
@Utils.timeIt
def _evalSmoothzDeriv(self, m):
if self.mrefInSmooth == True:
r = self.Wz * ( self.mapping * ( m - self.mref ) )
return r.T * ( self.Wz * self.mapping.deriv(m - self.mref) )
elif self.mrefInSmooth == False:
r = self.Wz * ( self.mapping * m )
return r.T * ( self.Wz * self.mapping.deriv(m) )
@Utils.timeIt
def _evalSmoothz2Deriv(self, m, v=None):
if self.mrefInSmooth == True:
rDeriv = self.Wz * ( self.mapping.deriv( m - self.mref ) )
elif self.mrefInSmooth == False:
rDeriv = self.Wz * ( self.mapping.deriv(m) )
if v is not None:
return rDeriv.T * ( rDeriv * v )
return rDeriv.T * rDeriv
@Utils.timeIt
def _evalSmoothDeriv(self, m):
deriv = self._evalSmoothxDeriv(m)
if self.regmesh.dim > 1:
deriv += self._evalSmoothyDeriv(m)
if self.regmesh.dim > 2:
deriv += self._evalSmoothzDeriv(m)
return deriv
@Utils.timeIt
def _evalSmooth2Deriv(self, m, v=None):
deriv = self._evalSmoothx2Deriv(m, v)
if self.regmesh.dim > 1:
deriv += self._evalSmoothy2Deriv(m, v)
if self.regmesh.dim > 2:
deriv += self._evalSmoothz2Deriv(m, v)
return deriv
@Utils.timeIt
def eval(self, m):
return self._evalSmall(m) + self._evalSmooth(m)
@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})}
"""
return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)
@Utils.timeIt
def eval2Deriv(self, m, v=None):
return self._evalSmall2Deriv(m, v) + self._evalSmooth2Deriv(m, v)
class Tikhonov(Simple):
"""
L2 Tikhonov regularization with both smallness and smoothness (first order
derivative) contributions.
@@ -493,56 +726,131 @@ class Tikhonov(BaseRegularization):
self._Wzz = Utils.sdiag((self.regmesh.vol*self.alpha_zz)**0.5)*self.regmesh.faceDiffz*self.regmesh.cellDiffz
return self._Wzz
@property
def Wsmooth(self):
def Wsmooth2(self):
"""Full smoothness regularization matrix W"""
if getattr(self, '_Wsmooth', None) is None:
wlist = (self.Wx, self.Wxx)
wlist = (self.Wxx)
if self.regmesh.dim > 1:
wlist += (self.Wy, self.Wyy)
wlist += (self.Wyy)
if self.regmesh.dim > 2:
wlist += (self.Wz, self.Wzz)
wlist += (self.Wzz)
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 _evalSmall(self, m):
r = self.Wsmall * ( self.mapping * (m - self.mref) )
return 0.5 * r.dot(r)
@Utils.timeIt
def _evalSmooth(self, m):
def _evalSmoothxx(self, m):
if self.mrefInSmooth == True:
r = self.Wsmooth * ( self.mapping * (m - self.mref) )
r = self.Wxx * ( self.mapping * (m - self.mref) )
elif self.mrefInSmooth == False:
r = self.Wsmooth * ( self.mapping * (m) )
r = self.Wxx * ( self.mapping * (m) )
return 0.5 * r.dot(r)
@Utils.timeIt
def _evalSmoothyy(self, m):
if self.mrefInSmooth == True:
r = self.Wyy * ( self.mapping * (m - self.mref) )
elif self.mrefInSmooth == False:
r = self.Wyy * ( self.mapping * (m) )
return 0.5 * r.dot(r)
@Utils.timeIt
def _evalSmoothzz(self, m):
if self.mrefInSmooth == True:
r = self.Wzz * ( self.mapping * (m - self.mref) )
elif self.mrefInSmooth == False:
r = self.Wzz * ( self.mapping * (m) )
return 0.5 * r.dot(r)
@Utils.timeIt
def _evalSmooth2(self, m):
phiSmooth2 = self._evalSmoothxx(m)
if self.regmesh.dim > 1:
phiSmooth2 += self._evalSmoothyy(m)
if self.regmesh.dim > 2:
phiSmooth2 += self._evalSmoothzz(m)
return phiSmooth2
@Utils.timeIt
def _evalSmoothxxDeriv(self, m):
if self.mrefInSmooth == True:
r = self.Wxx * ( self.mapping * ( m - self.mref ) )
return r.T * ( self.Wxx * self.mapping.deriv(m - self.mref) )
elif self.mrefInSmooth == False:
r = self.Wxx * ( self.mapping * m )
return r.T * ( self.Wxx * self.mapping.deriv(m) )
@Utils.timeIt
def _evalSmoothyyDeriv(self, m):
if self.mrefInSmooth == True:
r = self.Wyy * ( self.mapping * ( m - self.mref ) )
return r.T * ( self.Wyy * self.mapping.deriv(m - self.mref) )
elif self.mrefInSmooth == False:
r = self.Wyy * ( self.mapping * m )
return r.T * ( self.Wyy * self.mapping.deriv(m) )
@Utils.timeIt
def _evalSmoothzzDeriv(self, m):
if self.mrefInSmooth == True:
r = self.Wzz * ( self.mapping * ( m - self.mref ) )
return r.T * ( self.Wzz * self.mapping.deriv(m - self.mref) )
elif self.mrefInSmooth == False:
r = self.Wzz * ( self.mapping * m )
return r.T * ( self.Wzz * self.mapping.deriv(m) )
@Utils.timeIt
def _evalSmoothxx2Deriv(self, m, v=None):
if self.mrefInSmooth == True:
rDeriv = self.Wxx * ( self.mapping.deriv( m - self.mref ) )
elif self.mrefInSmooth == False:
rDeriv = self.Wxx * self.mapping.deriv(m)
if v is not None:
return rDeriv.T * (rDeriv * v)
return rDeriv.T * rDeriv
@Utils.timeIt
def _evalSmoothyy2Deriv(self, m, v=None):
if self.mrefInSmooth == True:
rDeriv = self.Wyy * ( self.mapping.deriv( m - self.mref ) )
elif self.mrefInSmooth == False:
rDeriv = self.Wyy * self.mapping.deriv(m)
if v is not None:
return rDeriv.T * (rDeriv * v)
return rDeriv.T * rDeriv
@Utils.timeIt
def _evalSmoothzz2Deriv(self, m, v=None):
if self.mrefInSmooth == True:
rDeriv = self.Wzz * ( self.mapping.deriv( m - self.mref ) )
elif self.mrefInSmooth == False:
rDeriv = self.Wzz * self.mapping.deriv(m)
if v is not None:
return rDeriv.T * (rDeriv * v)
return rDeriv.T * rDeriv
@Utils.timeIt
def _evalSmoothDeriv2(self, m):
deriv = self._evalSmoothxxDeriv(m)
if self.regmesh.dim > 1:
deriv += self._evalSmoothyyDeriv(m)
if self.regmesh.dim > 2:
deriv += self._evalSmoothzzDeriv(m)
return deriv
@Utils.timeIt
def _evalSmooth2Deriv2(self, m, v=None):
deriv = self._evalSmoothxx2Deriv(m, v)
if self.regmesh.dim > 1:
deriv += self._evalSmoothyy2Deriv(m, v)
if self.regmesh.dim > 2:
deriv += self._evalSmoothzz2Deriv(m, v)
return deriv
@Utils.timeIt
def eval(self, m):
return self._evalSmall(m) + self._evalSmooth(m)
@Utils.timeIt
def _evalSmallDeriv(self,m):
r = self.Wsmall * ( self.mapping * (m - self.mref) )
return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) )
@Utils.timeIt
def _evalSmoothDeriv(self,m):
if self.mrefInSmooth == True:
r = self.Wsmooth * ( self.mapping * ( m - self.mref ) )
return r.T * ( self.Wsmooth * self.mapping.deriv(m - self.mref) )
elif self.mrefInSmooth == False:
r = self.Wsmooth * ( self.mapping * m )
return r.T * ( self.Wsmooth * self.mapping.deriv(m) )
return self._evalSmall(m) + self._evalSmooth(m) + self._evalSmooth2(m)
@Utils.timeIt
def evalDeriv(self, m):
@@ -560,184 +868,134 @@ class Tikhonov(BaseRegularization):
R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
"""
return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)
return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + self._evalSmoothDeriv2(m)
def eval2Deriv(self, m, v=None):
"""
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})}
"""
return self._evalSmall2Deriv(m, v) + self._evalSmooth2Deriv(m, v) + self._evalSmooth2Deriv2(m, v)
class Simple(Tikhonov):
class Sparse(Simple):
"""
Simple regularization that does not include length scales in the derivatives.
The regularization is:
.. math::
R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top R^\\top R W(m-m_\\text{ref})}
where the IRLS weight
.. math::
R = \eta TO FINISH LATER!!!
So the derivative is straight forward:
.. math::
R(m) = \mathbf{W^\\top R^\\top R W (m-m_\\text{ref})}
The IRLS weights are recomputed after each beta solves.
It is strongly recommended to do a few Gauss-Newton iterations
before updating.
"""
mrefInSmooth = False #: SMOOTH and SMOOTH_MOD_DIF options
alpha_s = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "Smallness weight")
alpha_x = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
alpha_y = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
alpha_z = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
wght = 1.
# set default values
eps_p = 1e-1 # Threshold value for the model norm
eps_q = 1e-1 # Threshold value for the model gradient norm
curModel = None # Requires model to compute the weights
l2model = None
gamma = 1. # Model norm scaling to smooth out convergence
norms = [0., 2., 2., 2.] # Values for norm on (m, dmdx, dmdy, dmdz)
cell_weights = 1. # Consider overwriting with sensitivity weights
def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
if isinstance(self.wght,float):
self.wght = np.ones(self.regmesh.nC) * self.wght
if isinstance(self.cell_weights,float):
self.cell_weights = np.ones(self.regmesh.nC) * self.cell_weights
@property
def Wsmall(self):
"""Regularization matrix Wsmall"""
if getattr(self,'_Wsmall', None) is None:
self._Wsmall = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5)
if getattr(self, 'curModel', None) is None:
self.Rs = Utils.speye(self.regmesh.nC)
else:
f_m = self.mapping * (self.curModel - self.reg.mref)
self.rs = self.R(f_m , self.eps_p, self.norms[0])
self.Rs = Utils.sdiag( self.rs )
self._Wsmall = Utils.sdiag((self.alpha_s*self.gamma*self.cell_weights)**0.5)*self.Rs
return self._Wsmall
@property
def Wx(self):
"""Regularization matrix Wx"""
if getattr(self, '_Wx', None) is None:
self._Wx = Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.regmesh.cellDiffxStencil
if getattr(self,'_Wx', None) is None:
if getattr(self, 'curModel', None) is None:
self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
else:
f_m = self.regmesh.cellDiffxStencil * (self.mapping * self.curModel)
self.rx = self.R( f_m , self.eps_q, self.norms[1])
self.Rx = Utils.sdiag( self.rx )
self._Wx = Utils.sdiag(( self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.cell_weights))**0.5)*self.Rx*self.regmesh.cellDiffxStencil
return self._Wx
@property
def Wy(self):
"""Regularization matrix Wy"""
if getattr(self, '_Wy', None) is None:
self._Wy = Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol * self.alpha_y*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.regmesh.cellDiffyStencil
if getattr(self,'_Wy', None) is None:
if getattr(self, 'curModel', None) is None:
self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
else:
f_m = self.regmesh.cellDiffyStencil * (self.mapping * self.curModel)
self.ry = self.R( f_m , self.eps_q, self.norms[2])
self.Ry = Utils.sdiag( self.ry )
self._Wy = Utils.sdiag((self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.cell_weights))**0.5)*self.Ry*self.regmesh.cellDiffyStencil
return self._Wy
@property
def Wz(self):
"""Regularization matrix Wz"""
if getattr(self, '_Wz', None) is None:
self._Wz = Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.regmesh.cellDiffzStencil
if getattr(self,'_Wz', None) is None:
if getattr(self, 'curModel', None) is None:
self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
else:
f_m = self.regmesh.cellDiffzStencil * (self.mapping * self.curModel)
self.rz = self.R( f_m , self.eps_q, self.norms[3])
self.Rz = Utils.sdiag( self.rz )
self._Wz = Utils.sdiag((self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
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 _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)
class Sparse(Simple):
# set default values
eps_p = 1e-1
eps_q = 1e-1
curModel = None # use a model to compute the weights
gamma = 1.
norms = [0., 2., 2., 2.]
wght = 1.
def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
if isinstance(self.wght,float):
self.wght = np.ones(self.regmesh.nC) * self.wght
@property
def Wsmall(self):
"""Regularization matrix Wsmall"""
if getattr(self, 'curModel', None) is None:
self.Rs = Utils.speye(self.regmesh.nC)
else:
f_m = self.curModel - self.reg.mref
self.rs = self.R(f_m , self.eps_p, self.norms[0])
#print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
self.Rs = Utils.sdiag( self.rs )
return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.wght)**0.5)*self.Rs
@property
def Wx(self):
"""Regularization matrix Wx"""
if getattr(self, 'curModel', None) is None:
self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
else:
f_m = self.regmesh.cellDiffxStencil * self.curModel
self.rx = self.R( f_m , self.eps_q, self.norms[1])
self.Rx = Utils.sdiag( self.rx )
return Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.Rx*self.regmesh.cellDiffxStencil
@property
def Wy(self):
"""Regularization matrix Wy"""
if getattr(self, 'curModel', None) is None:
self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
else:
f_m = self.regmesh.cellDiffyStencil * self.curModel
self.ry = self.R( f_m , self.eps_q, self.norms[2])
self.Ry = Utils.sdiag( self.ry )
return Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.Ry*self.regmesh.cellDiffyStencil
@property
def Wz(self):
"""Regularization matrix Wz"""
if getattr(self, 'curModel', None) is None:
self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
else:
f_m = self.regmesh.cellDiffzStencil * self.curModel
self.rz = self.R( f_m , self.eps_q, self.norms[3])
self.Rz = Utils.sdiag( self.rz )
return Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
@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 sp.vstack(wlist)
@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
def R(self, f_m , eps, exponent):
# Eta scaling is important for mix-norms...do not mess with it
eta = (eps**(1.-exponent/2.))**0.5
r = eta / (f_m**2.+ eps**2.)**((1.-exponent/2.)/2.)
SimPEG/Utils/__init__.py
+1
View File
@@ -7,3 +7,4 @@ from CounterUtils import *
import ModelBuilder
import SolverUtils
from coordutils import *
from modelutils import *
SimPEG/Utils/modelutils.py
+63
View File
@@ -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)
docs/examples/DC_Forward_PseudoSection.rst
+2 -2
View File
@@ -20,8 +20,8 @@ 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)
surveyType = survey type 'pole-dipole' or 'dipole-dipole'
unitType = Data type "appResistivity" | "appConductivity" | "volt"
Created by @fourndo
docs/examples/Mesh_Basic_ForwardDC.rst
+25
View File
@@ -0,0 +1,25 @@
.. _examples_Mesh_Basic_ForwardDC:
.. --------------------------------- ..
.. ..
.. THIS FILE IS AUTO GENEREATED ..
.. ..
.. SimPEG/Examples/__init__.py ..
.. ..
.. --------------------------------- ..
Mesh: Basic Forward 2D DC Resistivity
=====================================
2D DC forward modeling example with Tensor and Curvilinear Meshes
.. plot::
from SimPEG import Examples
Examples.Mesh_Basic_ForwardDC.run()
.. literalinclude:: ../../SimPEG/Examples/Mesh_Basic_ForwardDC.py
:language: python
:linenos:
docs/examples/Forward_BasicDirectCurrent.rst → docs/examples/Utils_surface2ind_topo.rst
+8 -5
View File
@@ -1,4 +1,4 @@
.. _examples_Forward_BasicDirectCurrent:
.. _examples_Utils_surface2ind_topo:
.. --------------------------------- ..
.. ..
@@ -8,14 +8,17 @@
.. ..
.. --------------------------------- ..
Forward BasicDirectCurrent
==========================
Here we show how to use :code:`Utils.surface2ind_topo` to identify cells below
a topographic surface.
.. plot::
from SimPEG import Examples
Examples.Forward_BasicDirectCurrent.run()
Examples.Utils_surface2ind_topo.run()
.. literalinclude:: ../../SimPEG/Examples/Forward_BasicDirectCurrent.py
.. literalinclude:: ../../SimPEG/Examples/Utils_surface2ind_topo.py
:language: python
:linenos:
setup.py
+12 -6
View File
@@ -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
)
tests/base/test_PropMaps.py
+29
View File
@@ -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()
tests/em/fdem/forward/test_FDEM_analytics.py
+1 -1
View File
@@ -28,7 +28,7 @@ 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.eField(XYZ, orientation='x', real_or_imag='imag')
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])
tests/em/static/__init__.py
+12
View File
@@ -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)
tests/em/static/test_DC_2D_analytic.py
+69
View File
@@ -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()
tests/em/static/test_DC_2D_jvecjtvecadj.py
+127
View File
@@ -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()
tests/em/static/test_DC_analytic.py
+71
View File
@@ -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()
tests/em/static/test_DC_jvecjtvecadj.py
+127
View File
@@ -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()
tests/em/static/test_IP_fwd.py
+96
View File
@@ -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()
tests/em/static/test_IP_jvecjtvecadj.py
+126
View File
@@ -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()
tests/em/static/test_SIP_jvecjtvecadj.py
+232
View File
@@ -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()
tests/mesh/test_Mixed_boundaryPoisson.py
+411
View File
@@ -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()