diff --git a/SimPEG/DCIP/BaseDC.py b/SimPEG/DCIP/BaseDC.py
index e3d353d7..475c621e 100644
--- a/SimPEG/DCIP/BaseDC.py
+++ b/SimPEG/DCIP/BaseDC.py
@@ -200,11 +200,11 @@ class ProblemDC_CC(Problem.BaseProblem):
 
         return F
 
-    def Jvec(self, m, v, u=None):
+    def Jvec(self, m, v, f=None):
         """
             :param numpy.array m: model
             :param numpy.array v: vector to multiply
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: numpy.array
             :return: Jv
 
@@ -225,11 +225,10 @@ class ProblemDC_CC(Problem.BaseProblem):
         # Set current model; clear dependent property $\mathbf{A(m)}$
         self.curModel = m
         sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
-        if u is None:
+        if f is None:
             # Run forward simulation if $u$ not provided
-            u = self.fields(self.curModel)[self.survey.srcList, 'phi_sol']
-        else:
-            u = u[self.survey.srcList, 'phi_sol']
+            f = self.fields(self.curModel)
+        u = f[self.survey.srcList, 'phi_sol']
 
         D = self.mesh.faceDiv
         G = self.mesh.cellGrad
@@ -251,19 +250,18 @@ class ProblemDC_CC(Problem.BaseProblem):
         if self.Ainv is None:
             self.Ainv = self.Solver(dA_du, **self.solverOpts)
 
-        P        = self.survey.getP(self.mesh)
+        P     = self.survey.getP(self.mesh)
         Jv    = - P * mkvc( self.Ainv * dCdm_x_v )
         return Jv
 
-    def Jtvec(self, m, v, u=None):
+    def Jtvec(self, m, v, f=None):
 
         self.curModel = m
         sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
-        if u is None:
-            # Run forward simulation if $u$ not provided
-            u = self.fields(self.curModel)[self.survey.srcList, 'phi_sol']
-        else:
-            u = u[self.survey.srcList, 'phi_sol']
+        if f is None:
+            # Run forward simulation if $f$ not provided
+            f = self.fields(self.curModel)
+        u = f[self.survey.srcList, 'phi_sol']
 
         shp = u.shape
         P = self.survey.getP(self.mesh)
diff --git a/SimPEG/DCIP/BaseIP.py b/SimPEG/DCIP/BaseIP.py
index cec0ea2e..a18b5a47 100644
--- a/SimPEG/DCIP/BaseIP.py
+++ b/SimPEG/DCIP/BaseIP.py
@@ -14,12 +14,12 @@ class SurveyIP(SurveyDC):
         Survey.BaseSurvey.__init__(self, **kwargs)
         self._Ps = {}
 
-    def dpred(self, m, u=None):
+    def dpred(self, m, f=None):
         """
             Predicted data.
 
             .. math::
