Construct K (Curl-Curl) matrix during setup. Modify K constructor fun…

…ction to delete temporary arrays that had been stored in the mesh structure. Small change to explicit sensitivity matrix generator to improve code readability.

src/MaxwellTimeParam.jl

@@ -4,7 +4,7 @@ export MaxwellTimeParam, getMaxwellTimeParam
 # Inner and outer constructors are required to allow instantiating the type with
 # last three fields left uninitialized. Use getMaxwellTimeParam function
 # for instantiation. It validates input and sets things to default values.
-type MaxwellTimeParam{S<:Real,T<:AbstractSolver,U<:AbstractSolver} <: ForwardProbType
+mutable struct MaxwellTimeParam{S<:Real,T<:AbstractSolver,U<:AbstractSolver} <: ForwardProbType
     Mesh::AbstractMesh
     Sources::AbstractArray{S}
     Obs::AbstractArray{S}
@@ -232,7 +232,7 @@ function getMaxwellTimeParam{S<:Real}(Mesh::AbstractMesh,
     # to observation times
     ObsTimeMat = getObsTimeMatrix(ObsTimes,t0,dt,size(Obs,2),size(s,2),sourceType)
 
-    if sourceType == :InductiveLoopPotential
+    if in(sourceType,[:InductiveLoopPotential,:InductiveDiscreteWire])
         K = getMaxwellCurlCurlMatrix(Mesh,fill(pi*4e-7,Mesh.nc))
         s = -0.5*K*s
         return MaxwellTimeParam(Mesh,s,Obs,ObsTimeMat,t0,dt,wave,sourceType,storageLevel,

src/getFields.jl

@@ -7,6 +7,8 @@
 # 3) TRBDF-2 (timeIntegrationMethod=:TRBDF2)
 #    TRBDF-2 seems to be unstable. Only for experimentation.
 
+import JOcTree.MassMatrix
+
 function getFieldsBE{T,N}(K::SparseMatrixCSC{T,N},
                           Msig::SparseMatrixCSC{T,N},
                           s::AbstractArray{T},
@@ -315,12 +317,21 @@ end
 
 function getMaxwellCurlCurlMatrix(M::AbstractMesh,mu)
     Curl   = getCurlMatrix(M)
+    ftype  = eltype(Curl)
+    itype  = eltype(Curl.colptr)
     Mmu    = getFaceMassMatrix(M,1./mu)
     Nf,Qf, = getFaceConstraints(M)
     Ne,    = getEdgeConstraints(M)
     Curl   = Qf*Curl*Ne
     Mmu    = Nf'*Mmu*Nf
     K      = Curl'*Mmu*Curl
+    M.Curl = spzeros(ftype,itype,0,0)
+    M.Nf   = spzeros(ftype,itype,0,0)
+    M.Qf   = spzeros(ftype,itype,0,0)
+    M.activeFaces = zeros(itype,0)
+    #sM.Qe   = spzeros(ftype,itype,0,0)
+    M.activeEdges = zeros(itype,0)
+    M.Pf   = Dict{Int64,MassMatrix}()
     return K
 end
 

src/getSensMatVec.jl

@@ -57,13 +57,11 @@ function getSensMatVec(DsigDmz::MaxwellTimeModel,model::MaxwellTimeModel,
         end
         mod = model.invertSigma ? model.sigma : model.mu
         if isempty(param.Sens)
-            nt = length(param.dt)
-            ns = size(param.Sources,2)
-            nr = size(param.Obs,2)
-            J  = zeros(size(param.ObsTimes,1), length(mod))
+            ndata     = size(param.ObsTimes,1)
+            J         = zeros(ndata, length(mod))
             sensTFunc = sensitivityTFunctions[param.timeIntegrationMethod]
             for k=1:size(J,1)
-            	v        = zeros(size(param.ObsTimes,1))
+            	v        = zeros(ndata)
             	v[k]     = 1.0
             	JTvStruc = sensTFunc(v,model,param)
             	JTv      = model.invertSigma ? JTvStruc.sigma : JTvStruc.mu

src/getSensTMatVec.jl

@@ -18,13 +18,11 @@ function getSensTMatVec{T<:Real}(z::Vector{T},model::MaxwellTimeModel,
         end
         mod = model.invertSigma ? model.sigma : model.mu
         if isempty(param.Sens)
-            nt = length(param.dt)
-            ns = size(param.Sources,2)
-            nr = size(param.Obs,2)
-            J  = zeros(size(param.ObsTimes,1), length(mod))
+            ndata     = size(param.ObsTimes,1)
+            J         = zeros(ndata, length(mod))
             sensTFunc = sensitivityTFunctions[param.timeIntegrationMethod]
             for k=1:size(J,1)
-            	v        = zeros(size(param.ObsTimes,1))
+            	v        = zeros(ndata)
             	v[k]     = 1.0
             	JTvStruc = sensTFunc(v,model,param)
             	JTv      = model.invertSigma ? JTvStruc.sigma : JTvStruc.mu