update to SimPEG 0.23.0 (#78)

* replace SimPEG import with simpeg

* use default solver rather than pardiso

* update SimPEG version

* update the model builder calls, test against python 3.11, drop testing against python 3.8

* add empymod to environment, update setKwargs to set_kwargs

* only test on python 3.10 for now

Author: lheagy

.github/workflows/tests.yml

@@ -15,7 +15,7 @@ jobs:
       fail-fast: false
       matrix:
         os: ["ubuntu-latest"]
-        python-version: ["3.8", "3.10"]
+        python-version: ["3.10"]
         test-file: ["tests/test_dcip.py", "tests/test_em.py", "tests/test_gpr.py tests/test_seismic.py", "tests/test_inversion.py", "tests/test_gravity.py tests/test_mag.py"]
 
     steps:

environment-dev.yml

@@ -6,14 +6,15 @@ dependencies:
   - cvxopt
   - deepdish
   - discretize
+  - empymod
   - ipywidgets>=0.6.0
   - jupyter
   - matplotlib=3.4.3
   - Pillow
   - pip
   - pymatsolver
   - requests
-  - SimPEG==0.19.0
+  - SimPEG==0.23.0
   - pytest
   - flake8
   - black

geoscilabs/dcip/DCIP_overburden_PseudoSection.py

@@ -25,14 +25,16 @@
 
 from discretize import TensorMesh
 
-from SimPEG import maps, SolverLU, utils
-from SimPEG.utils import extract_core_mesh
-from SimPEG.electromagnetics.static import resistivity as DC
-from SimPEG.electromagnetics.static import induced_polarization as IP
-from pymatsolver import Pardiso
+from simpeg import maps, SolverLU, utils
+from simpeg.utils import extract_core_mesh
+from simpeg.electromagnetics.static import resistivity as DC
+from simpeg.electromagnetics.static import induced_polarization as IP
+from simpeg.utils.solver_utils import get_default_solver
 
 from ..base import widgetify
 
+Solver = get_default_solver()
+
 # Mesh, sigmaMap can be globals global
 npad = 12
 growrate = 2.0
@@ -179,13 +181,13 @@ def model_fields(A, B, mtrue, mhalf, mair, mover, whichprimary="overburden"):
         survey = DC.Survey([src])
         # Create three simulations so the fields object is accurate
         sim_primary = DC.Simulation2DCellCentered(
-            mesh, survey=survey, sigmaMap=mapping, solver=Pardiso
+            mesh, survey=survey, sigmaMap=mapping, solver=Solver
         )
         sim_total = DC.Simulation2DCellCentered(
-            mesh, survey=survey, sigmaMap=mapping, solver=Pardiso
+            mesh, survey=survey, sigmaMap=mapping, solver=Solver
         )
         sim_air = DC.Simulation2DCellCentered(
-            mesh, survey=survey, sigmaMap=mapping, solver=Pardiso
+            mesh, survey=survey, sigmaMap=mapping, solver=Solver
         )
 
         if whichprimary == "air":
@@ -302,7 +304,7 @@ def getSensitivity(survey, A, B, M, N, model):
 
     survey = DC.Survey([src])
     problem = DC.Simulation2DCellCentered(
-        mesh, sigmaMap=mapping, solver=Pardiso, survey=survey
+        mesh, sigmaMap=mapping, solver=Solver, survey=survey
     )
 
     J = problem.getJ(model)[0]
@@ -458,7 +460,7 @@ def DC2Dsimulation(mtrue, flag="PoleDipole", nmax=8):
 
     survey = DC.Survey(txList)
     simulation = DC.Simulation2DCellCentered(
-        mesh, sigmaMap=mapping, survey=survey, solver=Pardiso
+        mesh, sigmaMap=mapping, survey=survey, solver=Solver
     )
 
     return simulation, xzlocs
@@ -548,7 +550,7 @@ def IP2Dsimulation(miptrue, sigmadc, flag="PoleDipole", nmax=8):
         sigma=sigmadc,
         etaMap=maps.IdentityMap(mesh),
         survey=survey,
-        solver=Pardiso,
+        solver=Solver,
     )
 
     return simulation, xzlocs