-                d_\\text{pred} = Pu(m)
+                d_\\text{pred} = Pf(m)
         """
 
         return self.prob.forward(m)
@@ -143,10 +143,10 @@ class ProblemIP(Problem.BaseProblem):
         J_x_v    = - P * mkvc( self.Ainv * dCdm_x_v )
         return -J_x_v
 
-    def Jvec(self, m, v, u=None):
+    def Jvec(self, m, v, f=None):
         return self.forward(v)
 
-    def Jtvec(self, m, v, u=None):
+    def Jtvec(self, m, v, f=None):
 
         self.curModel = m
         # sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
diff --git a/SimPEG/DataMisfit.py b/SimPEG/DataMisfit.py
index 13574771..53728c4e 100644
--- a/SimPEG/DataMisfit.py
+++ b/SimPEG/DataMisfit.py
@@ -22,11 +22,11 @@ class BaseDataMisfit(object):
         Utils.setKwargs(self,**kwargs)
 
     @Utils.timeIt
-    def eval(self, m, u=None):
-        """eval(m, u=None)
+    def eval(self, m, f=None):
+        """eval(m, f=None)
 
             :param numpy.array m: geophysical model
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: float
             :return: data misfit
 
@@ -34,11 +34,11 @@ class BaseDataMisfit(object):
         raise NotImplementedError('This method should be overwritten.')
 
     @Utils.timeIt
-    def evalDeriv(self, m, u=None):
-        """evalDeriv(m, u=None)
+    def evalDeriv(self, m, f=None):
+        """evalDeriv(m, f=None)
 
             :param numpy.array m: geophysical model
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: numpy.array
             :return: data misfit derivative
 
@@ -47,12 +47,12 @@ class BaseDataMisfit(object):
 
 
     @Utils.timeIt
-    def eval2Deriv(self, m, v, u=None):
-        """eval2Deriv(m, v, u=None)
+    def eval2Deriv(self, m, v, f=None):
+        """eval2Deriv(m, v, f=None)
 
             :param numpy.array m: geophysical model
             :param numpy.array v: vector to multiply
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: numpy.array
             :return: data misfit derivative
 
@@ -108,24 +108,20 @@ class l2_DataMisfit(BaseDataMisfit):
         self._Wd = value
 
     @Utils.timeIt
-    def eval(self, m, u=None):
-        "eval(m, u=None)"
-        prob   = self.prob
-        survey = self.survey
-        R = self.Wd * survey.residual(m, u)
+    def eval(self, m, f=None):
+        "eval(m, f=None)"
+        if f is None: f = self.prob.fields(m)
+        R = self.Wd * self.survey.residual(m, f)
         return 0.5*np.vdot(R, R)
 
     @Utils.timeIt
-    def evalDeriv(self, m, u=None):
-        "evalDeriv(m, u=None)"
-        prob   = self.prob
-        survey = self.survey
-        if u is None: u = prob.fields(m)
-        return prob.Jtvec(m, self.Wd * (self.Wd * survey.residual(m, u)), u)
+    def evalDeriv(self, m, f=None):
+        "evalDeriv(m, f=None)"
+        if f is None: f = self.prob.fields(m)
+        return self.prob.Jtvec(m, self.Wd * (self.Wd * self.survey.residual(m, f=f)), f=f)
 
     @Utils.timeIt
-    def eval2Deriv(self, m, v, u=None):
-        "eval2Deriv(m, v, u=None)"
-        prob   = self.prob
-        if u is None: u = prob.fields(m)
-        return prob.Jtvec_approx(m, self.Wd * (self.Wd * prob.Jvec_approx(m, v, u=u)), u)
+    def eval2Deriv(self, m, v, f=None):
+        "eval2Deriv(m, v, f=None)"
+        if f is None: f = self.prob.fields(m)
+        return self.prob.Jtvec_approx(m, self.Wd * (self.Wd * self.prob.Jvec_approx(m, v, f=f)), f=f)
diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py
index 2ed27c20..3d726712 100644
--- a/SimPEG/Directives.py
+++ b/SimPEG/Directives.py
@@ -123,10 +123,10 @@ class BetaEstimate_ByEig(InversionDirective):
         if self.debug: print 'Calculating the beta0 parameter.'
 
         m = self.invProb.curModel
-        u = self.invProb.getFields(m, store=True, deleteWarmstart=False)
+        f = self.invProb.getFields(m, store=True, deleteWarmstart=False)
 
         x0 = np.random.rand(*m.shape)
-        t = x0.dot(self.dmisfit.eval2Deriv(m,x0,u=u))
+        t = x0.dot(self.dmisfit.eval2Deriv(m,x0,f=f))
         b = x0.dot(self.reg.eval2Deriv(m, v=x0))
         self.beta0 = self.beta0_ratio*(t/b)
 
diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py
index 163f6914..1552a12c 100644
--- a/SimPEG/EM/FDEM/SurveyFDEM.py
+++ b/SimPEG/EM/FDEM/SurveyFDEM.py
@@ -67,7 +67,7 @@ class Rx(SimPEG.Survey.BaseRx):
         """Grid Location projection (e.g. Ex Fy ...)"""
         return u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1]
 
-    def eval(self, src, mesh, u):
+    def eval(self, src, mesh, f):
         """
         Project fields to recievers to get data.
 
@@ -80,27 +80,27 @@ class Rx(SimPEG.Survey.BaseRx):
         # projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0])
         # projGLoc += self.knownRxTypes[self.rxType][1]
 
