From 7a5fde8a3eac998db6159c7c7587b93f16f218b2 Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Thu, 27 Feb 2014 09:44:23 -0800
Subject: [PATCH] Minor changes to problem

---
 docs/api_Richards.rst                               |  4 ++--
 .../examples/richards/comparisonToCeiliaEtAl1990.py | 12 ++++++------
 simpegFLOW/Richards/RichardsProblem.py              | 13 +++++++------
 3 files changed, 15 insertions(+), 14 deletions(-)

diff --git a/docs/api_Richards.rst b/docs/api_Richards.rst
index d1d51809..b0a69b3b 100644
--- a/docs/api_Richards.rst
+++ b/docs/api_Richards.rst
@@ -35,8 +35,8 @@ Here we reproduce the results from Ceilia et al. (1990):
 
 .. plot:: examples/richards/comparisonToCeiliaEtAl1990.py
 
-Richards
-========
+Richards API
+============
 
 .. automodule:: simpegFLOW.Richards.Empirical
     :show-inheritance:
diff --git a/docs/examples/richards/comparisonToCeiliaEtAl1990.py b/docs/examples/richards/comparisonToCeiliaEtAl1990.py
index 22abb2f2..1d3c8451 100644
--- a/docs/examples/richards/comparisonToCeiliaEtAl1990.py
+++ b/docs/examples/richards/comparisonToCeiliaEtAl1990.py
@@ -16,12 +16,12 @@ def getFields(timeStep,method):
                                     doNewton=False, method=method)
     return prob.fields(params['Ks'])
 
-Hs_M10 = getFields(10, 'mixed')
-Hs_M30 = getFields(30, 'mixed')
-Hs_M120= getFields(120,'mixed')
-Hs_H10 = getFields(10, 'head')
-Hs_H30 = getFields(30, 'head')
-Hs_H120= getFields(120,'head')
+Hs_M10 = getFields(10., 'mixed')
+Hs_M30 = getFields(30., 'mixed')
+Hs_M120= getFields(120.,'mixed')
+Hs_H10 = getFields(10., 'head')
+Hs_H30 = getFields(30., 'head')
+Hs_H120= getFields(120.,'head')
 
 plt.figure(figsize=(13,5))
 plt.subplot(121)
diff --git a/simpegFLOW/Richards/RichardsProblem.py b/simpegFLOW/Richards/RichardsProblem.py
index caf0ffc8..2a37308b 100644
--- a/simpegFLOW/Richards/RichardsProblem.py
+++ b/simpegFLOW/Richards/RichardsProblem.py
@@ -1,5 +1,5 @@
 from SimPEG import *
-from BaseRichards import RichardsModel
+from Empirical import RichardsModel
 
 class RichardsData(Data.BaseData):
     """docstring for RichardsData"""
@@ -66,7 +66,7 @@ class RichardsProblem(Problem.BaseProblem):
     boundaryConditions = None
     initialConditions  = None
 
-    dataPair = RichardsData
+    dataPair  = RichardsData
     modelPair = RichardsModel
 
     def __init__(self, mesh, model, **kwargs):
@@ -117,7 +117,7 @@ class RichardsProblem(Problem.BaseProblem):
         DIV  = self.mesh.faceDiv
         GRAD = self.mesh.cellGrad
         BC   = self.mesh.cellGradBC
-        AV   = self.mesh.aveCC2F
+        AV   = self.mesh.aveF2CC.T
         if self.mesh.dim == 1:
             Dz = self.mesh.faceDivx
         elif self.mesh.dim == 2:
@@ -171,7 +171,7 @@ class RichardsProblem(Problem.BaseProblem):
         DIV  = self.mesh.faceDiv
         GRAD = self.mesh.cellGrad
         BC   = self.mesh.cellGradBC
-        AV   = self.mesh.aveCC2F
+        AV   = self.mesh.aveF2CC.T
         if self.mesh.dim == 1:
             Dz = self.mesh.faceDivx
         elif self.mesh.dim == 2:
@@ -188,7 +188,7 @@ class RichardsProblem(Problem.BaseProblem):
         K  = self.model.k(h, m)
         dK = self.model.kDerivU(h, m)
 
-        aveK = 1./(AV*(1./K));
+        aveK = 1./(AV*(1./K))
 
         RHS = DIV*Utils.sdiag(aveK)*(GRAD*h+BC*bc) + Dz*aveK
         if self.method == 'mixed':
@@ -221,7 +221,8 @@ class RichardsProblem(Problem.BaseProblem):
         B = np.array(sp.vstack(Bs).todense())
 
         Ainv = Solver(A)
-        J = Ainv.solve(B)
+        P = self.data.projectFieldsDeriv(u, m)
+        J = P * Ainv.solve(B)
         return J
 
     def Jvec(self, m, v, u=None):