Skip to content

Commit

Permalink
Updates to Mag problem for new property maps.
Browse files Browse the repository at this point in the history
  • Loading branch information
rowanc1 committed Oct 23, 2016
1 parent a577812 commit fed2566
Show file tree
Hide file tree
Showing 3 changed files with 48 additions and 16 deletions.
35 changes: 27 additions & 8 deletions SimPEG/PF/Magnetics.py
@@ -1,9 +1,15 @@
from __future__ import print_function

import numpy as np
import scipy.sparse as sp
from SimPEG import Utils, Problem, Solver
from . import BaseMag as MAG
from scipy.constants import mu_0

from SimPEG import Utils
from SimPEG import Problem
from SimPEG import Solver
from SimPEG import Props

from . import BaseMag as MAG
from .MagAnalytics import spheremodel, CongruousMagBC


Expand Down Expand Up @@ -286,8 +292,21 @@ class Problem3D_DiffSecondary(Problem.BaseProblem):
surveyPair = MAG.BaseMagSurvey
modelPair = MAG.BaseMagMap

def __init__(self, model, mapping=None, **kwargs):
Problem.BaseProblem.__init__(self, model, mapping=mapping, **kwargs)
_depreciate_main_map = 'muMap'

mu, muMap, muDeriv = Props.Invertible(
"Magnetic Permeability (H/m)",
default=mu_0
)

mui, muiMap, muiDeriv = Props.Invertible(
"Inverse Magnetic Permeability (m/H)"
)

Props.Reciprocal(mu, mui)

def __init__(self, mesh, **kwargs):
Problem.BaseProblem.__init__(self, mesh, **kwargs)

Pbc, Pin, self._Pout = \
self.mesh.getBCProjWF('neumann', discretization='CC')
Expand All @@ -306,7 +325,7 @@ def MfMui(self): return self._MfMui
def MfMu0(self): return self._MfMu0

def makeMassMatrices(self, m):
mu = self.mapping*m
mu = self.muMap*m
self._MfMui = self.mesh.getFaceInnerProduct(1./mu)/self.mesh.dim
# self._MfMui = self.mesh.getFaceInnerProduct(1./mu)
# TODO: this will break if tensor mu
Expand Down Expand Up @@ -334,7 +353,7 @@ def getRHS(self, m):
Dface = self.mesh.faceDiv
Mc = Utils.sdiag(self.mesh.vol)

mu = self.mapping*m
mu = self.muMap*m
chi = mu/mu_0-1

# Temporary fix
Expand Down Expand Up @@ -458,8 +477,8 @@ def Jvec(self, m, v, u=None):
u = self.fields(m)

B, u = u['B'], u['u']
mu = self.mapping*(m)
dmudm = self.mapping.deriv(m)
mu = self.muMap*(m)
dmudm = self.muDeriv
dchidmu = Utils.sdiag(1/mu_0*np.ones(self.mesh.nC))

vol = self.mesh.vol
Expand Down
3 changes: 1 addition & 2 deletions tests/pf/test_forward_PFproblem.py
@@ -1,6 +1,5 @@
import unittest
from SimPEG import Mesh, Utils, np, PF
import matplotlib.pyplot as plt


class MagFwdProblemTests(unittest.TestCase):
Expand All @@ -19,7 +18,7 @@ def setUp(self):
sph_ind = PF.MagAnalytics.spheremodel(M, 0., 0., 0., 100)
chi[sph_ind] = chiblk
model = PF.BaseMag.BaseMagMap(M)
prob = PF.Magnetics.Problem3D_DiffSecondary(M, mapping=model)
prob = PF.Magnetics.Problem3D_DiffSecondary(M, muMap=model)
self.prob = prob
self.M = M
self.chi = chi
Expand Down
26 changes: 20 additions & 6 deletions tests/pf/test_grav_inversion_linear.py
@@ -1,10 +1,21 @@
from __future__ import print_function
import unittest
from SimPEG import *
import numpy as np
from SimPEG import Mesh
from SimPEG import Utils
from SimPEG import Maps
from SimPEG import Regularization
from SimPEG import DataMisfit
from SimPEG import Optimization
from SimPEG import InvProblem
from SimPEG import Directives
from SimPEG import Inversion

import SimPEG.PF as PF

np.random.seed(43)


class GravInvLinProblemTest(unittest.TestCase):

def setUp(self):
Expand All @@ -28,7 +39,7 @@ def setUp(self):
zz = -np.exp((xx**2 + yy**2) / 75**2) + mesh.vectorNz[-1]

# Go from topo to actv cells
topo = np.c_[mkvc(xx), mkvc(yy), mkvc(zz)]
topo = np.c_[Utils.mkvc(xx), Utils.mkvc(yy), Utils.mkvc(zz)]
actv = Utils.surface2ind_topo(mesh, topo, 'N')
actv = np.asarray([inds for inds, elem in enumerate(actv, 1)
if elem], dtype=int) - 1
Expand All @@ -55,18 +66,21 @@ def setUp(self):
# Here a simple block in half-space
model = np.zeros((mesh.nCx, mesh.nCy, mesh.nCz))
model[(midx-2):(midx+2), (midy-2):(midy+2), -6:-2] = 0.5
model = mkvc(model)
model = Utils.mkvc(model)
self.model = model[actv]

# Create active map to go from reduce set to full
actvMap = Maps.InjectActiveCells(mesh, actv, ndv)

# Creat reduced identity map
# Create reduced identity map
idenMap = Maps.IdentityMap(nP=nC)

# Create the forward model operator
prob = PF.Gravity.GravityIntegral(mesh, mapping=idenMap,
actInd=actv)
prob = PF.Gravity.GravityIntegral(
mesh,
mapping=idenMap,
actInd=actv
)

# Pair the survey and problem
survey.pair(prob)
Expand Down

0 comments on commit fed2566

Please sign in to comment.