From 58b65e0acc440a62d6bd0db8827ec7e033e03715 Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Thu, 13 Feb 2014 18:33:12 -0800
Subject: [PATCH] working on adjoint

---
 simpegEM/TDEM/BaseTDEM.py   | 36 +++++++++++++++++++++++++++++++++---
 simpegEM/TDEM/FieldsTDEM.py |  2 +-
 simpegEM/TDEM/TDEM_b.py     |  2 +-
 3 files changed, 35 insertions(+), 5 deletions(-)

diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py
index 943c2611..d0e7f992 100644
--- a/simpegEM/TDEM/BaseTDEM.py
+++ b/simpegEM/TDEM/BaseTDEM.py
@@ -10,6 +10,8 @@ import numpy as np
 class MixinInitialFieldCalc(object):
     """docstring for MixinInitialFieldCalc"""
 
+    storeTheseFields = 'b'
+
     def getInitialFields(self):
         if self.data.txType == 'VMD_MVP':
             # Vertical magnetic dipole, magnetic vector potential
@@ -29,7 +31,7 @@ class MixinInitialFieldCalc(object):
             MVP = np.concatenate((MVPx, MVPy, MVPz))
 
         # Initialize field object
-        F = FieldsTDEM(self.mesh, 1, self.times.size, 'b')
+        F = FieldsTDEM(self.mesh, 1, self.times.size, store=self.storeTheseFields)
 
         # Set initial B
         F.b0 = self.mesh.edgeCurl*MVP
@@ -140,8 +142,14 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem):
         self.makeMassMatrices(m)
 
         F = self.getInitialFields()
-        #TODO: Split next code to forward and adjoint.
-        # fields would call forward
+
+        return self.forward(m, RHS, CalcFields, F=F)
+
+
+    def forward(self, m, RHS, CalcFields, F=None):
+        if F is None:
+            F = FieldsTDEM(self.mesh, self.data.nTx, self.nTimes, store=self.storeTheseFields)
+
         dtFact = None
         for tInd, t in enumerate(self.times):
             dt = self.getDt(tInd)
@@ -158,3 +166,25 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem):
             newFields = CalcFields(sol, self.solType, tInd)
             F.update(newFields, tInd)
         return F
+
+    def adjoint(self, m, RHS, CalcFields, F=None):
+        if F is None:
+            F = FieldsTDEM(self.mesh, self.data.nTx, self.nTimes, store=self.storeTheseFields)
+
+        dtFact = None
+        for tInd, t in reversed(list(enumerate(self.times))):
+            dt = self.getDt(tInd)
+            if dt!=dtFact:
+                dtFact = dt
+                A = self.getA(tInd)
+                # print 'Factoring...   (dt = ' + str(dt) + ')'
+                Asolve = Solver(A, options=self.solveOpts)
+                # print 'Done'
+            rhs = RHS(tInd, F)
+            sol = Asolve.solve(rhs)
+            if sol.ndim == 1:
+                sol.shape = (sol.size,1)
+            newFields = CalcFields(sol, self.solType, tInd)
+            F.update(newFields, tInd)
+        return F
+
diff --git a/simpegEM/TDEM/FieldsTDEM.py b/simpegEM/TDEM/FieldsTDEM.py
index e1bb482b..25234892 100644
--- a/simpegEM/TDEM/FieldsTDEM.py
+++ b/simpegEM/TDEM/FieldsTDEM.py
@@ -18,7 +18,7 @@ class FieldsTDEM(object):
     j = None #: Current density
     h = None #: Magnetic field
 
-    def __init__(self, mesh, nTx, nTimes, store):
+    def __init__(self, mesh, nTx, nTimes, store='b'):
 
         self.nTimes = nTimes #: Number of times
         self.nTx = nTx #: Number of transmitters
diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py
index 58d44339..a72dda53 100644
--- a/simpegEM/TDEM/TDEM_b.py
+++ b/simpegEM/TDEM/TDEM_b.py
@@ -40,7 +40,7 @@ class ProblemTDEM_b(ProblemBaseTDEM):
     # Derivatives
     ####################################################
 
-    def J(self, m, v, u=None):
+    def Jvec(self, m, v, u=None):
         if u is None:
             u = self.fields(m)
         p = self.G(m, v, u)