geoscilabs/dcip/DCLayers.py

@@ -9,11 +9,13 @@
 
 from discretize import TensorMesh
 
-from SimPEG import maps, utils
-from pymatsolver import Pardiso
+from simpeg import maps, utils
+from simpeg.utils.solver_utils import get_default_solver
 
 from ..base import widgetify
 
+Solver = get_default_solver()
+
 rcParams["font.size"] = 16
 
 # Mesh parameters
@@ -153,7 +155,7 @@ def rho_a(VM, VN, A, B, M, N):
 
 def solve_2D_potentials(rho1, rho2, h, A, B):
     """
-    Here we solve the 2D DC problem for potentials (using SimPEG Mesg Class)
+    Here we solve the 2D DC problem for potentials (using simpeg Mesg Class)
     """
     sigma = 1.0 / rho2 * np.ones(mesh.nC)
     sigma[mesh.gridCC[:, 1] >= -h] = 1.0 / rho1  # since the model is 2D
@@ -172,7 +174,7 @@ def solve_2D_potentials(rho1, rho2, h, A, B):
         * utils.sdiag(1.0 / (mesh.dim * mesh.aveF2CC.T * (1.0 / sigma)))
         * mesh.cell_gradient
     )
-    Ainv = Pardiso(A)
+    Ainv = Solver(A)
 
     V = Ainv * q
     return V

geoscilabs/dcip/DCWidgetPlate2_5D.py

@@ -12,14 +12,16 @@
 from matplotlib.path import Path
 import matplotlib.patches as patches
 
-from pymatsolver import Pardiso
+from simpeg.utils.solver_utils import get_default_solver
 
 from discretize import TensorMesh
-from SimPEG import maps, SolverLU, utils
-from SimPEG.utils import extract_core_mesh
-from SimPEG.electromagnetics.static import resistivity as DC
+from simpeg import maps, SolverLU, utils
+from simpeg.utils import extract_core_mesh
+from simpeg.electromagnetics.static import resistivity as DC
 from ..base import widgetify
 
+Solver = get_default_solver()
+
 # Mesh, mapping can be globals global
 npad = 15
 growrate = 2.0
@@ -88,10 +90,10 @@ def plate_fields(A, B, dx, dz, xc, zc, rotAng, sigplate, sighalf):
             src = DC.sources.Dipole([], np.r_[A, 0.0], np.r_[B, 0.0])
         survey = DC.survey.Survey([src])
         problem = DC.Simulation2DCellCentered(
-            mesh, survey=survey, sigmaMap=mapping, solver=Pardiso
+            mesh, survey=survey, sigmaMap=mapping, solver=Solver
         )
         problem_prim = DC.Simulation2DCellCentered(
-            mesh, survey=survey, sigmaMap=mapping, solver=Pardiso
+            mesh, survey=survey, sigmaMap=mapping, solver=Solver
         )
 
         total_field = problem.fields(mtrue)
@@ -272,7 +274,7 @@ def sumPlateCharges(xc, zc, dx, dz, rotAng, qSecondary):
     return qPosSum, qNegSum, qPosAvgLoc, qNegAvgLoc
 
 
-# The only thing we need to make it work is a 2.5D field object in SimPEG
+# The only thing we need to make it work is a 2.5D field object in simpeg
 
 
 def getSensitivity(survey, A, B, M, N, model):
@@ -292,7 +294,7 @@ def getSensitivity(survey, A, B, M, N, model):
 
     survey = DC.Survey([src])
     sim = DC.Simulation3DCellCentered(
-        mesh, sigmaMap=mapping, solver=Pardiso, survey=survey
+        mesh, sigmaMap=mapping, solver=Solver, survey=survey
     )
     J = sim.getJ(model)[0]
 

geoscilabs/dcip/DCWidgetPlate_2D.py

@@ -11,15 +11,17 @@
 from matplotlib.ticker import LogFormatter
 from matplotlib.path import Path
 import matplotlib.patches as patches
-from pymatsolver import Pardiso
+from simpeg.utils.solver_utils import get_default_solver
 
 from discretize import TensorMesh