-        P = self.getP(mesh, self.projGLoc(u))
-        u_part_complex = u[src, self.projField]
+        P = self.getP(mesh, self.projGLoc(f))
+        f_part_complex = f[src, self.projField]
         # get the real or imag component
         real_or_imag = self.projComp
-        u_part = getattr(u_part_complex, real_or_imag)
+        f_part = getattr(f_part_complex, real_or_imag)
 
-        return P*u_part
+        return P*f_part
 
-    def evalDeriv(self, src, mesh, u, v, adjoint=False):
+    def evalDeriv(self, src, mesh, f, v, adjoint=False):
         """
         Derivative of projected fields with respect to the inversion model times a vector.
 
         :param Source src: FDEM source
         :param Mesh mesh: mesh used
-        :param Fields u: fields object
+        :param Fields f: fields object
         :param numpy.ndarray v: vector to multiply
         :rtype: numpy.ndarray
         :return: fields projected to recievers
         """
 
-        P = self.getP(mesh, self.projGLoc(u))
+        P = self.getP(mesh, self.projGLoc(f))
 
         if not adjoint:
             Pv_complex = P * v
diff --git a/SimPEG/EM/TDEM/BaseTDEM.py b/SimPEG/EM/TDEM/BaseTDEM.py
index 0da22072..15fc19e3 100644
--- a/SimPEG/EM/TDEM/BaseTDEM.py
+++ b/SimPEG/EM/TDEM/BaseTDEM.py
@@ -108,11 +108,11 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
         Ainv.clean()
         return F
 
-    def Jvec(self, m, v, u=None):
+    def Jvec(self, m, v, f=None):
         """
             :param numpy.array m: Conductivity model
             :param numpy.ndarray v: vector (model object)
-            :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
+            :param simpegEM.TDEM.FieldsTDEM f: Fields resulting from m
             :rtype: numpy.ndarray
             :return: w (data object)
 
@@ -125,15 +125,15 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
         """
         if self.verbose: print '%s\nCalculating J(v)\n%s'%('*'*50,'*'*50)
         self.curModel = m
-        if u is None:
-            u = self.fields(m)
-        p = self.Gvec(m, v, u)
+        if f is None:
+            f = self.fields(m)
+        p = self.Gvec(m, v, f)
         y = self.solveAh(m, p)
-        Jv = self.survey.evalDeriv(u, v=y)
+        Jv = self.survey.evalDeriv(f, v=y)
         if self.verbose: print '%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50)
         return - mkvc(Jv)
 
-    def Jtvec(self, m, v, u=None):
+    def Jtvec(self, m, v, f=None):
         """
             :param numpy.array m: Conductivity model
             :param numpy.ndarray,SimPEG.Survey.Data v: vector (data object)
@@ -150,15 +150,15 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
         """
         if self.verbose: print '%s\nCalculating J^T(v)\n%s'%('*'*50,'*'*50)
         self.curModel = m
-        if u is None:
-            u = self.fields(m)
+        if f is None:
+            f = self.fields(m)
 
         if not isinstance(v, self.dataPair):
             v = self.dataPair(self.survey, v)
 
-        p = self.survey.evalDeriv(u, v=v, adjoint=True)
+        p = self.survey.evalDeriv(f, v=v, adjoint=True)
         y = self.solveAht(m, p)
-        w = self.Gtvec(m, y, u)
+        w = self.Gtvec(m, y, f)
         if self.verbose: print '%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50)
         return - mkvc(w)
 
diff --git a/SimPEG/FLOW/Richards/RichardsProblem.py b/SimPEG/FLOW/Richards/RichardsProblem.py
index 4dcabe60..2346f4da 100644
--- a/SimPEG/FLOW/Richards/RichardsProblem.py
+++ b/SimPEG/FLOW/Richards/RichardsProblem.py
@@ -45,19 +45,19 @@ class RichardsSurvey(Survey.BaseSurvey):
 
     @Utils.count
     @Utils.requires('prob')
