diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py
index 4cb63c41..8b876f45 100644
--- a/simpegEM/TDEM/BaseTDEM.py
+++ b/simpegEM/TDEM/BaseTDEM.py
@@ -15,21 +15,21 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
         BaseTimeProblem.__init__(self, mesh, mapping=mapping, **kwargs)
 
     surveyPair = SurveyTDEM
-    _FieldsTDEM_pair = FieldsTDEM  #: used for the forward calculation only
+    _FieldsForward_pair = FieldsTDEM  #: used for the forward calculation only
 
     def fields(self, m):
         if self.verbose: print '%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50)
         self.curModel = m
         # Create a fields storage object
-        F = self._FieldsTDEM_pair(self.mesh, self.survey)
+        F = self._FieldsForward_pair(self.mesh, self.survey)
         for tx in self.survey.txList:
             # Set the initial conditions
             F[tx,:,0] = tx.getInitialFields(self.mesh)
-        F = self.forward(m, self.getRHS, self.calcFields, F=F)
+        F = self.forward(m, self.getRHS, F=F)
         if self.verbose: print '%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50)
         return F
 
-    def forward(self, m, RHS, CalcFields, F=None):
+    def forward(self, m, RHS, F=None):
         self.curModel = m
         F = F or FieldsTDEM(self.mesh, self.survey)
 
@@ -50,11 +50,11 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
             if self.verbose: print '    Done...'
             if sol.ndim == 1:
                 sol.shape = (sol.size,1)
-            F[:,:,tInd+1] = CalcFields(sol, tInd)
+            F[:,self.solType,tInd+1] = sol
         Ainv.clean()
         return F
 
-    def adjoint(self, m, RHS, CalcFields, F=None):
+    def adjoint(self, m, RHS, F=None):
         self.curModel = m
         F = F or FieldsTDEM(self.mesh, self.survey)
 
@@ -75,7 +75,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
             if self.verbose: print '    Done...'
             if sol.ndim == 1:
                 sol.shape = (sol.size,1)
-            F[:,:,tInd+1] = CalcFields(sol, tInd)
+            F[:,self.solType,tInd+1] = sol
         Ainv.clean()
         return F
 
diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py
index bd553df3..4ceed5d9 100644
--- a/simpegEM/TDEM/TDEM_b.py
+++ b/simpegEM/TDEM/TDEM_b.py
@@ -14,10 +14,31 @@ class FieldsTDEM_e_from_b(FieldsTDEM):
         self.edgeCurlT = self.survey.prob.mesh.edgeCurl.T
         self.MfMui     = self.survey.prob.MfMui
 
-    def e_from_b(self, b, ind):
+    def e_from_b(self, b, txInd, timeInd):
         # TODO: implement non-zero js
         return self.MeSigmaI*(self.edgeCurlT*(self.MfMui*b))
 