-from SimPEG import maps, SolverLU, utils
-from SimPEG.utils import extract_core_mesh
-from SimPEG.electromagnetics.static import resistivity as DC
+from simpeg import maps, SolverLU, utils
+from simpeg.utils import extract_core_mesh
+from simpeg.electromagnetics.static import resistivity as DC
 
 from ..base import widgetify
 
+Solver = get_default_solver()
+
 # Mesh, mapping can be globals global
 npad = 15
 growrate = 2.0
@@ -91,10 +93,10 @@ def plate_fields(A, B, dx, dz, xc, zc, rotAng, sigplate, sighalf):
         survey = DC.survey.Survey([src])
 
         problem = DC.Simulation3DCellCentered(
-            mesh, survey=survey, sigmaMap=mapping, solver=Pardiso
+            mesh, survey=survey, sigmaMap=mapping, solver=Solver
         )
         problem_prim = DC.Simulation3DCellCentered(
-            mesh, survey=survey, sigmaMap=mapping, solver=Pardiso
+            mesh, survey=survey, sigmaMap=mapping, solver=Solver
         )
 
         primary_field = problem_prim.fields(mhalf)
@@ -293,7 +295,7 @@ def getSensitivity(survey, A, B, M, N, model):
 
     survey = DC.Survey([src])
     sim = DC.Simulation3DCellCentered(
-        mesh, sigmaMap=mapping, solver=Pardiso, survey=survey
+        mesh, sigmaMap=mapping, solver=Solver, survey=survey
     )
     J = sim.getJ(model)[0]
 

geoscilabs/dcip/DCWidgetResLayer2D.py

@@ -1,8 +1,8 @@
 from discretize import TensorMesh
-from SimPEG import maps, SolverLU, utils
-from SimPEG.utils import extract_core_mesh
+from simpeg import maps, SolverLU, utils
+from simpeg.utils import extract_core_mesh
 import numpy as np
-from SimPEG.electromagnetics.static import resistivity as DC
+from simpeg.electromagnetics.static import resistivity as DC
 import matplotlib
 import matplotlib.pyplot as plt
 import matplotlib.pylab as pylab
@@ -11,12 +11,14 @@
 import matplotlib.patches as patches
 from scipy.constants import epsilon_0
 import copy
-from pymatsolver import Pardiso
+from simpeg.utils.solver_utils import get_default_solver
 
 from ipywidgets import interact, IntSlider, FloatSlider, FloatText, ToggleButtons
 
 from ..base import widgetify
 
+Solver = get_default_solver()
+
 # Mesh, sigmaMap can be globals global
 npad = 15
 growrate = 2.0
