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

simpeg updates

Browse Source
This commit is contained in:
rowanc1
2014-05-18 23:18:24 -07:00
parent 96e129aae2
commit 26079333bb
3 changed files with 121 additions and 40 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegFLOW/Richards/Empirical.py
+73 -10
View File
@@ -1,6 +1,69 @@
from SimPEG import Mesh, Maps, Utils, np
class NonLinearMap(object):
"""
SimPEG NonLinearMap
"""
__metaclass__ = Utils.SimPEGMetaClass
counter = None #: A SimPEG.Utils.Counter object
mesh = None #: A SimPEG Mesh
def __init__(self, mesh):
self.mesh = mesh
def _transform(self, u, m):
"""
:param numpy.array u: fields
:param numpy.array m: model
:rtype: numpy.array
:return: transformed model
The *transform* changes the model into the physical property.
"""
return m
def derivU(self, u, m):
"""
:param numpy.array u: fields
:param numpy.array m: model
:rtype: scipy.csr_matrix
:return: derivative of transformed model
The *transform* changes the model into the physical property.
The *transformDerivU* provides the derivative of the *transform* with respect to the fields.
"""
raise NotImplementedError('The transformDerivU is not implemented.')
def derivM(self, u, m):
"""
:param numpy.array u: fields
:param numpy.array m: model
:rtype: scipy.csr_matrix
:return: derivative of transformed model
The *transform* changes the model into the physical property.
The *transformDerivU* provides the derivative of the *transform* with respect to the model.
"""
raise NotImplementedError('The transformDerivM is not implemented.')
@property
def nP(self):
"""Number of parameters in the model."""
return self.mesh.nC
def example(self):
raise NotImplementedError('The example is not implemented.')
def test(self, m=None):
raise NotImplementedError('The test is not implemented.')
class RichardsMap(object):
"""docstring for RichardsMap"""
@@ -18,8 +81,8 @@ class RichardsMap(object):
def __init__(self, mesh, thetaModel, kModel):
self.mesh = mesh
assert isinstance(thetaModel, Maps.NonLinearMap)
assert isinstance(kModel, Maps.NonLinearMap)
assert isinstance(thetaModel, NonLinearMap)
assert isinstance(kModel, NonLinearMap)
self._thetaModel = thetaModel
self._kModel = kModel
@@ -94,7 +157,7 @@ class HaverkampParams(object):
'gamma':4.74}
class _haverkamp_theta(Maps.NonLinearMap):
class _haverkamp_theta(NonLinearMap):
theta_s = 0.430
theta_r = 0.078
@@ -102,7 +165,7 @@ class _haverkamp_theta(Maps.NonLinearMap):
beta = 3.960
def __init__(self, mesh, **kwargs):
Maps.NonLinearMap.__init__(self, mesh)
NonLinearMap.__init__(self, mesh)
Utils.setKwargs(self, **kwargs)
def setModel(self, m):
@@ -131,14 +194,14 @@ class _haverkamp_theta(Maps.NonLinearMap):
return g
class _haverkamp_k(Maps.NonLinearMap):
class _haverkamp_k(NonLinearMap):
A = 1.175e+06
gamma = 4.74
Ks = np.log(24.96)
def __init__(self, mesh, **kwargs):
Maps.NonLinearMap.__init__(self, mesh)
NonLinearMap.__init__(self, mesh)
Utils.setKwargs(self, **kwargs)
def setModel(self, m):
@@ -193,7 +256,7 @@ class Haverkamp(RichardsMap):
class _vangenuchten_theta(Maps.NonLinearMap):
class _vangenuchten_theta(NonLinearMap):
theta_s = 0.430
theta_r = 0.078
@@ -201,7 +264,7 @@ class _vangenuchten_theta(Maps.NonLinearMap):
n = 1.560
def __init__(self, mesh, **kwargs):
Maps.NonLinearMap.__init__(self, mesh)
NonLinearMap.__init__(self, mesh)
Utils.setKwargs(self, **kwargs)
def setModel(self, m):
@@ -229,7 +292,7 @@ class _vangenuchten_theta(Maps.NonLinearMap):
return g
class _vangenuchten_k(Maps.NonLinearMap):
class _vangenuchten_k(NonLinearMap):
I = 0.500
alpha = 0.036
@@ -237,7 +300,7 @@ class _vangenuchten_k(Maps.NonLinearMap):
Ks = np.log(24.96)
def __init__(self, mesh, **kwargs):
Maps.NonLinearMap.__init__(self, mesh)
NonLinearMap.__init__(self, mesh)
Utils.setKwargs(self, **kwargs)
def setModel(self, m):
simpegFLOW/Richards/RichardsProblem.py
+37 -28
View File
@@ -7,14 +7,16 @@ class RichardsRx(Survey.BaseTimeRx):
knownRxTypes = ['saturation','pressureHead']
def projectFields(self, u, m, mapping, mesh, timeMesh):
def projectFields(self, U, m, mapping, mesh, timeMesh):
if self.rxType == 'saturation':
u = mapping.theta(u, m)
if self.rxType == 'pressureHead':
u = np.concatenate(U)
elif self.rxType == 'saturation':
u = np.concatenate([mapping.theta(ui, m) for ui in U])
return self.getP(mesh, timeMesh) * u
def projectFieldsDeriv(self, u, m, mapping, mesh, timeMesh):
def projectFieldsDeriv(self, U, m, mapping, mesh, timeMesh):
P = self.getP(mesh, timeMesh)
if self.rxType == 'pressureHead':
@@ -23,7 +25,7 @@ class RichardsRx(Survey.BaseTimeRx):
#TODO: if m is a parameter in the theta
# distribution, we may need to do
# some more chain rule here.
dT = mapping.thetaDerivU(u, m)
dT = sp.block_diag([mapping.thetaDerivU(ui, m) for ui in U])
return P*dT
@@ -58,12 +60,9 @@ class RichardsSurvey(Survey.BaseSurvey):
@Utils.requires('prob')
def projectFields(self, U, m):
u = np.concatenate(U)
Ds = range(len(self.rxList))
for ii, rx in enumerate(self.rxList):
Ds[ii] = rx.projectFields(u, m,
Ds[ii] = rx.projectFields(U, m,
self.prob.mapping,
self.prob.mesh,
self.prob.timeMesh)
@@ -73,12 +72,9 @@ class RichardsSurvey(Survey.BaseSurvey):
@Utils.requires('prob')
def projectFieldsDeriv(self, U, m):
"""The Derivative with respect to the fields."""
u = np.concatenate(U)
Ds = range(len(self.rxList))
for ii, rx in enumerate(self.rxList):
Ds[ii] = rx.projectFieldsDeriv(u, m,
Ds[ii] = rx.projectFieldsDeriv(U, m,
self.prob.mapping,
self.prob.mesh,
self.prob.timeMesh)
@@ -100,6 +96,14 @@ class RichardsProblem(Problem.BaseTimeProblem):
def __init__(self, mesh, mapping=None, **kwargs):
Problem.BaseTimeProblem.__init__(self, mesh, mapping=mapping, **kwargs)
def getBoundaryConditions(self, ii):
if type(self.boundaryConditions) is np.ndarray:
return self.boundaryConditions
time = self.timeMesh.vectorCCx[ii]
return self.boundaryConditions(time)
@property
def method(self):
"""Method must be either 'mixed' or 'head'. See notes in Celia et al., 1990."""
@@ -127,10 +131,11 @@ class RichardsProblem(Problem.BaseTimeProblem):
u = range(self.nT+1)
u[0] = self.initialConditions
for ii, dt in enumerate(self.timeSteps):
u[ii+1] = self.rootFinder.root(lambda hn1m, return_g=True: self.getResidual(m, u[ii], hn1m, dt, return_g=return_g), u[ii])
bc = self.getBoundaryConditions(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])
return u
def diagsJacobian(self, m, hn, hn1, dt):
def diagsJacobian(self, m, hn, hn1, dt, bc):
DIV = self.mesh.faceDiv
GRAD = self.mesh.cellGrad
@@ -143,8 +148,6 @@ class RichardsProblem(Problem.BaseTimeProblem):
elif self.mesh.dim == 3:
Dz = sp.hstack((Utils.spzeros(self.mesh.nC,self.mesh.vnF[0]+self.mesh.vnF[1]), self.mesh.faceDivz),format='csr')
bc = self.boundaryConditions
dT = self.mapping.thetaDerivU(hn, m)
dT1 = self.mapping.thetaDerivU(hn1, m)
K1 = self.mapping.k(hn1, m)
@@ -181,7 +184,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
return Asub, Adiag, B
def getResidual(self, m, hn, h, dt, return_g=True):
def getResidual(self, m, hn, h, dt, bc, return_g=True):
"""
Where h is the proposed value for the next time iterate (h_{n+1})
"""
@@ -196,8 +199,6 @@ class RichardsProblem(Problem.BaseTimeProblem):
elif self.mesh.dim == 3:
Dz = sp.hstack((Utils.spzeros(self.mesh.nC,self.mesh.vnF[0]+self.mesh.vnF[1]), self.mesh.faceDivz),format='csr')
bc = self.boundaryConditions
T = self.mapping.theta(h, m)
dT = self.mapping.thetaDerivU(h, m)
Tn = self.mapping.theta(hn, m)
@@ -228,7 +229,9 @@ class RichardsProblem(Problem.BaseTimeProblem):
nn = len(u)-1
Asubs, Adiags, Bs = range(nn), range(nn), range(nn)
for ii in range(nn):
Asubs[ii], Adiags[ii], Bs[ii] = self.diagsJacobian(m, u[ii],u[ii+1])
dt = self.timeSteps[ii]
bc = self.getBoundaryConditions(ii)
Asubs[ii], Adiags[ii], Bs[ii] = self.diagsJacobian(m, u[ii], u[ii+1], dt, bc)
Ad = sp.block_diag(Adiags)
zRight = Utils.spzeros((len(Asubs)-1)*Asubs[0].shape[0],Adiags[0].shape[1])
zTop = Utils.spzeros(Adiags[0].shape[0], len(Adiags)*Adiags[0].shape[1])
@@ -238,7 +241,10 @@ class RichardsProblem(Problem.BaseTimeProblem):
Ainv = self.Solver(A, **self.solverOpts)
P = self.survey.projectFieldsDeriv(u, m)
J = P * Ainv.solve(B)
AinvB = Ainv * B
z = np.zeros((self.mesh.nC, B.shape[1]))
zAinvB = np.vstack((z, AinvB))
J = P * zAinvB
return J
def Jvec(self, m, v, u=None):
@@ -248,14 +254,16 @@ class RichardsProblem(Problem.BaseTimeProblem):
JvC = range(len(u)-1) # Cell to hold each row of the long vector.
# This is done via forward substitution.
temp, Adiag, B = self.diagsJacobian(m, u[0], u[1], self.timeSteps[0])
bc = self.getBoundaryConditions(0)
temp, Adiag, B = self.diagsJacobian(m, u[0], u[1], self.timeSteps[0], bc)
Adiaginv = self.Solver(Adiag, **self.solverOpts)
JvC[0] = Adiaginv.solve(B*v)
JvC[0] = Adiaginv * (B*v)
for ii in range(1,len(u)-1):
Asub, Adiag, B = self.diagsJacobian(m, u[ii], u[ii+1], self.timeSteps[ii])
bc = self.getBoundaryConditions(ii)
Asub, Adiag, B = self.diagsJacobian(m, u[ii], u[ii+1], self.timeSteps[ii], bc)
Adiaginv = self.Solver(Adiag, **self.solverOpts)
JvC[ii] = Adiaginv.solve(B*v - Asub*JvC[ii-1])
JvC[ii] = Adiaginv * (B*v - Asub*JvC[ii-1])
P = self.survey.projectFieldsDeriv(u, m)
return P * np.concatenate([np.zeros(self.mesh.nC)] + JvC)
@@ -271,11 +279,12 @@ class RichardsProblem(Problem.BaseTimeProblem):
minus = 0
BJtv = 0
for ii in range(len(u)-1,0,-1):
Asub, Adiag, B = self.diagsJacobian(m, u[ii-1], u[ii], self.timeSteps[ii-1])
bc = self.getBoundaryConditions(ii-1)
Asub, Adiag, B = self.diagsJacobian(m, u[ii-1], u[ii], self.timeSteps[ii-1], bc)
#select the correct part of v
vpart = range((ii)*Adiag.shape[0], (ii+1)*Adiag.shape[0])
AdiaginvT = self.Solver(Adiag.T, **self.solverOpts)
JTvC = AdiaginvT.solve(PTv[vpart] - minus)
JTvC = AdiaginvT * (PTv[vpart] - minus)
minus = Asub.T*JTvC # this is now the super diagonal.
BJtv = BJtv + B.T*JTvC
simpegFLOW/Tests/test_Richards.py
+11 -2
View File
@@ -93,13 +93,13 @@ class RichardsTests1D(unittest.TestCase):
def test_Richards_getResidual_Newton(self):
self.prob.doNewton = True
m = self.Ks
passed = checkDerivative(lambda hn1: self.prob.getResidual(m, self.h0, hn1, self.prob.timeSteps[0]), self.h0, plotIt=False)
passed = checkDerivative(lambda hn1: self.prob.getResidual(m, self.h0, hn1, self.prob.timeSteps[0], self.prob.boundaryConditions), self.h0, plotIt=False)
self.assertTrue(passed,True)
def test_Richards_getResidual_Picard(self):
self.prob.doNewton = False
m = self.Ks
passed = checkDerivative(lambda hn1: self.prob.getResidual(m, self.h0, hn1, self.prob.timeSteps[0]), self.h0, plotIt=False, expectedOrder=1)
passed = checkDerivative(lambda hn1: self.prob.getResidual(m, self.h0, hn1, self.prob.timeSteps[0], self.prob.boundaryConditions), self.h0, plotIt=False, expectedOrder=1)
self.assertTrue(passed,True)
def test_Adjoint(self):
@@ -122,6 +122,15 @@ class RichardsTests1D(unittest.TestCase):
self.assertTrue(passed,True)
def test_Sensitivity_full(self):
mTrue = self.Ks*np.ones(self.M.nC)
J = self.prob.Jfull(mTrue)
derChk = lambda m: [self.survey.dpred(m), J]
print 'Testing Richards Derivative'
passed = checkDerivative(derChk, mTrue, num=4, plotIt=False)
self.assertTrue(passed,True)
# class RichardsTests2D(object):
# def setUp(self):