From dd671e127765dfa547bbda531ff610b7f804129a Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Wed, 19 Mar 2014 11:24:30 -0700
Subject: [PATCH] Test with an analytic and some documentation.

---
 docs/api_Utils.rst                    |  7 +++
 simpegEM/FDEM/FDEM.py                 |  6 +--
 simpegEM/Tests/test_FDEM_analytics.py | 77 +++++++++++++++++++++++++++
 simpegEM/Utils/Ana/FEM.py             | 27 ++++++++--
 4 files changed, 109 insertions(+), 8 deletions(-)
 create mode 100644 simpegEM/Tests/test_FDEM_analytics.py

diff --git a/docs/api_Utils.rst b/docs/api_Utils.rst
index ef96eaf3..46a8e64a 100644
--- a/docs/api_Utils.rst
+++ b/docs/api_Utils.rst
@@ -11,6 +11,13 @@ Analytic Functions
     :inherited-members:
 
 
+.. automodule:: simpegEM.Utils.Ana.FEM
+    :show-inheritance:
+    :members:
+    :undoc-members:
+    :inherited-members:
+
+
 Sources
 *******
 
diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py
index d37bf038..80a82f0c 100644
--- a/simpegEM/FDEM/FDEM.py
+++ b/simpegEM/FDEM/FDEM.py
@@ -27,7 +27,7 @@ class BaseProblemFDEM(Problem.BaseProblem):
     dataPair = DataFDEM
 
     Solver = SimpegSolver
-    solverOpts = {'doDirect':True, 'options':{'factorize':False, 'backend':'scipy'}}
+    solverOpts = {}
 
     ####################################################
     # Mass Matrices
@@ -110,8 +110,6 @@ class ProblemFDEM_e(BaseProblemFDEM):
             solver = self.Solver(A, **self.solverOpts)
             e = solver.solve(rhs)
 
-            print np.linalg.norm(A*Utils.mkvc(e) - Utils.mkvc(rhs)) / np.linalg.norm(Utils.mkvc(rhs))
-
             F[freq, 'e'] = e
             b = -1./(1j*omega(freq))*self.mesh.edgeCurl*e
             F[freq, 'b'] = b
@@ -224,8 +222,6 @@ class ProblemFDEM_b(BaseProblemFDEM):
             #Note that we are solving for b_s
             b = solver.solve(rhs)
 
-            print np.linalg.norm(A*Utils.mkvc(b) - Utils.mkvc(rhs)) / np.linalg.norm(Utils.mkvc(rhs))
-
             F[freq, 'b'] = b
             e = self.MeSigmaI*self.mesh.edgeCurl.T*self.MfMui*b
             F[freq, 'e'] = e