@@ -94,10 +96,10 @@ def model_fields(A, B, zcLayer, dzLayer, xc, zc, r, sigLayer, sigTarget, sigHalf
         survey = DC.Survey([src])
 
         problem = DC.Simulation3DCellCentered(
-            mesh, sigmaMap=sigmaMap, solver=Pardiso, survey=survey
+            mesh, sigmaMap=sigmaMap, solver=Solver, survey=survey
         )
         problem_prim = DC.Simulation3DCellCentered(
-            mesh, sigmaMap=sigmaMap, solver=Pardiso, survey=survey
+            mesh, sigmaMap=sigmaMap, solver=Solver, survey=survey
         )
 
         primary_field = problem_prim.fields(mhalf)
@@ -335,7 +337,7 @@ def getSensitivity(survey, A, B, M, N, model):
 
     survey = DC.survey.Survey([src])
     problem = DC.Simulation3DCellCentered(
-        mesh, survey=survey, sigmaMap=sigmaMap, solver=Pardiso
+        mesh, survey=survey, sigmaMap=sigmaMap, solver=Solver
     )
     J = problem.getJ(model)
 

geoscilabs/dcip/DCWidgetResLayer2_5D.py

@@ -10,17 +10,19 @@
 from matplotlib.path import Path
 import matplotlib.patches as patches
 
-from SimPEG import maps, SolverLU, utils
-from SimPEG.utils import extract_core_mesh
-from SimPEG.electromagnetics.static import resistivity as DC
-from pymatsolver import Pardiso
+from simpeg import maps, SolverLU, utils
+from simpeg.utils import extract_core_mesh
+from simpeg.electromagnetics.static import resistivity as DC
+from simpeg.utils.solver_utils import get_default_solver
 
 from discretize import TensorMesh
 
 from ipywidgets import interact, IntSlider, FloatSlider, FloatText, ToggleButtons
 
 from ..base import widgetify
 
+Solver = get_default_solver()
+
 # Mesh, sigmaMap can be globals global
 npad = 15
 growrate = 2.0
@@ -99,11 +101,11 @@ def model_fields(A, B, zcLayer, dzLayer, xc, zc, r, sigLayer, sigTarget, sigHalf
             src = DC.sources.Dipole([], np.r_[A, 0.0], np.r_[B, 0.0])
         survey = DC.Survey([src])
         sim = DC.Simulation2DCellCentered(
-            mesh, survey=survey, sigmaMap=mapping, solver=Pardiso
+            mesh, survey=survey, sigmaMap=mapping, solver=Solver
         )
         total_field = sim.fields(mtrue)
         sim_prim = DC.Simulation2DCellCentered(
-            mesh, survey=survey, sigmaMap=mapping, solver=Pardiso
+            mesh, survey=survey, sigmaMap=mapping, solver=Solver
         )
         primary_field = sim_prim.fields(mhalf)
 
@@ -331,7 +333,7 @@ def sumCylinderCharges(xc, zc, r, qSecondary):
     return qPosSum, qNegSum, qPosAvgLoc, qNegAvgLoc
 
 
-# The only thing we need to make it work is a 2.5D field object in SimPEG
+# The only thing we need to make it work is a 2.5D field object in simpeg
 
 
 def getSensitivity(survey, A, B, M, N, model):
@@ -347,7 +349,7 @@ def getSensitivity(survey, A, B, M, N, model):
 
     Src = DC.Survey([src])
     sim = DC.Simulation2DCellCentered(
-        mesh, survey=Src, sigmaMap=mapping, solver=Pardiso
+        mesh, survey=Src, sigmaMap=mapping, solver=Solver
     )
     J = sim.getJ(model)
 

geoscilabs/dcip/DCWidget_Overburden_2_5D.py

@@ -21,14 +21,16 @@
 )
 
 from discretize import TensorMesh
-from SimPEG import maps, utils
-from SimPEG.utils import extract_core_mesh
-from SimPEG.electromagnetics.static import resistivity as DC
+from simpeg import maps, utils
+from simpeg.utils import extract_core_mesh
+from simpeg.electromagnetics.static import resistivity as DC
 
-from pymatsolver import Pardiso
+from simpeg.utils.solver_utils import get_default_solver
 
 from ..base import widgetify
 
+Solver = get_default_solver()
+
 # Mesh, sigmaMap can be globals global
 npad = 12
 growrate = 2.0
@@ -165,13 +167,13 @@ def model_fields(A, B, mtrue, mhalf, mair, mover, whichprimary="air"):
             src = DC.sources.Dipole([], np.r_[A, surfaceA], np.r_[B, surfaceB])
         survey = DC.survey.Survey([src])
         problem = DC.Simulation2DCellCentered(
-            mesh, survey=survey, sigmaMap=mapping, solver=Pardiso
+            mesh, survey=survey, sigmaMap=mapping, solver=Solver
         )
         problem_prim = DC.Simulation2DCellCentered(
-            mesh, survey=survey, sigmaMap=mapping, solver=Pardiso
+            mesh, survey=survey, sigmaMap=mapping, solver=Solver
         )
         problem_air = DC.Simulation2DCellCentered(
-            mesh, survey=survey, sigmaMap=mapping, solver=Pardiso
+            mesh, survey=survey, sigmaMap=mapping, solver=Solver
         )
 
         if whichprimary == "air":
@@ -307,7 +309,7 @@ def getSensitivity(survey, A, B, M, N, model):
 
     Src = DC.Survey([src])
     sim = DC.Simulation2DCellCentered(
-        mesh, survey=Src, sigmaMap=mapping, solver=Pardiso
+        mesh, survey=Src, sigmaMap=mapping, solver=Solver
     )
     J = sim.getJ(model)[0]