-    def dpred(self, m, u=None):
+    def dpred(self, m, f=None):
         """
             Create the projected data from a model.
-            The field, u, (if provided) will be used for the predicted data
+            The field, f, (if provided) will be used for the predicted data
             instead of recalculating the fields (which may be expensive!).
 
             .. math::
-                d_\\text{pred} = P(u(m), m)
+                d_\\text{pred} = P(f(m), m)
 
             Where P is a projection of the fields onto the data space.
         """
-        if u is None: u = self.prob.fields(m)
-        return Utils.mkvc(self.eval(u, m))
+        if f is None: f = self.prob.fields(m)
+        return Utils.mkvc(self.eval(f, m))
 
     @Utils.requires('prob')
     def eval(self, U, m):
@@ -233,16 +233,16 @@ class RichardsProblem(Problem.BaseTimeProblem):
         return r, J
 
     @Utils.timeIt
-    def Jfull(self, m, u=None):
-        if u is None:
-            u = self.fields(m)
+    def Jfull(self, m, f=None):
+        if f is None:
+            f = self.fields(m)
 
-        nn = len(u)-1
+        nn = len(f)-1
         Asubs, Adiags, Bs = range(nn), range(nn), range(nn)
         for ii in range(nn):
             dt = self.timeSteps[ii]
-            bc = self.getBoundaryConditions(ii, u[ii])
-            Asubs[ii], Adiags[ii], Bs[ii] = self.diagsJacobian(m, u[ii], u[ii+1], dt, bc)
+            bc = self.getBoundaryConditions(ii, f[ii])
+            Asubs[ii], Adiags[ii], Bs[ii] = self.diagsJacobian(m, f[ii], f[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])
@@ -251,7 +251,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
         B = np.array(sp.vstack(Bs).todense())
 
         Ainv = self.Solver(A, **self.solverOpts)
-        P = self.survey.evalDeriv(u, m)
+        P = self.survey.evalDeriv(f, m)
         AinvB = Ainv * B
         z = np.zeros((self.mesh.nC, B.shape[1]))
         zAinvB = np.vstack((z, AinvB))
@@ -259,41 +259,41 @@ class RichardsProblem(Problem.BaseTimeProblem):
         return J
 
     @Utils.timeIt
-    def Jvec(self, m, v, u=None):
-        if u is None:
-            u = self.fields(m)
+    def Jvec(self, m, v, f=None):
+        if f is None:
+            f = self.fields(m)
 
-        JvC = range(len(u)-1) # Cell to hold each row of the long vector.
+        JvC = range(len(f)-1) # Cell to hold each row of the long vector.
 
         # This is done via forward substitution.
-        bc = self.getBoundaryConditions(0, u[0])
-        temp, Adiag, B = self.diagsJacobian(m, u[0], u[1], self.timeSteps[0], bc)
+        bc = self.getBoundaryConditions(0, f[0])
+        temp, Adiag, B = self.diagsJacobian(m, f[0], f[1], self.timeSteps[0], bc)
         Adiaginv = self.Solver(Adiag, **self.solverOpts)
         JvC[0] = Adiaginv * (B*v)
 
-        for ii in range(1,len(u)-1):
-            bc = self.getBoundaryConditions(ii, u[ii])
-            Asub, Adiag, B = self.diagsJacobian(m, u[ii], u[ii+1], self.timeSteps[ii], bc)
+        for ii in range(1,len(f)-1):
+            bc = self.getBoundaryConditions(ii, f[ii])
+            Asub, Adiag, B = self.diagsJacobian(m, f[ii], f[ii+1], self.timeSteps[ii], bc)
             Adiaginv = self.Solver(Adiag, **self.solverOpts)
             JvC[ii] = Adiaginv * (B*v - Asub*JvC[ii-1])
 
-        P = self.survey.evalDeriv(u, m)
+        P = self.survey.evalDeriv(f, m)
         return P * np.concatenate([np.zeros(self.mesh.nC)] + JvC)
 
     @Utils.timeIt
-    def Jtvec(self, m, v, u=None):
-        if u is None:
-            u = self.field(m)
+    def Jtvec(self, m, v, f=None):
+        if f is None:
+            f = self.field(m)
 
-        P = self.survey.evalDeriv(u, m)
+        P = self.survey.evalDeriv(f, m)
         PTv = P.T*v
 
         # This is done via backward substitution.
         minus = 0
         BJtv = 0
