Skip to content

Commit

Permalink
Merge pull request #14 from simpeg-research/simulation
Browse files Browse the repository at this point in the history
Simulation
  • Loading branch information
lheagy authored Jul 6, 2020
2 parents e68e4c8 + 98d5a4a commit ba55d58
Show file tree
Hide file tree
Showing 13 changed files with 216 additions and 250 deletions.
2 changes: 1 addition & 1 deletion .bumpversion.cfg
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
[bumpversion]
current_version = 0.0.7
current_version = 0.1.0
files = casingSimulations/info.py setup.py docs/conf.py

16 changes: 2 additions & 14 deletions .travis.yml
Original file line number Diff line number Diff line change
@@ -1,33 +1,21 @@
language: python
python:
- 3.6
- "3.7-dev"
matrix:
allow_failures:
- python: "3.7-dev"

env:
- TEST_DIR=tests

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

Expand Down
4 changes: 2 additions & 2 deletions casingSimulations/info.py
Original file line number Diff line number Diff line change
@@ -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'
25 changes: 14 additions & 11 deletions casingSimulations/mesh.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
import inspect

import properties
from SimPEG import Utils
from SimPEG.utils import setKwargs

import discretize
from discretize import utils
Expand Down Expand Up @@ -48,7 +48,7 @@ class BaseMeshGenerator(BaseCasing):
)

def __init__(self, **kwargs):
Utils.setKwargs(self, **kwargs)
setKwargs(self, **kwargs)

@property
def mesh(self):
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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):
Expand All @@ -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])

Expand Down
9 changes: 5 additions & 4 deletions casingSimulations/model.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
26 changes: 13 additions & 13 deletions casingSimulations/physics.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
from SimPEG import Utils
from SimPEG.utils import ndgrid, mkvc, sdiag
import discretize

import numpy as np
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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')
Expand Down Expand Up @@ -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()
)

Expand Down Expand Up @@ -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')
Expand Down Expand Up @@ -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()
)

Expand Down Expand Up @@ -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')
Expand Down Expand Up @@ -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()
)

Expand Down
43 changes: 25 additions & 18 deletions casingSimulations/run.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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,
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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,
Expand All @@ -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
Expand Down Expand Up @@ -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
Loading

0 comments on commit ba55d58

Please sign in to comment.