Compare commits

...
17 Commits
Author SHA1 Message Date
Lindsey Heagy 8f1df12252 copy docs badges from readme to docs index 2016-07-21 22:30:49 -07:00
Lindsey Heagy 002554209b listing and remove wildcard inputs from tes_Mixed_boundaryPoisson 2016-07-21 22:28:20 -07:00
Lindsey Heagy 76f140f89b pep8 linting 2016-07-21 21:58:42 -07:00
Lindsey Heagy 5ce26a97ae remove wildcard imports from test_Fields 2016-07-21 21:39:11 -07:00
Lindsey Heagy a087c6e5b1 fix broken link 2016-07-21 21:34:17 -07:00
Lindsey Heagy 44ec0623d6 fixed formatting changes from quantified code, test links in docs 2016-07-21 14:38:05 -07:00
Lindsey Heagy 57d0c84b41 use @property for nodal and cell centered curvi grids 2016-07-21 13:19:23 -07:00
Lindsey Heagy a7befd6783 fix typo in cell grad 2016-07-21 12:38:21 -07:00
Lindsey Heagy 7acf46ddb8 fix typo in faceDivz 2016-07-21 12:32:49 -07:00
Lindsey Heagy b54d2494f5 add codacy to readme 2016-07-21 12:20:58 -07:00
Cody 45f5906554 Migrated % string formating 2016-07-21 12:08:19 -07:00
Lindsey Heagy 8093288391 fix double negative 2016-07-21 11:59:59 -07:00
Lindsey Heagy 7457912d84 avoid double import of numpy 2016-07-21 11:58:52 -07:00
Lindsey Heagy 910d4214b3 remove extra import of unittest 2016-07-21 11:56:49 -07:00
Lindsey Heagy 20e80ed983 use @property decorator in DiffOperators.py 2016-07-21 11:55:14 -07:00
Lindsey Heagy 3e8580f28d update readme 2016-07-21 09:59:43 -07:00
Lindsey Heagy 3ad2c98d43 use @property decorator in curve mesh (see: https://www.quantifiedcode.com/app/project/933aa3decf444538aa432c8817169b6d?groups=code_patterns%3A3bECxdfc%3Af0&tab=basics) 2016-07-21 09:58:14 -07:00
57 changed files with 1233 additions and 994 deletions
+7 -4
View File
@@ -21,10 +21,6 @@ SimPEG
:target: https://travis-ci.org/simpeg/simpeg :target: https://travis-ci.org/simpeg/simpeg
:alt: Travis CI build status :alt: Travis CI build status
.. image:: https://img.shields.io/coveralls/simpeg/simpeg.svg
:target: https://coveralls.io/r/simpeg/simpeg?branch=master
:alt: Coverage status
.. image:: http://img.shields.io/badge/GITTER-JOIN_CHAT-brightgreen.svg?style=flat-square .. image:: http://img.shields.io/badge/GITTER-JOIN_CHAT-brightgreen.svg?style=flat-square
:alt: gitter chat room at https://gitter.im/simpeg/simpeg :alt: gitter chat room at https://gitter.im/simpeg/simpeg
:target: https://gitter.im/simpeg/simpeg :target: https://gitter.im/simpeg/simpeg
@@ -32,6 +28,13 @@ SimPEG
.. image:: https://codecov.io/gh/simpeg/simpeg/branch/master/graph/badge.svg .. image:: https://codecov.io/gh/simpeg/simpeg/branch/master/graph/badge.svg
   :target: https://codecov.io/gh/simpeg/simpeg    :target: https://codecov.io/gh/simpeg/simpeg
.. image:: https://www.quantifiedcode.com/api/v1/project/933aa3decf444538aa432c8817169b6d/badge.svg
:target: https://www.quantifiedcode.com/app/project/933aa3decf444538aa432c8817169b6d
:alt: Code issues
.. image:: https://api.codacy.com/project/badge/Grade/4fc959a5294a418fa21fc7bc3b3aa078
:target: https://www.codacy.com/app/lindseyheagy/simpeg?utm_source=github.com&utm_medium=referral&utm_content=simpeg/simpeg&utm_campaign=Badge_Grade
:alt: codacy
Simulation and Parameter Estimation in Geophysics - A python package for simulation and gradient based parameter estimation in the context of geophysical applications. Simulation and Parameter Estimation in Geophysics - A python package for simulation and gradient based parameter estimation in the context of geophysical applications.
+7 -7
View File
@@ -476,7 +476,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
fid.write('! ' + surveyType + ' FORMAT\n') fid.write('! ' + surveyType + ' FORMAT\n')
if iptype!=0: if iptype!=0:
fid.write('IPTYPE=%i\n'%iptype) fid.write('IPTYPE={0:d}\n'.format(iptype))
else: else:
fid.write('! ' + stype + ' FORMAT\n') fid.write('! ' + stype + ' FORMAT\n')
@@ -512,7 +512,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
if surveyType == 'SURFACE': if surveyType == 'SURFACE':
fid.writelines("%f " % ii for ii in mkvc(tx[0,:])) fid.writelines("{0:f} ".format(ii) for ii in mkvc(tx[0,:]))
M = M[:,0] M = M[:,0]
N = N[:,0] N = N[:,0]
@@ -521,7 +521,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
# Flip sign for z-elevation to depth # Flip sign for z-elevation to depth
tx[2::2,:] = -tx[2::2,:] tx[2::2,:] = -tx[2::2,:]
fid.writelines("%e " % ii for ii in mkvc(tx[::2,:])) fid.writelines("{0:e} ".format(ii) for ii in mkvc(tx[::2,:]))
M = M[:,0::2] M = M[:,0::2]
N = N[:,0::2] N = N[:,0::2]
@@ -529,22 +529,22 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
M[:,1::2] = -M[:,1::2] M[:,1::2] = -M[:,1::2]
N[:,1::2] = -N[:,1::2] N[:,1::2] = -N[:,1::2]
fid.write('%i\n'% nD) fid.write('{0:d}\n'.format(nD))
np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%f',delimiter=' ',newline='\n') np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%f',delimiter=' ',newline='\n')
if dim=='3D': if dim=='3D':
if surveyType == 'SURFACE': if surveyType == 'SURFACE':
fid.writelines("%e " % ii for ii in mkvc(tx[0:2,:])) fid.writelines("{0:e} ".format(ii) for ii in mkvc(tx[0:2,:]))
M = M[:,0:2] M = M[:,0:2]
N = N[:,0:2] N = N[:,0:2]
if surveyType == 'GENERAL': if surveyType == 'GENERAL':
fid.writelines("%e " % ii for ii in mkvc(tx[0:3,:])) fid.writelines("{0:e} ".format(ii) for ii in mkvc(tx[0:3,:]))
fid.write('%i\n'% nD) fid.write('{0:d}\n'.format(nD))
np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n') 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') fid.write('\n')
+14 -14
View File
@@ -15,7 +15,7 @@ class InversionDirective(object):
@inversion.setter @inversion.setter
def inversion(self, i): def inversion(self, i):
if getattr(self,'_inversion',None) is not None: if getattr(self,'_inversion',None) is not None:
print 'Warning: InversionDirective %s has switched to a new inversion.' % self.__name__ print 'Warning: InversionDirective {0!s} has switched to a new inversion.'.format(self.__name__)
self._inversion = i self._inversion = i
@property @property
@@ -47,7 +47,7 @@ class DirectiveList(object):
def __init__(self, *directives, **kwargs): def __init__(self, *directives, **kwargs):
self.dList = [] self.dList = []
for d in directives: for d in directives:
assert isinstance(d, InversionDirective), 'All directives must be InversionDirectives not %s' % d.__name__ assert isinstance(d, InversionDirective), 'All directives must be InversionDirectives not {0!s}'.format(d.__name__)
self.dList.append(d) self.dList.append(d)
Utils.setKwargs(self, **kwargs) Utils.setKwargs(self, **kwargs)
@@ -68,7 +68,7 @@ class DirectiveList(object):
def inversion(self, i): def inversion(self, i):
if self.inversion is i: return if self.inversion is i: return
if getattr(self,'_inversion',None) is not None: if getattr(self,'_inversion',None) is not None:
print 'Warning: %s has switched to a new inversion.' % self.__name__ print 'Warning: {0!s} has switched to a new inversion.'.format(self.__name__)
for d in self.dList: for d in self.dList:
d.inversion = i d.inversion = i
self._inversion = i self._inversion = i
@@ -79,7 +79,7 @@ class DirectiveList(object):
return return
directives = ['initialize', 'endIter', 'finish'] directives = ['initialize', 'endIter', 'finish']
assert ruleType in directives, 'Directive type must be in ["%s"]' % '", "'.join(directives) assert ruleType in directives, 'Directive type must be in ["{0!s}"]'.format('", "'.join(directives))
for r in self.dList: for r in self.dList:
getattr(r, ruleType)() getattr(r, ruleType)()
@@ -141,7 +141,7 @@ class BetaSchedule(InversionDirective):
def endIter(self): def endIter(self):
if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0: 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 if self.debug: print 'BetaSchedule is cooling Beta. Iteration: {0:d}'.format(self.opt.iter)
self.invProb.beta /= self.coolingFactor self.invProb.beta /= self.coolingFactor
@@ -181,7 +181,7 @@ class SaveEveryIteration(InversionDirective):
def fileName(self): def fileName(self):
if getattr(self, '_fileName', None) is None: if getattr(self, '_fileName', None) is None:
from datetime import datetime from datetime import datetime
self._fileName = '%s-%s'%(self.name, datetime.now().strftime('%Y-%m-%d-%H-%M')) self._fileName = '{0!s}-{1!s}'.format(self.name, datetime.now().strftime('%Y-%m-%d-%H-%M'))
return self._fileName return self._fileName
@fileName.setter @fileName.setter
def fileName(self, value): def fileName(self, value):
@@ -192,31 +192,31 @@ class SaveModelEveryIteration(SaveEveryIteration):
"""SaveModelEveryIteration""" """SaveModelEveryIteration"""
def initialize(self): def initialize(self):
print "SimPEG.SaveModelEveryIteration will save your models as: '###-%s.npy'"%self.fileName print "SimPEG.SaveModelEveryIteration will save your models as: '###-{0!s}.npy'".format(self.fileName)
def endIter(self): def endIter(self):
np.save('%03d-%s' % (self.opt.iter, self.fileName), self.opt.xc) np.save('{0:03d}-{1!s}'.format(self.opt.iter, self.fileName), self.opt.xc)
class SaveOutputEveryIteration(SaveEveryIteration): class SaveOutputEveryIteration(SaveEveryIteration):
"""SaveModelEveryIteration""" """SaveModelEveryIteration"""
def initialize(self): def initialize(self):
print "SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-%s.txt'"%self.fileName print "SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-{0!s}.txt'".format(self.fileName)
f = open(self.fileName+'.txt', 'w') f = open(self.fileName+'.txt', 'w')
f.write(" # beta phi_d phi_m f\n") f.write(" # beta phi_d phi_m f\n")
f.close() f.close()
def endIter(self): def endIter(self):
f = open(self.fileName+'.txt', 'a') f = open(self.fileName+'.txt', 'a')
f.write(' %3d %1.4e %1.4e %1.4e %1.4e\n'%(self.opt.iter, self.invProb.beta, self.invProb.phi_d, self.invProb.phi_m, self.opt.f)) f.write(' {0:3d} {1:1.4e} {2:1.4e} {3:1.4e} {4:1.4e}\n'.format(self.opt.iter, self.invProb.beta, self.invProb.phi_d, self.invProb.phi_m, self.opt.f))
f.close() f.close()
class SaveOutputDictEveryIteration(SaveEveryIteration): class SaveOutputDictEveryIteration(SaveEveryIteration):
"""SaveOutputDictEveryIteration""" """SaveOutputDictEveryIteration"""
def initialize(self): def initialize(self):
print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-{0!s}.npz'".format(self.fileName)
def endIter(self): def endIter(self):
# Save the data. # Save the data.
@@ -328,7 +328,7 @@ class Update_IRLS(InversionDirective):
# Beta Schedule # Beta Schedule
if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0: 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 if self.debug: print 'BetaSchedule is cooling Beta. Iteration: {0:d}'.format(self.opt.iter)
self.invProb.beta /= self.coolingFactor self.invProb.beta /= self.coolingFactor
@@ -340,11 +340,11 @@ class Update_IRLS(InversionDirective):
phim_new = self.reg.eval(self.invProb.curModel) phim_new = self.reg.eval(self.invProb.curModel)
self.f_change = np.abs(self.f_old - phim_new) / self.f_old self.f_change = np.abs(self.f_old - phim_new) / self.f_old
print "Regularization decrease: %6.3e" % (self.f_change) print "Regularization decrease: {0:6.3e}".format((self.f_change))
# Check for maximum number of IRLS cycles # Check for maximum number of IRLS cycles
if self.IRLSiter == self.maxIRLSiter: if self.IRLSiter == self.maxIRLSiter:
print "Reach maximum number of IRLS cycles: %i" % self.maxIRLSiter print "Reach maximum number of IRLS cycles: {0:d}".format(self.maxIRLSiter)
self.opt.stopNextIteration = True self.opt.stopNextIteration = True
return return
+8 -8
View File
@@ -42,7 +42,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: total electric field :return: total electric field
""" """
if getattr(self, '_ePrimary', None) is None or getattr(self, '_eSecondary', None) is None: if getattr(self, '_ePrimary', None) is None or getattr(self, '_eSecondary', None) is None:
raise NotImplementedError ('Getting e from %s is not implemented' %self.knownFields.keys()[0]) raise NotImplementedError ('Getting e from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
return self._ePrimary(solution,srcList) + self._eSecondary(solution,srcList) return self._ePrimary(solution,srcList) + self._eSecondary(solution,srcList)
@@ -56,7 +56,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: total magnetic flux density :return: total magnetic flux density
""" """
if getattr(self, '_bPrimary', None) is None or getattr(self, '_bSecondary', None) is None: if getattr(self, '_bPrimary', None) is None or getattr(self, '_bSecondary', None) is None:
raise NotImplementedError ('Getting b from %s is not implemented' %self.knownFields.keys()[0]) raise NotImplementedError ('Getting b from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
return self._bPrimary(solution, srcList) + self._bSecondary(solution, srcList) return self._bPrimary(solution, srcList) + self._bSecondary(solution, srcList)
@@ -70,7 +70,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: total magnetic field :return: total magnetic field
""" """
if getattr(self, '_hPrimary', None) is None or getattr(self, '_hSecondary', None) is None: if getattr(self, '_hPrimary', None) is None or getattr(self, '_hSecondary', None) is None:
raise NotImplementedError ('Getting h from %s is not implemented' %self.knownFields.keys()[0]) raise NotImplementedError ('Getting h from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
return self._hPrimary(solution, srcList) + self._hSecondary(solution, srcList) return self._hPrimary(solution, srcList) + self._hSecondary(solution, srcList)
@@ -84,7 +84,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: total current density :return: total current density
""" """
if getattr(self, '_jPrimary', None) is None or getattr(self, '_jSecondary', None) is None: if getattr(self, '_jPrimary', None) is None or getattr(self, '_jSecondary', None) is None:
raise NotImplementedError ('Getting j from %s is not implemented' %self.knownFields.keys()[0]) raise NotImplementedError ('Getting j from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
return self._jPrimary(solution, srcList) + self._jSecondary(solution, srcList) return self._jPrimary(solution, srcList) + self._jSecondary(solution, srcList)
@@ -100,7 +100,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: derivative times a vector (or tuple for adjoint) :return: derivative times a vector (or tuple for adjoint)
""" """
if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None: 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]) raise NotImplementedError ('Getting eDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
if adjoint: if adjoint:
return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint) return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
@@ -118,7 +118,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: derivative times a vector (or tuple for adjoint) :return: derivative times a vector (or tuple for adjoint)
""" """
if getattr(self, '_bDeriv_u', None) is None or getattr(self, '_bDeriv_m', None) is None: if getattr(self, '_bDeriv_u', None) is None or getattr(self, '_bDeriv_m', None) is None:
raise NotImplementedError ('Getting bDerivs from %s is not implemented' %self.knownFields.keys()[0]) raise NotImplementedError ('Getting bDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
if adjoint: if adjoint:
return self._bDeriv_u(src, v, adjoint), self._bDeriv_m(src, v, adjoint) return self._bDeriv_u(src, v, adjoint), self._bDeriv_m(src, v, adjoint)
@@ -136,7 +136,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: derivative times a vector (or tuple for adjoint) :return: derivative times a vector (or tuple for adjoint)
""" """
if getattr(self, '_hDeriv_u', None) is None or getattr(self, '_hDeriv_m', None) is None: if getattr(self, '_hDeriv_u', None) is None or getattr(self, '_hDeriv_m', None) is None:
raise NotImplementedError ('Getting hDerivs from %s is not implemented' %self.knownFields.keys()[0]) raise NotImplementedError ('Getting hDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
if adjoint: if adjoint:
return self._hDeriv_u(src, v, adjoint), self._hDeriv_m(src, v, adjoint) return self._hDeriv_u(src, v, adjoint), self._hDeriv_m(src, v, adjoint)
@@ -154,7 +154,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: derivative times a vector (or tuple for adjoint) :return: derivative times a vector (or tuple for adjoint)
""" """
if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None: 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]) raise NotImplementedError ('Getting jDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
if adjoint: if adjoint:
return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint) return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
+2 -2
View File
@@ -11,8 +11,8 @@ class BaseRx(SimPEG.Survey.BaseRx):
""" """
def __init__(self, locs, orientation=None, component=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(orientation in ['x','y','z']), "Orientation {0!s} not known. Orientation must be in 'x', 'y', 'z'. Arbitrary orientations have not yet been implemented.".format(orientation)
assert(component in ['real', 'imag']), "'component' must be 'real' or 'imag', not %s"%component assert(component in ['real', 'imag']), "'component' must be 'real' or 'imag', not {0!s}".format(component)
self.projComp = orientation self.projComp = orientation
self.component = component self.component = component
+3 -3
View File
@@ -9,7 +9,7 @@ class Fields(SimPEG.Problem.Fields):
def _phiDeriv(self, src, du_dm_v, v, adjoint=False): 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: 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]) raise NotImplementedError ('Getting phiDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
if adjoint: if adjoint:
return self._phiDeriv_u(src, v, adjoint=adjoint), self._phiDeriv_m(src, v, adjoint=adjoint) return self._phiDeriv_u(src, v, adjoint=adjoint), self._phiDeriv_m(src, v, adjoint=adjoint)
@@ -18,7 +18,7 @@ class Fields(SimPEG.Problem.Fields):
def _eDeriv(self, src, du_dm_v, v, adjoint=False): 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: 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]) raise NotImplementedError ('Getting eDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
if adjoint: if adjoint:
return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint) return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
@@ -26,7 +26,7 @@ class Fields(SimPEG.Problem.Fields):
def _jDeriv(self, src, du_dm_v, v, adjoint=False): 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: 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]) raise NotImplementedError ('Getting jDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
if adjoint: if adjoint:
return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint) return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
+3 -3
View File
@@ -32,7 +32,7 @@ class Fields_ky(SimPEG.Problem.TimeFields):
def _phiDeriv(self,kyInd, src, du_dm_v, v, adjoint=False): 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: 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]) raise NotImplementedError ('Getting phiDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
if adjoint: if adjoint:
return self._phiDeriv_u(kyInd, src, v, adjoint=adjoint), self._phiDeriv_m(kyInd, src, v, adjoint=adjoint) return self._phiDeriv_u(kyInd, src, v, adjoint=adjoint), self._phiDeriv_m(kyInd, src, v, adjoint=adjoint)
@@ -41,7 +41,7 @@ class Fields_ky(SimPEG.Problem.TimeFields):
def _eDeriv(self,kyInd, src, du_dm_v, v, adjoint=False): 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: 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]) raise NotImplementedError ('Getting eDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
if adjoint: if adjoint:
return self._eDeriv_u(kyInd, src, v, adjoint), self._eDeriv_m(kyInd, src, v, adjoint) return self._eDeriv_u(kyInd, src, v, adjoint), self._eDeriv_m(kyInd, src, v, adjoint)
@@ -49,7 +49,7 @@ class Fields_ky(SimPEG.Problem.TimeFields):
def _jDeriv(self,kyInd, src, du_dm_v, v, adjoint=False): 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: 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]) raise NotImplementedError ('Getting jDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
if adjoint: if adjoint:
return self._jDeriv_u(kyInd, src, v, adjoint), self._jDeriv_m(kyInd, src, v, adjoint) return self._jDeriv_u(kyInd, src, v, adjoint), self._jDeriv_m(kyInd, src, v, adjoint)
+2 -2
View File
@@ -46,7 +46,7 @@ class BaseDCProblem(BaseEMProblem):
du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v ) du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v )
for rx in src.rxList: for rx in src.rxList:
df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) df_dmFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False) df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v) Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
return Utils.mkvc(Jv) return Utils.mkvc(Jv)
@@ -69,7 +69,7 @@ class BaseDCProblem(BaseEMProblem):
u_src = f[src, self._solutionType] u_src = f[src, self._solutionType]
for rx in src.rxList: for rx in src.rxList:
PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m 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!s}Deriv'.format(rx.projField), None)
df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True) df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
ATinvdf_duT = self.Ainv * df_duT ATinvdf_duT = self.Ainv * df_duT
+2 -2
View File
@@ -60,7 +60,7 @@ class BaseDCProblem_2D(BaseEMProblem):
dRHS_dm_v = self.getRHSDeriv(ky, src, v) dRHS_dm_v = self.getRHSDeriv(ky, src, v)
du_dm_v = self.Ainv[iky] * ( - dA_dm_v + dRHS_dm_v ) du_dm_v = self.Ainv[iky] * ( - dA_dm_v + dRHS_dm_v )
for rx in src.rxList: for rx in src.rxList:
df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) df_dmFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
df_dm_v = df_dmFun(iky, src, du_dm_v, v, adjoint=False) df_dm_v = df_dmFun(iky, src, du_dm_v, v, adjoint=False)
# Trapezoidal intergration # Trapezoidal intergration
Jv1_temp = 1./np.pi*rx.evalDeriv(ky, src, self.mesh, f, df_dm_v) Jv1_temp = 1./np.pi*rx.evalDeriv(ky, src, self.mesh, f, df_dm_v)
@@ -101,7 +101,7 @@ class BaseDCProblem_2D(BaseEMProblem):
ky = self.kys[iky] ky = self.kys[iky]
AT = self.getA(ky) AT = self.getA(ky)
PTv = rx.evalDeriv(ky, src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m 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_duTFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
df_duT, df_dmT = df_duTFun(iky, src, None, PTv, adjoint=True) df_duT, df_dmT = df_duTFun(iky, src, None, PTv, adjoint=True)
ATinvdf_duT = self.Ainv[iky] * df_duT ATinvdf_duT = self.Ainv[iky] * df_duT
+2 -2
View File
@@ -56,7 +56,7 @@ class BaseIPProblem(BaseEMProblem):
du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v ) du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v )
for rx in src.rxList: for rx in src.rxList:
df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) df_dmFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False) df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v) Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
# Conductivity (d u / d log sigma) # Conductivity (d u / d log sigma)
@@ -83,7 +83,7 @@ class BaseIPProblem(BaseEMProblem):
u_src = f[src, self._solutionType] u_src = f[src, self._solutionType]
for rx in src.rxList: for rx in src.rxList:
PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m 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!s}Deriv'.format(rx.projField), None)
df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True) df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
ATinvdf_duT = self.Ainv * df_duT ATinvdf_duT = self.Ainv * df_duT
dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True) dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True)
+3 -3
View File
@@ -83,7 +83,7 @@ class BaseSIPProblem(BaseEMProblem):
for rx in src.rxList: for rx in src.rxList:
timeindex = rx.getTimeP(self.survey.times) timeindex = rx.getTimeP(self.survey.times)
if timeindex[tind]: if timeindex[tind]:
df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) df_dmFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False) 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) Jv[src, rx, t] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
@@ -122,7 +122,7 @@ class BaseSIPProblem(BaseEMProblem):
for rx in src.rxList: for rx in src.rxList:
timeindex = rx.getTimeP(self.survey.times) timeindex = rx.getTimeP(self.survey.times)
if timeindex[tind]: if timeindex[tind]:
df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) df_dmFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
df_dm_v0 = df_dmFun(src, du_dm_v0, v0, adjoint=False) df_dm_v0 = df_dmFun(src, du_dm_v0, v0, adjoint=False)
df_dm_v1 = df_dmFun(src, du_dm_v1, v1, 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_v0)
@@ -153,7 +153,7 @@ class BaseSIPProblem(BaseEMProblem):
timeindex = rx.getTimeP(self.survey.times) timeindex = rx.getTimeP(self.survey.times)
if timeindex[tind]: if timeindex[tind]:
PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx, t], adjoint=True) # wrt f, need possibility wrt m 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_duTFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True) df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
ATinvdf_duT = self.Ainv * df_duT ATinvdf_duT = self.Ainv * df_duT
dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True) dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True)
+10 -10
View File
@@ -47,7 +47,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
self.waveformType = "GENERAL" self.waveformType = "GENERAL"
def fields(self, m): def fields(self, m):
if self.verbose: print '%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50) if self.verbose: print '{0!s}\nCalculating fields(m)\n{1!s}'.format('*'*50, '*'*50)
self.curModel = m self.curModel = m
# Create a fields storage object # Create a fields storage object
F = self._FieldsForward_pair(self.mesh, self.survey) F = self._FieldsForward_pair(self.mesh, self.survey)
@@ -55,7 +55,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
# Set the initial conditions # Set the initial conditions
F[src,:,0] = src.getInitialFields(self.mesh) F[src,:,0] = src.getInitialFields(self.mesh)
F = self.forward(m, self.getRHS, F=F) F = self.forward(m, self.getRHS, F=F)
if self.verbose: print '%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50) if self.verbose: print '{0!s}\nDone calculating fields(m)\n{1!s}'.format('*'*50, '*'*50)
return F return F
def forward(self, m, RHS, F=None): def forward(self, m, RHS, F=None):
@@ -70,11 +70,11 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
if Ainv is not None: if Ainv is not None:
Ainv.clean() Ainv.clean()
A = self.getA(tInd) A = self.getA(tInd)
if self.verbose: print 'Factoring... (dt = %e)'%dt if self.verbose: print 'Factoring... (dt = {0:e})'.format(dt)
Ainv = self.Solver(A, **self.solverOpts) Ainv = self.Solver(A, **self.solverOpts)
if self.verbose: print 'Done' if self.verbose: print 'Done'
rhs = RHS(tInd, F) rhs = RHS(tInd, F)
if self.verbose: print ' Solving... (tInd = %d)'%tInd if self.verbose: print ' Solving... (tInd = {0:d})'.format(tInd)
sol = Ainv * rhs sol = Ainv * rhs
if self.verbose: print ' Done...' if self.verbose: print ' Done...'
if sol.ndim == 1: if sol.ndim == 1:
@@ -95,11 +95,11 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
if Ainv is not None: if Ainv is not None:
Ainv.clean() Ainv.clean()
A = self.getA(tInd) A = self.getA(tInd)
if self.verbose: print 'Factoring (Adjoint)... (dt = %e)'%dt if self.verbose: print 'Factoring (Adjoint)... (dt = {0:e})'.format(dt)
Ainv = self.Solver(A, **self.solverOpts) Ainv = self.Solver(A, **self.solverOpts)
if self.verbose: print 'Done' if self.verbose: print 'Done'
rhs = RHS(tInd, F) rhs = RHS(tInd, F)
if self.verbose: print ' Solving (Adjoint)... (tInd = %d)'%tInd if self.verbose: print ' Solving (Adjoint)... (tInd = {0:d})'.format(tInd)
sol = Ainv * rhs sol = Ainv * rhs
if self.verbose: print ' Done...' if self.verbose: print ' Done...'
if sol.ndim == 1: if sol.ndim == 1:
@@ -123,14 +123,14 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
* Compute \\\(\\\\vec{w} = -\\\mathbf{Q} \\\\vec{y}\\\) * Compute \\\(\\\\vec{w} = -\\\mathbf{Q} \\\\vec{y}\\\)
""" """
if self.verbose: print '%s\nCalculating J(v)\n%s'%('*'*50,'*'*50) if self.verbose: print '{0!s}\nCalculating J(v)\n{1!s}'.format('*'*50, '*'*50)
self.curModel = m self.curModel = m
if f is None: if f is None:
f = self.fields(m) f = self.fields(m)
p = self.Gvec(m, v, f) p = self.Gvec(m, v, f)
y = self.solveAh(m, p) y = self.solveAh(m, p)
Jv = self.survey.evalDeriv(f, v=y) Jv = self.survey.evalDeriv(f, v=y)
if self.verbose: print '%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50) if self.verbose: print '{0!s}\nDone calculating J(v)\n{1!s}'.format('*'*50, '*'*50)
return - mkvc(Jv) return - mkvc(Jv)
def Jtvec(self, m, v, f=None): def Jtvec(self, m, v, f=None):
@@ -148,7 +148,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
* Compute \\\(\\\\vec{w} = -\\\mathbf{G}^\\\\top y\\\) * Compute \\\(\\\\vec{w} = -\\\mathbf{G}^\\\\top y\\\)
""" """
if self.verbose: print '%s\nCalculating J^T(v)\n%s'%('*'*50,'*'*50) if self.verbose: print '{0!s}\nCalculating J^T(v)\n{1!s}'.format('*'*50, '*'*50)
self.curModel = m self.curModel = m
if f is None: if f is None:
f = self.fields(m) f = self.fields(m)
@@ -159,6 +159,6 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
p = self.survey.evalDeriv(f, v=v, adjoint=True) p = self.survey.evalDeriv(f, v=v, adjoint=True)
y = self.solveAht(m, p) y = self.solveAht(m, p)
w = self.Gtvec(m, y, f) w = self.Gtvec(m, y, f)
if self.verbose: print '%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50) if self.verbose: print '{0!s}\nDone calculating J^T(v)\n{1!s}'.format('*'*50, '*'*50)
return - mkvc(w) return - mkvc(w)
+2 -2
View File
@@ -58,7 +58,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
Src.append(EM.FDEM.Src.RawVec([rx0], freq, mesh.getEdgeInnerProduct()*S_m, S_e)) Src.append(EM.FDEM.Src.RawVec([rx0], freq, mesh.getEdgeInnerProduct()*S_m, S_e))
if verbose: if verbose:
print ' Fetching %s problem' % (fdemType) print ' Fetching {0!s} problem'.format((fdemType))
if fdemType == 'e': if fdemType == 'e':
survey = EM.FDEM.Survey(Src) survey = EM.FDEM.Survey(Src)
@@ -94,7 +94,7 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose) prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose)
mesh = prb1.mesh mesh = prb1.mesh
print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp) print 'Cross Checking Forward: {0!s}, {1!s} formulations - {2!s}'.format(fdemType1, fdemType2, comp)
logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY) logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
mu = np.ones(mesh.nC)*MU mu = np.ones(mesh.nC)*MU
@@ -110,7 +110,7 @@ def run(plotIt=True):
# Mesh # Mesh
mesh = Mesh.CylMesh([hx,1.,hz], [0.,0.,-np.sum(hz[:npadzu+ncz-nza])]) mesh = Mesh.CylMesh([hx,1.,hz], [0.,0.,-np.sum(hz[:npadzu+ncz-nza])])
print 'Mesh Extent xmax: %f,: zmin: %f, zmax: %f'%(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max()) print 'Mesh Extent xmax: {0:f},: zmin: {1:f}, zmax: {2:f}'.format(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max())
print 'Number of cells', mesh.nC print 'Number of cells', mesh.nC
if plotIt is True: if plotIt is True:
@@ -98,7 +98,7 @@ def run(plotIt=True, n=60):
ii = int(ii) ii = int(ii)
out = M.plotImage(PHIS[ii][1],ax=ax) out = M.plotImage(PHIS[ii][1],ax=ax)
ax.axis('off') ax.axis('off')
ax.set_title('Elapsed Time: %4.1f'%PHIS[ii][0]) ax.set_title('Elapsed Time: {0:4.1f}'.format(PHIS[ii][0]))
plt.show() plt.show()
if __name__ == '__main__': if __name__ == '__main__':
+3 -3
View File
@@ -29,15 +29,15 @@ def run(plotIt=True, n=60):
axes[0].set_ylim([-1,17]) axes[0].set_ylim([-1,17])
for ii, loc in zip(range(M.nC),M.gridCC): for ii, loc in zip(range(M.nC),M.gridCC):
axes[0].text(loc[0]+0.2,loc[1],'%d'%ii, color='r') axes[0].text(loc[0]+0.2,loc[1],'{0:d}'.format(ii), color='r')
axes[0].plot(M.gridFx[:,0],M.gridFx[:,1], 'g>') axes[0].plot(M.gridFx[:,0],M.gridFx[:,1], 'g>')
for ii, loc in zip(range(M.nFx),M.gridFx): for ii, loc in zip(range(M.nFx),M.gridFx):
axes[0].text(loc[0]+0.2,loc[1],'%d'%ii, color='g') axes[0].text(loc[0]+0.2,loc[1],'{0:d}'.format(ii), color='g')
axes[0].plot(M.gridFy[:,0],M.gridFy[:,1], 'm^') axes[0].plot(M.gridFy[:,0],M.gridFy[:,1], 'm^')
for ii, loc in zip(range(M.nFy),M.gridFy): for ii, loc in zip(range(M.nFy),M.gridFy):
axes[0].text(loc[0]+0.2,loc[1]+0.2,'%d'%(ii+M.nFx), color='m') axes[0].text(loc[0]+0.2,loc[1]+0.2,'{0:d}'.format((ii+M.nFx)), color='m')
axes[1].spy(M.faceDiv) axes[1].spy(M.faceDiv)
axes[1].set_title('Face Divergence') axes[1].set_title('Face Divergence')
+8 -8
View File
@@ -59,7 +59,7 @@ if __name__ == '__main__':
if line == "##### AUTOIMPORTS #####\n": if line == "##### AUTOIMPORTS #####\n":
inimports = not inimports inimports = not inimports
if inimports: if inimports:
out += '\n'.join(["import %s"%_ for _ in exfiles]) out += '\n'.join(["import {0!s}".format(_) for _ in exfiles])
out += '\n\n__examples__ = ["' + '", "'.join(exfiles)+ '"]\n' out += '\n\n__examples__ = ["' + '", "'.join(exfiles)+ '"]\n'
out += '\n##### AUTOIMPORTS #####\n' out += '\n##### AUTOIMPORTS #####\n'
f.close() f.close()
@@ -76,11 +76,11 @@ if __name__ == '__main__':
docstr = runFunction.__doc__ docstr = runFunction.__doc__
if docstr is None: if docstr is None:
doc = '%s\n%s'%(name.replace('_',' '),'='*len(name)) doc = '{0!s}\n{1!s}'.format(name.replace('_',' '), '='*len(name))
else: else:
doc = '\n'.join([_[8:].rstrip() for _ in docstr.split('\n')]) doc = '\n'.join([_[8:].rstrip() for _ in docstr.split('\n')])
out = """.. _examples_%s: out = """.. _examples_{0!s}:
.. --------------------------------- .. .. --------------------------------- ..
.. .. .. ..
@@ -90,21 +90,21 @@ if __name__ == '__main__':
.. .. .. ..
.. --------------------------------- .. .. --------------------------------- ..
%s {1!s}
.. plot:: .. plot::
from SimPEG import Examples from SimPEG import Examples
Examples.%s.run() Examples.{2!s}.run()
.. literalinclude:: ../../../SimPEG/Examples/%s.py .. literalinclude:: ../../../SimPEG/Examples/{3!s}.py
:language: python :language: python
:linenos: :linenos:
"""%(name,doc,name,name) """.format(name, doc, name, name)
rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'content', 'examples', name + '.rst'])) rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'content', 'examples', name + '.rst']))
print 'Creating: %s.rst'%name print 'Creating: {0!s}.rst'.format(name)
f = open(rst, 'w') f = open(rst, 'w')
f.write(out) f.write(out)
f.close() f.close()
+1 -1
View File
@@ -116,7 +116,7 @@ class RichardsMap(object):
ax.semilogx(self.k(h, m), h) ax.semilogx(self.k(h, m), h)
def _assertMatchesPair(self, pair): def _assertMatchesPair(self, pair):
assert isinstance(self, pair), "Mapping object must be an instance of a %s class."%(pair.__name__) assert isinstance(self, pair), "Mapping object must be an instance of a {0!s} class.".format((pair.__name__))
+1 -1
View File
@@ -140,7 +140,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
for ii, dt in enumerate(self.timeSteps): for ii, dt in enumerate(self.timeSteps):
bc = self.getBoundaryConditions(ii, u[ii]) bc = self.getBoundaryConditions(ii, u[ii])
u[ii+1] = self.rootFinder.root(lambda hn1m, return_g=True: self.getResidual(m, u[ii], hn1m, dt, bc, return_g=return_g), u[ii]) u[ii+1] = self.rootFinder.root(lambda hn1m, return_g=True: self.getResidual(m, u[ii], hn1m, dt, bc, return_g=return_g), u[ii])
if self.debug: print "Solving Fields (%4d/%d - %3.1f%% Done) %d Iterations, %4.2f seconds"%(ii+1, self.nT, 100.0*(ii+1)/self.nT, self.rootFinder.iter, time.time() - tic) if self.debug: print "Solving Fields ({0:4d}/{1:d} - {2:3.1f}% Done) {3:d} Iterations, {4:4.2f} seconds".format(ii+1, self.nT, 100.0*(ii+1)/self.nT, self.rootFinder.iter, time.time() - tic)
return u return u
@Utils.timeIt @Utils.timeIt
+4 -4
View File
@@ -37,7 +37,7 @@ class Fields(object):
for f in self.knownFields: for f in self.knownFields:
loc =self.knownFields[f] loc =self.knownFields[f]
sz += np.array(self._storageShape(loc)).prod()*8.0/(1024**2) sz += np.array(self._storageShape(loc)).prod()*8.0/(1024**2)
return "%e MB"%sz return "{0:e} MB".format(sz)
def _storageShape(self, loc): def _storageShape(self, loc):
nSrc = self.survey.nSrc nSrc = self.survey.nSrc
@@ -84,12 +84,12 @@ class Fields(object):
return return
if accessType=='set' and name not in self.knownFields: if accessType=='set' and name not in self.knownFields:
if name in self.aliasFields: if name in self.aliasFields:
raise KeyError("Invalid field name (%s) for setter, you can't set an aliased property"%name) raise KeyError("Invalid field name ({0!s}) for setter, you can't set an aliased property".format(name))
else: else:
raise KeyError('Invalid field name (%s) for setter'%name) raise KeyError('Invalid field name ({0!s}) for setter'.format(name))
elif accessType=='get' and (name not in self.knownFields and name not in self.aliasFields): elif accessType=='get' and (name not in self.knownFields and name not in self.aliasFields):
raise KeyError('Invalid field name (%s) for getter'%name) raise KeyError('Invalid field name ({0!s}) for getter'.format(name))
return name return name
def _indexAndNameFromKey(self, key, accessType): def _indexAndNameFromKey(self, key, accessType):
+7 -7
View File
@@ -101,7 +101,7 @@ class IdentityMap(object):
:return: passed the test? :return: passed the test?
""" """
print 'Testing %s' % str(self) print 'Testing {0!s}'.format(str(self))
if m is None: if m is None:
m = abs(np.random.rand(self.nP)) m = abs(np.random.rand(self.nP))
if 'plotIt' not in kwargs: if 'plotIt' not in kwargs:
@@ -111,21 +111,21 @@ class IdentityMap(object):
def _assertMatchesPair(self, pair): def _assertMatchesPair(self, pair):
assert (isinstance(self, pair) or assert (isinstance(self, pair) or
isinstance(self, ComboMap) and isinstance(self.maps[0], pair) isinstance(self, ComboMap) and isinstance(self.maps[0], pair)
), "Mapping object must be an instance of a %s class."%(pair.__name__) ), "Mapping object must be an instance of a {0!s} class.".format((pair.__name__))
def __mul__(self, val): def __mul__(self, val):
if isinstance(val, IdentityMap): if isinstance(val, IdentityMap):
if not (self.shape[1] == '*' or val.shape[0] == '*') and not self.shape[1] == val.shape[0]: if not (self.shape[1] == '*' or val.shape[0] == '*') and not self.shape[1] == val.shape[0]:
raise ValueError('Dimension mismatch in %s and %s.' % (str(self), str(val))) raise ValueError('Dimension mismatch in {0!s} and {1!s}.'.format(str(self), str(val)))
return ComboMap([self, val]) return ComboMap([self, val])
elif isinstance(val, np.ndarray): elif isinstance(val, np.ndarray):
if not self.shape[1] == '*' and not self.shape[1] == val.shape[0]: if not self.shape[1] == '*' and not self.shape[1] == val.shape[0]:
raise ValueError('Dimension mismatch in %s and np.ndarray%s.' % (str(self), str(val.shape))) raise ValueError('Dimension mismatch in {0!s} and np.ndarray{1!s}.'.format(str(self), str(val.shape)))
return self._transform(val) return self._transform(val)
raise Exception('Unrecognized data type to multiply. Try a map or a numpy.ndarray!') raise Exception('Unrecognized data type to multiply. Try a map or a numpy.ndarray!')
def __str__(self): def __str__(self):
return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0], self.shape[1]) return "{0!s}({1!s},{2!s})".format(self.__class__.__name__, self.shape[0], self.shape[1])
class ComboMap(IdentityMap): class ComboMap(IdentityMap):
@@ -140,7 +140,7 @@ class ComboMap(IdentityMap):
if ii > 0 and not (self.shape[1] == '*' or m.shape[0] == '*') and not self.shape[1] == m.shape[0]: if ii > 0 and not (self.shape[1] == '*' or m.shape[0] == '*') and not self.shape[1] == m.shape[0]:
prev = self.maps[-1] prev = self.maps[-1]
errArgs = (prev.__class__.__name__, prev.shape[0], prev.shape[1], m.__class__.__name__, m.shape[0], m.shape[1]) errArgs = (prev.__class__.__name__, prev.shape[0], prev.shape[1], m.__class__.__name__, m.shape[0], m.shape[1])
raise ValueError('Dimension mismatch in map[%s] (%s, %s) and map[%s] (%s, %s).' % errArgs) raise ValueError('Dimension mismatch in map[{0!s}] ({1!s}, {2!s}) and map[{3!s}] ({4!s}, {5!s}).'.format(*errArgs))
if isinstance(m, ComboMap): if isinstance(m, ComboMap):
self.maps += m.maps self.maps += m.maps
@@ -173,7 +173,7 @@ class ComboMap(IdentityMap):
return deriv return deriv
def __str__(self): def __str__(self):
return 'ComboMap[%s](%s,%s)' % (' * '.join([m.__str__() for m in self.maps]), self.shape[0], self.shape[1]) return 'ComboMap[{0!s}]({1!s},{2!s})'.format(' * '.join([m.__str__() for m in self.maps]), self.shape[0], self.shape[1])
class ExpMap(IdentityMap): class ExpMap(IdentityMap):
+1 -1
View File
@@ -522,7 +522,7 @@ class BaseRectangularMesh(BaseMesh):
assert xType in outType, 'You cannot change type of components.' assert xType in outType, 'You cannot change type of components.'
if type(x) == list: if type(x) == list:
for i, xi in enumerate(x): for i, xi in enumerate(x):
assert isinstance(x, np.ndarray), "x[%i] must be a numpy array" % i assert isinstance(x, np.ndarray), "x[{0:d}] must be a numpy array".format(i)
assert xi.size == x[0].size, "Number of elements in list must not change." assert xi.size == x[0].size, "Number of elements in list must not change."
x_array = np.ones((x.size, len(x))) x_array = np.ones((x.size, len(x)))
+166 -134
View File
@@ -4,14 +4,28 @@ from DiffOperators import DiffOperators
from InnerProducts import InnerProducts from InnerProducts import InnerProducts
from View import CurvView from View import CurvView
# Some helper functions. # Some helper functions.
length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5 def length2D(x):
length3D = lambda x: (x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)**0.5 return (x[:, 0]**2 + x[:, 1]**2)**0.5
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, CurvView): def length3D(x):
return (x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)**0.5
def normalize2D(x):
return x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2))
def normalize3D(x):
return x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2))
# Curvi Mesh
class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts,
CurvView):
""" """
CurvilinearMesh is a mesh class that deals with curvilinear meshes. CurvilinearMesh is a mesh class that deals with curvilinear meshes.
@@ -31,12 +45,16 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie
_meshType = 'Curv' _meshType = 'Curv'
def __init__(self, nodes): def __init__(self, nodes):
assert type(nodes) == list, "'nodes' variable must be a list of np.ndarray" assert type(nodes) == list, ("'nodes' variable must be a list of "
"np.ndarray")
assert len(nodes) > 1, "len(node) must be greater than 1" assert len(nodes) > 1, "len(node) must be greater than 1"
for i, nodes_i in enumerate(nodes): for i, nodes_i in enumerate(nodes):
assert isinstance(nodes_i, np.ndarray), ("nodes[%i] is not a numpy array." % i) assert isinstance(nodes_i, np.ndarray), ("nodes[{0:d}] is not a"
assert nodes_i.shape == nodes[0].shape, ("nodes[%i] is not the same shape as nodes[0]" % i) "numpy array.".format(i))
assert nodes_i.shape == nodes[0].shape, ("nodes[{0:d}] is not the "
"same shape as nodes[0]"
.format(i))
assert len(nodes[0].shape) == len(nodes), "Dimension mismatch" assert len(nodes[0].shape) == len(nodes), "Dimension mismatch"
assert len(nodes[0].shape) > 1, "Not worth using Curv for a 1D mesh." assert len(nodes[0].shape) > 1, "Not worth using Curv for a 1D mesh."
@@ -48,80 +66,79 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie
for i, node_i in enumerate(nodes): for i, node_i in enumerate(nodes):
self._gridN[:, i] = Utils.mkvc(node_i.astype(float)) self._gridN[:, i] = Utils.mkvc(node_i.astype(float))
def gridCC(): @property
doc = "Cell-centered grid." def gridCC(self):
"""
def fget(self): Cell-centered grid
if self._gridCC is None: """
self._gridCC = np.concatenate([self.aveN2CC*self.gridN[:,i] for i in range(self.dim)]).reshape((-1,self.dim), order='F') if getattr(self, '_gridCC', None) is None:
self._gridCC = np.concatenate([self.aveN2CC*self.gridN[:, i]
for i in range(self.dim)]).reshape(
(-1, self.dim), order='F')
return self._gridCC return self._gridCC
return locals()
_gridCC = None # Store grid by default
gridCC = property(**gridCC())
def gridN(): @property
doc = "Nodal grid." def gridN(self):
"""
def fget(self): Nodal grid.
if self._gridN is None: """
if getattr(self, '_gridN', None) is None:
raise Exception("Someone deleted this. I blame you.") raise Exception("Someone deleted this. I blame you.")
return self._gridN return self._gridN
return locals()
_gridN = None # Store grid by default
gridN = property(**gridN())
def gridFx(): @property
doc = "Face staggered grid in the x direction." def gridFx(self):
"""
Face staggered grid in the x direction.
"""
def fget(self): if getattr(self, '_gridFx', None) is None:
if self._gridFx is None:
N = self.r(self.gridN, 'N', 'N', 'M') N = self.r(self.gridN, 'N', 'N', 'M')
if self.dim == 2: if self.dim == 2:
XY = [Utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N] XY = [Utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N]
self._gridFx = np.c_[XY[0], XY[1]] self._gridFx = np.c_[XY[0], XY[1]]
elif self.dim == 3: elif self.dim == 3:
XYZ = [Utils.mkvc(0.25 * (n[:, :-1, :-1] + n[:, :-1, 1:] + n[:, 1:, :-1] + n[:, 1:, 1:])) for n in N] XYZ = [Utils.mkvc(0.25 * (n[:, :-1, :-1] + n[:, :-1, 1:] +
n[:, 1:, :-1] + n[:, 1:, 1:])) for n in N]
self._gridFx = np.c_[XYZ[0], XYZ[1], XYZ[2]] self._gridFx = np.c_[XYZ[0], XYZ[1], XYZ[2]]
return self._gridFx return self._gridFx
return locals()
_gridFx = None # Store grid by default
gridFx = property(**gridFx())
def gridFy(): @property
doc = "Face staggered grid in the y direction." def gridFy(self):
"""
Face staggered grid in the y direction.
"""
def fget(self): if getattr(self, '_gridFy', None) is None:
if self._gridFy is None:
N = self.r(self.gridN, 'N', 'N', 'M') N = self.r(self.gridN, 'N', 'N', 'M')
if self.dim == 2: if self.dim == 2:
XY = [Utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N] XY = [Utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N]
self._gridFy = np.c_[XY[0], XY[1]] self._gridFy = np.c_[XY[0], XY[1]]
elif self.dim == 3: elif self.dim == 3:
XYZ = [Utils.mkvc(0.25 * (n[:-1, :, :-1] + n[:-1, :, 1:] + n[1:, :, :-1] + n[1:, :, 1:])) for n in N] XYZ = [Utils.mkvc(0.25 * (n[:-1, :, :-1] + n[:-1, :, 1:] +
n[1:, :, :-1] + n[1:, :, 1:])) for n in N]
self._gridFy = np.c_[XYZ[0], XYZ[1], XYZ[2]] self._gridFy = np.c_[XYZ[0], XYZ[1], XYZ[2]]
return self._gridFy return self._gridFy
return locals()
_gridFy = None # Store grid by default
gridFy = property(**gridFy())
def gridFz(): @property
doc = "Face staggered grid in the z direction." def gridFz(self):
"""
Face staggered grid in the y direction.
"""
def fget(self): if getattr(self, '_gridFz', None) is None:
if self._gridFz is None and self.dim == 3:
N = self.r(self.gridN, 'N', 'N', 'M') N = self.r(self.gridN, 'N', 'N', 'M')
XYZ = [Utils.mkvc(0.25 * (n[:-1, :-1, :] + n[:-1, 1:, :] + n[1:, :-1, :] + n[1:, 1:, :])) for n in N] XYZ = [Utils.mkvc(0.25 * (n[:-1, :-1, :] + n[:-1, 1:, :] +
n[1:, :-1, :] + n[1:, 1:, :])) for n in N]
self._gridFz = np.c_[XYZ[0], XYZ[1], XYZ[2]] self._gridFz = np.c_[XYZ[0], XYZ[1], XYZ[2]]
return self._gridFz return self._gridFz
return locals()
_gridFz = None # Store grid by default
gridFz = property(**gridFz())
def gridEx(): @property
doc = "Edge staggered grid in the x direction." def gridEx(self):
"""
def fget(self): Edge staggered grid in the x direction.
if self._gridEx is None: """
if getattr(self, '_gridEx', None) is None:
N = self.r(self.gridN, 'N', 'N', 'M') N = self.r(self.gridN, 'N', 'N', 'M')
if self.dim == 2: if self.dim == 2:
XY = [Utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N] XY = [Utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N]
@@ -130,15 +147,13 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie
XYZ = [Utils.mkvc(0.5 * (n[:-1, :, :] + n[1:, :, :])) for n in N] XYZ = [Utils.mkvc(0.5 * (n[:-1, :, :] + n[1:, :, :])) for n in N]
self._gridEx = np.c_[XYZ[0], XYZ[1], XYZ[2]] self._gridEx = np.c_[XYZ[0], XYZ[1], XYZ[2]]
return self._gridEx return self._gridEx
return locals()
_gridEx = None # Store grid by default
gridEx = property(**gridEx())
def gridEy(): @property
doc = "Edge staggered grid in the y direction." def gridEy(self):
"""
def fget(self): Edge staggered grid in the y direction.
if self._gridEy is None: """
if getattr(self, '_gridEy', None) is None:
N = self.r(self.gridN, 'N', 'N', 'M') N = self.r(self.gridN, 'N', 'N', 'M')
if self.dim == 2: if self.dim == 2:
XY = [Utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N] XY = [Utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N]
@@ -147,22 +162,17 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie
XYZ = [Utils.mkvc(0.5 * (n[:, :-1, :] + n[:, 1:, :])) for n in N] XYZ = [Utils.mkvc(0.5 * (n[:, :-1, :] + n[:, 1:, :])) for n in N]
self._gridEy = np.c_[XYZ[0], XYZ[1], XYZ[2]] self._gridEy = np.c_[XYZ[0], XYZ[1], XYZ[2]]
return self._gridEy return self._gridEy
return locals()
_gridEy = None # Store grid by default
gridEy = property(**gridEy())
def gridEz(): @property
doc = "Edge staggered grid in the z direction." def gridEz(self):
"""
def fget(self): Edge staggered grid in the z direction.
if self._gridEz is None and self.dim == 3: """
if getattr(self, '_gridEz', None) is None and self.dim == 3:
N = self.r(self.gridN, 'N', 'N', 'M') N = self.r(self.gridN, 'N', 'N', 'M')
XYZ = [Utils.mkvc(0.5 * (n[:, :, :-1] + n[:, :, 1:])) for n in N] XYZ = [Utils.mkvc(0.5 * (n[:, :, :-1] + n[:, :, 1:])) for n in N]
self._gridEz = np.c_[XYZ[0], XYZ[1], XYZ[2]] self._gridEz = np.c_[XYZ[0], XYZ[1], XYZ[2]]
return self._gridEz return self._gridEz
return locals()
_gridEz = None # Store grid by default
gridEz = property(**gridEz())
# --------------- Geometries --------------------- # --------------- Geometries ---------------------
# #
@@ -194,19 +204,25 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie
# | / | / # | / | /
# D -------------- C # D -------------- C
# node(i+1,j,k) node(i+1,j+1,k) # node(i+1,j,k) node(i+1,j+1,k)
def vol():
doc = "Construct cell volumes of the 3D model as 1d array."
def fget(self): @property
if(self._vol is None): def vol(self):
"""
Construct cell volumes of the 3D model as 1d array
"""
if getattr(self, '_vol', None) is None:
if self.dim == 2: if self.dim == 2:
A, B, C, D = Utils.indexCube('ABCD', self.vnC+1) A, B, C, D = Utils.indexCube('ABCD', self.vnC+1)
normal, area = Utils.faceInfo(np.c_[self.gridN, np.zeros((self.nN, 1))], A, B, C, D) normal, area = Utils.faceInfo(np.c_[self.gridN, np.zeros(
(self.nN, 1))], A, B, C, D)
self._vol = area self._vol = area
elif self.dim == 3: elif self.dim == 3:
# Each polyhedron can be decomposed into 5 tetrahedrons # Each polyhedron can be decomposed into 5 tetrahedrons
# However, this presents a choice so we may as well divide in two ways and average. # However, this presents a choice so we may as well divide in
A, B, C, D, E, F, G, H = Utils.indexCube('ABCDEFGH', self.vnC+1) # two ways and average.
A, B, C, D, E, F, G, H = Utils.indexCube('ABCDEFGH', self.vnC +
1)
vol1 = (Utils.volTetra(self.gridN, A, B, D, E) + # cutted edge top vol1 = (Utils.volTetra(self.gridN, A, B, D, E) + # cutted edge top
Utils.volTetra(self.gridN, B, E, F, G) + # cutted edge top Utils.volTetra(self.gridN, B, E, F, G) + # cutted edge top
@@ -222,50 +238,60 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie
self._vol = (vol1 + vol2)/2 self._vol = (vol1 + vol2)/2
return self._vol return self._vol
return locals()
_vol = None
vol = property(**vol())
def area(): @property
doc = "Face areas." def area(self):
if (getattr(self, '_area', None) is None or
def fget(self): getattr(self, '_normals', None) is None):
if(self._area is None or self._normals is None):
# Compute areas of cell faces # Compute areas of cell faces
if(self.dim == 2): if(self.dim == 2):
xy = self.gridN xy = self.gridN
A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy])) A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx,
self.nCy]))
edge1 = xy[B, :] - xy[A, :] edge1 = xy[B, :] - xy[A, :]
normal1 = np.c_[edge1[:, 1], -edge1[:, 0]] normal1 = np.c_[edge1[:, 1], -edge1[:, 0]]
area1 = length2D(edge1) area1 = length2D(edge1)
A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy])) A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx,
# Note that we are doing A-D to make sure the normal points the right way. self.nNy]))
# Think about it. Look at the picture. Normal points towards C iff you do this. # Note that we are doing A-D to make sure the normal points the
# right way.
# Think about it. Look at the picture. Normal points towards C
# iff you do this.
edge2 = xy[A, :] - xy[D, :] edge2 = xy[A, :] - xy[D, :]
normal2 = np.c_[edge2[:, 1], -edge2[:, 0]] normal2 = np.c_[edge2[:, 1], -edge2[:, 0]]
area2 = length2D(edge2) area2 = length2D(edge2)
self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2)] self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2)]
self._normals = [normalize2D(normal1), normalize2D(normal2)] self._normals = [normalize2D(normal1), normalize2D(normal2)]
elif(self.dim == 3): elif(self.dim == 3):
A, E, F, B = Utils.indexCube('AEFB', self.vnC+1, np.array([self.nNx, self.nCy, self.nCz])) A, E, F, B = Utils.indexCube('AEFB', self.vnC+1, np.array(
normal1, area1 = Utils.faceInfo(self.gridN, A, E, F, B, average=False, normalizeNormals=False) [self.nNx, self.nCy, self.nCz]))
normal1, area1 = Utils.faceInfo(self.gridN, A, E, F, B,
average=False,
normalizeNormals=False)
A, D, H, E = Utils.indexCube('ADHE', self.vnC+1, np.array([self.nCx, self.nNy, self.nCz])) A, D, H, E = Utils.indexCube('ADHE', self.vnC+1, np.array(
normal2, area2 = Utils.faceInfo(self.gridN, A, D, H, E, average=False, normalizeNormals=False) [self.nCx, self.nNy, self.nCz]))
normal2, area2 = Utils.faceInfo(self.gridN, A, D, H, E,
average=False,
normalizeNormals=False)
A, B, C, D = Utils.indexCube('ABCD', self.vnC+1, np.array([self.nCx, self.nCy, self.nNz])) A, B, C, D = Utils.indexCube('ABCD', self.vnC+1, np.array(
normal3, area3 = Utils.faceInfo(self.gridN, A, B, C, D, average=False, normalizeNormals=False) [self.nCx, self.nCy, self.nNz]))
normal3, area3 = Utils.faceInfo(self.gridN, A, B, C, D,
average=False,
normalizeNormals=False)
self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2), Utils.mkvc(area3)] self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2),
Utils.mkvc(area3)]
self._normals = [normal1, normal2, normal3] self._normals = [normal1, normal2, normal3]
return self._area return self._area
return locals()
_area = None
area = property(**area())
def normals(): @property
doc = """Face normals: calling this will average def normals(self):
"""
Face normals: calling this will average
the computed normals so that there is one the computed normals so that there is one
per face. This is especially relevant in per face. This is especially relevant in
3D, as there are up to 4 different normals 3D, as there are up to 4 different normals
@@ -276,8 +302,7 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie
NyX, NyY, NyZ = M.r(M.normals, 'F', 'Fy', 'M') NyX, NyY, NyZ = M.r(M.normals, 'F', 'Fy', 'M')
""" """
def fget(self): if getattr(self, '_normals', None) is None:
if(self._normals is None):
self.area # calling .area will create the face normals self.area # calling .area will create the face normals
if self.dim == 2: if self.dim == 2:
return normalize2D(np.r_[self._normals[0], self._normals[1]]) return normalize2D(np.r_[self._normals[0], self._normals[1]])
@@ -286,48 +311,55 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie
normal2 = (self._normals[1][0] + self._normals[1][1] + self._normals[1][2] + self._normals[1][3])/4 normal2 = (self._normals[1][0] + self._normals[1][1] + self._normals[1][2] + self._normals[1][3])/4
normal3 = (self._normals[2][0] + self._normals[2][1] + self._normals[2][2] + self._normals[2][3])/4 normal3 = (self._normals[2][0] + self._normals[2][1] + self._normals[2][2] + self._normals[2][3])/4
return normalize3D(np.r_[normal1, normal2, normal3]) return normalize3D(np.r_[normal1, normal2, normal3])
return locals()
_normals = None
normals = property(**normals())
def edge(): @property
doc = "Edge legnths." def edge(self):
"""
def fget(self): Edge lengths
if(self._edge is None or self._tangents is None): """
if getattr(self, '_edge', None) is None:
if(self.dim == 2): if(self.dim == 2):
xy = self.gridN xy = self.gridN
A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy])) A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx,
self.nNy]))
edge1 = xy[D, :] - xy[A, :] edge1 = xy[D, :] - xy[A, :]
A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy])) A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx,
self.nCy]))
edge2 = xy[B, :] - xy[A, :] edge2 = xy[B, :] - xy[A, :]
self._edge = np.r_[Utils.mkvc(length2D(edge1)), Utils.mkvc(length2D(edge2))] self._edge = np.r_[Utils.mkvc(length2D(edge1)),
self._tangents = np.r_[edge1, edge2]/np.c_[self._edge, self._edge] Utils.mkvc(length2D(edge2))]
self._tangents = np.r_[edge1, edge2]/np.c_[self._edge,
self._edge]
elif(self.dim == 3): elif(self.dim == 3):
xyz = self.gridN xyz = self.gridN
A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx, self.nNy, self.nNz])) A, D = Utils.indexCube('AD', self.vnC+1, np.array([self.nCx,
self.nNy,
self.nNz]))
edge1 = xyz[D, :] - xyz[A, :] edge1 = xyz[D, :] - xyz[A, :]
A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx, self.nCy, self.nNz])) A, B = Utils.indexCube('AB', self.vnC+1, np.array([self.nNx,
self.nCy,
self.nNz]))
edge2 = xyz[B, :] - xyz[A, :] edge2 = xyz[B, :] - xyz[A, :]
A, E = Utils.indexCube('AE', self.vnC+1, np.array([self.nNx, self.nNy, self.nCz])) A, E = Utils.indexCube('AE', self.vnC+1, np.array([self.nNx,
self.nNy,
self.nCz]))
edge3 = xyz[E, :] - xyz[A, :] edge3 = xyz[E, :] - xyz[A, :]
self._edge = np.r_[Utils.mkvc(length3D(edge1)), Utils.mkvc(length3D(edge2)), Utils.mkvc(length3D(edge3))] self._edge = np.r_[Utils.mkvc(length3D(edge1)),
self._tangents = np.r_[edge1, edge2, edge3]/np.c_[self._edge, self._edge, self._edge] Utils.mkvc(length3D(edge2)),
Utils.mkvc(length3D(edge3))]
self._tangents = (np.r_[edge1, edge2, edge3] /
np.c_[self._edge, self._edge, self._edge])
return self._edge
return self._edge return self._edge
return locals()
_edge = None
edge = property(**edge())
def tangents(): @property
doc = "Edge tangents." def tangents(self):
"""
def fget(self): Edge tangents
if(self._tangents is None): """
if getattr(self, '_tangents', None) is None:
self.edge # calling .edge will create the tangents self.edge # calling .edge will create the tangents
return self._tangents return self._tangents
return locals()
_tangents = None
tangents = property(**tangents())
+222 -159
View File
@@ -18,13 +18,15 @@ def checkBC(bc):
for bc_i in bc: for bc_i in bc:
assert type(bc_i) is str, "each bc must be a string" assert type(bc_i) is str, "each bc must be a string"
assert bc_i in ['dirichlet', 'neumann'], "each bc must be either, 'dirichlet' or 'neumann'" assert bc_i in ['dirichlet', 'neumann'], ("each bc must be either,"
"'dirichlet' or 'neumann'")
return bc return bc
def ddxCellGrad(n, bc): def ddxCellGrad(n, bc):
""" """
Create 1D derivative operator from cell-centers to nodes this means we go from n to n+1 Create 1D derivative operator from cell-centers to nodes this means we
go from n to n+1
For Cell-Centered **Dirichlet**, use a ghost point:: For Cell-Centered **Dirichlet**, use a ghost point::
@@ -52,7 +54,8 @@ def ddxCellGrad(n, bc):
""" """
bc = checkBC(bc) bc = checkBC(bc)
D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n, format="csr") D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n,
format="csr")
# Set the first side # Set the first side
if(bc[0] == 'dirichlet'): if(bc[0] == 'dirichlet'):
D[0, 0] = 2 D[0, 0] = 2
@@ -65,10 +68,11 @@ def ddxCellGrad(n, bc):
D[-1, -1] = 0 D[-1, -1] = 0
return D return D
def ddxCellGradBC(n, bc): def ddxCellGradBC(n, bc):
""" """
Create 1D derivative operator from cell-centers to nodes this means we
Create 1D derivative operator from cell-centers to nodes this means we go from n to n+1 go from n to n+1
For Cell-Centered **Dirichlet**, use a ghost point:: For Cell-Centered **Dirichlet**, use a ghost point::
@@ -121,14 +125,16 @@ class DiffOperators(object):
Class creates the differential operators that you need! Class creates the differential operators that you need!
""" """
def __init__(self): def __init__(self):
raise Exception('DiffOperators is a base class providing differential operators on meshes and cannot run on its own. Inherit to your favorite Mesh class.') raise Exception('DiffOperators is a base class providing differential'
'operators on meshes and cannot run on its own.'
'Inherit to your favorite Mesh class.')
def faceDiv(): @property
doc = "Construct divergence operator (face-stg to cell-centres)." def faceDiv(self):
"""
def fget(self): Construct divergence operator (face-stg to cell-centres).
if(self._faceDiv is None): """
# The number of cell centers in each direction if getattr(self, '_faceDiv', None) is None:
n = self.vnC n = self.vnC
# Compute faceDivergence operator on faces # Compute faceDivergence operator on faces
if(self.dim == 1): if(self.dim == 1):
@@ -146,17 +152,15 @@ class DiffOperators(object):
S = self.area S = self.area
V = self.vol V = self.vol
self._faceDiv = sdiag(1/V)*D*sdiag(S) self._faceDiv = sdiag(1/V)*D*sdiag(S)
return self._faceDiv return self._faceDiv
return locals()
_faceDiv = None
faceDiv = property(**faceDiv())
def faceDivx(): @property
doc = "Construct divergence operator in the x component (face-stg to cell-centres)." def faceDivx(self):
"""
def fget(self): Construct divergence operator in the x component (face-stg to
if(self._faceDivx is None): cell-centres).
"""
if getattr(self, '_faceDivx', None) is None:
# The number of cell centers in each direction # The number of cell centers in each direction
n = self.vnC n = self.vnC
# Compute faceDivergence operator on faces # Compute faceDivergence operator on faces
@@ -172,16 +176,12 @@ class DiffOperators(object):
self._faceDivx = sdiag(1/V)*D1*sdiag(S) self._faceDivx = sdiag(1/V)*D1*sdiag(S)
return self._faceDivx return self._faceDivx
return locals()
_faceDivx = None
faceDivx = property(**faceDivx())
def faceDivy(): @property
doc = "Construct divergence operator in the y component (face-stg to cell-centres)." def faceDivy(self):
if(self.dim < 2):
def fget(self): return None
if(self.dim < 2): return None if getattr(self, '_faceDivy', None) is None:
if(self._faceDivy is None):
# The number of cell centers in each direction # The number of cell centers in each direction
n = self.vnC n = self.vnC
# Compute faceDivergence operator on faces # Compute faceDivergence operator on faces
@@ -193,18 +193,17 @@ class DiffOperators(object):
S = self.r(self.area, 'F', 'Fy', 'V') S = self.r(self.area, 'F', 'Fy', 'V')
V = self.vol V = self.vol
self._faceDivy = sdiag(1/V)*D2*sdiag(S) self._faceDivy = sdiag(1/V)*D2*sdiag(S)
return self._faceDivy return self._faceDivy
return locals()
_faceDivy = None
faceDivy = property(**faceDivy())
def faceDivz(): @property
doc = "Construct divergence operator in the z component (face-stg to cell-centres)." def faceDivz(self):
"""
def fget(self): Construct divergence operator in the z component (face-stg to
if(self.dim < 3): return None cell-centres).
if(self._faceDivz is None): """
if(self.dim < 3):
return None
if getattr(self, '_faceDivz', None) is None:
# The number of cell centers in each direction # The number of cell centers in each direction
n = self.vnC n = self.vnC
# Compute faceDivergence operator on faces # Compute faceDivergence operator on faces
@@ -213,17 +212,14 @@ class DiffOperators(object):
S = self.r(self.area, 'F', 'Fz', 'V') S = self.r(self.area, 'F', 'Fz', 'V')
V = self.vol V = self.vol
self._faceDivz = sdiag(1/V)*D3*sdiag(S) self._faceDivz = sdiag(1/V)*D3*sdiag(S)
return self._faceDivz return self._faceDivz
return locals()
_faceDivz = None
faceDivz = property(**faceDivz())
def nodalGrad(): @property
doc = "Construct gradient operator (nodes to edges)." def nodalGrad(self):
"""
def fget(self): Construct gradient operator (nodes to edges).
if(self._nodalGrad is None): """
if getattr(self, '_nodalGrad', None) is None:
# The number of cell centers in each direction # The number of cell centers in each direction
n = self.vnC n = self.vnC
# Compute divergence operator on faces # Compute divergence operator on faces
@@ -242,15 +238,13 @@ class DiffOperators(object):
L = self.edge L = self.edge
self._nodalGrad = sdiag(1/L)*G self._nodalGrad = sdiag(1/L)*G
return self._nodalGrad return self._nodalGrad
return locals()
_nodalGrad = None
nodalGrad = property(**nodalGrad())
def nodalLaplacian(): @property
doc = "Construct laplacian operator (nodes to edges)." def nodalLaplacian(self):
"""
def fget(self): Construct laplacian operator (nodes to edges).
if(self._nodalLaplacian is None): """
if getattr(self, '_nodalLaplacian', None) is None:
print 'Warning: Laplacian has not been tested rigorously.' print 'Warning: Laplacian has not been tested rigorously.'
# The number of cell centers in each direction # The number of cell centers in each direction
n = self.vnC n = self.vnC
@@ -274,22 +268,23 @@ class DiffOperators(object):
L = L1 + L2 + L3 L = L1 + L2 + L3
self._nodalLaplacian = L self._nodalLaplacian = L
return self._nodalLaplacian return self._nodalLaplacian
return locals()
_nodalLaplacian = None
nodalLaplacian = property(**nodalLaplacian())
def setCellGradBC(self, BC): def setCellGradBC(self, BC):
""" """
Function that sets the boundary conditions for cell-centred derivative operators. Function that sets the boundary conditions for cell-centred derivative
operators.
Examples:: Examples::
# Neumann in all directions
BC = 'neumann'
BC = 'neumann' # Neumann in all directions # 3D, Dirichlet in y Neumann else
BC = ['neumann', 'dirichlet', 'neumann'] # 3D, Dirichlet in y Neumann else BC = ['neumann', 'dirichlet', 'neumann']
BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet'] # 3D, Neumann in x on bottom of domain,
# Dirichlet else
# 3D, Neumann in x on bottom of domain, Dirichlet else
BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet']
""" """
if(type(BC) is str): if(type(BC) is str):
BC = [BC]*self.dim BC = [BC]*self.dim
if(type(BC) is list): if(type(BC) is list):
@@ -323,26 +318,24 @@ class DiffOperators(object):
G = sp.vstack((G1, G2, G3), format="csr") G = sp.vstack((G1, G2, G3), format="csr")
return G return G
def cellGrad(): @property
doc = "The cell centered Gradient, takes you to cell faces." def cellGrad(self):
"""
def fget(self): The cell centered Gradient, takes you to cell faces.
if(self._cellGrad is None): """
if getattr(self, '_cellGrad', None) is None:
G = self._cellGradStencil() G = self._cellGradStencil()
# Compute areas of cell faces & volumes S = self.area # Compute areas of cell faces & volumes
S = self.area
V = self.aveCC2F*self.vol # Average volume between adjacent cells V = self.aveCC2F*self.vol # Average volume between adjacent cells
self._cellGrad = sdiag(S/V)*G self._cellGrad = sdiag(S/V)*G
return self._cellGrad return self._cellGrad
return locals()
_cellGrad = None
cellGrad = property(**cellGrad())
def cellGradBC(): @property
doc = "The cell centered Gradient boundary condition matrix" def cellGradBC(self):
"""
def fget(self): The cell centered Gradient boundary condition matrix
if(self._cellGradBC is None): """
if getattr(self, '_cellGradBC', None) is None:
BC = self.setCellGradBC(self._cellGradBC_list) BC = self.setCellGradBC(self._cellGradBC_list)
n = self.vnC n = self.vnC
if(self.dim == 1): if(self.dim == 1):
@@ -361,9 +354,33 @@ class DiffOperators(object):
V = self.aveCC2F*self.vol # Average volume between adjacent cells V = self.aveCC2F*self.vol # Average volume between adjacent cells
self._cellGradBC = sdiag(S/V)*G self._cellGradBC = sdiag(S/V)*G
return self._cellGradBC return self._cellGradBC
return locals()
_cellGradBC = None # def cellGradBC():
cellGradBC = property(**cellGradBC()) # doc = "The cell centered Gradient boundary condition matrix"
# def fget(self):
# if(self._cellGradBC is None):
# BC = self.setCellGradBC(self._cellGradBC_list)
# n = self.vnC
# if(self.dim == 1):
# G = ddxCellGradBC(n[0], BC[0])
# elif(self.dim == 2):
# G1 = sp.kron(speye(n[1]), ddxCellGradBC(n[0], BC[0]))
# G2 = sp.kron(ddxCellGradBC(n[1], BC[1]), speye(n[0]))
# G = sp.block_diag((G1, G2), format="csr")
# elif(self.dim == 3):
# G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGradBC(n[0], BC[0]))
# G2 = kron3(speye(n[2]), ddxCellGradBC(n[1], BC[1]), speye(n[0]))
# G3 = kron3(ddxCellGradBC(n[2], BC[2]), speye(n[1]), speye(n[0]))
# G = sp.block_diag((G1, G2, G3), format="csr")
# # Compute areas of cell faces & volumes
# S = self.area
# V = self.aveCC2F*self.vol # Average volume between adjacent cells
# self._cellGradBC = sdiag(S/V)*G
# return self._cellGradBC
# return locals()
# _cellGradBC = None
# cellGradBC = property(**cellGradBC())
def _cellGradxStencil(self): def _cellGradxStencil(self):
BC = ['neumann', 'neumann'] BC = ['neumann', 'neumann']
@@ -376,11 +393,12 @@ class DiffOperators(object):
G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC)) G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC))
return G1 return G1
@property
def cellGradx(): def cellGradx(self):
doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." """
Cell centered Gradient in the x dimension. Has neumann boundary
def fget(self): conditions.
"""
if getattr(self, '_cellGradx', None) is None: if getattr(self, '_cellGradx', None) is None:
G1 = self._cellGradxStencil() G1 = self._cellGradxStencil()
# Compute areas of cell faces & volumes # Compute areas of cell faces & volumes
@@ -388,8 +406,6 @@ class DiffOperators(object):
L = self.r(self.area/V, 'F','Fx', 'V') L = self.r(self.area/V, 'F','Fx', 'V')
self._cellGradx = sdiag(L)*G1 self._cellGradx = sdiag(L)*G1
return self._cellGradx return self._cellGradx
return locals()
cellGradx = property(**cellGradx())
def _cellGradyStencil(self): def _cellGradyStencil(self):
if self.dim < 2: return None if self.dim < 2: return None
@@ -401,10 +417,10 @@ class DiffOperators(object):
G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0])) G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0]))
return G2 return G2
def cellGrady(): @property
doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." def cellGrady(self):
def fget(self): if self.dim < 2:
if self.dim < 2: return None return None
if getattr(self, '_cellGrady', None) is None: if getattr(self, '_cellGrady', None) is None:
G2 = self._cellGradyStencil() G2 = self._cellGradyStencil()
# Compute areas of cell faces & volumes # Compute areas of cell faces & volumes
@@ -412,8 +428,6 @@ class DiffOperators(object):
L = self.r(self.area/V, 'F', 'Fy', 'V') L = self.r(self.area/V, 'F', 'Fy', 'V')
self._cellGrady = sdiag(L)*G2 self._cellGrady = sdiag(L)*G2
return self._cellGrady return self._cellGrady
return locals()
cellGrady = property(**cellGrady())
def _cellGradzStencil(self): def _cellGradzStencil(self):
if self.dim < 3: return None if self.dim < 3: return None
@@ -422,10 +436,14 @@ class DiffOperators(object):
G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0])) G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0]))
return G3 return G3
def cellGradz(): @property
doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." def cellGradz(self):
def fget(self): """
if self.dim < 3: return None Cell centered Gradient in the x dimension. Has neumann boundary
conditions.
"""
if self.dim < 3:
return None
if getattr(self, '_cellGradz', None) is None: if getattr(self, '_cellGradz', None) is None:
G3 = self._cellGradzStencil() G3 = self._cellGradzStencil()
# Compute areas of cell faces & volumes # Compute areas of cell faces & volumes
@@ -433,23 +451,18 @@ class DiffOperators(object):
L = self.r(self.area/V, 'F', 'Fz', 'V') L = self.r(self.area/V, 'F', 'Fz', 'V')
self._cellGradz = sdiag(L)*G3 self._cellGradz = sdiag(L)*G3
return self._cellGradz return self._cellGradz
return locals()
cellGradz = property(**cellGradz())
def edgeCurl(): @property
doc = "Construct the 3D curl operator." def edgeCurl(self):
"""
def fget(self): Construct the 3D curl operator.
if(self._edgeCurl is None): """
if getattr(self, '_edgeCurl', None) is None:
assert self.dim > 1, "Edge Curl only programed for 2 or 3D." assert self.dim > 1, "Edge Curl only programed for 2 or 3D."
# The number of cell centers in each direction
n = self.vnC
# Compute lengths of cell edges n = self.vnC # The number of cell centers in each direction
L = self.edge L = self.edge # Compute lengths of cell edges
S = self.area # Compute areas of cell faces
# Compute areas of cell faces
S = self.area
# Compute divergence operator on faces # Compute divergence operator on faces
if self.dim == 2: if self.dim == 2:
@@ -477,11 +490,7 @@ class DiffOperators(object):
sp.hstack((-D21, D12, O3))), format="csr") sp.hstack((-D21, D12, O3))), format="csr")
self._edgeCurl = sdiag(1/S)*(C*sdiag(L)) self._edgeCurl = sdiag(1/S)*(C*sdiag(L))
return self._edgeCurl return self._edgeCurl
return locals()
_edgeCurl = None
edgeCurl = property(**edgeCurl())
def getBCProjWF(self, BC, discretization='CC'): def getBCProjWF(self, BC, discretization='CC'):
""" """
@@ -489,16 +498,19 @@ class DiffOperators(object):
The weak form boundary condition projection matrices. The weak form boundary condition projection matrices.
Examples:: Examples::
# Neumann in all directions
BC = 'neumann'
BC = 'neumann' # Neumann in all directions # 3D, Dirichlet in y Neumann else
BC = ['neumann', 'dirichlet', 'neumann'] # 3D, Dirichlet in y Neumann else BC = ['neumann', 'dirichlet', 'neumann']
BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet'] # 3D, Neumann in x on bottom of domain,
# Dirichlet else
# 3D, Neumann in x on bottom of domain, Dirichlet else
BC = [['neumann', 'dirichlet'], 'dirichlet', 'dirichlet']
""" """
if discretization is not 'CC': if discretization is not 'CC':
raise NotImplementedError('Boundary conditions only implemented for CC discretization.') raise NotImplementedError('Boundary conditions only implemented'
'for CC discretization.')
if(type(BC) is str): if(type(BC) is str):
BC = [BC for _ in self.vnC] # Repeat the str self.dim times BC = [BC for _ in self.vnC] # Repeat the str self.dim times
@@ -510,7 +522,6 @@ class DiffOperators(object):
for i, bc_i in enumerate(BC): for i, bc_i in enumerate(BC):
BC[i] = checkBC(bc_i) BC[i] = checkBC(bc_i)
def projDirichlet(n, bc): def projDirichlet(n, bc):
bc = checkBC(bc) bc = checkBC(bc)
ij = ([0, n], [0, 1]) ij = ([0, n], [0, 1])
@@ -550,6 +561,7 @@ class DiffOperators(object):
Pin = projNeumannIn(n[0], BC[0]) Pin = projNeumannIn(n[0], BC[0])
Pout = projNeumannOut(n[0], BC[0]) Pout = projNeumannOut(n[0], BC[0])
elif(self.dim == 2): elif(self.dim == 2):
Pbc1 = sp.kron(speye(n[1]), projDirichlet(n[0], BC[0])) Pbc1 = sp.kron(speye(n[1]), projDirichlet(n[0], BC[0]))
Pbc2 = sp.kron(projDirichlet(n[1], BC[1]), speye(n[0])) Pbc2 = sp.kron(projDirichlet(n[1], BC[1]), speye(n[0]))
@@ -564,12 +576,14 @@ class DiffOperators(object):
P1 = sp.kron(speye(n[1]), projNeumannOut(n[0], BC[0])) P1 = sp.kron(speye(n[1]), projNeumannOut(n[0], BC[0]))
P2 = sp.kron(projNeumannOut(n[1], BC[1]), speye(n[0])) P2 = sp.kron(projNeumannOut(n[1], BC[1]), speye(n[0]))
Pout = sp.block_diag((P1, P2), format="csr") Pout = sp.block_diag((P1, P2), format="csr")
elif(self.dim == 3): elif(self.dim == 3):
Pbc1 = kron3(speye(n[2]), speye(n[1]), projDirichlet(n[0], BC[0])) 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])) 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])) Pbc3 = kron3(projDirichlet(n[2], BC[2]), speye(n[1]), speye(n[0]))
Pbc = sp.block_diag((Pbc1, Pbc2, Pbc3), format="csr") Pbc = sp.block_diag((Pbc1, Pbc2, Pbc3), format="csr")
indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3]), (indF[4] | indF[5])] indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3]), (indF[4] |
indF[5])]
Pbc = Pbc*sdiag(self.area[indF]) Pbc = Pbc*sdiag(self.area[indF])
P1 = kron3(speye(n[2]), speye(n[1]), projNeumannIn(n[0], BC[0])) P1 = kron3(speye(n[2]), speye(n[1]), projNeumannIn(n[0], BC[0]))
@@ -586,15 +600,13 @@ class DiffOperators(object):
def getBCProjWF_simple(self, discretization='CC'): def getBCProjWF_simple(self, discretization='CC'):
""" """
The weak form boundary condition projection matrices The weak form boundary condition projection matrices
when mixed boundary condition is used when mixed boundary condition is used
""" """
if discretization is not 'CC': if discretization is not 'CC':
raise NotImplementedError('Boundary conditions only implemented for CC discretization.') raise NotImplementedError('Boundary conditions only implemented'
'for CC discretization.')
def projBC(n): def projBC(n):
ij = ([0, n], [0, 1]) ij = ([0, n], [0, 1])
@@ -613,9 +625,11 @@ class DiffOperators(object):
vals[1] = 1 vals[1] = 1
return sp.csr_matrix((vals, ij), shape=(n+1, 2)) return sp.csr_matrix((vals, ij), shape=(n+1, 2))
BC = [['dirichlet','dirichlet'],['dirichlet','dirichlet'],['dirichlet','dirichlet']] BC = [['dirichlet', 'dirichlet'], ['dirichlet', 'dirichlet'],
['dirichlet', 'dirichlet']]
n = self.vnC n = self.vnC
indF = self.faceBoundaryInd indF = self.faceBoundaryInd
if(self.dim == 1): if(self.dim == 1):
Pbc = projDirichlet(n[0], BC[0]) Pbc = projDirichlet(n[0], BC[0])
B = projBC(n[0]) B = projBC(n[0])
@@ -653,9 +667,11 @@ class DiffOperators(object):
if(self.dim == 1): if(self.dim == 1):
return self.aveFx2CC return self.aveFx2CC
elif(self.dim == 2): elif(self.dim == 2):
return (0.5)*sp.hstack((self.aveFx2CC, self.aveFy2CC), format="csr") return (0.5)*sp.hstack((self.aveFx2CC, self.aveFy2CC),
format="csr")
elif(self.dim == 3): elif(self.dim == 3):
return (1./3.)*sp.hstack((self.aveFx2CC, self.aveFy2CC, self.aveFz2CC), format="csr") return (1./3.)*sp.hstack((self.aveFx2CC, self.aveFy2CC,
self.aveFz2CC), format="csr")
@property @property
def aveF2CCV(self): def aveF2CCV(self):
@@ -665,11 +681,16 @@ class DiffOperators(object):
elif(self.dim == 2): elif(self.dim == 2):
return sp.block_diag((self.aveFx2CC, self.aveFy2CC), format="csr") return sp.block_diag((self.aveFx2CC, self.aveFy2CC), format="csr")
elif(self.dim == 3): elif(self.dim == 3):
return sp.block_diag((self.aveFx2CC, self.aveFy2CC, self.aveFz2CC), format="csr") return sp.block_diag((self.aveFx2CC, self.aveFy2CC, self.aveFz2CC),
format="csr")
@property @property
def aveFx2CC(self): def aveFx2CC(self):
"Construct the averaging operator on cell faces in the x direction to cell centers." """
Construct the averaging operator on cell faces in the x direction to
cell centers.
"""
if getattr(self, '_aveFx2CC', None) is None: if getattr(self, '_aveFx2CC', None) is None:
n = self.vnC n = self.vnC
if(self.dim == 1): if(self.dim == 1):
@@ -682,8 +703,12 @@ class DiffOperators(object):
@property @property
def aveFy2CC(self): def aveFy2CC(self):
"Construct the averaging operator on cell faces in the y direction to cell centers." """
if self.dim < 2: return None Construct the averaging operator on cell faces in the y direction to
cell centers.
"""
if self.dim < 2:
return None
if getattr(self, '_aveFy2CC', None) is None: if getattr(self, '_aveFy2CC', None) is None:
n = self.vnC n = self.vnC
if(self.dim == 2): if(self.dim == 2):
@@ -694,7 +719,10 @@ class DiffOperators(object):
@property @property
def aveFz2CC(self): def aveFz2CC(self):
"Construct the averaging operator on cell faces in the z direction to cell centers." """
Construct the averaging operator on cell faces in the z direction to
cell centers.
"""
if self.dim < 3: return None if self.dim < 3: return None
if getattr(self, '_aveFz2CC', None) is None: if getattr(self, '_aveFz2CC', None) is None:
n = self.vnC n = self.vnC
@@ -711,12 +739,18 @@ class DiffOperators(object):
if(self.dim == 1): if(self.dim == 1):
self._aveCC2F = avExtrap(n[0]) self._aveCC2F = avExtrap(n[0])
elif(self.dim == 2): elif(self.dim == 2):
self._aveCC2F = sp.vstack((sp.kron(speye(n[1]), avExtrap(n[0])), self._aveCC2F = sp.vstack((sp.kron(speye(n[1]),
sp.kron(avExtrap(n[1]), speye(n[0]))), format="csr") avExtrap(n[0])),
sp.kron(avExtrap(n[1]),
speye(n[0]))), format="csr")
elif(self.dim == 3): elif(self.dim == 3):
self._aveCC2F = sp.vstack((kron3(speye(n[2]), speye(n[1]), avExtrap(n[0])), self._aveCC2F = sp.vstack((kron3(speye(n[2]), speye(n[1]),
kron3(speye(n[2]), avExtrap(n[1]), speye(n[0])), avExtrap(n[0])),
kron3(avExtrap(n[2]), speye(n[1]), speye(n[0]))), format="csr") kron3(speye(n[2]), avExtrap(n[1]),
speye(n[0])),
kron3(avExtrap(n[2]), speye(n[1]),
speye(n[0]))),
format="csr")
return self._aveCC2F return self._aveCC2F
@property @property
@@ -727,7 +761,8 @@ class DiffOperators(object):
elif(self.dim == 2): elif(self.dim == 2):
return 0.5*sp.hstack((self.aveEx2CC, self.aveEy2CC), format="csr") return 0.5*sp.hstack((self.aveEx2CC, self.aveEy2CC), format="csr")
elif(self.dim == 3): elif(self.dim == 3):
return (1./3)*sp.hstack((self.aveEx2CC, self.aveEy2CC, self.aveEz2CC), format="csr") return (1./3)*sp.hstack((self.aveEx2CC, self.aveEy2CC,
self.aveEz2CC), format="csr")
@property @property
def aveE2CCV(self): def aveE2CCV(self):
@@ -737,11 +772,15 @@ class DiffOperators(object):
elif(self.dim == 2): elif(self.dim == 2):
return sp.block_diag((self.aveEx2CC, self.aveEy2CC), format="csr") return sp.block_diag((self.aveEx2CC, self.aveEy2CC), format="csr")
elif(self.dim == 3): elif(self.dim == 3):
return sp.block_diag((self.aveEx2CC, self.aveEy2CC, self.aveEz2CC), format="csr") return sp.block_diag((self.aveEx2CC, self.aveEy2CC, self.aveEz2CC),
format="csr")
@property @property
def aveEx2CC(self): def aveEx2CC(self):
"Construct the averaging operator on cell edges in the x direction to cell centers." """
Construct the averaging operator on cell edges in the x direction to
cell centers.
"""
if getattr(self, '_aveEx2CC', None) is None: if getattr(self, '_aveEx2CC', None) is None:
# The number of cell centers in each direction # The number of cell centers in each direction
n = self.vnC n = self.vnC
@@ -755,8 +794,12 @@ class DiffOperators(object):
@property @property
def aveEy2CC(self): def aveEy2CC(self):
"Construct the averaging operator on cell edges in the y direction to cell centers." """
if self.dim < 2: return None Construct the averaging operator on cell edges in the y direction to
cell centers.
"""
if self.dim < 2:
return None
if getattr(self, '_aveEy2CC', None) is None: if getattr(self, '_aveEy2CC', None) is None:
# The number of cell centers in each direction # The number of cell centers in each direction
n = self.vnC n = self.vnC
@@ -768,8 +811,12 @@ class DiffOperators(object):
@property @property
def aveEz2CC(self): def aveEz2CC(self):
"Construct the averaging operator on cell edges in the z direction to cell centers." """
if self.dim < 3: return None Construct the averaging operator on cell edges in the z direction to
cell centers.
"""
if self.dim < 3:
return None
if getattr(self, '_aveEz2CC', None) is None: if getattr(self, '_aveEz2CC', None) is None:
# The number of cell centers in each direction # The number of cell centers in each direction
n = self.vnC n = self.vnC
@@ -793,7 +840,10 @@ class DiffOperators(object):
@property @property
def aveN2E(self): def aveN2E(self):
"Construct the averaging operator on cell nodes to cell edges, keeping each dimension separate." """
Construct the averaging operator on cell nodes to cell edges, keeping
each dimension separate.
"""
if getattr(self, '_aveN2E', None) is None: if getattr(self, '_aveN2E', None) is None:
# The number of cell centers in each direction # The number of cell centers in each direction
@@ -802,16 +852,24 @@ class DiffOperators(object):
self._aveN2E = av(n[0]) self._aveN2E = av(n[0])
elif(self.dim == 2): elif(self.dim == 2):
self._aveN2E = sp.vstack((sp.kron(speye(n[1]+1), av(n[0])), self._aveN2E = sp.vstack((sp.kron(speye(n[1]+1), av(n[0])),
sp.kron(av(n[1]), speye(n[0]+1))), format="csr") sp.kron(av(n[1]), speye(n[0]+1))),
format="csr")
elif(self.dim == 3): elif(self.dim == 3):
self._aveN2E = sp.vstack((kron3(speye(n[2]+1), speye(n[1]+1), av(n[0])), self._aveN2E = sp.vstack((kron3(speye(n[2]+1), speye(n[1]+1),
kron3(speye(n[2]+1), av(n[1]), speye(n[0]+1)), av(n[0])),
kron3(av(n[2]), speye(n[1]+1), speye(n[0]+1))), format="csr") kron3(speye(n[2]+1), av(n[1]),
speye(n[0]+1)),
kron3(av(n[2]), speye(n[1]+1),
speye(n[0]+1))),
format="csr")
return self._aveN2E return self._aveN2E
@property @property
def aveN2F(self): def aveN2F(self):
"Construct the averaging operator on cell nodes to cell faces, keeping each dimension separate." """
Construct the averaging operator on cell nodes to cell faces, keeping
each dimension separate.
"""
if getattr(self, '_aveN2F', None) is None: if getattr(self, '_aveN2F', None) is None:
# The number of cell centers in each direction # The number of cell centers in each direction
n = self.vnC n = self.vnC
@@ -819,9 +877,14 @@ class DiffOperators(object):
self._aveN2F = av(n[0]) self._aveN2F = av(n[0])
elif(self.dim == 2): elif(self.dim == 2):
self._aveN2F = sp.vstack((sp.kron(av(n[1]), speye(n[0]+1)), self._aveN2F = sp.vstack((sp.kron(av(n[1]), speye(n[0]+1)),
sp.kron(speye(n[1]+1), av(n[0]))), format="csr") sp.kron(speye(n[1]+1), av(n[0]))),
format="csr")
elif(self.dim == 3): elif(self.dim == 3):
self._aveN2F = sp.vstack((kron3(av(n[2]), av(n[1]), speye(n[0]+1)), self._aveN2F = sp.vstack((kron3(av(n[2]), av(n[1]),
kron3(av(n[2]), speye(n[1]+1), av(n[0])), speye(n[0]+1)),
kron3(speye(n[2]+1), av(n[1]), av(n[0]))), format="csr") kron3(av(n[2]), speye(n[1]+1),
av(n[0])),
kron3(speye(n[2]+1), av(n[1]),
av(n[0]))),
format="csr")
return self._aveN2F return self._aveN2F
+1 -1
View File
@@ -421,7 +421,7 @@ class InnerProducts(object):
def _getEdgePx(M): def _getEdgePx(M):
"""Returns a function for creating projection matrices""" """Returns a function for creating projection matrices"""
def Px(xEdge): def Px(xEdge):
assert xEdge == 'eX0', 'xEdge = %s, not eX0' % xEdge assert xEdge == 'eX0', 'xEdge = {0!s}, not eX0'.format(xEdge)
return sp.identity(M.nC) return sp.identity(M.nC)
return Px return Px
+2 -2
View File
@@ -198,11 +198,11 @@ class TensorMeshIO(object):
""" """
assert mesh.dim == 3 assert mesh.dim == 3
s = '' s = ''
s += '%i %i %i\n' %tuple(mesh.vnC) s += '{0:d} {1:d} {2:d}\n'.format(*tuple(mesh.vnC))
origin = mesh.x0 + np.array([0,0,mesh.hz.sum()]) # Have to it in the same operation or use mesh.x0.copy(), otherwise the mesh.x0 is updated. origin = mesh.x0 + np.array([0,0,mesh.hz.sum()]) # Have to it in the same operation or use mesh.x0.copy(), otherwise the mesh.x0 is updated.
origin.dtype = float origin.dtype = float
s += '%.2f %.2f %.2f\n' %tuple(origin) s += '{0:.2f} {1:.2f} {2:.2f}\n'.format(*tuple(origin))
s += ('%.2f '*mesh.nCx+'\n')%tuple(mesh.hx) s += ('%.2f '*mesh.nCx+'\n')%tuple(mesh.hx)
s += ('%.2f '*mesh.nCy+'\n')%tuple(mesh.hy) s += ('%.2f '*mesh.nCy+'\n')%tuple(mesh.hy)
s += ('%.2f '*mesh.nCz+'\n')%tuple(mesh.hz[::-1]) s += ('%.2f '*mesh.nCz+'\n')%tuple(mesh.hz[::-1])
+4 -4
View File
@@ -23,8 +23,8 @@ class BaseTensorMesh(BaseMesh):
h_i = self._unitDimensions[i] * np.ones(int(h_i))/int(h_i) h_i = self._unitDimensions[i] * np.ones(int(h_i))/int(h_i)
elif type(h_i) is list: elif type(h_i) is list:
h_i = Utils.meshTensor(h_i) h_i = Utils.meshTensor(h_i)
assert isinstance(h_i, np.ndarray), ("h[%i] is not a numpy array." % i) assert isinstance(h_i, np.ndarray), ("h[{0:d}] is not a numpy array.".format(i))
assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i) assert len(h_i.shape) == 1, ("h[{0:d}] must be a 1D numpy array.".format(i))
h[i] = h_i[:] # make a copy. h[i] = h_i[:] # make a copy.
x0 = np.zeros(len(h)) x0 = np.zeros(len(h))
@@ -41,7 +41,7 @@ class BaseTensorMesh(BaseMesh):
elif x_i == 'N': elif x_i == 'N':
x0[i] = -h_i.sum() x0[i] = -h_i.sum()
else: else:
raise Exception("x0[%i] must be a scalar or '0' to be zero, 'C' to center, or 'N' to be negative." % i) raise Exception("x0[{0:d}] must be a scalar or '0' to be zero, 'C' to center, or 'N' to be negative.".format(i))
if isinstance(self, BaseRectangularMesh): if isinstance(self, BaseRectangularMesh):
BaseRectangularMesh.__init__(self, np.array([x.size for x in h]), x0) BaseRectangularMesh.__init__(self, np.array([x.size for x in h]), x0)
@@ -239,7 +239,7 @@ class BaseTensorMesh(BaseMesh):
'CCVz' -> z-component of vector field defined on cell centers 'CCVz' -> z-component of vector field defined on cell centers
""" """
if self._meshType == 'CYL' and self.isSymmetric and locType in ['Ex','Ez','Fy']: if self._meshType == 'CYL' and self.isSymmetric and locType in ['Ex','Ez','Fy']:
raise Exception('Symmetric CylMesh does not support %s interpolation, as this variable does not exist.' % locType) raise Exception('Symmetric CylMesh does not support {0!s} interpolation, as this variable does not exist.'.format(locType))
loc = Utils.asArray_N_x_Dim(loc, self.dim) loc = Utils.asArray_N_x_Dim(loc, self.dim)
+4 -4
View File
@@ -177,7 +177,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
return l return l
def __str__(self): def __str__(self):
outStr = ' ---- %sTreeMesh ---- '%('Oc' if self.dim == 3 else 'Quad') outStr = ' ---- {0!s}TreeMesh ---- '.format(('Oc' if self.dim == 3 else 'Quad'))
def printH(hx, outStr=''): def printH(hx, outStr=''):
i = -1 i = -1
while True: while True:
@@ -213,7 +213,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
outStr += printH(self.hy, outStr='\n hy:') outStr += printH(self.hy, outStr='\n hy:')
outStr += printH(self.hz, outStr='\n hz:') outStr += printH(self.hz, outStr='\n hz:')
outStr += '\n nC: {0:d}'.format(self.nC) outStr += '\n nC: {0:d}'.format(self.nC)
outStr += '\n Fill: %2.2f%%'%(self.fill*100) outStr += '\n Fill: {0:2.2f}%'.format((self.fill*100))
return outStr return outStr
@property @property
@@ -2210,7 +2210,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
ax.set_xlabel('y' if normal == 'X' else 'x') ax.set_xlabel('y' if normal == 'X' else 'x')
ax.set_ylabel('y' if normal == 'Z' else 'z') ax.set_ylabel('y' if normal == 'Z' else 'z')
ax.set_title('Slice %d, %s = %4.2f' % (ind,normal,indLoc)) ax.set_title('Slice {0:d}, {1!s} = {2:4.2f}'.format(ind, normal, indLoc))
if grid: if grid:
_ = antiNormalInd _ = antiNormalInd
@@ -2240,7 +2240,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
if key < 0 : #Handle negative indices if key < 0 : #Handle negative indices
key += len( self ) key += len( self )
if key >= len( self ) : if key >= len( self ) :
raise IndexError, "The index (%d) is out of range."%key raise IndexError, "The index ({0:d}) is out of range.".format(key)
self._numberCells() # no-op if numbered self._numberCells() # no-op if numbered
index = self._i2cc[key] index = self._i2cc[key]
+8 -8
View File
@@ -171,7 +171,7 @@ class TensorView(object):
iz = ix + iy*nX iz = ix + iy*nX
if iz < self.nCz: if iz < self.nCz:
ax.text((ix+1)*(self.vectorNx[-1]-self.x0[0])-pad,(iy)*(self.vectorNy[-1]-self.x0[1])+pad, ax.text((ix+1)*(self.vectorNx[-1]-self.x0[0])-pad,(iy)*(self.vectorNy[-1]-self.x0[1])+pad,
'#%i'%iz,color=annotationColor,verticalalignment='bottom',horizontalalignment='right',size='x-large') '#{0:.0f}'.format(iz),color=annotationColor,verticalalignment='bottom',horizontalalignment='right',size='x-large')
ax.set_title(vType) ax.set_title(vType)
if showIt: plt.show() if showIt: plt.show()
@@ -221,10 +221,10 @@ class TensorView(object):
vTypeOpts = ['CC', 'CCv','N','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 # Some user error checking
assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts) assert vType in vTypeOpts, "vType must be in ['{0!s}']".format("','".join(vTypeOpts))
assert self.dim == 3, 'Must be a 3D mesh. Use plotImage.' assert self.dim == 3, 'Must be a 3D mesh. Use plotImage.'
assert view in viewOpts, "view must be in ['%s']" % "','".join(viewOpts) assert view in viewOpts, "view must be in ['{0!s}']".format("','".join(viewOpts))
assert normal in normalOpts, "normal must be in ['%s']" % "','".join(normalOpts) assert normal in normalOpts, "normal must be in ['{0!s}']".format("','".join(normalOpts))
assert type(grid) is bool, 'grid must be a boolean' assert type(grid) is bool, 'grid must be a boolean'
szSliceDim = getattr(self, 'nC'+normal.lower()) #: Size of the sliced dimension szSliceDim = getattr(self, 'nC'+normal.lower()) #: Size of the sliced dimension
@@ -295,7 +295,7 @@ class TensorView(object):
ax.set_xlabel('y' if normal == 'X' else 'x') ax.set_xlabel('y' if normal == 'X' else 'x')
ax.set_ylabel('y' if normal == 'Z' else 'z') ax.set_ylabel('y' if normal == 'Z' else 'z')
ax.set_title('Slice %d' % ind) ax.set_title('Slice {0:.0f}'.format(ind))
return out return out
@@ -316,11 +316,11 @@ class TensorView(object):
vTypeOptsV = ['CCv','F','E'] vTypeOptsV = ['CCv','F','E']
vTypeOpts = vTypeOptsCC + vTypeOptsV vTypeOpts = vTypeOptsCC + vTypeOptsV
if view == 'vec': if view == 'vec':
assert vType in vTypeOptsV, "vType must be in ['%s'] when view='vec'" % "','".join(vTypeOptsV) assert vType in vTypeOptsV, "vType must be in ['{0!s}'] when view='vec'".format("','".join(vTypeOptsV))
assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts) assert vType in vTypeOpts, "vType must be in ['{0!s}']".format("','".join(vTypeOpts))
viewOpts = ['real','imag','abs','vec'] viewOpts = ['real','imag','abs','vec']
assert view in viewOpts, "view must be in ['%s']" % "','".join(viewOpts) assert view in viewOpts, "view must be in ['{0!s}']".format("','".join(viewOpts))
if ax is None: if ax is None:
+3 -3
View File
@@ -121,7 +121,7 @@ class Minimize(object):
@callback.setter @callback.setter
def callback(self, value): def callback(self, value):
if self.callback is not None: if self.callback is not None:
print 'The callback on the %s Optimization was replaced.' % self.__name__ print 'The callback on the {0!s} Optimization was replaced.'.format(self.__name__)
self._callback = value self._callback = value
@@ -855,7 +855,7 @@ class NewtonRoot(object):
if self.comments and self.doLS: print '\tLinesearch:\n' if self.comments and self.doLS: print '\tLinesearch:\n'
# Enter Linesearch # Enter Linesearch
while True and self.doLS: while True and self.doLS:
if self.comments: print '\t\tResid: %e\n'%norm(rt) if self.comments: print '\t\tResid: {0:e}\n'.format(norm(rt))
if norm(rt) <= norm(r) or norm(rt) < self.tol: if norm(rt) <= norm(r) or norm(rt) < self.tol:
break break
@@ -873,7 +873,7 @@ class NewtonRoot(object):
if norm(rt) < self.tol: if norm(rt) < self.tol:
break break
if self.iter > self.maxIter: if self.iter > self.maxIter:
print 'NewtonRoot stopped by maxIters (%d). norm: %4.4e' % (self.maxIter, norm(rt)) print 'NewtonRoot stopped by maxIters ({0:d}). norm: {1:4.4e}'.format(self.maxIter, norm(rt))
break break
return x return x
+1 -1
View File
@@ -49,7 +49,7 @@ class BaseProblem(object):
def pair(self, d): def pair(self, d):
"""Bind a survey to this problem instance using pointers.""" """Bind a survey to this problem instance using pointers."""
assert isinstance(d, self.surveyPair), "Data object must be an instance of a %s class."%(self.surveyPair.__name__) assert isinstance(d, self.surveyPair), "Data object must be an instance of a {0!s} class.".format((self.surveyPair.__name__))
if d.ispaired: if d.ispaired:
raise Exception("The survey object is already paired to a problem. Use survey.unpair()") raise Exception("The survey object is already paired to a problem. Use survey.unpair()")
self._survey = d self._survey = d
+28 -28
View File
@@ -19,85 +19,85 @@ class Property(object):
return getattr(self, '_propertyLink', None) return getattr(self, '_propertyLink', None)
@propertyLink.setter @propertyLink.setter
def propertyLink(self, value): def propertyLink(self, value):
assert type(value) is tuple and len(value) == 2 and type(value[0]) is str and issubclass(value[1], Maps.IdentityMap), 'Use format: ("%s", Maps.ReciprocalMap)'%self.name assert type(value) is tuple and len(value) == 2 and type(value[0]) is str and issubclass(value[1], Maps.IdentityMap), 'Use format: ("{0!s}", Maps.ReciprocalMap)'.format(self.name)
self._propertyLink = value self._propertyLink = value
def _getMapProperty(self): def _getMapProperty(self):
prop = self prop = self
def fget(self): def fget(self):
return getattr(self, '_%sMap'%prop.name, None) return getattr(self, '_{0!s}Map'.format(prop.name), None)
def fset(self, val): def fset(self, val):
if prop.propertyLink is not None: if prop.propertyLink is not None:
linkName, linkMap = prop.propertyLink linkName, linkMap = prop.propertyLink
assert getattr(self, '%sMap'%linkName, None) is None, 'Cannot set both sides of a linked property.' assert getattr(self, '{0!s}Map'.format(linkName), None) is None, 'Cannot set both sides of a linked property.'
# TODO: Check if the mapping can be correct # TODO: Check if the mapping can be correct
setattr(self, '_%sMap'%prop.name, val) setattr(self, '_{0!s}Map'.format(prop.name), val)
return property(fget=fget, fset=fset, doc=prop.doc) return property(fget=fget, fset=fset, doc=prop.doc)
def _getIndexProperty(self): def _getIndexProperty(self):
prop = self prop = self
def fget(self): def fget(self):
return getattr(self, '_%sIndex'%prop.name, slice(None)) return getattr(self, '_{0!s}Index'.format(prop.name), slice(None))
def fset(self, val): def fset(self, val):
setattr(self, '_%sIndex'%prop.name, val) setattr(self, '_{0!s}Index'.format(prop.name), val)
return property(fget=fget, fset=fset, doc=prop.doc) return property(fget=fget, fset=fset, doc=prop.doc)
def _getProperty(self): def _getProperty(self):
prop = self prop = self
def fget(self): def fget(self):
mapping = getattr(self, '%sMap'%prop.name) mapping = getattr(self, '{0!s}Map'.format(prop.name))
if mapping is None and prop.propertyLink is None: if mapping is None and prop.propertyLink is None:
return prop.defaultVal return prop.defaultVal
if mapping is None and prop.propertyLink is not None: if mapping is None and prop.propertyLink is not None:
linkName, linkMapClass = prop.propertyLink linkName, linkMapClass = prop.propertyLink
linkMap = linkMapClass(None) linkMap = linkMapClass(None)
if getattr(self, '%sMap'%linkName, None) is None: if getattr(self, '{0!s}Map'.format(linkName), None) is None:
return prop.defaultVal return prop.defaultVal
m = getattr(self, '%s'%linkName) m = getattr(self, '{0!s}'.format(linkName))
return linkMap * m return linkMap * m
m = getattr(self, '%sModel'%prop.name) m = getattr(self, '{0!s}Model'.format(prop.name))
return mapping * m return mapping * m
return property(fget=fget) return property(fget=fget)
def _getModelDerivProperty(self): def _getModelDerivProperty(self):
prop = self prop = self
def fget(self): def fget(self):
mapping = getattr(self, '%sMap'%prop.name) mapping = getattr(self, '{0!s}Map'.format(prop.name))
if mapping is None and prop.propertyLink is None: if mapping is None and prop.propertyLink is None:
return None return None
if mapping is None and prop.propertyLink is not None: if mapping is None and prop.propertyLink is not None:
linkName, linkMapClass = prop.propertyLink linkName, linkMapClass = prop.propertyLink
linkedMap = getattr(self, '%sMap'%linkName) linkedMap = getattr(self, '{0!s}Map'.format(linkName))
if linkedMap is None: if linkedMap is None:
return None return None
linkMap = linkMapClass(None) * linkedMap linkMap = linkMapClass(None) * linkedMap
m = getattr(self, '%sModel'%linkName) m = getattr(self, '{0!s}Model'.format(linkName))
return linkMap.deriv( m ) return linkMap.deriv( m )
m = getattr(self, '%sModel'%prop.name) m = getattr(self, '{0!s}Model'.format(prop.name))
return mapping.deriv( m ) return mapping.deriv( m )
return property(fget=fget) return property(fget=fget)
def _getModelProperty(self): def _getModelProperty(self):
prop = self prop = self
def fget(self): def fget(self):
mapping = getattr(self, '%sMap'%prop.name) mapping = getattr(self, '{0!s}Map'.format(prop.name))
if mapping is None: if mapping is None:
return None return None
index = getattr(self.propMap, '%sIndex'%prop.name) index = getattr(self.propMap, '{0!s}Index'.format(prop.name))
return self.vector[index] return self.vector[index]
return property(fget=fget) return property(fget=fget)
def _getModelProjProperty(self): def _getModelProjProperty(self):
prop = self prop = self
def fget(self): def fget(self):
mapping = getattr(self, '%sMap'%prop.name) mapping = getattr(self, '{0!s}Map'.format(prop.name))
if mapping is None: if mapping is None:
return None return None
inds = getattr(self.propMap, '%sIndex'%prop.name) inds = getattr(self.propMap, '{0!s}Index'.format(prop.name))
if type(inds) is slice: if type(inds) is slice:
inds = range(*inds.indices(self.nP)) inds = range(*inds.indices(self.nP))
nI, nP = len(inds),self.nP nI, nP = len(inds),self.nP
@@ -107,7 +107,7 @@ class Property(object):
def _getModelMapProperty(self): def _getModelMapProperty(self):
prop = self prop = self
def fget(self): def fget(self):
return getattr(self.propMap, '_%sMap'%prop.name, None) return getattr(self.propMap, '_{0!s}Map'.format(prop.name), None)
return property(fget=fget) return property(fget=fget)
@@ -123,7 +123,7 @@ class PropModel(object):
inds = [] inds = []
if getattr(self, '_nP', None) is None: if getattr(self, '_nP', None) is None:
for name in self.propMap._properties: for name in self.propMap._properties:
index = getattr(self.propMap, '%sIndex'%name, None) index = getattr(self.propMap, '{0!s}Index'.format(name), None)
if index is not None: if index is not None:
if type(index) is slice: if type(index) is slice:
inds += range(*index.indices(len(self.vector))) inds += range(*index.indices(len(self.vector)))
@@ -163,9 +163,9 @@ class _PropMapMetaClass(type):
if prop.defaultInvProp: if prop.defaultInvProp:
defaultInvProps += [p] defaultInvProps += [p]
if prop.propertyLink is not None: if prop.propertyLink is not None:
assert prop.propertyLink[0] in _properties, "You can only link to things that exist: '%s' is trying to link to '%s'"%(prop.name, prop.propertyLink[0]) assert prop.propertyLink[0] in _properties, "You can only link to things that exist: '{0!s}' is trying to link to '{1!s}'".format(prop.name, prop.propertyLink[0])
if len(defaultInvProps) > 1: if len(defaultInvProps) > 1:
raise Exception('You have more than one default inversion property: %s' % defaultInvProps) raise Exception('You have more than one default inversion property: {0!s}'.format(defaultInvProps))
newClass = super(_PropMapMetaClass, cls).__new__(cls, name, bases, attrs) newClass = super(_PropMapMetaClass, cls).__new__(cls, name, bases, attrs)
@@ -223,7 +223,7 @@ class PropMap(object):
type(m[0]) is str and type(m[0]) is str and
m[0] in self._properties and m[0] in self._properties and
isinstance(m[1], Maps.IdentityMap) isinstance(m[1], Maps.IdentityMap)
for m in maps]), "Use signature: [%s]" % (', '.join(["('%s', %sMap)"%(p,p) for p in self._properties])) for m in maps]), "Use signature: [{0!s}]".format((', '.join(["('{0!s}', {1!s}Map)".format(p, p) for p in self._properties])))
if slices is None: if slices is None:
slices = dict() slices = dict()
else: else:
@@ -236,8 +236,8 @@ class PropMap(object):
nP = 0 nP = 0
for name, mapping in maps: for name, mapping in maps:
setattr(self, '%sMap'%name, mapping) setattr(self, '{0!s}Map'.format(name), mapping)
setattr(self, '%sIndex'%name, slices.get(name, slice(nP, nP + mapping.nP))) setattr(self, '{0!s}Index'.format(name), slices.get(name, slice(nP, nP + mapping.nP)))
nP += mapping.nP nP += mapping.nP
self.nP = nP self.nP = nP
@@ -250,12 +250,12 @@ class PropMap(object):
def clearMaps(self): def clearMaps(self):
for name in self._properties: for name in self._properties:
setattr(self, '%sMap'%name, None) setattr(self, '{0!s}Map'.format(name), None)
setattr(self, '%sIndex'%name, None) setattr(self, '{0!s}Index'.format(name), None)
def __call__(self, vec): def __call__(self, vec):
return self.PropModel(self, vec) return self.PropModel(self, vec)
def __contains__(self, val): def __contains__(self, val):
activeMaps = [name for name in self._properties if getattr(self, '%sMap'%name) is not None] activeMaps = [name for name in self._properties if getattr(self, '{0!s}Map'.format(name)) is not None]
return val in activeMaps return val in activeMaps
+6 -6
View File
@@ -26,7 +26,7 @@ class BaseRx(object):
def rxType(self, value): def rxType(self, value):
known = self.knownRxTypes known = self.knownRxTypes
if known is not None: if known is not None:
assert value in known, "rxType must be in ['%s']" % ("', '".join(known)) assert value in known, "rxType must be in ['{0!s}']".format(("', '".join(known)))
self._rxType = value self._rxType = value
@property @property
@@ -125,7 +125,7 @@ class BaseSrc(object):
def __init__(self, rxList, **kwargs): def __init__(self, rxList, **kwargs):
assert type(rxList) is list, 'rxList must be a list' assert type(rxList) is list, 'rxList must be a list'
for rx in rxList: for rx in rxList:
assert isinstance(rx, self.rxPair), 'rxList must be a %s'%self.rxPair.__name__ assert isinstance(rx, self.rxPair), 'rxList must be a {0!s}'.format(self.rxPair.__name__)
assert len(set(rxList)) == len(rxList), 'The rxList must be unique' assert len(set(rxList)) == len(rxList), 'The rxList must be unique'
self.uid = str(uuid.uuid4()) self.uid = str(uuid.uuid4())
self.rxList = rxList self.rxList = rxList
@@ -227,7 +227,7 @@ class BaseSurvey(object):
@srcList.setter @srcList.setter
def srcList(self, value): def srcList(self, value):
assert type(value) is list, 'srcList must be a list' assert type(value) is list, 'srcList must be a list'
assert np.all([isinstance(src, self.srcPair) for src in value]), 'All sources must be instances of %s' % self.srcPair.__name__ assert np.all([isinstance(src, self.srcPair) for src in value]), 'All sources must be instances of {0!s}'.format(self.srcPair.__name__)
assert len(set(value)) == len(value), 'The srcList must be unique' assert len(set(value)) == len(value), 'The srcList must be unique'
self._srcList = value self._srcList = value
self._sourceOrder = dict() self._sourceOrder = dict()
@@ -238,10 +238,10 @@ class BaseSurvey(object):
sources = [sources] sources = [sources]
for src in sources: for src in sources:
if getattr(src,'uid',None) is None: if getattr(src,'uid',None) is None:
raise KeyError('Source does not have a uid: %s'%str(src)) raise KeyError('Source does not have a uid: {0!s}'.format(str(src)))
inds = map(lambda src: self._sourceOrder.get(src.uid, None), sources) inds = map(lambda src: self._sourceOrder.get(src.uid, None), sources)
if None in inds: if None in inds:
raise KeyError('Some of the sources specified are not in this survey. %s'%str(inds)) raise KeyError('Some of the sources specified are not in this survey. {0!s}'.format(str(inds)))
return inds return inds
@property @property
@@ -263,7 +263,7 @@ class BaseSurvey(object):
def pair(self, p): def pair(self, p):
"""Bind a problem to this survey instance using pointers""" """Bind a problem to this survey instance using pointers"""
assert hasattr(p, 'surveyPair'), "Problem must have an attribute 'surveyPair'." assert hasattr(p, 'surveyPair'), "Problem must have an attribute 'surveyPair'."
assert isinstance(self, p.surveyPair), "Problem requires survey object must be an instance of a %s class."%(p.surveyPair.__name__) assert isinstance(self, p.surveyPair), "Problem requires survey object must be an instance of a {0!s} class.".format((p.surveyPair.__name__))
if p.ispaired: if p.ispaired:
raise Exception("The problem object is already paired to a survey. Use prob.unpair()") raise Exception("The problem object is already paired to a survey. Use prob.unpair()")
self._prob = p self._prob = p
+8 -9
View File
@@ -4,7 +4,6 @@ from SimPEG.Utils import mkvc, sdiag, diagEst
from SimPEG import Utils from SimPEG import Utils
from SimPEG.Mesh import TensorMesh, CurvilinearMesh, CylMesh from SimPEG.Mesh import TensorMesh, CurvilinearMesh, CylMesh
from SimPEG.Mesh.TreeMesh import TreeMesh as Tree from SimPEG.Mesh.TreeMesh import TreeMesh as Tree
import numpy as np
import scipy.sparse as sp import scipy.sparse as sp
import unittest import unittest
import inspect import inspect
@@ -200,10 +199,10 @@ class OrderTest(unittest.TestCase):
print '_____________________________________________' print '_____________________________________________'
print ' h | error | e(i-1)/e(i) | order' print ' h | error | e(i-1)/e(i) | order'
print '~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~' print '~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~'
print '%4i | %8.2e |' % (nc, err) print '{0:4d} | {1:8.2e} |'.format(nc, err)
else: else:
order.append(np.log(err/err_old)/np.log(max_h/max_h_old)) order.append(np.log(err/err_old)/np.log(max_h/max_h_old))
print '%4i | %8.2e | %6.4f | %6.4f' % (nc, err, err_old/err, order[-1]) print '{0:4d} | {1:8.2e} | {2:6.4f} | {3:6.4f}'.format(nc, err, err_old/err, order[-1])
err_old = err err_old = err
max_h_old = max_h max_h_old = max_h
print '---------------------------------------------' print '---------------------------------------------'
@@ -258,8 +257,8 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
Tests.checkDerivative(simplePass, np.random.randn(5)) Tests.checkDerivative(simplePass, np.random.randn(5))
""" """
print "%s checkDerivative %s" % ('='*20, '='*20) print "{0!s} checkDerivative {1!s}".format('='*20, '='*20)
print "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n%s" % ('-'*57) print "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n{0!s}".format(('-'*57))
f0, J0 = fctn(x0) f0, J0 = fctn(x0)
@@ -290,7 +289,7 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
order0 = np.log10(E0[:-1]/E0[1:]) order0 = np.log10(E0[:-1]/E0[1:])
order1 = np.log10(E1[:-1]/E1[1:]) order1 = np.log10(E1[:-1]/E1[1:])
print " %d %1.2e %1.3e %1.3e %1.3f" % (i, h[i], E0[i], E1[i], np.nan if i == 0 else order1[i-1]) print " {0:d} {1:1.2e} {2:1.3e} {3:1.3e} {4:1.3f}".format(i, h[i], E0[i], E1[i], np.nan if i == 0 else order1[i-1])
# Ensure we are about precision # Ensure we are about precision
order0 = order0[E0[1:] > eps] order0 = order0[E0[1:] > eps]
@@ -302,10 +301,10 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
passTest = belowTol or correctOrder passTest = belowTol or correctOrder
if passTest: if passTest:
print "%s PASS! %s" % ('='*25, '='*25) print "{0!s} PASS! {1!s}".format('='*25, '='*25)
print happiness[np.random.randint(len(happiness))]+'\n' print happiness[np.random.randint(len(happiness))]+'\n'
else: else:
print "%s\n%s FAIL! %s\n%s" % ('*'*57, '<'*25, '>'*25, '*'*57) print "{0!s}\n{1!s} FAIL! {2!s}\n{3!s}".format('*'*57, '<'*25, '>'*25, '*'*57)
print sadness[np.random.randint(len(sadness))]+'\n' print sadness[np.random.randint(len(sadness))]+'\n'
@@ -314,7 +313,7 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
ax = ax or plt.subplot(111) ax = ax or plt.subplot(111)
ax.loglog(h, E0, 'b') ax.loglog(h, E0, 'b')
ax.loglog(h, E1, 'g--') ax.loglog(h, E1, 'g--')
ax.set_title('Check Derivative - %s' % ('PASSED :)' if passTest else 'FAILED :(')) ax.set_title('Check Derivative - {0!s}'.format(('PASSED :)' if passTest else 'FAILED :(')))
ax.set_xlabel('h') ax.set_xlabel('h')
ax.set_ylabel('Error') ax.set_ylabel('Error')
leg = ax.legend(['$\mathcal{O}(h)$', '$\mathcal{O}(h^2)$'], loc='best', leg = ax.legend(['$\mathcal{O}(h)$', '$\mathcal{O}(h^2)$'], loc='best',
+1 -1
View File
@@ -8,7 +8,7 @@ def _checkAccuracy(A, b, X, accuracyTol):
if nrm_b > 0: if nrm_b > 0:
nrm /= nrm_b nrm /= nrm_b
if nrm > accuracyTol: if nrm > accuracyTol:
msg = '### SolverWarning ###: Accuracy on solve is above tolerance: %e > %e' % (nrm, accuracyTol) msg = '### SolverWarning ###: Accuracy on solve is above tolerance: {0:e} > {1:e}'.format(nrm, accuracyTol)
print msg print msg
warnings.warn(msg, RuntimeWarning) warnings.warn(msg, RuntimeWarning)
+15 -15
View File
@@ -32,7 +32,7 @@ def memProfileWrapper(towrap, *funNames):
if hasattr(towrap,f): if hasattr(towrap,f):
attrs[f] = profile(getattr(towrap,f)) attrs[f] = profile(getattr(towrap,f))
else: else:
print '%s not found in %s Class' % (f, towrap.__name__) print '{0!s} not found in {1!s} Class'.format(f, towrap.__name__)
return type(towrap.__name__ + 'MemProfileWrap', (towrap,), attrs) return type(towrap.__name__ + 'MemProfileWrap', (towrap,), attrs)
@@ -65,7 +65,7 @@ def setKwargs(obj, ignore=None, **kwargs):
if hasattr(obj, attr): if hasattr(obj, attr):
setattr(obj, attr, kwargs[attr]) setattr(obj, attr, kwargs[attr])
else: else:
raise Exception('%s attr is not recognized' % attr) raise Exception('{0!s} attr is not recognized'.format(attr))
hook(obj,hook, silent=True) hook(obj,hook, silent=True)
hook(obj,setKwargs, silent=True) hook(obj,setKwargs, silent=True)
@@ -74,7 +74,7 @@ def printTitles(obj, printers, name='Print Titles', pad=''):
titles = '' titles = ''
widths = 0 widths = 0
for printer in printers: for printer in printers:
titles += ('{:^%i}'%printer['width']).format(printer['title']) + '' titles += ('{{:^{0:d}}}'.format(printer['width'])).format(printer['title']) + ''
widths += printer['width'] widths += printer['width']
print pad + "{0} {1} {0}".format('='*((widths-1-len(name))/2), name) print pad + "{0} {1} {0}".format('='*((widths-1-len(name))/2), name)
print pad + titles print pad + titles
@@ -83,7 +83,7 @@ def printTitles(obj, printers, name='Print Titles', pad=''):
def printLine(obj, printers, pad=''): def printLine(obj, printers, pad=''):
values = '' values = ''
for printer in printers: for printer in printers:
values += ('{:^%i}'%printer['width']).format(printer['format'] % printer['value'](obj)) values += ('{{:^{0:d}}}'.format(printer['width'])).format(printer['format'] % printer['value'](obj))
print pad + values print pad + values
def checkStoppers(obj, stoppers): def checkStoppers(obj, stoppers):
@@ -104,12 +104,12 @@ def checkStoppers(obj, stoppers):
return (len(optimal)>0 and all(optimal)) | (len(critical)>0 and any(critical)) return (len(optimal)>0 and all(optimal)) | (len(critical)>0 and any(critical))
def printStoppers(obj, stoppers, pad='', stop='STOP!', done='DONE!'): def printStoppers(obj, stoppers, pad='', stop='STOP!', done='DONE!'):
print pad + "%s%s%s" % ('-'*25,stop,'-'*25) print pad + "{0!s}{1!s}{2!s}".format('-'*25, stop, '-'*25)
for stopper in stoppers: for stopper in stoppers:
l = stopper['left'](obj) l = stopper['left'](obj)
r = stopper['right'](obj) r = stopper['right'](obj)
print pad + stopper['str'] % (l<=r,l,r) print pad + stopper['str'] % (l<=r,l,r)
print pad + "%s%s%s" % ('-'*25,done,'-'*25) print pad + "{0!s}{1!s}{2!s}".format('-'*25, done, '-'*25)
def callHooks(match, mainFirst=False): def callHooks(match, mainFirst=False):
""" """
@@ -144,14 +144,14 @@ def callHooks(match, mainFirst=False):
extra = """ extra = """
If you have things that also need to run in the method %s, you can create a method:: If you have things that also need to run in the method {0!s}, you can create a method::
def _%s*(self, ... ): def _{1!s}*(self, ... ):
pass pass
Where the * can be any string. If present, _%s* will be called at the start of the default %s call. Where the * can be any string. If present, _{2!s}* will be called at the start of the default {3!s} call.
You may also completely overwrite this function. You may also completely overwrite this function.
""" % (match, match, match, match) """.format(match, match, match, match)
doc = wrapper.__doc__ doc = wrapper.__doc__
wrapper.__doc__ = ('' if doc is None else doc) + extra wrapper.__doc__ = ('' if doc is None else doc) + extra
return wrapper return wrapper
@@ -186,7 +186,7 @@ def asArray_N_x_Dim(pts, dim):
elif len(pts.shape) == 1: elif len(pts.shape) == 1:
pts = pts[:,np.newaxis] pts = pts[:,np.newaxis]
assert pts.shape[1] == dim, "pts must be a column vector of shape (nPts, %d) not (%d, %d)" % ((dim,)+pts.shape) assert pts.shape[1] == dim, "pts must be a column vector of shape (nPts, {0:d}) not ({1:d}, {2:d})".format(*((dim,)+pts.shape))
return pts return pts
@@ -207,17 +207,17 @@ def requires(var):
.. note:: .. note::
To use survey.%s(), SimPEG requires that a problem be bound to the survey. To use survey.{0!s}(), SimPEG requires that a problem be bound to the survey.
If a problem has not been bound, an Exception will be raised. If a problem has not been bound, an Exception will be raised.
To bind a problem to the Data object:: To bind a problem to the Data object::
survey.pair(myProblem) survey.pair(myProblem)
""" % f.__name__ """.format(f.__name__)
else: else:
extra = """ extra = """
To use *%s* method, SimPEG requires that the %s be specified. To use *{0!s}* method, SimPEG requires that the {1!s} be specified.
""" % (f.__name__, var) """.format(f.__name__, var)
@wraps(f) @wraps(f)
def requiresVarWrapper(self,*args,**kwargs): def requiresVarWrapper(self,*args,**kwargs):
if getattr(self, var, None) is None: if getattr(self, var, None) is None:
+1 -1
View File
@@ -80,7 +80,7 @@ def indexCube(nodes, gridSize, n=None):
# Make sure that we choose from the possible nodes. # Make sure that we choose from the possible nodes.
possibleNodes = 'ABCD' if gridSize.size == 2 else 'ABCDEFGH' possibleNodes = 'ABCD' if gridSize.size == 2 else 'ABCDEFGH'
for node in nodes: for node in nodes:
assert node in possibleNodes, "Nodes must be chosen from: '%s'" % possibleNodes assert node in possibleNodes, "Nodes must be chosen from: '{0!s}'".format(possibleNodes)
dim = gridSize.size dim = gridSize.size
if n is None: if n is None:
n = gridSize - 1 n = gridSize - 1
+1 -1
View File
@@ -278,7 +278,7 @@ class TensorType(object):
else: else:
raise Exception('Unexpected shape of tensor') raise Exception('Unexpected shape of tensor')
def __str__(self): def __str__(self):
return 'TensorType[%i]: %s' % (self._tt, self._tts) return 'TensorType[{0:d}]: {1!s}'.format(self._tt, self._tts)
def __eq__(self, v): return self._tt == v def __eq__(self, v): return self._tt == v
def __le__(self, v): return self._tt <= v def __le__(self, v): return self._tt <= v
def __ge__(self, v): return self._tt >= v def __ge__(self, v): return self._tt >= v
+2 -2
View File
@@ -26,7 +26,7 @@ def surface2ind_topo(mesh, topo, gridLoc='CC'):
gridTopo = Ftopo(XY).reshape(mesh.vnN[:2], order='F') gridTopo = Ftopo(XY).reshape(mesh.vnN[:2], order='F')
if mesh._meshType not in ['TENSOR', 'CYL', 'BASETENSOR']: if mesh._meshType not in ['TENSOR', 'CYL', 'BASETENSOR']:
raise NotImplementedError('Nodal surface2ind_topo not implemented for %s mesh'%mesh._meshType) raise NotImplementedError('Nodal surface2ind_topo not implemented for {0!s} mesh'.format(mesh._meshType))
Nz = mesh.vectorNz[1:] # TODO: this will only work for tensor meshes Nz = mesh.vectorNz[1:] # TODO: this will only work for tensor meshes
actind = np.array([False]*mesh.nC).reshape(mesh.vnC, order='F') actind = np.array([False]*mesh.nC).reshape(mesh.vnC, order='F')
@@ -47,7 +47,7 @@ def surface2ind_topo(mesh, topo, gridLoc='CC'):
gridTopo = Ftopo(mesh.vectorNx) gridTopo = Ftopo(mesh.vectorNx)
if mesh._meshType not in ['TENSOR', 'CYL', 'BASETENSOR']: if mesh._meshType not in ['TENSOR', 'CYL', 'BASETENSOR']:
raise NotImplementedError('Nodal surface2ind_topo not implemented for %s mesh'%mesh._meshType) raise NotImplementedError('Nodal surface2ind_topo not implemented for {0!s} mesh'.format(mesh._meshType))
Ny = mesh.vectorNy[1:] # TODO: this will only work for tensor meshes Ny = mesh.vectorNy[1:] # TODO: this will only work for tensor meshes
actind = np.array([False]*mesh.nC).reshape(mesh.vnC, order='F') actind = np.array([False]*mesh.nC).reshape(mesh.vnC, order='F')
+1 -1
View File
@@ -266,7 +266,7 @@ def _supress_nonlocal_image_warn(self, msg, node):
from docutils.utils import get_source_line from docutils.utils import get_source_line
if not msg.startswith('nonlocal image URI found:'): if not msg.startswith('nonlocal image URI found:'):
self._warnfunc(msg, '%s:%s' % get_source_line(node)) self._warnfunc(msg, '{0!s}:{1!s}'.format(*get_source_line(node)))
supress_nonlocal_image_warn() supress_nonlocal_image_warn()
+1 -1
View File
@@ -66,7 +66,7 @@ Numpy and Matlab
Lessons in Python Lessons in Python
----------------- -----------------
* `Software Carpentry <http://software-carpentry.org/v4/python/index.html>`_ * `Software Carpentry <http://swcarpentry.github.io/python-novice-inflammation/>`_
* `Introduction to NumPy and Matplotlib <http://www.youtube.com/watch?v=3Fp1zn5ao2M>`_ * `Introduction to NumPy and Matplotlib <http://www.youtube.com/watch?v=3Fp1zn5ao2M>`_
Editing Python Editing Python
+10 -1
View File
@@ -47,7 +47,16 @@ direct current (DC) resistivity and induced polarization (IP) geophysical proble
DC resistivity survey DC resistivity survey
===================== =====================
Electrical resistivity of subsurface materials is measured by causing an electrical current to flow in the earth between one pair of electrodes while the voltage across a second pair of electrodes is measured. The result is an "apparent" resistivity which is a value representing the weighted average resistivity over a volume of the earth. Variations in this measurement are caused by variations in the soil, rock, and pore fluid electrical resistivity. Surveys require contact with the ground, so they can be labour intensive. Results are sometimes interpreted directly, but more commonly, 1D, 2D or 3D models are estimated using inversion procedures (`GPG <http://www.eos.ubc.ca/courses/eosc350/content/>`_). Electrical resistivity of subsurface materials is measured by causing an
electrical current to flow in the earth between one pair of electrodes while
the voltage across a second pair of electrodes is measured. The result is an
"apparent" resistivity which is a value representing the weighted average
resistivity over a volume of the earth. Variations in this measurement are
caused by variations in the soil, rock, and pore fluid electrical resistivity.
Surveys require contact with the ground, so they can be labour intensive.
Results are sometimes interpreted directly, but more commonly, 1D, 2D or 3D
models are estimated using inversion procedures (`GPG
<http://gpg.geosci.xyz>`_).
Background Background
+14 -6
View File
@@ -16,16 +16,24 @@ SimPEG Documentation
:target: https://github.com/simpeg/simpeg/blob/master/LICENSE :target: https://github.com/simpeg/simpeg/blob/master/LICENSE
:alt: BSD 3 clause license. :alt: BSD 3 clause license.
.. image:: https://img.shields.io/travis/simpeg/simpeg.svg .. image:: https://api.travis-ci.org/simpeg/simpeg.svg?branch=master
:target: https://travis-ci.org/simpeg/simpeg?branch=master :target: https://travis-ci.org/simpeg/simpeg
:alt: Travis CI build status :alt: Travis CI build status
.. image:: https://img.shields.io/coveralls/simpeg/simpeg.svg .. image:: http://img.shields.io/badge/GITTER-JOIN_CHAT-brightgreen.svg?style=flat-square
:target: https://coveralls.io/r/simpeg/simpeg?branch=master :alt: gitter chat room at https://gitter.im/simpeg/simpeg
:alt: Coverage status :target: https://gitter.im/simpeg/simpeg
.. image:: https://codecov.io/gh/simpeg/simpeg/branch/master/graph/badge.svg .. image:: https://codecov.io/gh/simpeg/simpeg/branch/master/graph/badge.svg
:target: https://codecov.io/gh/simpeg/simpeg    :target: https://codecov.io/gh/simpeg/simpeg
.. image:: https://www.quantifiedcode.com/api/v1/project/933aa3decf444538aa432c8817169b6d/badge.svg
:target: https://www.quantifiedcode.com/app/project/933aa3decf444538aa432c8817169b6d
:alt: Code issues
.. image:: https://api.codacy.com/project/badge/Grade/4fc959a5294a418fa21fc7bc3b3aa078
:target: https://www.codacy.com/app/lindseyheagy/simpeg?utm_source=github.com&amp;utm_medium=referral&amp;utm_content=simpeg/simpeg&amp;utm_campaign=Badge_Grade
:alt: codacy
Simulation and Parameter Estimation in Geophysics - A python package for simulation and gradient based parameter estimation in the context of geophysical applications. Simulation and Parameter Estimation in Geophysics - A python package for simulation and gradient based parameter estimation in the context of geophysical applications.
+1 -1
View File
@@ -54,7 +54,7 @@ class Images(webapp2.RequestHandler):
class Redirect(webapp2.RequestHandler): class Redirect(webapp2.RequestHandler):
def get(self): def get(self):
path = str(self.request.path).split(os.path.sep)[3:] path = str(self.request.path).split(os.path.sep)[3:]
self.redirect(('/%s'%os.path.sep.join(path)), permanent=True) self.redirect(('/{0!s}'.format(os.path.sep.join(path))), permanent=True)
class MainPage(webapp2.RequestHandler): class MainPage(webapp2.RequestHandler):
+69 -30
View File
@@ -1,14 +1,16 @@
import unittest import unittest
from SimPEG import * from SimPEG import Mesh, Problem, Fields, Survey, Utils
import numpy as np
class FieldsTest(unittest.TestCase): class FieldsTest(unittest.TestCase):
def setUp(self): def setUp(self):
mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30]) mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10, 11, 12]],
[0, 0, -30])
x = np.linspace(5, 10, 3) x = np.linspace(5, 10, 3)
XYZ = Utils.ndgrid(x, x, np.r_[0.]) XYZ = Utils.ndgrid(x, x, np.r_[0.])
srcLoc = np.r_[0,0,0.] srcLoc = np.r_[0., 0., 0.]
rxList0 = Survey.BaseRx(XYZ, 'exi') rxList0 = Survey.BaseRx(XYZ, 'exi')
Src0 = Survey.BaseSrc([rxList0], loc=srcLoc) Src0 = Survey.BaseSrc([rxList0], loc=srcLoc)
rxList1 = Survey.BaseRx(XYZ, 'bxi') rxList1 = Survey.BaseRx(XYZ, 'bxi')
@@ -21,7 +23,10 @@ class FieldsTest(unittest.TestCase):
srcList = [Src0, Src1, Src2, Src3, Src4] srcList = [Src0, Src1, Src2, Src3, Src4]
survey = Survey.BaseSurvey(srcList=srcList) survey = Survey.BaseSurvey(srcList=srcList)
self.D = Survey.Data(survey) self.D = Survey.Data(survey)
self.F = Problem.Fields(mesh, survey, knownFields={'phi':'CC','e':'E','b':'F'}, dtype={"phi":float,"e":complex,"b":complex}) self.F = Problem.Fields(mesh, survey, knownFields={'phi': 'CC',
'e': 'E', 'b': 'F'},
dtype={"phi": float, "e": complex,
"b": complex})
self.Src0 = Src0 self.Src0 = Src0
self.Src1 = Src1 self.Src1 = Src1
self.mesh = mesh self.mesh = mesh
@@ -38,16 +43,18 @@ class FieldsTest(unittest.TestCase):
self.assertTrue('e' in F) self.assertTrue('e' in F)
def test_overlappingFields(self): def test_overlappingFields(self):
self.assertRaises(AssertionError, Problem.Fields, self.F.mesh, self.F.survey, self.assertRaises(AssertionError, Problem.Fields, self.F.mesh,
knownFields={'b':'F'}, self.F.survey, knownFields={'b': 'F'},
aliasFields={'b': ['b', (lambda F, b, ind: b)]}) aliasFields={'b': ['b', (lambda F, b, ind: b)]})
def test_SetGet(self): def test_SetGet(self):
F = self.F F = self.F
nSrc = F.survey.nSrc nSrc = F.survey.nSrc
e = np.random.rand(F.mesh.nE, nSrc) + np.random.rand(F.mesh.nE, nSrc)*1j e = (np.random.rand(F.mesh.nE, nSrc) +
np.random.rand(F.mesh.nE, nSrc)*1j)
F[:, 'e'] = e F[:, 'e'] = e
b = np.random.rand(F.mesh.nF, nSrc) + np.random.rand(F.mesh.nF, nSrc)*1j b = (np.random.rand(F.mesh.nF, nSrc) +
np.random.rand(F.mesh.nF, nSrc)*1j)
F[:, 'b'] = b F[:, 'b'] = b
self.assertTrue(np.all(F[:, 'e'] == e)) self.assertTrue(np.all(F[:, 'e'] == e))
@@ -86,12 +93,16 @@ class FieldsTest(unittest.TestCase):
def test_assertions(self): def test_assertions(self):
freq = [self.Src0, self.Src1] freq = [self.Src0, self.Src1]
bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nSrc) bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nSrc)
def fun(): self.F[freq, 'b'] = bWrongSize def fun(): self.F[freq, 'b'] = bWrongSize
self.assertRaises(ValueError, fun) self.assertRaises(ValueError, fun)
def fun(): self.F[-999.] def fun(): self.F[-999.]
self.assertRaises(KeyError, fun) self.assertRaises(KeyError, fun)
def fun(): self.F['notRight'] def fun(): self.F['notRight']
self.assertRaises(KeyError, fun) self.assertRaises(KeyError, fun)
def fun(): self.F[freq, 'notThere'] def fun(): self.F[freq, 'notThere']
self.assertRaises(KeyError, fun) self.assertRaises(KeyError, fun)
@@ -99,7 +110,8 @@ class FieldsTest(unittest.TestCase):
class FieldsTest_Alias(unittest.TestCase): class FieldsTest_Alias(unittest.TestCase):
def setUp(self): def setUp(self):
mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30]) mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10, 11, 12]],
[0, 0, -30])
x = np.linspace(5, 10, 3) x = np.linspace(5, 10, 3)
XYZ = Utils.ndgrid(x, x, np.r_[0.]) XYZ = Utils.ndgrid(x, x, np.r_[0.])
srcLoc = np.r_[0, 0, 0.] srcLoc = np.r_[0, 0, 0.]
@@ -114,7 +126,11 @@ class FieldsTest_Alias(unittest.TestCase):
Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3], loc=srcLoc) Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3], loc=srcLoc)
srcList = [Src0, Src1, Src2, Src3, Src4] srcList = [Src0, Src1, Src2, Src3, Src4]
survey = Survey.BaseSurvey(srcList=srcList) survey = Survey.BaseSurvey(srcList=srcList)
self.F = Problem.Fields(mesh, survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',(lambda e, ind: self.F.mesh.edgeCurl * e)]}) self.F = Problem.Fields(mesh, survey, knownFields={'e': 'E'},
aliasFields={'b': ['e', 'F',
(lambda e, ind:
self.F.mesh.edgeCurl *
e)]})
self.Src0 = Src0 self.Src0 = Src0
self.Src1 = Src1 self.Src1 = Src1
self.mesh = mesh self.mesh = mesh
@@ -149,18 +165,20 @@ class FieldsTest_Alias(unittest.TestCase):
def alias(e, ind): def alias(e, ind):
self.assertTrue(ind[0] is self.Src0) self.assertTrue(ind[0] is self.Src0)
return self.F.mesh.edgeCurl * e return self.F.mesh.edgeCurl * e
F = Problem.Fields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]}) F = Problem.Fields(self.F.mesh, self.F.survey, knownFields={'e': 'E'},
aliasFields={'b': ['e', 'F', alias]})
e = np.random.rand(F.mesh.nE, 1) e = np.random.rand(F.mesh.nE, 1)
F[self.Src0, 'e'] = e F[self.Src0, 'e'] = e
F[self.Src0, 'b'] F[self.Src0, 'b']
def alias(e, ind): def alias(e, ind):
self.assertTrue(type(ind) is list) self.assertTrue(type(ind) is list)
self.assertTrue(ind[0] is self.Src0) self.assertTrue(ind[0] is self.Src0)
self.assertTrue(ind[1] is self.Src1) self.assertTrue(ind[1] is self.Src1)
return self.F.mesh.edgeCurl * e return self.F.mesh.edgeCurl * e
F = Problem.Fields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]})
F = Problem.Fields(self.F.mesh, self.F.survey, knownFields={'e': 'E'},
aliasFields={'b': ['e', 'F', alias]})
e = np.random.rand(F.mesh.nE, 2) e = np.random.rand(F.mesh.nE, 2)
F[[self.Src0, self.Src1], 'e'] = e F[[self.Src0, self.Src1], 'e'] = e
F[[self.Src0, self.Src1], 'b'] F[[self.Src0, self.Src1], 'b']
@@ -169,7 +187,8 @@ class FieldsTest_Alias(unittest.TestCase):
class FieldsTest_Time(unittest.TestCase): class FieldsTest_Time(unittest.TestCase):
def setUp(self): def setUp(self):
mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30]) mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10, 11, 12]],
[0, 0, -30])
x = np.linspace(5, 10, 3) x = np.linspace(5, 10, 3)
XYZ = Utils.ndgrid(x, x, np.r_[0.]) XYZ = Utils.ndgrid(x, x, np.r_[0.])
srcLoc = np.r_[0, 0, 0.] srcLoc = np.r_[0, 0, 0.]
@@ -186,7 +205,9 @@ class FieldsTest_Time(unittest.TestCase):
survey = Survey.BaseSurvey(srcList=srcList) survey = Survey.BaseSurvey(srcList=srcList)
prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10., 3), (20., 2)]) prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10., 3), (20., 2)])
survey.pair(prob) survey.pair(prob)
self.F = Problem.TimeFields(mesh, survey, knownFields={'phi':'CC','e':'E','b':'F'}) self.F = Problem.TimeFields(mesh, survey, knownFields={'phi': 'CC',
'e': 'E',
'b': 'F'})
self.Src0 = Src0 self.Src0 = Src0
self.Src1 = Src1 self.Src1 = Src1
self.mesh = mesh self.mesh = mesh
@@ -230,7 +251,8 @@ class FieldsTest_Time(unittest.TestCase):
b = np.random.rand(F.mesh.nF, 1, nT) b = np.random.rand(F.mesh.nF, 1, nT)
F[self.Src0, 'b', 0] = b[:, :, 0] F[self.Src0, 'b', 0] = b[:, :, 0]
self.assertTrue(np.all(F[self.Src0, 'b', 0] == Utils.mkvc(b[:,0,0],2))) self.assertTrue(np.all(F[self.Src0, 'b', 0] == Utils.mkvc(b[:, 0, 0],
2)))
phi = np.random.rand(F.mesh.nC, 2, nT) phi = np.random.rand(F.mesh.nC, 2, nT)
F[[self.Src0, self.Src1], 'phi'] = phi F[[self.Src0, self.Src1], 'phi'] = phi
@@ -246,11 +268,14 @@ class FieldsTest_Time(unittest.TestCase):
self.assertTrue(F[self.Src0, 'b'].shape == (F.mesh.nF, nT)) self.assertTrue(F[self.Src0, 'b'].shape == (F.mesh.nF, nT))
self.assertTrue(np.all(F[self.Src0, 'b'] == b[:, 0, :])) self.assertTrue(np.all(F[self.Src0, 'b'] == b[:, 0, :]))
self.assertTrue(np.all(F[self.Src1, 'b'] == b[:, 1, :])) self.assertTrue(np.all(F[self.Src1, 'b'] == b[:, 1, :]))
self.assertTrue(np.all(F[self.Src0,'b',1] == Utils.mkvc(b[:,0,1],2))) self.assertTrue(np.all(F[self.Src0, 'b', 1] ==
self.assertTrue(np.all(F[self.Src1,'b',1] == Utils.mkvc(b[:,1,1],2))) Utils.mkvc(b[:, 0, 1], 2)))
self.assertTrue(np.all(F[self.Src0,'b',4] == Utils.mkvc(b[:,0,4],2))) self.assertTrue(np.all(F[self.Src1, 'b', 1] ==
self.assertTrue(np.all(F[self.Src1,'b',4] == Utils.mkvc(b[:,1,4],2))) Utils.mkvc(b[:, 1, 1], 2)))
self.assertTrue(np.all(F[self.Src0, 'b', 4] ==
Utils.mkvc(b[:, 0, 4], 2)))
self.assertTrue(np.all(F[self.Src1, 'b', 4] ==
Utils.mkvc(b[:, 1, 4], 2)))
b = np.random.rand(F.mesh.nF, 2, nT) b = np.random.rand(F.mesh.nF, 2, nT)
F[[self.Src0, self.Src1], 'b', 0] = b[:, :, 0] F[[self.Src0, self.Src1], 'b', 0] = b[:, :, 0]
@@ -258,12 +283,16 @@ class FieldsTest_Time(unittest.TestCase):
def test_assertions(self): def test_assertions(self):
freq = [self.Src0, self.Src1] freq = [self.Src0, self.Src1]
bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nSrc) bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nSrc)
def fun(): self.F[freq, 'b'] = bWrongSize def fun(): self.F[freq, 'b'] = bWrongSize
self.assertRaises(ValueError, fun) self.assertRaises(ValueError, fun)
def fun(): self.F[-999.] def fun(): self.F[-999.]
self.assertRaises(KeyError, fun) self.assertRaises(KeyError, fun)
def fun(): self.F['notRight'] def fun(): self.F['notRight']
self.assertRaises(KeyError, fun) self.assertRaises(KeyError, fun)
def fun(): self.F[freq, 'notThere'] def fun(): self.F[freq, 'notThere']
self.assertRaises(KeyError, fun) self.assertRaises(KeyError, fun)
@@ -271,7 +300,8 @@ class FieldsTest_Time(unittest.TestCase):
class FieldsTest_Time_Aliased(unittest.TestCase): class FieldsTest_Time_Aliased(unittest.TestCase):
def setUp(self): def setUp(self):
mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30]) mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10, 11, 12]],
[0, 0, -30])
x = np.linspace(5, 10, 3) x = np.linspace(5, 10, 3)
XYZ = Utils.ndgrid(x, x, np.r_[0.]) XYZ = Utils.ndgrid(x, x, np.r_[0.])
srcLoc = np.r_[0, 0, 0.] srcLoc = np.r_[0, 0, 0.]
@@ -288,9 +318,11 @@ class FieldsTest_Time_Aliased(unittest.TestCase):
survey = Survey.BaseSurvey(srcList=srcList) survey = Survey.BaseSurvey(srcList=srcList)
prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10., 3), (20., 2)]) prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10., 3), (20., 2)])
survey.pair(prob) survey.pair(prob)
def alias(b, srcInd, timeInd): def alias(b, srcInd, timeInd):
return self.F.mesh.edgeCurl.T * b + timeInd return self.F.mesh.edgeCurl.T * b + timeInd
self.F = Problem.TimeFields(mesh, survey, knownFields={'b':'F'}, aliasFields={'e':['b','E',alias]}) self.F = Problem.TimeFields(mesh, survey, knownFields={'b': 'F'},
aliasFields={'e': ['b', 'E', alias]})
self.Src0 = Src0 self.Src0 = Src0
self.Src1 = Src1 self.Src1 = Src1
self.mesh = mesh self.mesh = mesh
@@ -307,7 +339,6 @@ class FieldsTest_Time_Aliased(unittest.TestCase):
self.assertTrue('e' in F) self.assertTrue('e' in F)
self.assertTrue('b' in F) self.assertTrue('b' in F)
def test_simpleAlias(self): def test_simpleAlias(self):
F = self.F F = self.F
nSrc = F.survey.nSrc nSrc = F.survey.nSrc
@@ -325,7 +356,8 @@ class FieldsTest_Time_Aliased(unittest.TestCase):
self.assertTrue(np.all(F[self.Src0, 'e', :] == e[:, 0, :])) self.assertTrue(np.all(F[self.Src0, 'e', :] == e[:, 0, :]))
self.assertTrue(np.all(F[self.Src1, 'e', :] == e[:, 1, :])) self.assertTrue(np.all(F[self.Src1, 'e', :] == e[:, 1, :]))
for t in range(nT): for t in range(nT):
self.assertTrue(np.all(F[self.Src1, 'e', t] == Utils.mkvc(e[:,1,t],2) )) self.assertTrue(np.all(F[self.Src1, 'e', t] ==
Utils.mkvc(e[:, 1, t], 2)))
b = np.random.rand(F.mesh.nF, nT) b = np.random.rand(F.mesh.nF, nT)
F[self.Src0, 'b', :] = b F[self.Src0, 'b', :] = b
@@ -341,34 +373,41 @@ class FieldsTest_Time_Aliased(unittest.TestCase):
def test_aliasFunction(self): def test_aliasFunction(self):
nT = self.F.survey.prob.nT + 1 nT = self.F.survey.prob.nT + 1
count = [0] count = [0]
def alias(e, srcInd, timeInd): def alias(e, srcInd, timeInd):
count[0] += 1 count[0] += 1
self.assertTrue(srcInd[0] is self.Src0) self.assertTrue(srcInd[0] is self.Src0)
return self.F.mesh.edgeCurl * e return self.F.mesh.edgeCurl * e
F = Problem.TimeFields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]}) F = Problem.TimeFields(self.F.mesh, self.F.survey,
knownFields={'e': 'E'},
aliasFields={'b': ['e', 'F', alias]})
e = np.random.rand(F.mesh.nE, 1, nT) e = np.random.rand(F.mesh.nE, 1, nT)
F[self.Src0, 'e', :] = e F[self.Src0, 'e', :] = e
F[self.Src0, 'b', :] F[self.Src0, 'b', :]
self.assertTrue(count[0] == nT) # ensure that this is called for every time separately. # ensure that this is called for every time separately.
self.assertTrue(count[0] == nT)
e = np.random.rand(F.mesh.nE, 1, 1) e = np.random.rand(F.mesh.nE, 1, 1)
F[self.Src0, 'e', 1] = e F[self.Src0, 'e', 1] = e
count[0] = 0 count[0] = 0
F[self.Src0, 'b', 1] F[self.Src0, 'b', 1]
self.assertTrue(count[0] == 1) # ensure that this is called only once. self.assertTrue(count[0] == 1) # ensure that this is called only once.
def alias(e, srcInd, timeInd): def alias(e, srcInd, timeInd):
count[0] += 1 count[0] += 1
self.assertTrue(type(srcInd) is list) self.assertTrue(type(srcInd) is list)
self.assertTrue(srcInd[0] is self.Src0) self.assertTrue(srcInd[0] is self.Src0)
self.assertTrue(srcInd[1] is self.Src1) self.assertTrue(srcInd[1] is self.Src1)
return self.F.mesh.edgeCurl * e return self.F.mesh.edgeCurl * e
F = Problem.TimeFields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]}) F = Problem.TimeFields(self.F.mesh, self.F.survey,
knownFields={'e': 'E'},
aliasFields={'b': ['e', 'F', alias]})
e = np.random.rand(F.mesh.nE, 2, nT) e = np.random.rand(F.mesh.nE, 2, nT)
F[[self.Src0, self.Src1], 'e', :] = e F[[self.Src0, self.Src1], 'e', :] = e
count[0] = 0 count[0] = 0
F[[self.Src0, self.Src1], 'b', :] F[[self.Src0, self.Src1], 'b', :]
self.assertTrue(count[0] == nT) # ensure that this is called for every time separately.
# ensure that this is called for every time separately.
self.assertTrue(count[0] == nT)
e = np.random.rand(F.mesh.nE, 2, 1) e = np.random.rand(F.mesh.nE, 2, 1)
F[[self.Src0, self.Src1], 'e', 1] = e F[[self.Src0, self.Src1], 'e', 1] = e
count[0] = 0 count[0] = 0
+5 -5
View File
@@ -28,14 +28,14 @@ class RegularizationTests(unittest.TestCase):
for i, mesh in enumerate(self.meshlist): for i, mesh in enumerate(self.meshlist):
print 'Testing %iD'%mesh.dim print 'Testing {0:d}D'.format(mesh.dim)
mapping = r.mapPair(mesh) mapping = r.mapPair(mesh)
reg = r(mesh, mapping=mapping) reg = r(mesh, mapping=mapping)
m = np.random.rand(mapping.nP) m = np.random.rand(mapping.nP)
reg.mref = np.ones_like(m)*np.mean(m) reg.mref = np.ones_like(m)*np.mean(m)
print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref) print 'Check: phi_m (mref) = {0:f}'.format(reg.eval(reg.mref))
passed = reg.eval(reg.mref) < TOL passed = reg.eval(reg.mref) < TOL
self.assertTrue(passed) self.assertTrue(passed)
@@ -56,7 +56,7 @@ class RegularizationTests(unittest.TestCase):
for i, mesh in enumerate(self.meshlist): for i, mesh in enumerate(self.meshlist):
print 'Testing Active Cells %iD'%(mesh.dim) print 'Testing Active Cells {0:d}D'.format((mesh.dim))
if mesh.dim == 1: if mesh.dim == 1:
indActive = Utils.mkvc(mesh.gridCC <= 0.8) indActive = Utils.mkvc(mesh.gridCC <= 0.8)
@@ -70,7 +70,7 @@ class RegularizationTests(unittest.TestCase):
m = np.random.rand(mesh.nC)[indAct] m = np.random.rand(mesh.nC)[indAct]
reg.mref = np.ones_like(m)*np.mean(m) reg.mref = np.ones_like(m)*np.mean(m)
print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref) print 'Check: phi_m (mref) = {0:f}'.format(reg.eval(reg.mref))
passed = reg.eval(reg.mref) < TOL passed = reg.eval(reg.mref) < TOL
self.assertTrue(passed) self.assertTrue(passed)
@@ -87,7 +87,7 @@ class RegularizationTests(unittest.TestCase):
for i, mesh in enumerate(self.meshlist): for i, mesh in enumerate(self.meshlist):
print 'Testing %iD'%mesh.dim print 'Testing {0:d}D'.format(mesh.dim)
# mapping = r.mapPair(mesh) # mapping = r.mapPair(mesh)
# reg = r(mesh, mapping=mapping) # reg = r(mesh, mapping=mapping)
+11 -11
View File
@@ -14,9 +14,9 @@ class Doc_Test(unittest.TestCase):
html_path = os.path.sep.join(self.path_to_docs.split(os.path.sep) + ['_build']+['html']) html_path = os.path.sep.join(self.path_to_docs.split(os.path.sep) + ['_build']+['html'])
check = subprocess.call(["sphinx-build", "-nW", "-b", "html", "-d", check = subprocess.call(["sphinx-build", "-nW", "-b", "html", "-d",
"%s"%(doctrees_path) , "{0!s}".format((doctrees_path)) ,
"%s"%(self.path_to_docs), "{0!s}".format((self.path_to_docs)),
"%s"%(html_path)]) "{0!s}".format((html_path))])
assert check == 0 assert check == 0
# def test_latex(self): # def test_latex(self):
@@ -29,15 +29,15 @@ class Doc_Test(unittest.TestCase):
# "%s"%(latex_path)]) # "%s"%(latex_path)])
# assert check == 0 # assert check == 0
# def test_linkcheck(self): def test_linkcheck(self):
# doctrees_path = os.path.sep.join(self.path_to_docs.split(os.path.sep) + ['_build']+['doctrees']) doctrees_path = os.path.sep.join(self.path_to_docs.split(os.path.sep) + ['_build']+['doctrees'])
# link_path = os.path.sep.join(self.path_to_docs.split(os.path.sep) + ['_build']) link_path = os.path.sep.join(self.path_to_docs.split(os.path.sep) + ['_build'])
# check = subprocess.call(["sphinx-build", "-nW", "-b", "linkcheck", "-d", check = subprocess.call(["sphinx-build", "-nW", "-b", "linkcheck", "-d",
# "%s"%(doctrees_path), "%s"%(doctrees_path),
# "%s"%(self.path_to_docs), "%s"%(self.path_to_docs),
# "%s"%(link_path)]) "%s"%(link_path)])
# assert check == 0 assert check == 0
if __name__ == '__main__': if __name__ == '__main__':
unittest.main() unittest.main()
@@ -21,7 +21,7 @@ SrcList = ['RawVec', 'MagDipole'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVe
def adjointTest(fdemType, comp): def adjointTest(fdemType, comp):
prb = getFDEMProblem(fdemType, comp, SrcList, freq) prb = getFDEMProblem(fdemType, comp, SrcList, freq)
print 'Adjoint %s formulation - %s' % (fdemType, comp) print 'Adjoint {0!s} formulation - {1!s}'.format(fdemType, comp)
m = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY) m = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
mu = np.ones(prb.mesh.nC)*MU mu = np.ones(prb.mesh.nC)*MU
@@ -21,7 +21,7 @@ SrcList = ['RawVec', 'MagDipole'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVe
def adjointTest(fdemType, comp): def adjointTest(fdemType, comp):
prb = getFDEMProblem(fdemType, comp, SrcList, freq) prb = getFDEMProblem(fdemType, comp, SrcList, freq)
print 'Adjoint %s formulation - %s' % (fdemType, comp) print 'Adjoint {0!s} formulation - {1!s}'.format(fdemType, comp)
m = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY) m = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
mu = np.ones(prb.mesh.nC)*MU mu = np.ones(prb.mesh.nC)*MU
@@ -26,7 +26,7 @@ SrcType = ['MagDipole', 'RawVec'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVe
def derivTest(fdemType, comp): def derivTest(fdemType, comp):
prb = getFDEMProblem(fdemType, comp, SrcType, freq) prb = getFDEMProblem(fdemType, comp, SrcType, freq)
print '%s formulation - %s' % (fdemType, comp) print '{0!s} formulation - {1!s}'.format(fdemType, comp)
x0 = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY) x0 = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
mu = np.log(np.ones(prb.mesh.nC)*MU) mu = np.log(np.ones(prb.mesh.nC)*MU)
+1 -1
View File
@@ -55,7 +55,7 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=None,
if showIt == True: if showIt == True:
plt.loglog(rx.times[bz_calc>0], bz_calc[bz_calc>0], 'r', rx.times[bz_calc<0], -bz_calc[bz_calc<0], 'r--') plt.loglog(rx.times[bz_calc>0], bz_calc[bz_calc>0], 'r', rx.times[bz_calc<0], -bz_calc[bz_calc<0], 'r--')
plt.loglog(rx.times, abs(bz_ana), 'b*') plt.loglog(rx.times, abs(bz_ana), 'b*')
plt.title('sig_half = %e'%sig_half) plt.title('sig_half = {0:e}'.format(sig_half))
plt.show() plt.show()
return log10diff return log10diff
+1 -1
View File
@@ -25,7 +25,7 @@ class compareInitFiles(unittest.TestCase):
def get(test): def get(test):
def test_func(self): def test_func(self):
print '\nTesting %s.run(plotIt=False)\n'%test print '\nTesting {0!s}.run(plotIt=False)\n'.format(test)
getattr(Examples, test).run(plotIt=False) getattr(Examples, test).run(plotIt=False)
self.assertTrue(True) self.assertTrue(True)
return test_func return test_func
+3 -3
View File
@@ -121,7 +121,7 @@ class RichardsTests1D(unittest.TestCase):
tol = TOL*(10**int(np.log10(np.abs(zJv)))) tol = TOL*(10**int(np.log10(np.abs(zJv))))
passed = np.abs(vJz - zJv) < tol passed = np.abs(vJz - zJv) < tol
print 'Richards Adjoint Test - PressureHead' print 'Richards Adjoint Test - PressureHead'
print '%4.4e === %4.4e, diff=%4.4e < %4.e'%(vJz, zJv,np.abs(vJz - zJv),tol) print '{0:4.4e} === {1:4.4e}, diff={2:4.4e} < {3:4e}'.format(vJz, zJv, np.abs(vJz - zJv), tol)
self.assertTrue(passed,True) self.assertTrue(passed,True)
def test_Sensitivity(self): def test_Sensitivity(self):
@@ -193,7 +193,7 @@ class RichardsTests2D(unittest.TestCase):
tol = TOL*(10**int(np.log10(np.abs(zJv)))) tol = TOL*(10**int(np.log10(np.abs(zJv))))
passed = np.abs(vJz - zJv) < tol passed = np.abs(vJz - zJv) < tol
print '2D: Richards Adjoint Test - PressureHead' print '2D: Richards Adjoint Test - PressureHead'
print '%4.4e === %4.4e, diff=%4.4e < %4.e'%(vJz, zJv,np.abs(vJz - zJv),tol) print '{0:4.4e} === {1:4.4e}, diff={2:4.4e} < {3:4e}'.format(vJz, zJv, np.abs(vJz - zJv), tol)
self.assertTrue(passed,True) self.assertTrue(passed,True)
def test_Sensitivity(self): def test_Sensitivity(self):
@@ -265,7 +265,7 @@ class RichardsTests3D(unittest.TestCase):
tol = TOL*(10**int(np.log10(np.abs(zJv)))) tol = TOL*(10**int(np.log10(np.abs(zJv))))
passed = np.abs(vJz - zJv) < tol passed = np.abs(vJz - zJv) < tol
print '3D: Richards Adjoint Test - PressureHead' print '3D: Richards Adjoint Test - PressureHead'
print '%4.4e === %4.4e, diff=%4.4e < %4.e'%(vJz, zJv,np.abs(vJz - zJv),tol) print '{0:4.4e} === {1:4.4e}, diff={2:4.4e} < {3:4e}'.format(vJz, zJv, np.abs(vJz - zJv), tol)
self.assertTrue(passed,True) self.assertTrue(passed,True)
def test_Sensitivity(self): def test_Sensitivity(self):
+66 -30
View File
@@ -10,7 +10,9 @@ class BasicCurvTests(unittest.TestCase):
a = np.array([1, 1, 1]) a = np.array([1, 1, 1])
b = np.array([1, 2]) b = np.array([1, 2])
c = np.array([1, 4]) c = np.array([1, 4])
gridIt = lambda h: [np.cumsum(np.r_[0, x]) for x in h]
def gridIt(h): return [np.cumsum(np.r_[0, x]) for x in h]
X, Y = ndgrid(gridIt([a, b]), vector=False) X, Y = ndgrid(gridIt([a, b]), vector=False)
self.TM2 = TensorMesh([a, b]) self.TM2 = TensorMesh([a, b])
self.Curv2 = CurvilinearMesh([X, Y]) self.Curv2 = CurvilinearMesh([X, Y])
@@ -19,7 +21,10 @@ class BasicCurvTests(unittest.TestCase):
self.Curv3 = CurvilinearMesh([X, Y, Z]) self.Curv3 = CurvilinearMesh([X, Y, Z])
def test_area_3D(self): def test_area_3D(self):
test_area = np.array([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2]) test_area = np.array([1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4, 8, 8, 8, 8,
1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4,
4, 4, 1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2, 1, 1,
1, 2, 2, 2])
self.assertTrue(np.all(self.Curv3.area == test_area)) self.assertTrue(np.all(self.Curv3.area == test_area))
def test_vol_3D(self): def test_vol_3D(self):
@@ -33,54 +38,85 @@ class BasicCurvTests(unittest.TestCase):
self.assertTrue(t1) self.assertTrue(t1)
def test_edge_3D(self): def test_edge_3D(self):
test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]) test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2,
2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2,
2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4,
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4])
t1 = np.all(self.Curv3.edge == test_edge) t1 = np.all(self.Curv3.edge == test_edge)
self.assertTrue(t1) self.assertTrue(t1)
def test_edge_2D(self): def test_edge_2D(self):
test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2]) test_edge = np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2,
2])
t1 = np.all(self.Curv2.edge == test_edge) t1 = np.all(self.Curv2.edge == test_edge)
self.assertTrue(t1) self.assertTrue(t1)
def test_tangents(self): def test_tangents(self):
T = self.Curv2.tangents T = self.Curv2.tangents
self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.Curv2.nEx))) self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[0] ==
self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.Curv2.nEx))) np.ones(self.Curv2.nEx)))
self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.Curv2.nEy))) self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ex', 'V')[1] ==
self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.Curv2.nEy))) np.zeros(self.Curv2.nEx)))
self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[0] ==
np.zeros(self.Curv2.nEy)))
self.assertTrue(np.all(self.Curv2.r(T, 'E', 'Ey', 'V')[1] ==
np.ones(self.Curv2.nEy)))
T = self.Curv3.tangents T = self.Curv3.tangents
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.Curv3.nEx))) self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[0] ==
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.Curv3.nEx))) np.ones(self.Curv3.nEx)))
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.Curv3.nEx))) self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[1] ==
np.zeros(self.Curv3.nEx)))
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ex', 'V')[2] ==
np.zeros(self.Curv3.nEx)))
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.Curv3.nEy))) self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[0] ==
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.Curv3.nEy))) np.zeros(self.Curv3.nEy)))
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.Curv3.nEy))) self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[1] ==
np.ones(self.Curv3.nEy)))
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ey', 'V')[2] ==
np.zeros(self.Curv3.nEy)))
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.Curv3.nEz))) self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[0] ==
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.Curv3.nEz))) np.zeros(self.Curv3.nEz)))
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.Curv3.nEz))) self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[1] ==
np.zeros(self.Curv3.nEz)))
self.assertTrue(np.all(self.Curv3.r(T, 'E', 'Ez', 'V')[2] ==
np.ones(self.Curv3.nEz)))
def test_normals(self): def test_normals(self):
N = self.Curv2.normals N = self.Curv2.normals
self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.Curv2.nFx))) self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[0] ==
self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.Curv2.nFx))) np.ones(self.Curv2.nFx)))
self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.Curv2.nFy))) self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fx', 'V')[1] ==
self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.Curv2.nFy))) np.zeros(self.Curv2.nFx)))
self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[0] ==
np.zeros(self.Curv2.nFy)))
self.assertTrue(np.all(self.Curv2.r(N, 'F', 'Fy', 'V')[1] ==
np.ones(self.Curv2.nFy)))
N = self.Curv3.normals N = self.Curv3.normals
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.Curv3.nFx))) self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[0] ==
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.Curv3.nFx))) np.ones(self.Curv3.nFx)))
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.Curv3.nFx))) self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[1] ==
np.zeros(self.Curv3.nFx)))
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fx', 'V')[2] ==
np.zeros(self.Curv3.nFx)))
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.Curv3.nFy))) self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[0] ==
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.Curv3.nFy))) np.zeros(self.Curv3.nFy)))
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.Curv3.nFy))) self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[1] ==
np.ones(self.Curv3.nFy)))
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fy', 'V')[2] ==
np.zeros(self.Curv3.nFy)))
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.Curv3.nFz))) self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[0] ==
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.Curv3.nFz))) np.zeros(self.Curv3.nFz)))
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.Curv3.nFz))) self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[1] ==
np.zeros(self.Curv3.nFz)))
self.assertTrue(np.all(self.Curv3.r(N, 'F', 'Fz', 'V')[2] ==
np.ones(self.Curv3.nFz)))
def test_grid(self): def test_grid(self):
self.assertTrue(np.all(self.Curv2.gridCC == self.TM2.gridCC)) self.assertTrue(np.all(self.Curv2.gridCC == self.TM2.gridCC))
+109 -57
View File
@@ -2,12 +2,12 @@ import numpy as np
import scipy.sparse as sp import scipy.sparse as sp
import unittest import unittest
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from SimPEG import * from SimPEG import Mesh, Tests, Utils, Solver
MESHTYPES = ['uniformTensorMesh'] MESHTYPES = ['uniformTensorMesh']
def getxBCyBC_CC(mesh, alpha, beta, gamma): def getxBCyBC_CC(mesh, alpha, beta, gamma):
# def getxBCyBC(mesh, alpha, beta, gamma):
""" """
This is a subfunction generating mixed-boundary condition: This is a subfunction generating mixed-boundary condition:
@@ -17,13 +17,15 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
\rho \vec{j} = -\nabla \phi \phi \rho \vec{j} = -\nabla \phi \phi
\alpha \phi + \beta \frac{\partial \phi}{\partial r} = \gamma \ at \ r = \partial \Omega \alpha \phi + \beta \frac{\partial \phi}{\partial r} = \gamma \ at \ r
= \partial \Omega
xBC = f_1(\alpha, \beta, \gamma) xBC = f_1(\alpha, \beta, \gamma)
yBC = f(\alpha, \beta, \gamma) yBC = f(\alpha, \beta, \gamma)
Computes xBC and yBC for cell-centered discretizations Computes xBC and yBC for cell-centered discretizations
""" """
if mesh.dim == 1: # 1D if mesh.dim == 1: # 1D
if (len(alpha) != 2 or len(beta) != 2 or len(gamma) != 2): if (len(alpha) != 2 or len(beta) != 2 or len(gamma) != 2):
raise Exception("Lenght of list, alpha should be 2") raise Exception("Lenght of list, alpha should be 2")
@@ -65,8 +67,10 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
# h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0] # 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_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_xm = mesh.hx[0]*np.ones_like(alpha_xm)
h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp) h_xp = mesh.hx[-1]*np.ones_like(alpha_xp)
h_ym = mesh.hy[0]*np.ones_like(alpha_ym)
h_yp = mesh.hy[-1]*np.ones_like(alpha_yp)
a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) 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) b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
@@ -87,8 +91,10 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
yBC_ym = 0.5*(1.-b_ym) yBC_ym = 0.5*(1.-b_ym)
yBC_yp = 0.5*(1.-1./b_yp) yBC_yp = 0.5*(1.-1./b_yp)
sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]]) sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm],
sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]]) 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_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy] xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
@@ -116,9 +122,12 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
# h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1] # 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_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_xm = mesh.hx[0]*np.ones_like(alpha_xm)
h_ym, h_yp = mesh.hy[0]*np.ones_like(alpha_ym), mesh.hy[-1]*np.ones_like(alpha_yp) h_xp = mesh.hx[-1]*np.ones_like(alpha_xp)
h_zm, h_zp = mesh.hz[0]*np.ones_like(alpha_zm), mesh.hz[-1]*np.ones_like(alpha_zp) h_ym = mesh.hy[0]*np.ones_like(alpha_ym)
h_yp = mesh.hy[-1]*np.ones_like(alpha_yp)
h_zm = mesh.hz[0]*np.ones_like(alpha_zm)
h_zp = mesh.hz[-1]*np.ones_like(alpha_zp)
a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm) 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) b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
@@ -148,9 +157,12 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
yBC_zm = 0.5*(1.-b_zm) yBC_zm = 0.5*(1.-b_zm)
yBC_zp = 0.5*(1.-1./b_zp) yBC_zp = 0.5*(1.-1./b_zp)
sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]]) sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm],
sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]]) np.arange(mesh.nFx)[fxp]])
sortindsfz = np.argsort(np.r_[np.arange(mesh.nFz)[fzm], np.arange(mesh.nFz)[fzp]]) 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_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy] xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
@@ -165,6 +177,7 @@ def getxBCyBC_CC(mesh, alpha, beta, gamma):
return xBC, yBC return xBC, yBC
class Test1D_InhomogeneousMixed(Tests.OrderTest): class Test1D_InhomogeneousMixed(Tests.OrderTest):
name = "1D - Mixed" name = "1D - Mixed"
meshTypes = MESHTYPES meshTypes = MESHTYPES
@@ -174,10 +187,13 @@ class Test1D_InhomogeneousMixed(Tests.OrderTest):
def getError(self): def getError(self):
# Test function # Test function
phi_fun = lambda x: np.cos(np.pi*x) def phi_fun(x): return np.cos(np.pi*x)
j_fun = lambda x: np.pi*np.sin(np.pi*x)
phi_deriv = lambda x: -j_fun(x) def j_fun(x): return np.pi*np.sin(np.pi*x)
q_fun = lambda x: (np.pi**2)*np.cos(np.pi*x)
def phi_deriv(x): return -j_fun(x)
def q_fun(x): return (np.pi**2)*np.cos(np.pi*x)
xc_ana = phi_fun(self.M.gridCC) xc_ana = phi_fun(self.M.gridCC)
q_ana = q_fun(self.M.gridCC) q_ana = q_fun(self.M.gridCC)
@@ -199,7 +215,6 @@ class Test1D_InhomogeneousMixed(Tests.OrderTest):
gamma = alpha*phi_bc + beta*phi_deriv_bc gamma = alpha*phi_bc + beta*phi_deriv_bc
x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma) x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma)
sigma = np.ones(self.M.nC) sigma = np.ones(self.M.nC)
Mfrho = self.M.getFaceInnerProduct(1./sigma) Mfrho = self.M.getFaceInnerProduct(1./sigma)
MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True) MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True)
@@ -222,13 +237,13 @@ class Test1D_InhomogeneousMixed(Tests.OrderTest):
NotImplementedError NotImplementedError
return err return err
def test_order(self): def test_order(self):
print "==== Testing Mixed boudary conduction for CC-problem ====" print "==== Testing Mixed boudary conduction for CC-problem ===="
self.name = "1D" self.name = "1D"
self.myTest = 'xc' self.myTest = 'xc'
self.orderTest() self.orderTest()
class Test2D_InhomogeneousMixed(Tests.OrderTest): class Test2D_InhomogeneousMixed(Tests.OrderTest):
name = "2D - Mixed" name = "2D - Mixed"
meshTypes = MESHTYPES meshTypes = MESHTYPES
@@ -238,12 +253,23 @@ class Test2D_InhomogeneousMixed(Tests.OrderTest):
def getError(self): def getError(self):
# Test function # Test function
phi_fun = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1]) def phi_fun(x):
j_funX = lambda x: +np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1]) return np.cos(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) def j_funX(x):
phideriv_funY = lambda x: -j_funY(x) return +np.pi*np.sin(np.pi*x[:, 0])*np.cos(np.pi*x[:, 1])
q_fun = lambda x: +2*(np.pi**2)*phi_fun(x)
def j_funY(x):
return +np.pi*np.cos(np.pi*x[:, 0])*np.sin(np.pi*x[:, 1])
def phideriv_funX(x):
return -j_funX(x)
def phideriv_funY(x):
return -j_funY(x)
def q_fun(x):
return +2*(np.pi**2)*phi_fun(x)
xc_ana = phi_fun(self.M.gridCC) xc_ana = phi_fun(self.M.gridCC)
q_ana = q_fun(self.M.gridCC) q_ana = q_fun(self.M.gridCC)
@@ -259,18 +285,26 @@ class Test2D_InhomogeneousMixed(Tests.OrderTest):
gBFyp = self.M.gridFy[fyp, :] gBFyp = self.M.gridFy[fyp, :]
# Setup Mixed B.C (alpha, beta, gamma) # Setup Mixed B.C (alpha, beta, gamma)
alpha_xm, alpha_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) alpha_xm = np.ones_like(gBFxm[:, 0])
beta_xm, beta_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) alpha_xp = np.ones_like(gBFxp[:, 0])
alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) beta_xm = np.ones_like(gBFxm[:, 0])
beta_ym, beta_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) beta_xp = np.ones_like(gBFxp[:, 0])
alpha_ym = np.ones_like(gBFym[:, 1])
alpha_yp = np.ones_like(gBFyp[:, 1])
beta_ym = np.ones_like(gBFym[:, 1])
beta_yp = np.ones_like(gBFyp[:, 1])
phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp) 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_ym, phi_bc_yp = phi_fun(gBFym), phi_fun(gBFyp)
phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp) phiderivX_bc_xm = phideriv_funX(gBFxm)
phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp) phiderivX_bc_xp = phideriv_funX(gBFxp)
phiderivY_bc_ym = phideriv_funY(gBFym)
phiderivY_bc_yp = phideriv_funY(gBFyp)
def gamma_fun(alpha, beta, phi, phi_deriv):
return alpha*phi + beta*phi_deriv
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_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_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_ym = gamma_fun(alpha_ym, beta_ym, phi_bc_ym, phiderivY_bc_ym)
@@ -282,7 +316,6 @@ class Test2D_InhomogeneousMixed(Tests.OrderTest):
x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma) x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma)
sigma = np.ones(self.M.nC) sigma = np.ones(self.M.nC)
Mfrho = self.M.getFaceInnerProduct(1./sigma) Mfrho = self.M.getFaceInnerProduct(1./sigma)
MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True) MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True)
@@ -303,13 +336,13 @@ class Test2D_InhomogeneousMixed(Tests.OrderTest):
NotImplementedError NotImplementedError
return err return err
def test_order(self): def test_order(self):
print "==== Testing Mixed boudary conduction for CC-problem ====" print "==== Testing Mixed boudary conduction for CC-problem ===="
self.name = "2D" self.name = "2D"
self.myTest = 'xc' self.myTest = 'xc'
self.orderTest() self.orderTest()
class Test3D_InhomogeneousMixed(Tests.OrderTest): class Test3D_InhomogeneousMixed(Tests.OrderTest):
name = "3D - Mixed" name = "3D - Mixed"
meshTypes = MESHTYPES meshTypes = MESHTYPES
@@ -319,16 +352,29 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest):
def getError(self): def getError(self):
# Test function # Test function
phi_fun = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.cos(np.pi*x[:,2]) def phi_fun(x):
j_funX = lambda x: +np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.cos(np.pi*x[:,2]) return (np.cos(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])*np.cos(np.pi*x[:,2]) 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) def j_funX(x):
phideriv_funY = lambda x: -j_funY(x) return (np.pi*np.sin(np.pi*x[:, 0])*np.cos(np.pi*x[:, 1]) *
phideriv_funZ = lambda x: -j_funZ(x) np.cos(np.pi*x[:, 2]))
q_fun = lambda x: 3*(np.pi**2)*phi_fun(x) def j_funY(x):
return (np.pi*np.cos(np.pi*x[:, 0])*np.sin(np.pi*x[:, 1]) *
np.cos(np.pi*x[:, 2]))
def j_funZ(x):
return (np.pi*np.cos(np.pi*x[:, 0])*np.cos(np.pi*x[:, 1]) *
np.sin(np.pi*x[:, 2]))
def phideriv_funX(x): return -j_funX(x)
def phideriv_funY(x): return -j_funY(x)
def phideriv_funZ(x): return -j_funZ(x)
def q_fun(x): return 3*(np.pi**2)*phi_fun(x)
xc_ana = phi_fun(self.M.gridCC) xc_ana = phi_fun(self.M.gridCC)
q_ana = q_fun(self.M.gridCC) q_ana = q_fun(self.M.gridCC)
@@ -346,23 +392,33 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest):
gBFzp = self.M.gridFz[fzp, :] gBFzp = self.M.gridFz[fzp, :]
# Setup Mixed B.C (alpha, beta, gamma) # Setup Mixed B.C (alpha, beta, gamma)
alpha_xm, alpha_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) alpha_xm = np.ones_like(gBFxm[:, 0])
beta_xm, beta_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) alpha_xp = np.ones_like(gBFxp[:, 0])
alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) beta_xm = np.ones_like(gBFxm[:, 0])
beta_ym, beta_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) beta_xp = np.ones_like(gBFxp[:, 0])
alpha_zm, alpha_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2]) alpha_ym = np.ones_like(gBFym[:, 1])
beta_zm, beta_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2]) alpha_yp = np.ones_like(gBFyp[:, 1])
beta_ym = np.ones_like(gBFym[:, 1])
beta_yp = np.ones_like(gBFyp[:, 1])
alpha_zm = np.ones_like(gBFzm[:, 2])
alpha_zp = np.ones_like(gBFzp[:, 2])
beta_zm = np.ones_like(gBFzm[:, 2])
beta_zp = np.ones_like(gBFzp[:, 2])
phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp) 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_ym, phi_bc_yp = phi_fun(gBFym), phi_fun(gBFyp)
phi_bc_zm, phi_bc_zp = phi_fun(gBFzm), phi_fun(gBFzp) phi_bc_zm, phi_bc_zp = phi_fun(gBFzm), phi_fun(gBFzp)
phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp) phiderivX_bc_xm = phideriv_funX(gBFxm)
phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp) phiderivX_bc_xp = phideriv_funX(gBFxp)
phiderivY_bc_zm, phiderivY_bc_zp = phideriv_funZ(gBFzm), phideriv_funZ(gBFzp) phiderivY_bc_ym = phideriv_funY(gBFym)
phiderivY_bc_yp = phideriv_funY(gBFyp)
phiderivY_bc_zm = phideriv_funZ(gBFzm)
phiderivY_bc_zp = phideriv_funZ(gBFzp)
def gamma_fun(alpha, beta, phi, phi_deriv):
return alpha*phi + beta*phi_deriv
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_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_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_ym = gamma_fun(alpha_ym, beta_ym, phi_bc_ym, phiderivY_bc_ym)
@@ -376,7 +432,6 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest):
x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma) x_BC, y_BC = getxBCyBC_CC(self.M, alpha, beta, gamma)
sigma = np.ones(self.M.nC) sigma = np.ones(self.M.nC)
Mfrho = self.M.getFaceInnerProduct(1./sigma) Mfrho = self.M.getFaceInnerProduct(1./sigma)
MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True) MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True)
@@ -398,14 +453,11 @@ class Test3D_InhomogeneousMixed(Tests.OrderTest):
NotImplementedError NotImplementedError
return err return err
def test_order(self): def test_order(self):
print "==== Testing Mixed boudary conduction for CC-problem ====" print "==== Testing Mixed boudary conduction for CC-problem ===="
self.name = "3D" self.name = "3D"
self.myTest = 'xc' self.myTest = 'xc'
self.orderTest() self.orderTest()
if __name__ == '__main__': if __name__ == '__main__':
unittest.main() unittest.main()
-2
View File
@@ -2,8 +2,6 @@ import numpy as np
import unittest import unittest
from SimPEG.Utils import mkvc from SimPEG.Utils import mkvc
from SimPEG import Mesh, Tests from SimPEG import Mesh, Tests
import unittest
MESHTYPES = ['uniformTensorMesh', 'randomTensorMesh'] MESHTYPES = ['uniformTensorMesh', 'randomTensorMesh']
TOLERANCES = [0.9, 0.5, 0.5] TOLERANCES = [0.9, 0.5, 0.5]
+1 -1
View File
@@ -52,7 +52,7 @@ class Tests(unittest.TestCase):
assert o <= -1 assert o <= -1
assert not (o > -1) assert not (o > -1)
assert o >= -1 assert o >= -1
assert -(-o)*o == -o assert -1.*(-o)*o == -o
o = Identity() o = Identity()
assert +o == o assert +o == o
assert -o == -o assert -o == -o