mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-25 13:30:06 +08:00
Merge branch 'master' into em/ZeroIdentity
This commit is contained in:
+5
-1
@@ -5,7 +5,11 @@ python:
|
||||
sudo: false
|
||||
|
||||
env:
|
||||
- TEST_DIR=tests/em
|
||||
- TEST_DIR=tests/em/examples
|
||||
- TEST_DIR=tests/em/fdem/forward
|
||||
- TEST_DIR=tests/em/fdem/inverse/derivs
|
||||
- TEST_DIR=tests/em/fdem/inverse/adjoint
|
||||
- TEST_DIR=tests/em/tdem
|
||||
- TEST_DIR=tests/mesh
|
||||
- TEST_DIR=tests/flow
|
||||
- TEST_DIR=tests/utils
|
||||
|
||||
@@ -0,0 +1,75 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
from SimPEG import EM
|
||||
import sys
|
||||
from scipy.constants import mu_0
|
||||
|
||||
def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False):
|
||||
cs = 5.
|
||||
ncx, ncy, ncz = 6, 6, 6
|
||||
npad = 3
|
||||
hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)]
|
||||
hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
|
||||
hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
|
||||
mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C'])
|
||||
|
||||
mapping = Maps.ExpMap(mesh)
|
||||
|
||||
x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the source
|
||||
XYZ = Utils.ndgrid(x,x,np.r_[0.])
|
||||
Rx0 = EM.FDEM.RxFDEM(XYZ, comp)
|
||||
|
||||
Src = []
|
||||
|
||||
for SrcType in SrcList:
|
||||
if SrcType is 'MagDipole':
|
||||
Src.append(EM.FDEM.SrcFDEM_MagDipole([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
|
||||
elif SrcType is 'MagDipole_Bfield':
|
||||
Src.append(EM.FDEM.SrcFDEM_MagDipole_Bfield([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
|
||||
elif SrcType is 'CircularLoop':
|
||||
Src.append(EM.FDEM.SrcFDEM_CircularLoop([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
|
||||
elif SrcType is 'RawVec':
|
||||
if fdemType is 'e' or fdemType is 'b':
|
||||
S_m = np.zeros(mesh.nF)
|
||||
S_e = np.zeros(mesh.nE)
|
||||
S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1.
|
||||
S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1.
|
||||
Src.append(EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e))
|
||||
|
||||
elif fdemType is 'h' or fdemType is 'j':
|
||||
S_m = np.zeros(mesh.nE)
|
||||
S_e = np.zeros(mesh.nF)
|
||||
S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1.
|
||||
S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1.
|
||||
Src.append(EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e))
|
||||
|
||||
if verbose:
|
||||
print ' Fetching %s problem' % (fdemType)
|
||||
|
||||
if fdemType == 'e':
|
||||
survey = EM.FDEM.SurveyFDEM(Src)
|
||||
prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping)
|
||||
|
||||
elif fdemType == 'b':
|
||||
survey = EM.FDEM.SurveyFDEM(Src)
|
||||
prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping)
|
||||
|
||||
elif fdemType == 'j':
|
||||
survey = EM.FDEM.SurveyFDEM(Src)
|
||||
prb = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping)
|
||||
|
||||
elif fdemType == 'h':
|
||||
survey = EM.FDEM.SurveyFDEM(Src)
|
||||
prb = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping)
|
||||
|
||||
else:
|
||||
raise NotImplementedError()
|
||||
prb.pair(survey)
|
||||
|
||||
try:
|
||||
from pymatsolver import MumpsSolver
|
||||
prb.Solver = MumpsSolver
|
||||
except ImportError, e:
|
||||
pass
|
||||
|
||||
return prb
|
||||
@@ -7,4 +7,4 @@ from ipythonutils import easyAnimate as animate
|
||||
from CounterUtils import *
|
||||
import ModelBuilder
|
||||
import SolverUtils
|
||||
|
||||
from coordutils import *
|
||||
|
||||
@@ -0,0 +1,62 @@
|
||||
import numpy as np
|
||||
from SimPEG.Utils import mkvc
|
||||
|
||||
def rotationMatrixFromNormals(v0,v1,tol=1e-20):
|
||||
"""
|
||||
Performs the minimum number of rotations to define a rotation from the direction indicated by the vector n0 to the direction indicated by n1.
|
||||
The axis of rotation is n0 x n1
|
||||
https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula
|
||||
|
||||
:param numpy.array v0: vector of length 3
|
||||
:param numpy.array v1: vector of length 3
|
||||
:param tol = 1e-20: tolerance. If the norm of the cross product between the two vectors is below this, no rotation is performed
|
||||
:rtype: numpy.array, 3x3
|
||||
:return: rotation matrix which rotates the frame so that n0 is aligned with n1
|
||||
|
||||
"""
|
||||
|
||||
# ensure both n0, n1 are vectors of length 1
|
||||
assert len(v0) == 3, "Length of n0 should be 3"
|
||||
assert len(v1) == 3, "Length of n1 should be 3"
|
||||
|
||||
# ensure both are true normals
|
||||
n0 = v0*1./np.linalg.norm(v0)
|
||||
n1 = v1*1./np.linalg.norm(v1)
|
||||
|
||||
n0dotn1 = n0.dot(n1)
|
||||
|
||||
# define the rotation axis, which is the cross product of the two vectors
|
||||
rotAx = np.cross(n0,n1)
|
||||
|
||||
if np.linalg.norm(rotAx) < tol:
|
||||
return np.eye(3,dtype=float)
|
||||
|
||||
rotAx *= 1./np.linalg.norm(rotAx)
|
||||
|
||||
cosT = n0dotn1/(np.linalg.norm(n0)*np.linalg.norm(n1))
|
||||
sinT = np.sqrt(1.-n0dotn1**2)
|
||||
|
||||
ux = np.array([[0., -rotAx[2], rotAx[1]], [rotAx[2], 0., -rotAx[0]], [-rotAx[1], rotAx[0], 0.]],dtype=float)
|
||||
|
||||
return np.eye(3,dtype=float) + sinT*ux + (1.-cosT)*(ux.dot(ux))
|
||||
|
||||
|
||||
def rotatePointsFromNormals(XYZ,n0,n1,x0=np.r_[0.,0.,0.]):
|
||||
"""
|
||||
rotates a grid so that the vector n0 is aligned with the vector n1
|
||||
|
||||
:param numpy.array n0: vector of length 3, should have norm 1
|
||||
:param numpy.array n1: vector of length 3, should have norm 1
|
||||
:param numpy.array x0: vector of length 3, point about which we perform the rotation
|
||||
:rtype: numpy.array, 3x3
|
||||
:return: rotation matrix which rotates the frame so that n0 is aligned with n1
|
||||
"""
|
||||
|
||||
R = rotationMatrixFromNormals(n0, n1)
|
||||
|
||||
assert XYZ.shape[1] == 3, "Grid XYZ should be 3 wide"
|
||||
assert len(x0) == 3, "x0 should have length 3"
|
||||
|
||||
X0 = np.ones([XYZ.shape[0],1])*mkvc(x0)
|
||||
|
||||
return (XYZ - X0).dot(R.T) + X0 # equivalent to (R*(XYZ - X0)).T + X0
|
||||
@@ -0,0 +1,11 @@
|
||||
if __name__ == '__main__':
|
||||
import os
|
||||
import glob
|
||||
import unittest
|
||||
test_file_strings = glob.glob('test_*.py')
|
||||
module_strings = [str[0:len(str)-3] for str in test_file_strings]
|
||||
suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
|
||||
in module_strings]
|
||||
testSuite = unittest.TestSuite(suites)
|
||||
|
||||
unittest.TextTestRunner(verbosity=2).run(testSuite)
|
||||
@@ -0,0 +1,127 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
from SimPEG import EM
|
||||
import sys
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.EM.Utils.testingUtils import getFDEMProblem
|
||||
|
||||
testEB = True
|
||||
testHJ = True
|
||||
|
||||
verbose = False
|
||||
|
||||
TOL = 1e-5
|
||||
FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
|
||||
CONDUCTIVITY = 1e1
|
||||
MU = mu_0
|
||||
freq = 1e-1
|
||||
addrandoms = True
|
||||
|
||||
SrcList = ['RawVec', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop']
|
||||
|
||||
|
||||
def crossCheckTest(fdemType, comp):
|
||||
|
||||
l2norm = lambda r: np.sqrt(r.dot(r))
|
||||
|
||||
prb1 = getFDEMProblem(fdemType, comp, SrcList, freq, verbose)
|
||||
mesh = prb1.mesh
|
||||
print 'Cross Checking Forward: %s formulation - %s' % (fdemType, comp)
|
||||
m = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
|
||||
mu = np.log(np.ones(mesh.nC)*MU)
|
||||
|
||||
if addrandoms is True:
|
||||
m = m + np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1
|
||||
mu = mu + np.random.randn(mesh.nC)*MU*1e-1
|
||||
|
||||
# prb1.PropMap.PropModel.mu = mu
|
||||
# prb1.PropMap.PropModel.mui = 1./mu
|
||||
survey1 = prb1.survey
|
||||
d1 = survey1.dpred(m)
|
||||
|
||||
if verbose:
|
||||
print ' Problem 1 solved'
|
||||
|
||||
if fdemType == 'e':
|
||||
prb2 = getFDEMProblem('b', comp, SrcList, freq, verbose)
|
||||
elif fdemType == 'b':
|
||||
prb2 = getFDEMProblem('e', comp, SrcList, freq, verbose)
|
||||
elif fdemType == 'j':
|
||||
prb2 = getFDEMProblem('h', comp, SrcList, freq, verbose)
|
||||
elif fdemType == 'h':
|
||||
prb2 = getFDEMProblem('j', comp, SrcList, freq, verbose)
|
||||
else:
|
||||
raise NotImplementedError()
|
||||
|
||||
# prb2.mu = mu
|
||||
survey2 = prb2.survey
|
||||
d2 = survey2.dpred(m)
|
||||
|
||||
if verbose:
|
||||
print ' Problem 2 solved'
|
||||
|
||||
r = d2-d1
|
||||
l2r = l2norm(r)
|
||||
|
||||
tol = np.max([TOL*(10**int(np.log10(l2norm(d1)))),FLR])
|
||||
print l2norm(d1), l2norm(d2), l2r , tol, l2r < tol
|
||||
return l2r < tol
|
||||
|
||||
|
||||
class FDEM_CrossCheck(unittest.TestCase):
|
||||
if testEB:
|
||||
def test_EB_CrossCheck_exr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'exr'))
|
||||
def test_EB_CrossCheck_eyr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'eyr'))
|
||||
def test_EB_CrossCheck_ezr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'ezr'))
|
||||
def test_EB_CrossCheck_exi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'exi'))
|
||||
def test_EB_CrossCheck_eyi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'eyi'))
|
||||
def test_EB_CrossCheck_ezi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'ezi'))
|
||||
|
||||
def test_EB_CrossCheck_bxr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bxr'))
|
||||
def test_EB_CrossCheck_byr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'byr'))
|
||||
def test_EB_CrossCheck_bzr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bzr'))
|
||||
def test_EB_CrossCheck_bxi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bxi'))
|
||||
def test_EB_CrossCheck_byi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'byi'))
|
||||
def test_EB_CrossCheck_bzi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bzi'))
|
||||
|
||||
if testHJ:
|
||||
def test_HJ_CrossCheck_jxr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jxr'))
|
||||
def test_HJ_CrossCheck_jyr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jyr'))
|
||||
def test_HJ_CrossCheck_jzr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jzr'))
|
||||
def test_HJ_CrossCheck_jxi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jxi'))
|
||||
def test_HJ_CrossCheck_jyi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jyi'))
|
||||
def test_HJ_CrossCheck_jzi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jzi'))
|
||||
|
||||
def test_HJ_CrossCheck_hxr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hxr'))
|
||||
def test_HJ_CrossCheck_hyr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hyr'))
|
||||
def test_HJ_CrossCheck_hzr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hzr'))
|
||||
def test_HJ_CrossCheck_hxi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hxi'))
|
||||
def test_HJ_CrossCheck_hyi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hyi'))
|
||||
def test_HJ_CrossCheck_hzi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hzi'))
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,11 @@
|
||||
if __name__ == '__main__':
|
||||
import os
|
||||
import glob
|
||||
import unittest
|
||||
test_file_strings = glob.glob('test_*.py')
|
||||
module_strings = [str[0:len(str)-3] for str in test_file_strings]
|
||||
suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
|
||||
in module_strings]
|
||||
testSuite = unittest.TestSuite(suites)
|
||||
|
||||
unittest.TextTestRunner(verbosity=2).run(testSuite)
|
||||
@@ -0,0 +1,156 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
from SimPEG import EM
|
||||
import sys
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.EM.Utils.testingUtils import getFDEMProblem
|
||||
|
||||
testEB = True
|
||||
testHJ = True
|
||||
|
||||
verbose = False
|
||||
|
||||
TOL = 1e-5
|
||||
FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
|
||||
CONDUCTIVITY = 1e1
|
||||
MU = mu_0
|
||||
freq = 1e-1
|
||||
addrandoms = True
|
||||
|
||||
SrcType = 'RawVec' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
|
||||
|
||||
def adjointTest(fdemType, comp):
|
||||
prb = getFDEMProblem(fdemType, comp, [SrcType], freq)
|
||||
print 'Adjoint %s formulation - %s' % (fdemType, comp)
|
||||
|
||||
m = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY)
|
||||
mu = np.ones(prb.mesh.nC)*MU
|
||||
|
||||
if addrandoms is True:
|
||||
m = m + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1
|
||||
mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1
|
||||
|
||||
survey = prb.survey
|
||||
# prb.PropMap.PropModel.mu = mu
|
||||
# prb.PropMap.PropModel.mui = 1./mu
|
||||
u = prb.fields(m)
|
||||
|
||||
v = np.random.rand(survey.nD)
|
||||
w = np.random.rand(prb.mesh.nC)
|
||||
|
||||
vJw = v.dot(prb.Jvec(m, w, u))
|
||||
wJtv = w.dot(prb.Jtvec(m, v, u))
|
||||
tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR])
|
||||
print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
|
||||
return np.abs(vJw - wJtv) < tol
|
||||
|
||||
class FDEM_AdjointTests(unittest.TestCase):
|
||||
if testEB:
|
||||
def test_Jtvec_adjointTest_exr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'exr'))
|
||||
def test_Jtvec_adjointTest_eyr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'eyr'))
|
||||
def test_Jtvec_adjointTest_ezr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'ezr'))
|
||||
def test_Jtvec_adjointTest_exi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'exi'))
|
||||
def test_Jtvec_adjointTest_eyi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'eyi'))
|
||||
def test_Jtvec_adjointTest_ezi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'ezi'))
|
||||
|
||||
def test_Jtvec_adjointTest_bxr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'bxr'))
|
||||
def test_Jtvec_adjointTest_byr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'byr'))
|
||||
def test_Jtvec_adjointTest_bzr_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'bzr'))
|
||||
def test_Jtvec_adjointTest_bxi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'bxi'))
|
||||
def test_Jtvec_adjointTest_byi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'byi'))
|
||||
def test_Jtvec_adjointTest_bzi_Eform(self):
|
||||
self.assertTrue(adjointTest('e', 'bzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_exr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'exr'))
|
||||
def test_Jtvec_adjointTest_eyr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'eyr'))
|
||||
def test_Jtvec_adjointTest_ezr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'ezr'))
|
||||
def test_Jtvec_adjointTest_exi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'exi'))
|
||||
def test_Jtvec_adjointTest_eyi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'eyi'))
|
||||
def test_Jtvec_adjointTest_ezi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'ezi'))
|
||||
def test_Jtvec_adjointTest_bxr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'bxr'))
|
||||
def test_Jtvec_adjointTest_byr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'byr'))
|
||||
def test_Jtvec_adjointTest_bzr_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'bzr'))
|
||||
def test_Jtvec_adjointTest_bxi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'bxi'))
|
||||
def test_Jtvec_adjointTest_byi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'byi'))
|
||||
def test_Jtvec_adjointTest_bzi_Bform(self):
|
||||
self.assertTrue(adjointTest('b', 'bzi'))
|
||||
|
||||
|
||||
if testHJ:
|
||||
def test_Jtvec_adjointTest_jxr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jxr'))
|
||||
def test_Jtvec_adjointTest_jyr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jyr'))
|
||||
def test_Jtvec_adjointTest_jzr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jzr'))
|
||||
def test_Jtvec_adjointTest_jxi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jxi'))
|
||||
def test_Jtvec_adjointTest_jyi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jyi'))
|
||||
def test_Jtvec_adjointTest_jzi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'jzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_hxr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hxr'))
|
||||
def test_Jtvec_adjointTest_hyr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hyr'))
|
||||
def test_Jtvec_adjointTest_hzr_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hzr'))
|
||||
def test_Jtvec_adjointTest_hxi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hxi'))
|
||||
def test_Jtvec_adjointTest_hyi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hyi'))
|
||||
def test_Jtvec_adjointTest_hzi_Jform(self):
|
||||
self.assertTrue(adjointTest('j', 'hzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_hxr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hxr'))
|
||||
def test_Jtvec_adjointTest_hyr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hyr'))
|
||||
def test_Jtvec_adjointTest_hzr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hzr'))
|
||||
def test_Jtvec_adjointTest_hxi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hxi'))
|
||||
def test_Jtvec_adjointTest_hyi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hyi'))
|
||||
def test_Jtvec_adjointTest_hzi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'hzi'))
|
||||
|
||||
def test_Jtvec_adjointTest_hxr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jxr'))
|
||||
def test_Jtvec_adjointTest_hyr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jyr'))
|
||||
def test_Jtvec_adjointTest_hzr_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jzr'))
|
||||
def test_Jtvec_adjointTest_hxi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jxi'))
|
||||
def test_Jtvec_adjointTest_hyi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jyi'))
|
||||
def test_Jtvec_adjointTest_hzi_Hform(self):
|
||||
self.assertTrue(adjointTest('h', 'jzi'))
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,11 @@
|
||||
if __name__ == '__main__':
|
||||
import os
|
||||
import glob
|
||||
import unittest
|
||||
test_file_strings = glob.glob('test_*.py')
|
||||
module_strings = [str[0:len(str)-3] for str in test_file_strings]
|
||||
suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
|
||||
in module_strings]
|
||||
testSuite = unittest.TestSuite(suites)
|
||||
|
||||
unittest.TextTestRunner(verbosity=2).run(testSuite)
|
||||
@@ -5,7 +5,6 @@ import sys
|
||||
from scipy.constants import mu_0
|
||||
|
||||
testDerivs = True
|
||||
testCrossCheck = True
|
||||
testAdjoint = True
|
||||
testEB = True
|
||||
testHJ = True
|
||||
@@ -149,57 +148,6 @@ def derivTest(fdemType, comp):
|
||||
return Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR)
|
||||
|
||||
|
||||
def crossCheckTest(fdemType, comp):
|
||||
|
||||
l2norm = lambda r: np.sqrt(r.dot(r))
|
||||
|
||||
prb1 = getProblem(fdemType, comp)
|
||||
mesh = prb1.mesh
|
||||
print 'Cross Checking Forward: %s formulation - %s' % (fdemType, comp)
|
||||
m = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
|
||||
mu = np.log(np.ones(mesh.nC)*MU)
|
||||
|
||||
if addrandoms is True:
|
||||
m = m + np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1
|
||||
mu = mu + np.random.randn(mesh.nC)*MU*1e-1
|
||||
|
||||
# prb1.PropMap.PropModel.mu = mu
|
||||
# prb1.PropMap.PropModel.mui = 1./mu
|
||||
survey1 = prb1.survey
|
||||
d1 = survey1.dpred(m)
|
||||
|
||||
if verbose:
|
||||
print ' Problem 1 solved'
|
||||
|
||||
if fdemType == 'e':
|
||||
prb2 = getProblem('b', comp)
|
||||
elif fdemType == 'b':
|
||||
prb2 = getProblem('e', comp)
|
||||
elif fdemType == 'j':
|
||||
prb2 = getProblem('h', comp)
|
||||
elif fdemType == 'h':
|
||||
prb2 = getProblem('j', comp)
|
||||
else:
|
||||
raise NotImplementedError()
|
||||
|
||||
# prb2.mu = mu
|
||||
survey2 = prb2.survey
|
||||
d2 = survey2.dpred(m)
|
||||
|
||||
if verbose:
|
||||
print ' Problem 2 solved'
|
||||
|
||||
r = d2-d1
|
||||
l2r = l2norm(r)
|
||||
|
||||
tol = np.max([TOL*(10**int(np.log10(l2norm(d1)))),FLR])
|
||||
print l2norm(d1), l2norm(d2), l2r , tol, l2r < tol
|
||||
return l2r < tol
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
class FDEM_DerivTests(unittest.TestCase):
|
||||
|
||||
if testDerivs:
|
||||
@@ -418,61 +366,5 @@ class FDEM_DerivTests(unittest.TestCase):
|
||||
self.assertTrue(adjointTest('h', 'jzi'))
|
||||
|
||||
|
||||
if testCrossCheck:
|
||||
if testEB:
|
||||
def test_EB_CrossCheck_exr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'exr'))
|
||||
def test_EB_CrossCheck_eyr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'eyr'))
|
||||
def test_EB_CrossCheck_ezr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'ezr'))
|
||||
def test_EB_CrossCheck_exi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'exi'))
|
||||
def test_EB_CrossCheck_eyi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'eyi'))
|
||||
def test_EB_CrossCheck_ezi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'ezi'))
|
||||
|
||||
def test_EB_CrossCheck_bxr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bxr'))
|
||||
def test_EB_CrossCheck_byr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'byr'))
|
||||
def test_EB_CrossCheck_bzr_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bzr'))
|
||||
def test_EB_CrossCheck_bxi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bxi'))
|
||||
def test_EB_CrossCheck_byi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'byi'))
|
||||
def test_EB_CrossCheck_bzi_Eform(self):
|
||||
self.assertTrue(crossCheckTest('e', 'bzi'))
|
||||
|
||||
if testHJ:
|
||||
def test_HJ_CrossCheck_jxr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jxr'))
|
||||
def test_HJ_CrossCheck_jyr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jyr'))
|
||||
def test_HJ_CrossCheck_jzr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jzr'))
|
||||
def test_HJ_CrossCheck_jxi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jxi'))
|
||||
def test_HJ_CrossCheck_jyi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jyi'))
|
||||
def test_HJ_CrossCheck_jzi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'jzi'))
|
||||
|
||||
def test_HJ_CrossCheck_hxr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hxr'))
|
||||
def test_HJ_CrossCheck_hyr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hyr'))
|
||||
def test_HJ_CrossCheck_hzr_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hzr'))
|
||||
def test_HJ_CrossCheck_hxi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hxi'))
|
||||
def test_HJ_CrossCheck_hyi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hyi'))
|
||||
def test_HJ_CrossCheck_hzi_Jform(self):
|
||||
self.assertTrue(crossCheckTest('j', 'hzi'))
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -0,0 +1,11 @@
|
||||
if __name__ == '__main__':
|
||||
import os
|
||||
import glob
|
||||
import unittest
|
||||
test_file_strings = glob.glob('test_*.py')
|
||||
module_strings = [str[0:len(str)-3] for str in test_file_strings]
|
||||
suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
|
||||
in module_strings]
|
||||
testSuite = unittest.TestSuite(suites)
|
||||
|
||||
unittest.TextTestRunner(verbosity=2).run(testSuite)
|
||||
@@ -0,0 +1,45 @@
|
||||
import unittest, os
|
||||
import numpy as np
|
||||
from SimPEG import Utils
|
||||
|
||||
tol = 1e-15
|
||||
|
||||
class coorUtilsTest(unittest.TestCase):
|
||||
|
||||
def test_rotationMatrixFromNormals(self):
|
||||
v0 = np.random.rand(3)
|
||||
v0 *= 1./np.linalg.norm(v0)
|
||||
v1 = np.random.rand(3)
|
||||
v1 *= 1./np.linalg.norm(v1)
|
||||
Rf = Utils.coordutils.rotationMatrixFromNormals(v0,v1)
|
||||
Ri = Utils.coordutils.rotationMatrixFromNormals(v1,v0)
|
||||
|
||||
self.assertTrue(np.linalg.norm(Utils.mkvc(Rf.dot(v0) - v1)) < tol)
|
||||
self.assertTrue(np.linalg.norm(Utils.mkvc(Ri.dot(v1) - v0)) < tol)
|
||||
|
||||
def test_rotatePointsFromNormals(self):
|
||||
v0 = np.random.rand(3)
|
||||
v0*= 1./np.linalg.norm(v0)
|
||||
v1 = np.random.rand(3)
|
||||
v1*= 1./np.linalg.norm(v1)
|
||||
|
||||
v2 = Utils.mkvc(Utils.coordutils.rotatePointsFromNormals(Utils.mkvc(v0,2).T,v0,v1))
|
||||
|
||||
self.assertTrue(np.linalg.norm(v2-v1) < tol)
|
||||
|
||||
def test_rotateMatrixFromNormals(self):
|
||||
n0 = np.random.rand(3)
|
||||
n0*= 1./np.linalg.norm(n0)
|
||||
n1 = np.random.rand(3)
|
||||
n1*= 1./np.linalg.norm(n1)
|
||||
|
||||
scale = np.random.rand(100,1)
|
||||
XYZ0 = scale * n0
|
||||
XYZ1 = scale * n1
|
||||
|
||||
XYZ2 = Utils.coordutils.rotatePointsFromNormals(XYZ0,n0,n1)
|
||||
self.assertTrue(np.linalg.norm(Utils.mkvc(XYZ1) - Utils.mkvc(XYZ2))/np.linalg.norm(Utils.mkvc(XYZ1)) < tol)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
|
||||
Reference in New Issue
Block a user