-        for ii in range(len(u)-1,0,-1):
-            bc = self.getBoundaryConditions(ii-1, u[ii-1])
-            Asub, Adiag, B = self.diagsJacobian(m, u[ii-1], u[ii], self.timeSteps[ii-1], bc)
+        for ii in range(len(f)-1,0,-1):
+            bc = self.getBoundaryConditions(ii-1, f[ii-1])
+            Asub, Adiag, B = self.diagsJacobian(m, f[ii-1], f[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)
diff --git a/SimPEG/InvProblem.py b/SimPEG/InvProblem.py
index 0296bf4b..fd6c48c3 100644
--- a/SimPEG/InvProblem.py
+++ b/SimPEG/InvProblem.py
@@ -82,23 +82,23 @@ class BaseInvProblem(object):
         self._warmstart = value
 
     def getFields(self, m, store=False, deleteWarmstart=True):
-        u = None
+        f = None
 
         for mtest, u_ofmtest in self.warmstart:
             if m is mtest:
-                u = u_ofmtest
+                f = u_ofmtest
                 if self.debug: print 'InvProb is Warm Starting!'
                 break
 
-        if u is None:
-            u = self.prob.fields(m)
+        if f is None:
+            f = self.prob.fields(m)
 
         if deleteWarmstart:
             self.warmstart = []
         if store:
-            self.warmstart += [(m,u)]
+            self.warmstart += [(m,f)]
 
-        return u
+        return f
 
     @Utils.timeIt
     def evalFunction(self, m, return_g=True, return_H=True):
@@ -109,21 +109,21 @@ class BaseInvProblem(object):
         gc.collect()
 
         # Store fields if doing a line-search
-        u = self.getFields(m, store=(return_g==False and return_H==False))
+        f = self.getFields(m, store=(return_g==False and return_H==False))
 
-        phi_d = self.dmisfit.eval(m, u=u)
+        phi_d = self.dmisfit.eval(m, f=f)
         phi_m = self.reg.eval(m)
 
-        self.dpred = self.survey.dpred(m, u=u)  # This is a cheap matrix vector calculation.
+        self.dpred = self.survey.dpred(m, f=f)  # This is a cheap matrix vector calculation.
 
         self.phi_d, self.phi_d_last  = phi_d, self.phi_d
         self.phi_m, self.phi_m_last  = phi_m, self.phi_m
 
-        f = phi_d + self.beta * phi_m
+        phi = phi_d + self.beta * phi_m
 
-        out = (f,)
+        out = (phi,)
         if return_g:
-            phi_dDeriv = self.dmisfit.evalDeriv(m, u=u)
+            phi_dDeriv = self.dmisfit.evalDeriv(m, f=f)
             phi_mDeriv = self.reg.evalDeriv(m)
 
             g = phi_dDeriv + self.beta * phi_mDeriv
@@ -131,7 +131,7 @@ class BaseInvProblem(object):
 
         if return_H:
             def H_fun(v):
-                phi_d2Deriv = self.dmisfit.eval2Deriv(m, v, u=u)
+                phi_d2Deriv = self.dmisfit.eval2Deriv(m, v, f=f)
                 phi_m2Deriv = self.reg.eval2Deriv(m, v=v)
 
                 return phi_d2Deriv + self.beta * phi_m2Deriv
diff --git a/SimPEG/MT/BaseMT.py b/SimPEG/MT/BaseMT.py
index 36389430..c201dfb0 100644
--- a/SimPEG/MT/BaseMT.py
+++ b/SimPEG/MT/BaseMT.py
@@ -27,7 +27,7 @@ class BaseMTProblem(BaseFDEMProblem):
     # Might need to add more stuff here.
 
     ## NEED to clean up the Jvec and Jtvec to use Zero and Identities for None components.
-    def Jvec(self, m, v, u=None):
+    def Jvec(self, m, v, f=None):
         """
         Function to calculate the data sensitivities dD/dm times a vector.
 
@@ -39,8 +39,8 @@ class BaseMTProblem(BaseFDEMProblem):
         """
 
         # Calculate the fields
-        if u is None:
-           u = self.fields(m)
+        if f is None:
+           f= self.fields(m)
         # Set current model
         self.curModel = m
         # Initiate the Jv object
@@ -56,9 +56,9 @@ class BaseMTProblem(BaseFDEMProblem):
                 # We need fDeriv_m = df/du*du/dm + df/dm
                 # Construct du/dm, it requires a solve
                 # NOTE: need to account for the 2 polarizations in the derivatives.
-                u_src = u[src,:]
+                f_src = f[src,:]
                 # dA_dm and dRHS_dm should be of size nE,2, so that we can multiply by dA_duI. The 2 columns are each of the polarizations.
-                dA_dm = self.getADeriv_m(freq, u_src, v) # Size: nE,2 (u_px,u_py) in the columns.
+                dA_dm = self.getADeriv_m(freq, f_src, v) # Size: nE,2 (u_px,u_py) in the columns.
                 dRHS_dm = self.getRHSDeriv_m(freq, v) # Size: nE,2 (u_px,u_py) in the columns.
                 if dRHS_dm is None:
                     du_dm = dA_duI * ( -dA_dm )
@@ -68,13 +68,13 @@ class BaseMTProblem(BaseFDEMProblem):
                 for rx in src.rxList:
                     # Get the projection derivative
                     # v should be of size 2*nE (for 2 polarizations)
-                    PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, u, t) # wrt u, we don't have have PDeriv wrt m
+                    PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, f, t) # wrt u, we don't have have PDeriv wrt m
                     Jv[src, rx] = PDeriv_u(mkvc(du_dm))
             dA_duI.clean()
         # Return the vectorized sensitivities
         return mkvc(Jv)
 