diff --git a/simpegEM/Tests/test_FDEM_analytics.py b/simpegEM/Tests/test_FDEM_analytics.py
new file mode 100644
index 00000000..c4c29906
--- /dev/null
+++ b/simpegEM/Tests/test_FDEM_analytics.py
@@ -0,0 +1,77 @@
+import unittest
+from SimPEG import *
+import simpegEM as EM
+from scipy.constants import mu_0
+
+plotIt = False
+
+class FDEM_analyticTests(unittest.TestCase):
+
+    def setUp(self):
+
+        cs = 10.
+        ncx, ncy, ncz = 8, 8, 8
+        npad = 5
+        hx = Utils.meshTensors(((npad,cs), (ncx,cs), (npad,cs)))
+        hy = Utils.meshTensors(((npad,cs), (ncy,cs), (npad,cs)))
+        hz = Utils.meshTensors(((npad,cs), (ncz,cs), (npad,cs)))
+        mesh = Mesh.TensorMesh([hx,hy,hz], x0=[-hx.sum()/2.,-hy.sum()/2.,-hz.sum()/2.,])
+
+        model = Model.LogModel(mesh)
+
+        x = np.linspace(-10,10,5)
+        XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
+        rxList = EM.FDEM.RxListFDEM(XYZ, 'Ex')
+        Tx0 = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD', 1e2, rxList)
+
+        survey = EM.FDEM.SurveyFDEM([Tx0])
+
+        prb = EM.FDEM.ProblemFDEM_b(model)
+        prb.pair(survey)
+        prb.Solver = Utils.SolverUtils.DSolverWrap(sp.linalg.splu, checkAccuracy=False)
+
+        sig = 1e-1
+        sigma = np.ones(mesh.nC)*sig
+        sigma[mesh.gridCC[:,2] > 0] = 1e-8
+        m = np.log(sigma)
+
+        self.prb = prb
+        self.mesh = mesh
+        self.m = m
+        self.Tx0 = Tx0
+        self.sig = sig
+
+    def test_Transect(self):
+        print 'Testing Transect for analytic'
+
+        u = self.prb.fields(self.m)
+
+        bfz = self.mesh.r(u[self.Tx0, 'b'],'F','Fz','M')
+
+        x = np.linspace(-55,55,12)
+        XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
+
+        P = self.mesh.getInterpolationMat(XYZ, 'Fz')
+
+        an = EM.Utils.Ana.FEM.hzAnalyticDipoleF(x, self.Tx0.freq, self.sig)
+
+        diff = np.log10(np.abs(P*np.imag(u[self.Tx0, 'b']) - np.abs(mu_0*np.imag(an))))
+
+        if plotIt:
+            import matplotlib.pyplot as plt
+            plt.plot(x,np.log10(np.abs(P*np.imag(u[self.Tx0, 'b']))))
+            plt.plot(x,np.log10(np.abs(mu_0*np.imag(an))), 'r')
+            plt.plot(x,diff,'g')
+            plt.show()
+
+        # We want the difference to be an orderMag less
+        # than the analytic solution. Note that right at
+        # the source, both the analytic and the numerical
+        # solution will be poor. Use plotIt up top to see that...
+        orderMag = 1.6
+        passed = np.abs(np.mean(diff - np.log10(np.abs(mu_0*np.imag(an))))) > orderMag
+        self.assertTrue(passed)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/simpegEM/Utils/Ana/FEM.py b/simpegEM/Utils/Ana/FEM.py
index df99ece2..5a516ed6 100644
--- a/simpegEM/Utils/Ana/FEM.py
+++ b/simpegEM/Utils/Ana/FEM.py
@@ -2,9 +2,24 @@ import numpy as np
 from scipy.constants import mu_0, pi
 from scipy.special import erf
 
-def hzAnalyticDipoleF(r, freq, sigma):
+def hzAnalyticDipoleF(r, freq, sigma, secondary=True):
     """
     4.56 in Ward and Hohmann
+
+    .. plot::
+
+        import matplotlib.pyplot as plt
+        import simpegEM as EM
+        freq = np.logspace(-1, 6, 61)
+        test = EM.Utils.Ana.FEM.hzAnalyticDipoleF(100, freq, 0.001, secondary=False)
+        plt.loglog(freq, abs(test.real))
+        plt.loglog(freq, abs(test.imag))
+        plt.title('Response at $r$=100m')
+        plt.xlabel('Frequency')
+        plt.ylabel('Response')
+        plt.legend(('real','imag'))
+        plt.show()
+
     """
     r = np.abs(r)
     k = np.sqrt(-1j*2.*np.pi*freq*mu_0*sigma)
@@ -13,5 +28,11 @@ def hzAnalyticDipoleF(r, freq, sigma):
     front = m / (2. * np.pi * (k**2) * (r**5) )
     back = 9 - ( 9 + 9j * k * r - 4 * (k**2) * (r**2) - 1j * (k**3) * (r**3)) * np.exp(-1j*k*r)
     hz = front*back
-    hp =-1/(4*np.pi*r**3)
-    return hz-hp
+
+    if secondary:
+        hp =-1/(4*np.pi*r**3)
+        return hz-hp
+
+    return hz
+
+