diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py
index c91d6c59..eec4bab9 100644
--- a/SimPEG/EM/FDEM/FDEM.py
+++ b/SimPEG/EM/FDEM/FDEM.py
@@ -161,8 +161,8 @@ class BaseFDEMProblem(BaseEMProblem):
 
         for i, src in enumerate(Srcs):
             smi, sei = src.eval(self)
-            S_m[:,i] += smi
-            S_e[:,i] += sei
+            S_m[:,i] = S_m[:,i] + smi
+            S_e[:,i] = S_e[:,i] + sei
 
         return S_m, S_e
 
diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py
index fa6a6485..f90dbfef 100644
--- a/SimPEG/EM/FDEM/FieldsFDEM.py
+++ b/SimPEG/EM/FDEM/FieldsFDEM.py
@@ -33,7 +33,7 @@ class Fields_e(Fields):
         ePrimary = np.zeros_like(eSolution)
         for i, src in enumerate(srcList):
             ep = src.ePrimary(self.prob)
-            ePrimary[:,i] += ep
+            ePrimary[:,i] = ePrimary[:,i] + ep
         return ePrimary
 
     def _eSecondary(self, eSolution, srcList):
@@ -53,7 +53,7 @@ class Fields_e(Fields):
         bPrimary = np.zeros([self._edgeCurl.shape[0],eSolution.shape[1]],dtype = complex)
         for i, src in enumerate(srcList):
             bp = src.bPrimary(self.prob)
-            bPrimary[:,i] += bp
+            bPrimary[:,i] = bPrimary[:,i] + bp
         return bPrimary
 
     def _bSecondary(self, eSolution, srcList):
@@ -62,8 +62,7 @@ class Fields_e(Fields):
         for i, src in enumerate(srcList):
             b[:,i] *= - 1./(1j*omega(src.freq))
             S_m, _ = src.eval(self.prob)
-            if S_m is not None:
-                b[:,i] += 1./(1j*omega(src.freq)) * S_m
+            b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * S_m
         return b
 
     def _bSecondaryDeriv_u(self, src, v, adjoint = False):
@@ -75,9 +74,7 @@ class Fields_e(Fields):
     def _bSecondaryDeriv_m(self, src, v, adjoint = False):
         S_mDeriv, _ = src.evalDeriv(self.prob, adjoint)
         S_mDeriv = S_mDeriv(v)
-        # if S_mDeriv is not None:
         return 1./(1j * omega(src.freq)) * S_mDeriv
-        # return None
 
     def _b(self, eSolution, srcList):
         return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList)
@@ -117,8 +114,7 @@ class Fields_b(Fields):
         bPrimary = np.zeros_like(bSolution)
         for i, src in enumerate(srcList):
             bp = src.bPrimary(self.prob)
-            # if bp is not None:
-            bPrimary[:,i] += bp
+            bPrimary[:,i] = bPrimary[:,i] + bp
         return bPrimary
 
     def _bSecondary(self, bSolution, srcList):
@@ -138,15 +134,14 @@ class Fields_b(Fields):
         ePrimary = np.zeros([self._edgeCurl.shape[1],bSolution.shape[1]],dtype = complex)
         for i,src in enumerate(srcList):
             ep = src.ePrimary(self.prob)
-            if ep is not None:
-                ePrimary[:,i] += ep
+            ePrimary[:,i] = ePrimary[:,i] + ep
         return ePrimary
 
     def _eSecondary(self, bSolution, srcList):
         e = self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * bSolution))
         for i,src in enumerate(srcList):
             _,S_e = src.eval(self.prob)
-            e[:,i] += -self._MeSigmaI * S_e
+            e[:,i] = e[:,i]+ -self._MeSigmaI * S_e
         return e
 
     def _eSecondaryDeriv_u(self, src, v, adjoint=False):
@@ -164,7 +159,6 @@ class Fields_b(Fields):
             Me = Me.T
 
         w = self._edgeCurl.T * (self._MfMui * bSolution)
-        # if S_e is not None:
         w += -Utils.mkvc(Me * S_e,2)
 
         if not adjoint:
@@ -176,7 +170,7 @@ class Fields_b(Fields):
         Se_Deriv = S_eDeriv(v)
 
         if Se_Deriv is not None:
-            de_dm += -self._MeSigmaI * Se_Deriv
+            de_dm = de_dm-self._MeSigmaI * Se_Deriv
 
         return de_dm
 
@@ -217,7 +211,7 @@ class Fields_j(Fields):
         jPrimary = np.zeros_like(jSolution,dtype = complex)
         for i, src in enumerate(srcList):
             jp = src.jPrimary(self.prob)
-            jPrimary[:,i] += jp
+            jPrimary[:,i] = jPrimary[:,i] + jp
         return jPrimary
 
     def _jSecondary(self, jSolution, srcList):
@@ -237,7 +231,7 @@ class Fields_j(Fields):
         hPrimary = np.zeros([self._edgeCurl.shape[1],jSolution.shape[1]],dtype = complex)
         for i, src in enumerate(srcList):
             hp = src.hPrimary(self.prob)
-            hPrimary[:,i] += hp
+            hPrimary[:,i] = hPrimary[:,i] + hp
         return hPrimary
 
     def _hSecondary(self, jSolution, srcList):
@@ -245,7 +239,7 @@ class Fields_j(Fields):
         for i, src in enumerate(srcList):
             h[:,i] *= -1./(1j*omega(src.freq))
             S_m,_ = src.eval(self.prob)
-            h[:,i] += 1./(1j*omega(src.freq)) * self._MeMuI * (S_m)
+            h[:,i] = h[:,i]+ 1./(1j*omega(src.freq)) * self._MeMuI * (S_m)
         return h
 
     def _hSecondaryDeriv_u(self, src, v, adjoint=False):
@@ -271,10 +265,10 @@ class Fields_j(Fields):
 
         if not adjoint:
             S_mDeriv = S_mDeriv(v)
-            hDeriv_m += 1./(1j*omega(src.freq)) * MeMuI * (Me * S_mDeriv)
+            hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * MeMuI * (Me * S_mDeriv)
         elif adjoint:
             S_mDeriv = S_mDeriv(Me.T * (MeMuI.T * v))
-            hDeriv_m += 1./(1j*omega(src.freq)) * S_mDeriv
+            hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * S_mDeriv
         return hDeriv_m
 
 
@@ -313,7 +307,7 @@ class Fields_h(Fields):
         hPrimary = np.zeros_like(hSolution,dtype = complex)
         for i, src in enumerate(srcList):
             hp = src.hPrimary(self.prob)
-            hPrimary[:,i] += hp
+            hPrimary[:,i] = hPrimary[:,i] + hp
         return hPrimary
 
     def _hSecondary(self, hSolution, srcList):
@@ -333,14 +327,14 @@ class Fields_h(Fields):
         jPrimary = np.zeros([self._edgeCurl.shape[0], hSolution.shape[1]], dtype = complex)
         for i, src in enumerate(srcList):
             jp = src.jPrimary(self.prob)
-            jPrimary[:,i] += jp
+            jPrimary[:,i] = jPrimary[:,i] + jp
         return jPrimary
 
     def _jSecondary(self, hSolution, srcList):
         j = self._edgeCurl*hSolution
         for i, src in enumerate(srcList):
             _,S_e = src.eval(self.prob)
-            j[:,i] += -S_e
+            j[:,i] = j[:,i]+ -S_e
         return j
 
     def _jSecondaryDeriv_u(self, src, v, adjoint=False):