-    def Jtvec(self, m, v, u=None):
+    def Jtvec(self, m, v, f=None):
         """
         Function to calculate the transpose of the data sensitivities (dD/dm)^T times a vector.
 
@@ -85,8 +85,8 @@ class BaseMTProblem(BaseFDEMProblem):
             :return: Data sensitivities wrt m
         """
 
-        if u is None:
-            u = self.fields(m)
+        if f is None:
+            f = self.fields(m)
 
         self.curModel = m
 
@@ -103,15 +103,15 @@ class BaseMTProblem(BaseFDEMProblem):
 
             for src in self.survey.getSrcByFreq(freq):
                 ftype = self._fieldType + 'Solution'
-                u_src = u[src, :]
+                f_src = f[src, :]
 
                 for rx in src.rxList:
                     # Get the adjoint evalDeriv
                     # PTv needs to be nE,
-                    PTv = rx.evalDeriv(src, self.mesh, u, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
+                    PTv = rx.evalDeriv(src, self.mesh, f, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
                     # Get the
                     dA_duIT = ATinv * PTv
-                    dA_dmT = self.getADeriv_m(freq, u_src, mkvc(dA_duIT), adjoint=True)
+                    dA_dmT = self.getADeriv_m(freq, f_src, mkvc(dA_duIT), adjoint=True)
                     dRHS_dmT = self.getRHSDeriv_m(freq, mkvc(dA_duIT), adjoint=True)
                     # Make du_dmT
                     if dRHS_dmT is None:
@@ -129,4 +129,4 @@ class BaseMTProblem(BaseFDEMProblem):
                         raise Exception('Must be real or imag')
             # Clean the factorization, clear memory.
             ATinv.clean()
-        return Jtv
\ No newline at end of file
+        return Jtv
diff --git a/SimPEG/MT/SurveyMT.py b/SimPEG/MT/SurveyMT.py
index 4e4a8688..0ec91a0e 100644
--- a/SimPEG/MT/SurveyMT.py
+++ b/SimPEG/MT/SurveyMT.py
@@ -427,15 +427,15 @@ class Survey(SimPEGsurvey.BaseSurvey):
         assert freq in self._freqDict, "The requested frequency is not in this survey."
         return self._freqDict[freq]
 
-    def eval(self, u):
+    def eval(self, f):
         data = Data(self)
         for src in self.srcList:
             sys.stdout.flush()
             for rx in src.rxList:
-                data[src, rx] = rx.eval(src, self.mesh, u)
+                data[src, rx] = rx.eval(src, self.mesh, f)
         return data
 
-    def evalDeriv(self, u):
+    def evalDeriv(self, f):
         raise Exception('Use Transmitters to project fields deriv.')
 
 #################
diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py
index cd8a4aaa..f8520c5c 100644
--- a/SimPEG/Problem.py
+++ b/SimPEG/Problem.py
@@ -88,28 +88,28 @@ class BaseProblem(object):
         return self.survey is not None
 
     @Utils.timeIt
-    def Jvec(self, m, v, u=None):
-        """Jvec(m, v, u=None)
+    def Jvec(self, m, v, f=None):
+        """Jvec(m, v, f=None)
 
             Effect of J(m) on a vector v.
 
             :param numpy.array m: model
             :param numpy.array v: vector to multiply
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: numpy.array
             :return: Jv
         """
         raise NotImplementedError('J is not yet implemented.')
 
     @Utils.timeIt
-    def Jtvec(self, m, v, u=None):
-        """Jtvec(m, v, u=None)
+    def Jtvec(self, m, v, f=None):
+        """Jtvec(m, v, f=None)
 
             Effect of transpose of J(m) on a vector v.
 
             :param numpy.array m: model
             :param numpy.array v: vector to multiply
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: numpy.array
             :return: JTv
         """
@@ -117,32 +117,32 @@ class BaseProblem(object):
 
 
     @Utils.timeIt
-    def Jvec_approx(self, m, v, u=None):
-        """Jvec_approx(m, v, u=None)
+    def Jvec_approx(self, m, v, f=None):
+        """Jvec_approx(m, v, f=None)
 
             Approximate effect of J(m) on a vector v
 
             :param numpy.array m: model
             :param numpy.array v: vector to multiply
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: numpy.array
             :return: approxJv
         """
-        return self.Jvec(m, v, u)
+        return self.Jvec(m, v, f)
 
     @Utils.timeIt
-    def Jtvec_approx(self, m, v, u=None):
-        """Jtvec_approx(m, v, u=None)
+    def Jtvec_approx(self, m, v, f=None):
+        """Jtvec_approx(m, v, f=None)
 
             Approximate effect of transpose of J(m) on a vector v.
 
             :param numpy.array m: model
             :param numpy.array v: vector to multiply
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: numpy.array
             :return: JTv
         """
-        return self.Jtvec(m, v, u)
+        return self.Jtvec(m, v, f)
 
     def fields(self, m):
         """
@@ -224,9 +224,9 @@ class LinearProblem(BaseProblem):
     def fields(self, m):
         return self.G.dot(m)
 
-    def Jvec(self, m, v, u=None):
+    def Jvec(self, m, v, f=None):
         return self.G.dot(v)
 
-    def Jtvec(self, m, v, u=None):
+    def Jtvec(self, m, v, f=None):
         return self.G.T.dot(v)
 
diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py
index 37024028..fbc88276 100644
--- a/SimPEG/Survey.py
+++ b/SimPEG/Survey.py
@@ -295,38 +295,38 @@ class BaseSurvey(object):
 
     @Utils.count
     @Utils.requires('prob')
-    def dpred(self, m, u=None):
-        """dpred(m, u=None)
+    def dpred(self, m, f=None):
+        """dpred(m, f=None)
 
             Create the projected data from a model.
-            The field, u, (if provided) will be used for the predicted data
+            The fields, f, (if provided) will be used for the predicted data
             instead of recalculating the fields (which may be expensive!).
 
             .. math::
 
-                d_\\text{pred} = P(u(m))
+                d_\\text{pred} = P(f(m))
 
             Where P is a projection of the fields onto the data space.
         """
-        if u is None: u = self.prob.fields(m)
-        return Utils.mkvc(self.eval(u))
+        if f is None: f = self.prob.fields(m)
+        return Utils.mkvc(self.eval(f))
 
 
     @Utils.count
-    def eval(self, u):
-        """eval(u)
+    def eval(self, f):
+        """eval(f)
 
             This function projects the fields onto the data space.
 
             .. math::
 
-                d_\\text{pred} = \mathbf{P} u(m)
+                d_\\text{pred} = \mathbf{P} f(m)
         """
         raise NotImplemented('eval is not yet implemented.')
 
     @Utils.count
-    def evalDeriv(self, u):
-        """evalDeriv(u)
+    def evalDeriv(self, f):
+        """evalDeriv(f)
 
             This function s the derivative of projects the fields onto the data space.
 
@@ -337,11 +337,11 @@ class BaseSurvey(object):
         raise NotImplemented('eval is not yet implemented.')
 
     @Utils.count
-    def residual(self, m, u=None):
-        """residual(m, u=None)
+    def residual(self, m, f=None):
+        """residual(m, f=None)
 
             :param numpy.array m: geophysical model
-            :param numpy.array u: fields
+            :param numpy.array f: fields
             :rtype: numpy.array
             :return: data residual
 
@@ -352,14 +352,14 @@ class BaseSurvey(object):
                 \mu_\\text{data} = \mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}
 
         """
-        return Utils.mkvc(self.dpred(m, u=u) - self.dobs)
+        return Utils.mkvc(self.dpred(m, f=f) - self.dobs)
 
     @property
     def isSynthetic(self):
         "Check if the data is synthetic."
         return self.mtrue is not None
 
-    def makeSyntheticData(self, m, std=0.05, u=None, force=False):
+    def makeSyntheticData(self, m, std=0.05, f=None, force=False):
         """
             Make synthetic data given a model, and a standard deviation.
 
@@ -372,16 +372,16 @@ class BaseSurvey(object):
         if getattr(self, 'dobs', None) is not None and not force:
             raise Exception('Survey already has dobs. You can use force=True to override this exception.')
         self.mtrue = m
-        self.dtrue = self.dpred(m, u=u)
+        self.dtrue = self.dpred(m, f=f)
         noise = std*abs(self.dtrue)*np.random.randn(*self.dtrue.shape)
         self.dobs = self.dtrue+noise
         self.std = self.dobs*0 + std
         return self.dobs
 
 class LinearSurvey(BaseSurvey):
-    def eval(self, u):
-        return u
-    
+    def eval(self, f):
+        return f
+
     @property
     def nD(self):
         return self.prob.G.shape[0]
diff --git a/tests/flow/test_Richards.py b/tests/flow/test_Richards.py
index d63a6210..f67ec71d 100644
--- a/tests/flow/test_Richards.py
+++ b/tests/flow/test_Richards.py
@@ -116,8 +116,8 @@ class RichardsTests1D(unittest.TestCase):
         v = np.random.rand(self.survey.nD)
         z = np.random.rand(self.M.nC)
         Hs = self.prob.fields(self.Ks)
-        vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs))
-        zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs))
+        vJz = v.dot(self.prob.Jvec(self.Ks,z,f=Hs))
+        zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=Hs))
         tol = TOL*(10**int(np.log10(np.abs(zJv))))
         passed = np.abs(vJz - zJv) < tol
         print 'Richards Adjoint Test - PressureHead'
@@ -188,8 +188,8 @@ class RichardsTests2D(unittest.TestCase):
         v = np.random.rand(self.survey.nD)
         z = np.random.rand(self.M.nC)
         Hs = self.prob.fields(self.Ks)
-        vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs))
-        zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs))
+        vJz = v.dot(self.prob.Jvec(self.Ks,z,f=Hs))
+        zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=Hs))
         tol = TOL*(10**int(np.log10(np.abs(zJv))))
         passed = np.abs(vJz - zJv) < tol
         print '2D: Richards Adjoint Test - PressureHead'
@@ -260,8 +260,8 @@ class RichardsTests3D(unittest.TestCase):
         v = np.random.rand(self.survey.nD)
         z = np.random.rand(self.M.nC)
         Hs = self.prob.fields(self.Ks)
-        vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs))
-        zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs))
+        vJz = v.dot(self.prob.Jvec(self.Ks,z,f=Hs))
+        zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=Hs))
         tol = TOL*(10**int(np.log10(np.abs(zJv))))
         passed = np.abs(vJz - zJv) < tol
         print '3D: Richards Adjoint Test - PressureHead'