+class FieldsTDEM_e_from_b_Ah(FieldsTDEM):
+    """Fancy Field Storage for a TDEM survey.
+
+        This is used when solving Ahat and AhatT
+    """
+    knownFields = {'b': 'F'}
+    aliasFields = {'e': ['b','E','e_from_b']}
+    p = None
+
+    def startup(self):
+        self.MeSigmaI  = self.survey.prob.MeSigmaI
+        self.edgeCurlT = self.survey.prob.mesh.edgeCurl.T
+        self.MfMui     = self.survey.prob.MfMui
+
+    def e_from_b(self, y_b, txInd, tInd):
+        # TODO: implement non-zero js
+        y_e = self.MeSigmaI*(self.edgeCurlT*(self.MfMui*y_b))
+        if 'e' in self.p:
+            y_e = y_e - self.MeSigmaI*self.p[txInd,'e',tInd]
+        return y_e
+
 class ProblemTDEM_b(BaseTDEMProblem):
     """
         Time-Domain EM problem - B-formulation
@@ -36,7 +57,7 @@ class ProblemTDEM_b(BaseTDEMProblem):
     solType = 'b' #: Type of the solution, in this case the 'b' field
 
     surveyPair = SurveyTDEM
-    _FieldsTDEM_pair = FieldsTDEM_e_from_b  #: used for the forward calculation only
+    _FieldsForward_pair = FieldsTDEM_e_from_b     #: used for the forward calculation only
 
     ####################################################
     # Internal Methods
@@ -57,17 +78,6 @@ class ProblemTDEM_b(BaseTDEMProblem):
         RHS = (1.0/dt)*self.MfMui*B_n
         return RHS
 
-    def calcFields(self, sol, tInd):
-
-        if self.solType == 'b':
-            b = sol
-            # e = self.MeSigmaI*(self.mesh.edgeCurl.T*(self.MfMui*b))
-        else:
-            raise NotImplementedError('solType "%s" is not implemented in CalcFields.' % self.solType)
-
-        return {'b':b}
-
-
     ####################################################
     # Derivatives
     ####################################################
@@ -173,16 +183,9 @@ class ProblemTDEM_b(BaseTDEMProblem):
             dt = self.timeSteps[tInd]
             return rhs + 1.0/dt*self.MfMui*y[:,'b',tInd]
 
-        def AhCalcFields(sol, tInd):
-            y_b = sol
-            if self.survey.nTx == 1:
-                y_b = mkvc(y_b)
-            y_e = self.MeSigmaI*(self.mesh.edgeCurl.T*(self.MfMui*y_b))
-            if 'e' in p:
-                y_e = y_e - self.MeSigmaI*p[:,'e',tInd+1]
-            return {'b':y_b, 'e':y_e}
+        F = FieldsTDEM_e_from_b_Ah(self.mesh, self.survey, p=p)
 
-        return self.forward(m, AhRHS, AhCalcFields)
+        return self.forward(m, AhRHS, F)
 
     def solveAht(self, m, p):
         """
@@ -246,16 +249,9 @@ class ProblemTDEM_b(BaseTDEMProblem):
             dt = self.timeSteps[tInd+1]
             return rhs + 1.0/dt*self.MfMui*y[:,'b',tInd+2]
 
-        def AhtCalcFields(sol, tInd):
-            y_b = sol
-            if self.survey.nTx == 1:
-                y_b = mkvc(y_b)
-            y_e = self.MeSigmaI*(self.mesh.edgeCurl.T*(self.MfMui*y_b))
-            if 'e' in p:
-                y_e += - self.MeSigmaI*p[:,'e',tInd+1]
-            return {'b':y_b, 'e':y_e}
+        F = FieldsTDEM_e_from_b_Ah(self.mesh, self.survey, p=p)
 
-        return self.adjoint(m, AhtRHS, AhtCalcFields)
+        return self.adjoint(m, AhtRHS, F)
 
     ####################################################
     # Functions for tests
diff --git a/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py
index cd054d34..e17343e0 100644
--- a/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py
+++ b/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py
@@ -33,6 +33,12 @@ class TDEM_bDerivTests(unittest.TestCase):
         self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
         # self.prb.timeSteps = [(1e-05, 100)]
 
+        try:
+            from pymatsolver import MumpsSolver
+            self.prb.Solver = MumpsSolver
+        except ImportError, e:
+            self.prb.Solver  = SolverLU
+
         self.sigma = np.ones(mesh.nCz)*1e-8
         self.sigma[mesh.vectorCCz<0] = 1e-1
         self.sigma = np.log(self.sigma[active])
diff --git a/simpegEM/Tests/test_TDEM_combos.py b/simpegEM/Tests/test_TDEM_combos.py
index 286438ab..69f9e397 100644
--- a/simpegEM/Tests/test_TDEM_combos.py
+++ b/simpegEM/Tests/test_TDEM_combos.py
@@ -2,11 +2,6 @@ import unittest
 from SimPEG import *
 import simpegEM as EM
 
-try:
-    from pymatsolver import MumpsSolver
-except ImportError, e:
-    MumpsSolver = SolverLU
-
 plotIt = False
 
 def getProb(meshType='CYL',rxTypes='bx,bz',nTx=1):
@@ -35,7 +30,12 @@ def getProb(meshType='CYL',rxTypes='bx,bz',nTx=1):
     # prb.timeSteps = [1e-5]
     prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
     # prb.timeSteps = [(1e-05, 100)]
-    prb.Solver = MumpsSolver
+
+    try:
+        from pymatsolver import MumpsSolver
+        prb.Solver = MumpsSolver
+    except ImportError, e:
+        prb.Solver  = SolverLU
 
     sigma = np.ones(mesh.nCz)*1e-8
     sigma[mesh.vectorCCz<0] = 1e-1