diff --git a/.bumpversion.cfg b/.bumpversion.cfg index f5e90f2..e1793f3 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,4 +1,4 @@ [bumpversion] -current_version = 0.0.7 +current_version = 0.1.0 files = casingSimulations/info.py setup.py docs/conf.py diff --git a/.travis.yml b/.travis.yml index f425681..1c41eab 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,10 +1,6 @@ language: python python: - 3.6 -- "3.7-dev" -matrix: - allow_failures: - - python: "3.7-dev" env: - TEST_DIR=tests @@ -12,22 +8,14 @@ env: sudo: false before_install: - - if [ ${TRAVIS_PYTHON_VERSION:0:1} == "2" ]; then - wget http://repo.continuum.io/miniconda/Miniconda-3.8.3-Linux-x86_64.sh -O miniconda.sh; - else - wget http://repo.continuum.io/miniconda/Miniconda3-3.8.3-Linux-x86_64.sh -O miniconda.sh; - fi + - wget http://repo.continuum.io/miniconda/Miniconda3-3.8.3-Linux-x86_64.sh -O miniconda.sh; - chmod +x miniconda.sh - ./miniconda.sh -b -p $HOME/miniconda - export PATH=/home/travis/anaconda/bin:/home/travis/miniconda/bin:$PATH - conda update --yes conda install: - - if [ ${TRAVIS_PYTHON_VERSION:0:1} == "2" ]; then - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython jupyter mkl; - else - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython jupyter mkl; - fi + - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython jupyter mkl; - python setup.py install - pip install -r requirements_dev.txt diff --git a/casingSimulations/info.py b/casingSimulations/info.py index 1da67df..ab8a542 100644 --- a/casingSimulations/info.py +++ b/casingSimulations/info.py @@ -1,4 +1,4 @@ -__version__ = '0.0.7' +__version__ = '0.1.0' __author__ = 'Lindsey Heagy' __license__ = 'MIT' -__copyright__ = 'Copyright 2018 Lindsey Heagy' +__copyright__ = 'Copyright 2018-2019 Lindsey Heagy' diff --git a/casingSimulations/mesh.py b/casingSimulations/mesh.py index 01279ab..fb8e2e7 100644 --- a/casingSimulations/mesh.py +++ b/casingSimulations/mesh.py @@ -6,7 +6,7 @@ import inspect import properties -from SimPEG import Utils +from SimPEG.utils import setKwargs import discretize from discretize import utils @@ -48,7 +48,7 @@ class BaseMeshGenerator(BaseCasing): ) def __init__(self, **kwargs): - Utils.setKwargs(self, **kwargs) + setKwargs(self, **kwargs) @property def mesh(self): @@ -411,7 +411,7 @@ def hz(self): """ if getattr(self, '_hz', None) is None: - self._hz = Utils.meshTensor([ + self._hz = utils.meshTensor([ (self.csz, self.npadz, -self.pfz), (self.csz, self.ncz), (self.csz, self.npadz, self.pfz) @@ -531,7 +531,10 @@ def __init__(self, **kwargs): @property def ncx1(self): """number of cells with size csx1""" - return np.ceil(self.modelParameters.casing_b/self.csx1+2) + if getattr(self, '_ncx1', None) is None: + self._ncx1 = np.ceil(self.modelParameters.casing_b/self.csx1+2) + return self._ncx1 + @property def npadx1(self): @@ -546,22 +549,22 @@ def hx(self): if getattr(self, '_hx', None) is None: # finest uniform region - hx1a = Utils.meshTensor([(self.csx1, self.ncx1)]) + hx1a = utils.meshTensor([(self.csx1, self.ncx1)]) # pad to second uniform region - hx1b = Utils.meshTensor([(self.csx1, self.npadx1, self.pfx1)]) + hx1b = utils.meshTensor([(self.csx1, self.npadx1, self.pfx1)]) # scale padding so it matches cell size properly - dx1 = sum(hx1a)+sum(hx1b) - dx1 = np.floor(dx1/self.csx2) - hx1b *= (dx1*self.csx2 - sum(hx1a))/sum(hx1b) + dx1 = np.sum(hx1a)+np.sum(hx1b) + dx1 = 3 #np.floor(dx1/self.csx2) + hx1b *= (dx1*self.csx2 - np.sum(hx1a))/np.sum(hx1b) # second uniform chunk of mesh ncx2 = np.ceil((self.domain_x - dx1)/self.csx2) - hx2a = Utils.meshTensor([(self.csx2, ncx2)]) + hx2a = utils.meshTensor([(self.csx2, ncx2)]) # pad to infinity - hx2b = Utils.meshTensor([(self.csx2, self.npadx, self.pfx2)]) + hx2b = utils.meshTensor([(self.csx2, self.npadx, self.pfx2)]) self._hx = np.hstack([hx1a, hx1b, hx2a, hx2b]) diff --git a/casingSimulations/model.py b/casingSimulations/model.py index 45902c4..f34bac6 100644 --- a/casingSimulations/model.py +++ b/casingSimulations/model.py @@ -2,7 +2,8 @@ import properties import json import os -from SimPEG import Maps, Utils +from SimPEG import maps +from SimPEG.utils import setKwargs from scipy.constants import mu_0 import discretize @@ -132,7 +133,7 @@ class Wholespace(SurveyParametersMixin, BaseCasing): ) def __init__(self, filename=None, **kwargs): - Utils.setKwargs(self, **kwargs) + setKwargs(self, **kwargs) def __str__(self): return self.info @@ -858,10 +859,10 @@ def wires(self): """ wires to hook up maps to sigma, mu - :rtype: SimPEG.Maps.Wires + :rtype: SimPEG.maps.Wires """ if getattr(self, '_wires', None) is None: - self._wires = Maps.Wires( + self._wires = maps.Wires( ('sigma', self.mesh.nC), ('mu', self.mesh.nC) ) return self._wires diff --git a/casingSimulations/physics.py b/casingSimulations/physics.py index 2149731..701901d 100644 --- a/casingSimulations/physics.py +++ b/casingSimulations/physics.py @@ -1,4 +1,4 @@ -from SimPEG import Utils +from SimPEG.utils import ndgrid, mkvc, sdiag import discretize import numpy as np @@ -40,9 +40,9 @@ def casing_currents(j, mesh, model_parameters): (mesh.gridFz[:, 2] >= casing_z[0]) ) - jA = Utils.sdiag(mesh.area) * j + jA = sdiag(mesh.area) * j - jACasing = Utils.sdiag( + jACasing = sdiag( np.hstack([casing_faces_x, casing_faces_y, casing_faces_z]) ) * jA @@ -187,7 +187,7 @@ def plot_currents_over_freq( for i, f in enumerate(modelParameters.freqs): for srcind in srcinds: # src = survey.getSrcByFreq(survey.freqs[freqind])[srcind] - # j = Utils.mkvc(fields[mur][src, 'j'].copy()) + # j = mkvc(fields[mur][src, 'j'].copy()) Iind = i + srcind*len(modelParameters.freqs) @@ -360,7 +360,7 @@ def plot_j_over_mu_z( x_plt = np.r_[r] z_plt = np.linspace(xlim[0], xlim[1], int(xlim[1]-xlim[0])) - XYZ = Utils.ndgrid(x_plt, np.r_[0], z_plt) + XYZ = ndgrid(x_plt, np.r_[0], z_plt) Pfx = mesh.getInterpolationMat(XYZ, 'Fx') Pfz = mesh.getInterpolationMat(XYZ, 'Fz') @@ -391,10 +391,10 @@ def plot_j_over_mu_z( for i, mur in enumerate(modelParameters.muModels): for srcind in srcinds: src = survey.getSrcByFreq(survey.freqs[freqind])[srcind] - j = Utils.mkvc(fields[mur][src, 'j'].copy()) + j = mkvc(fields[mur][src, 'j'].copy()) if subtract is not None: - j = j - Utils.mkvc( + j = j - mkvc( fields[subtract][src, 'j'].copy() ) @@ -469,7 +469,7 @@ def plot_j_over_freq_z( x_plt = np.r_[r] z_plt = np.linspace(xlim[0], xlim[1], int(xlim[1]-xlim[0])) - XYZ = Utils.ndgrid(x_plt, np.r_[0], z_plt) + XYZ = ndgrid(x_plt, np.r_[0], z_plt) Pfx = mesh.getInterpolationMat(XYZ, 'Fx') Pfz = mesh.getInterpolationMat(XYZ, 'Fz') @@ -498,10 +498,10 @@ def plot_j_over_freq_z( for i, freq in enumerate(modelParameters.freqs): for srcind in srcinds: src = survey.getSrcByFreq(freq)[srcind] - j = Utils.mkvc(fields[mur][src, 'j'].copy()) + j = mkvc(fields[mur][src, 'j'].copy()) if subtract is not None: - j = j - Utils.mkvc( + j = j - mkvc( fields[subtract][src, 'j'].copy() ) @@ -576,7 +576,7 @@ def plot_j_over_mu_x( x_plt = np.linspace(xlim[0], xlim[1], xlim[1]) z_plt = np.r_[z] - XYZ = Utils.ndgrid(x_plt, np.r_[0], z_plt) + XYZ = ndgrid(x_plt, np.r_[0], z_plt) Pfx = mesh.getInterpolationMat(XYZ, 'Fx') Pfz = mesh.getInterpolationMat(XYZ, 'Fz') @@ -605,10 +605,10 @@ def plot_j_over_mu_x( for i, f in enumerate(modelParameters.freqs): for srcind in srcinds: src = survey.getSrcByFreq(survey.freqs[freqind])[srcind] - j = Utils.mkvc(fields[mur][src, 'j'].copy()) + j = mkvc(fields[mur][src, 'j'].copy()) if subtract is not None: - j = j - Utils.mkvc( + j = j - mkvc( fields[subtract][src, 'j'].copy() ) diff --git a/casingSimulations/run.py b/casingSimulations/run.py index 827462c..c40da3c 100644 --- a/casingSimulations/run.py +++ b/casingSimulations/run.py @@ -8,9 +8,11 @@ import discretize from discretize import utils import properties -from SimPEG.EM import FDEM, TDEM -from SimPEG import Utils, Maps -from SimPEG.EM.Static import DC +from SimPEG import maps +from SimPEG.electromagnetics import frequency_domain as fdem +from SimPEG.electromagnetics import time_domain as tdem +from SimPEG.electromagnetics import resistivity as dc +from SimPEG.utils import setKwargs try: from pymatsolver import Pardiso as Solver @@ -81,7 +83,7 @@ class BaseSimulation(BaseCasing): def __init__(self, **kwargs): # set keyword arguments - Utils.setKwargs(self, **kwargs) + setKwargs(self, **kwargs) # if the working directory does not exsist, create it if not os.path.isdir(self.directory): @@ -121,6 +123,11 @@ def prob(self): def survey(self): return self._survey + @property + def mesh(self): + return self.meshGenerator.mesh + + def write_py(self, includeDC=True, include2D=True): """ Write a python script for running the simulation @@ -222,7 +229,7 @@ class SimulationFDEM(BaseSimulation): choices=["e", "b", "h", "j"] ) - physics = "FDEM" + physics = "fdem" def __init__(self, **kwargs): super(SimulationFDEM, self).__init__(**kwargs) @@ -231,7 +238,7 @@ def __init__(self, **kwargs): def prob(self): if getattr(self, '_prob', None) is None: self._prob = getattr( - FDEM, 'Problem3D_{}'.format(self.formulation) + fdem, 'Problem3D_{}'.format(self.formulation) )( self.meshGenerator.mesh, sigmaMap=self.physprops.wires.sigma, @@ -241,12 +248,12 @@ def prob(self): ) if getattr(self, 'srcList') is not None: - self._survey = FDEM.Survey(self.srcList.srcList) + self._survey = fdem.Survey(self.srcList.srcList) elif getattr(self, 'src') is not None: - self._survey = FDEM.Survey(self.src.srcList) + self._survey = fdem.Survey(self.src.srcList) else: raise Exception("one of src, srcList must be set") - self._prob.pair(self._survey) + self._prob.survey = self._survey return self._prob @property @@ -269,7 +276,7 @@ class SimulationTDEM(BaseSimulation): choices=["e", "b", "h", "j"] ) - physics = "TDEM" + physics = "tdem" def __init__(self, **kwargs): super(SimulationTDEM, self).__init__(**kwargs) @@ -278,7 +285,7 @@ def __init__(self, **kwargs): def prob(self): if getattr(self, '_prob', None) is None: self._prob = getattr( - TDEM, 'Problem3D_{}'.format(self.formulation) + tdem, 'Problem3D_{}'.format(self.formulation) )( self.meshGenerator.mesh, timeSteps=self.modelParameters.timeSteps, @@ -288,9 +295,9 @@ def prob(self): verbose=self.verbose ) - self._survey = TDEM.Survey(self.srcList.srcList) + self._survey = tdem.Survey(self.srcList.srcList) - self._prob.pair(self._survey) + self._prob.survey = self._survey return self._prob @property @@ -324,22 +331,22 @@ class SimulationDC(BaseSimulation): default="phi" ) - physics = "DC" + physics = "dc" def __init__(self, **kwargs): super(SimulationDC, self).__init__(**kwargs) - self._prob = DC.Problem3D_CC( + self._prob = dc.Problem3D_CC( self.meshGenerator.mesh, sigmaMap=self.physprops.wires.sigma, bc_type='Dirichlet', Solver=Solver ) self._srcList = [ - DC.Src.Dipole([], self.src_a[i, :], self.src_b[i, :]) + dc.sources.Dipole([], self.src_a[i, :], self.src_b[i, :]) for i in range(self.src_a.shape[0]) ] # self._src = DC.Src.Dipole([], self.src_a, self.src_b) - self._survey = DC.Survey(self._srcList) + self._survey = dc.Survey(self._srcList) - self._prob.pair(self._survey) + self._prob.survey = self._survey diff --git a/casingSimulations/sources.py b/casingSimulations/sources.py index e77b011..1490270 100644 --- a/casingSimulations/sources.py +++ b/casingSimulations/sources.py @@ -6,8 +6,9 @@ import discretize from discretize.utils import closestPoints -from SimPEG import Utils -from SimPEG.EM import FDEM, TDEM +from SimPEG.utils import setKwargs +from SimPEG.electromagnetics import frequency_domain as fdem +from SimPEG.electromagnetics import time_domain as tdem from .base import LoadableInstance, BaseCasing from . import model @@ -36,8 +37,8 @@ class BaseCasingSrc(BaseCasing): ) physics = properties.StringChoice( - "FDEM or TDEM simulation?", - choices=["FDEM", "TDEM"], + "fdem or tdem simulation?", + choices=["fdem", "tdem"], required=False ) @@ -50,7 +51,7 @@ class BaseCasingSrc(BaseCasing): ) def __init__(self, **kwargs): - Utils.setKwargs(self, **kwargs) + setKwargs(self, **kwargs) if self.src_a is None: self.src_a = self.modelParameters.src_a @@ -115,13 +116,13 @@ def srcList(self): Source List """ if getattr(self, '_srcList', None) is None: - if self.physics == "FDEM": + if self.physics.lower() == "fdem": srcList = [ - FDEM.Src.RawVec_e([], f, self.s_e.astype("complex")) + fdem.sources.RawVec_e([], f, self.s_e.astype("complex")) for f in self.freqs ] - elif self.physics == "TDEM": - srcList = [TDEM.Src.RawVec_Grounded([], self.s_e)] + elif self.physics == "tdem": + srcList = [tdem.sources.RawVec_Grounded([], self.s_e)] self._srcList = srcList return self._srcList @@ -511,7 +512,7 @@ def surface_wire(self): ) surface_wirez = ( (mesh.gridFx[:, 2] > mesh.hz.min()) & - (mesh.gridFx[:, 2] <= 1.75*mesh.hz.min()) + (mesh.gridFx[:, 2] <= 2*mesh.hz.min()) ) self._surface_wire = surface_wirex & surface_wirez @@ -542,7 +543,7 @@ def surface_electrode(self): ) surface_electrodez = ( (mesh.gridFz[:, 2] >= src_b[2] - mesh.hz.min()) & - (mesh.gridFz[:, 2] < src_b[2] + 2*mesh.hz.min()) + (mesh.gridFz[:, 2] < 2*mesh.hz.min()) ) self._surface_electrode = surface_electrodex & surface_electrodez @@ -765,8 +766,8 @@ def positive_electrode(self): positive_electrodex = (mesh.gridFz[:, 0] == src_a[0]) positive_electrodez = ( - (mesh.gridFz[:, 2] >= src_a[2] - mesh.hz.min()) & - (mesh.gridFz[:, 2] <= src_a[2] + 2*mesh.hz.min()) + (mesh.gridFz[:, 2] >= src_a[2]) & + (mesh.gridFz[:, 2] < 1.5*mesh.hz.min()) ) self._positive_electrode = ( diff --git a/casingSimulations/utils.py b/casingSimulations/utils.py index 039700d..5a3aee1 100644 --- a/casingSimulations/utils.py +++ b/casingSimulations/utils.py @@ -3,7 +3,7 @@ import time import os import datetime -from discretize import utils +from discretize import utils, CylMesh import properties import casingSimulations @@ -86,6 +86,14 @@ def ccv3DthetaSlice(mesh3D, v3D, theta_ind=0): utils.mkvc(ccv_z[:, theta_ind, :], 2) ]) +def mesh2d_from_3d(mesh): + """ + create cylindrically symmetric mesh generator + """ + mesh2D = CylMesh( + [mesh.hx, 1., mesh.hz], x0=mesh.x0 + ) + return mesh2D def writeSimulationPy( modelParameters='ModelParameters.json', diff --git a/casingSimulations/view.py b/casingSimulations/view.py index 85b12ce..3fd04b9 100644 --- a/casingSimulations/view.py +++ b/casingSimulations/view.py @@ -4,12 +4,20 @@ import matplotlib.pyplot as plt from matplotlib.colors import LogNorm, SymLogNorm from scipy.spatial import cKDTree +from scipy.constants import mu_0 import ipywidgets import discretize import properties -from .utils import face3DthetaSlice +from .utils import face3DthetaSlice, mesh2d_from_3d + +from SimPEG.electromagnetics.static.resistivity.fields import FieldsDC +from SimPEG.electromagnetics.time_domain.fields import ( + FieldsTDEM, Fields3D_b, Fields3D_e, Fields3D_h, Fields3D_j +) +from SimPEG.electromagnetics.frequency_domain.fields import FieldsFDEM + # from .run import SimulationDC, SimulationFDEM, SimulationTDEM @@ -33,9 +41,7 @@ def plot_slice( pcolorOpts = {} # generate a 2D mesh for plotting slices - mesh2D = discretize.CylMesh( - [mesh.hx, 1., mesh.hz], x0=mesh.x0 - ) + mesh2D = mesh2d_from_3d(mesh) vplt = v.reshape(mesh.vnC, order='F') plotme = discretize.utils.mkvc(vplt[:, theta_ind, :]) @@ -199,8 +205,7 @@ def plotLinesFx( linestyle='-', alpha=None ): - - mesh2D = discretize.CylMesh([mesh.hx, 1., mesh.hz], x0=mesh.x0) + mesh2D = mesh2d_from_3d(mesh) if ax is None: fig, ax = plt.subplots(1, 1, figsize=(8, 6)) @@ -255,43 +260,37 @@ class FieldsViewer(properties.HasProperties): ) def __init__( - self, sim_dict, fields_dict, model_keys=None, + self, mesh, model_parameters_dict, survey_dict, fields_dict, model_keys=None, **kwargs ): super(FieldsViewer, self).__init__(**kwargs) - self.sim_dict = sim_dict + # self.sim_dict = sim_dict + self.model_parameters_dict = model_parameters_dict + self.mesh = mesh + self.survey_dict = survey_dict self.fields_dict = fields_dict self.model_keys = ( - model_keys if model_keys is not None else sorted(sim_dict.keys()) + model_keys if model_keys is not None else sorted(fields_dict.keys()) ) - if all( - sim.__class__.__name__ == "SimulationDC" - for sim in sim_dict.values() - ): + if all(isinstance(f, FieldsDC) for f in fields_dict.values()): self._physics = 'DC' self.fields_opts = [ 'sigma', 'e', 'j', 'phi', 'charge', 'charge_density' ] self.reim_opts = ["real"] - elif all( - sim.__class__.__name__ == "SimulationFDEM" - for sim in sim_dict.values() - ): + elif all(isinstance(f, FieldsFDEM) for f in fields_dict.values()): self._physics = 'FDEM' self.fields_opts = ['sigma', 'mur', 'e', 'j', 'h', 'b'] self.reim_opts = ["real", "imag"] - elif all( - sim.__class__.__name__ == "SimulationTDEM" - for sim in sim_dict.values() - ): + elif all(isinstance(f, FieldsTDEM) for f in fields_dict.values()): self._physics = 'TDEM' self.fields_opts = ['sigma', 'mur', 'e', 'j', 'dbdt', 'dhdt'] - if self.sim_dict[model_keys[0]].prob._fieldType=="j": + if isinstance(self.fields_dict[model_keys[0]], Fields3D_j): self.fields_opts += ["charge", "charge_density"] - elif self.sim_dict[model_keys[0]].prob._fieldType in ["b", "h"]: + elif isinstance(self.fields_dict[model_keys[0]], (Fields3D_b, Fields3D_h)): self.fields_opts += ['sigma', 'mur', 'h', 'b'] self.reim_opts = ["real"] @@ -307,13 +306,11 @@ def prim_sec_opts(self): else: return None - def _mesh2D(self, model_key): - if self.sim_dict[model_key].meshGenerator.mesh.isSymmetric: - return self.sim_dict[model_key].meshGenerator.mesh - return self.sim_dict[model_key].meshGenerator.create_2D_mesh().mesh - - def _mesh(self, model_key): - return self.sim_dict[model_key].meshGenerator.mesh + @property + def _mesh2D(self): + if self.mesh.isSymmetric: + return self.mesh + return mesh2d_from_3d(self.mesh) def _check_inputs( self, model_key, view, prim_sec, real_or_imag @@ -344,6 +341,19 @@ def run_assert(field, provided, allowed): run_assert('real_or_imag', real_or_imag,['real', 'imag']) + def fetch_field(self, model_key, view, src=None, time_ind=None): + + if view in ['sigma', 'mur']: + plotme = getattr(self.model_parameters_dict[model_key], view)(self.mesh) + else: + if self._physics == "TDEM": + plotme = self.fields_dict[model_key][src, view, time_ind] + else: + plotme = self.fields_dict[model_key][src, view] + + return plotme + + def plot_cross_section( self, ax=None, @@ -383,52 +393,32 @@ def plot_cross_section( if prim_sec is None and self.prim_sec_opts is not None: prim_sec = 'total' - self._check_inputs( - model_key, view, prim_sec, real_or_imag - ) + self._check_inputs(model_key, view, prim_sec, real_or_imag) # get background model if doing primsec if prim_sec == 'primary': model_key = self.primary_key # grab relevant parameters - src = self.sim_dict[model_key].survey.srcList[src_ind] - if view in ['sigma', 'mur']: - plotme = getattr(self.sim_dict[model_key].physprops, view) - else: - if self._physics == "TDEM": - plotme = self.fields_dict[model_key][src, view, time_ind] - else: - plotme = self.fields_dict[model_key][src, view] - - mesh = self._mesh(model_key) + src = self.survey_dict[model_key].source_list[src_ind] + plotme = self.fetch_field(model_key, view, src, time_ind) norm = None + if view == "sigma": + norm = LogNorm() if prim_sec in ['secondary', 'percent']: - if view in ['sigma', 'mur']: - background = getattr( - self.sim_dict[self.primary_key].physprops, model_key - ) - else: - prim_src = self.sim_dict[self.primary_key].survey.srcList[src_ind] - if self._physics == "TDEM": - background = self.fields_dict[self.primary_key][ - prim_src, view, time_ind - ] - else: - background = self.fields_dict[self.primary_key][ - prim_src, view - ] + prim_src = self.survey_dict[self.primary_key].source_list[src_ind] + background = self.fetch_field(self.primary_key, view, prim_src, time_ind) plotme = plotme - background if prim_sec == "percent": plotme = 100 * plotme / (np.absolute(background) + self.eps) - if not mesh.isSymmetric: + if not self.mesh.isSymmetric: theta_ind_mirror = ( - theta_ind+int(mesh.vnC[1]/2) - if theta_ind < int(mesh.vnC[1]/2) - else theta_ind-int(mesh.vnC[1]/2) + theta_ind+int(self.mesh.vnC[1]/2) + if theta_ind < int(self.mesh.vnC[1]/2) + else theta_ind-int(self.mesh.vnC[1]/2) ) else: theta_ind_mirror = 0 @@ -449,7 +439,7 @@ def plot_cross_section( clim = clim[1]*np.r_[-1., 1.] if clim is not None else None # if not mesh.isSymmetric: - plotme = plotme.reshape(mesh.vnC, order='F') + plotme = plotme.reshape(self.mesh.vnC, order='F') mirror_data = discretize.utils.mkvc( plotme[:, theta_ind_mirror, :] ) @@ -457,18 +447,13 @@ def plot_cross_section( elif view in ['j', 'e']: - if self.sim_dict[model_key].prob._formulation == "HJ": - plt_vec = face3DthetaSlice( - mesh, plotme, - theta_ind=theta_ind - ) - mirror_data = face3DthetaSlice( - mesh, plotme, theta_ind=theta_ind_mirror - ) + if len(plotme) == self.mesh.vnF.sum(): + plt_vec = face3DthetaSlice(self.mesh, plotme, theta_ind=theta_ind) + mirror_data = face3DthetaSlice(self.mesh, plotme, theta_ind=theta_ind_mirror) plot_type = "vec" - elif self.sim_dict[model_key].prob._formulation == "EB": - plotme = mesh.aveE2CC * plotme + elif len(plotme) == self.mesh.vnF.sum(): + plotme = self.mesh.aveE2CC * plotme mirror_data = -plotme plot_type = "scalar" norm = SymLogNorm( @@ -479,32 +464,24 @@ def plot_cross_section( elif view in ['h', 'b', 'dbdt', 'dhdt']: - if self.sim_dict[model_key].prob._formulation == "EB": - plt_vec = face3DthetaSlice( - mesh, plotme, theta_ind=theta_ind - ) - mirror_data = face3DthetaSlice( - mesh, plotme, theta_ind=theta_ind_mirror - ) + if len(plotme) == self.mesh.vnF.sum(): + plt_vec = face3DthetaSlice(self.mesh, plotme, theta_ind=theta_ind) + mirror_data = face3DthetaSlice(self.mesh, plotme, theta_ind=theta_ind_mirror) plot_type = "vec" - elif self.sim_dict[model_key].prob._formulation == "HJ": + else: plot_type = "scalar" - if len(mesh.hy) == 1: - plotme = mesh.aveE2CC * plotme + if len(self.mesh.hy) == 1: + plotme = self.mesh.aveE2CC * plotme else: - plotme = (mesh.aveE2CCV * plotme)[mesh.nC:2*mesh.nC] + plotme = (self.mesh.aveE2CCV * plotme)[self.mesh.nC:2*self.mesh.nC] - plotme = plotme.reshape(mesh.vnC, order="F") - mirror_data = discretize.utils.mkvc( - -plotme[:, theta_ind_mirror, :] - ) - plotme = discretize.utils.mkvc( - plotme[:, theta_ind, :] - ) + plotme = plotme.reshape(self.mesh.vnC, order="F") + mirror_data = discretize.utils.mkvc(-plotme[:, theta_ind_mirror, :]) + plotme = discretize.utils.mkvc(plotme[:, theta_ind, :]) norm = SymLogNorm( clim[0] if clim is not None else @@ -514,10 +491,10 @@ def plot_cross_section( if plot_type == "scalar": - out = self._mesh2D(model_key).plotImage( + out = self._mesh2D.plotImage( getattr(plotme, real_or_imag), ax=ax, pcolorOpts = { - 'cmap': 'bwr' if view in ['charge', 'charge_density'] else 'viridis', + 'cmap': 'RdBu_r' if view in ['charge', 'charge_density'] else 'viridis', 'norm': norm }, clim=clim, @@ -535,7 +512,7 @@ def plot_cross_section( elif plot_type == "vec": out = plotFace2D( - self._mesh2D(model_key), + self._mesh2D, plt_vec + self.eps, real_or_imag=real_or_imag, ax=ax, @@ -560,14 +537,14 @@ def plot_cross_section( cb.update_ticks() if show_mesh is True: - self._mesh2D(model_key).plotGrid(ax=ax) + self._mesh2D.plotGrid(ax=ax) title = "{} \n{} {}".format(model_key, prim_sec, view) if self._physics == "FDEM": - title += "\nf = {:1.1e} Hz".format(src.freq) + title += "\nf = {:1.1e} Hz".format(src.frequency) elif self._physics == "TDEM": title += "\n t = {:1.1e} s".format( - self.sim_dict[model_key].prob.times[time_ind] + self.fields_dict[model_key]._times[time_ind] ) ax.set_title( title, fontsize=13 @@ -577,7 +554,7 @@ def plot_cross_section( # plot outline of casing if casing_outline is True: - m = self.sim_dict[model_key].modelParameters + m = self.model_parameters_dict[model_key] factor = [-1, 1] [ ax.plot( @@ -602,9 +579,8 @@ def plot_cross_section( def _get_cKDTree(self, model_key, plan_mesh, k=10, theta_shift=None): # construct interpolation - mesh = self._mesh(model_key) - CCcart = mesh.cartesianGrid('CC', theta_shift=theta_shift) - tree = cKDTree(CCcart[:mesh.vnC[:2].prod(),:2]) + CCcart = self.mesh.cartesianGrid('CC', theta_shift=theta_shift) + tree = cKDTree(CCcart[:self.mesh.vnC[:2].prod(),:2]) d, ii = tree.query(plan_mesh.gridCC, k=k) weights = 1./d @@ -662,51 +638,33 @@ def plot_depth_slice( model_key = self.primary_key # grab relevant parameters + src = self.survey_dict[model_key].source_list[src_ind] + plotme = self.fetch_field(model_key, view, src, time_ind) norm = None - src = self.sim_dict[model_key].survey.srcList[src_ind] - if view in ['sigma', 'mur']: - plotme = getattr(self.sim_dict[model_key].physprops, view) - if view == "sigma": - norm = LogNorm() - else: - if self._physics == "TDEM": - plotme = self.fields_dict[model_key][src, view, time_ind] - else: - plotme = self.fields_dict[model_key][src, view] - mesh = self._mesh(model_key) + if view == "sigma": + norm = LogNorm() if prim_sec in ['secondary', 'percent']: - if view in ['sigma', 'mur']: - background = getattr( - self.sim_dict[self.primary_key].physprops, view - ) - else: - prim_src = self.sim_dict[self.primary_key].survey.srcList[src_ind] - if self._physics == "TDEM": - background = self.fields_dict[self.primary_key][ - prim_src, view, time_ind - ] - else: - background = self.fields_dict[self.primary_key][ - prim_src, view - ] + prim_src = self.survey_dict[self.primary_key].source_list[src_ind] + background = self.fetch_field(self.primary_key, view, prim_src, time_ind) plotme = plotme - background + if prim_sec == "percent": + plotme = 100 * plotme / (np.absolute(background) + self.eps) + # interpolate to cell centers if view in ['sigma', 'mur', 'phi', 'charge', 'charge_density']: plotme_cart = discretize.utils.mkvc(plotme) else: - if self.sim_dict[model_key].prob._formulation == 'HJ': - ave = mesh.aveF2CCV if view in ['e', 'j'] else mesh.aveE2CCV - elif self.sim_dict[model_key].prob._formulation == 'EB': - ave = mesh.aveF2CCV if view in ['h', 'b', 'dbdt', 'dhdt'] else mesh.aveE2CCV + if len(plotme) == self.mesh.vnE.sum(): + ave = self.mesh.aveE2CCV #if view in ['e', 'j'] else self.mesh.aveE2CCV + elif len(plotme) == self.mesh.vnF.sum(): + ave = self.mesh.aveF2CCV #if view in ['h', 'b', 'dbdt', 'dhdt'] else self.mesh.aveE2CCV plotme = ave * plotme - plotme = plotme.reshape(mesh.gridCC.shape, order='F') + plotme = plotme.reshape(self.mesh.gridCC.shape, order='F') if prim_sec == "percent": - background = (ave * background).reshape( - mesh.gridCC.shape, order='F' - ) + background = (ave * background).reshape(self.mesh.gridCC.shape, order='F') if denominator is None or denominator == "magnitude": den = np.outer(np.sqrt((background**2).sum(1)), np.ones(3)) elif denominator == "component": @@ -718,7 +676,7 @@ def plot_depth_slice( plotme = 100 * plotme / (den + self.eps) - gridCC = mesh.gridCC.copy() + gridCC = self.mesh.gridCC.copy() if theta_shift is not None: gridCC[:, 1] = gridCC[:, 1] - theta_shift @@ -739,10 +697,10 @@ def plot_depth_slice( # deal with vectors if view in ['e', 'b', 'h', 'j', 'dbdt', 'dhdt']: plotme_x = discretize.utils.mkvc( - (plotme_cart[:, 0]).reshape(mesh.vnC, order='F')[:, :, z_ind] + (plotme_cart[:, 0]).reshape(self.mesh.vnC, order='F')[:, :, z_ind] ) plotme_y = discretize.utils.mkvc( - (plotme_cart[:, 1]).reshape(mesh.vnC, order='F')[:, :, z_ind] + (plotme_cart[:, 1]).reshape(self.mesh.vnC, order='F')[:, :, z_ind] ) plotme = np.hstack([ @@ -782,7 +740,7 @@ def plot_depth_slice( out = plan_mesh.plotImage( getattr(plotme, real_or_imag), ax=ax, pcolorOpts = { - 'cmap': 'bwr' if view in ['charge', 'charge_density'] else 'viridis', + 'cmap': 'RdBu_r' if view in ['charge', 'charge_density'] else 'viridis', 'norm': norm }, clim=clim, @@ -795,13 +753,13 @@ def plot_depth_slice( ax.set_aspect(1) title = "{} \n{} {} \n z={:1.1e}m".format( - model_key, prim_sec, view, mesh.vectorCCz[z_ind] + model_key, prim_sec, view, self.mesh.vectorCCz[z_ind] ) if self._physics == "FDEM": - title += "\nf = {:1.1e} Hz".format(src.freq) + title += "\nf = {:1.1e} Hz".format(src.frequency) elif self._physics == "TDEM": title += "\n t = {:1.1e} s".format( - self.sim_dict[model_key].prob.times[time_ind] + self.fields_dict[model_key]._times[time_ind] ) ax.set_title(title) return out @@ -871,21 +829,21 @@ def widget_cross_section(self, ax=None, defaults={}, fixed={}, figwidth=5): widget_defaults = { "max_r": ( - 2*self.sim_dict[self.model_keys[0]].modelParameters.casing_b + 2*self.model_parameters_dict[self.model_keys[0]].casing_b if getattr( - self.sim_dict[self.model_keys[0]].modelParameters, + self.model_parameters_dict[self.model_keys[0]], 'casing_b', None ) is not None else - self.sim_dict[self.model_keys[0]].meshGenerator.domain_x + self.mesh.vectorNx.max() ), "min_depth": -10., "max_depth": ( - 1.25*self.sim_dict[self.model_keys[0]].modelParameters.casing_l + 1.25*self.model_parameters_dict[self.model_keys[0]].casing_l if getattr( - self.sim_dict[self.model_keys[0]].modelParameters, + self.model_parameters_dict[self.model_keys[0]], 'casing_b', None ) is not None else - self.sim_dict[self.model_keys[0]].meshGenerator.domain_z + - self.mesh.vectorNz.min() ), "clim_min": 0, "clim_max": 0, @@ -904,12 +862,12 @@ def widget_cross_section(self, ax=None, defaults={}, fixed={}, figwidth=5): fixed["ax"] = ax fixed["figwidth"] = figwidth - if not self.sim_dict[self.model_keys[0]].meshGenerator.mesh.isSymmetric: + if not self.mesh.isSymmetric: widget_defaults["theta_ind"]=0 else: fixed["theta_ind"] = 0 - if len(self.sim_dict[self.model_keys[0]].survey.srcList) == 1: + if len(self.survey_dict[self.model_keys[0]].source_list) == 1: fixed["src_ind"] = 0 else: widget_defaults["src_ind"] = 0 @@ -960,9 +918,9 @@ def widget_cross_section(self, ax=None, defaults={}, fixed={}, figwidth=5): for key, max_len in zip( ["theta_ind", "src_ind", "time_ind"], [ - self.sim_dict[self.model_keys[0]].meshGenerator.mesh.vnC[1] - 1, - len(self.sim_dict[self.model_keys[0]].survey.srcList) - 1, - self.sim_dict[self.model_keys[0]].prob.nT if + self.mesh.vnC[1] - 1, + len(self.survey_dict[self.model_keys[0]].source_list) - 1, + len(self.fields_dict[self.model_keys[0]]._times) if self._physics == "TDEM" else None ] ): @@ -1081,9 +1039,9 @@ def assign_cKDTree(): def widget_depth_slice(self, ax=None, defaults={}, fixed={}, figwidth=5): widget_defaults = { - "max_r": self.sim_dict[self.model_keys[0]].modelParameters.casing_l, + "max_r": self.model_parameters_dict[self.model_keys[0]].casing_l, "z_ind": int( - self.sim_dict[self.model_keys[0]].meshGenerator.mesh.vnC[2]/2. + self.mesh.vnC[2]/2. ), "clim_min": 0, "clim_max": 0, @@ -1106,7 +1064,7 @@ def widget_depth_slice(self, ax=None, defaults={}, fixed={}, figwidth=5): fixed["ax"] = ax fixed["figwidth"] = figwidth - if len(self.sim_dict[self.model_keys[0]].survey.srcList) == 1: + if len(self.survey_dict[self.model_keys[0]].source_list) == 1: fixed["src_ind"] = 0 else: widget_defaults["src_ind"] = 0 @@ -1157,9 +1115,9 @@ def widget_depth_slice(self, ax=None, defaults={}, fixed={}, figwidth=5): for key, max_len in zip( ["z_ind", "src_ind", "time_ind", "k"], [ - self.sim_dict[self.model_keys[0]].meshGenerator.mesh.vnC[2] - 1, - len(self.sim_dict[self.model_keys[0]].survey.srcList) - 1, - self.sim_dict[self.model_keys[0]].prob.nT if + self.mesh.vnC[2] - 1, + len(self.survey_dict[self.model_keys[0]].source_list) - 1, + len(self.fields_dict[self.model_keys[0]]._times) if self._physics == "TDEM" else None, 50 ] diff --git a/docs/conf.py b/docs/conf.py index 750c953..3e57d83 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -71,9 +71,9 @@ # built documents. # # The short X.Y version. -version = u'0.0.7' +version = u'0.1.0' # The full version, including alpha/beta/rc tags. -release = u'0.0.7' +release = u'0.1.0' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. @@ -148,7 +148,7 @@ # The name for this set of Sphinx documents. # " v documentation" by default. # -# html_title = u'EM Casing Research v0.0.7' +# html_title = u'EM Casing Research v0.1.0' # A shorter title for the navigation bar. Default is the same as html_title. # diff --git a/setup.py b/setup.py index 466641e..378ec50 100644 --- a/setup.py +++ b/setup.py @@ -28,18 +28,18 @@ setup( name="CasingSimulations", - version="0.0.7", + version="0.1.0", packages=find_packages(), install_requires=[ 'numpy>=1.7', 'scipy>=0.13', 'cython', 'pymatsolver>=0.1.1', - 'ipython', 'ipywidgets', 'jupyter', 'matplotlib', - 'properties[math]', + 'properties', + 'vectormath', 'SimPEG', ], @@ -49,8 +49,8 @@ long_description=LONG_DESCRIPTION, license="MIT", keywords="geophysics electromagnetics", - url="http://github.com/lheagy/casingResearch", - download_url="http://github.com/lheagy/casingResearch", + url="http://github.com/simpeg-research/casingResearch", + download_url="http://github.com/simpeg-research/casingResearch", classifiers=CLASSIFIERS, platforms=["Windows", "Linux", "Solaris", "Mac OS-X", "Unix"], use_2to3=False diff --git a/tests/test_cyl2D3D.py b/tests/test_cyl2D3D.py index 26ddfa7..a4bc749 100644 --- a/tests/test_cyl2D3D.py +++ b/tests/test_cyl2D3D.py @@ -7,8 +7,8 @@ from scipy.constants import mu_0 -from SimPEG.EM import FDEM -from SimPEG import Utils, Maps +from SimPEG.electromagnetics import frequency_domain as fdem +from SimPEG import utils, maps from pymatsolver import Pardiso @@ -110,11 +110,11 @@ def setUpClass(self): # create sources srcList2D = [ - FDEM.Src.RawVec_e(s_e=wire2D, freq=freq, rxList=[]) for freq in + fdem.sources.RawVec_e(s_e=wire2D, freq=freq, rxList=[]) for freq in modelParameters.freqs ] srcList3D = [ - FDEM.Src.RawVec_e(s_e=wire3D, freq=freq, rxList=[]) for freq in + fdem.sources.RawVec_e(s_e=wire3D, freq=freq, rxList=[]) for freq in modelParameters.freqs ] @@ -127,21 +127,21 @@ def setUpClass(self): ) # create the problems and surveys - prb2D = FDEM.Problem3D_h( + prb2D = fdem.Problem3D_h( mesh2D.mesh, sigmaMap=physprops2D.wires.sigma, muMap=physprops2D.wires.mu, Solver=Pardiso ) - prb3D = FDEM.Problem3D_h( + prb3D = fdem.Problem3D_h( mesh3D.mesh, sigmaMap=physprops3D.wires.sigma, muMap=physprops3D.wires.mu, Solver=Pardiso ) - survey2D = FDEM.Survey(srcList2D) - survey3D = FDEM.Survey(srcList3D) + survey2D = fdem.Survey(srcList2D) + survey3D = fdem.Survey(srcList3D) prb2D.pair(survey2D) prb3D.pair(survey3D)