simpeg/pymatsolver

Solve matrix equations in python.

SimPEG

Simulation and Parameter Estimation in Geophysics - A python package for simulation and gradient based parameter estimation in the context of geophysical applications.

The vision is to create a package for finite volume simulation with applications to geophysical imaging and subsurface flow. To enable the understanding of the many different components, this package has the following features:

• modular with respect to the spacial discretization, optimization routine, and geophysical problem
• built with the inverse problem in mind
• provides a framework for geophysical and hydrogeologic problems
• supports 1D, 2D and 3D problems
• designed for large-scale inversions

You are welcome to join our forum and engage with people who use and develop SimPEG at: https://simpeg.discourse.group/.

Weekly meetings are open to all. They are generally held on Wednesdays at 10:30am PDT. Please see the calendar (GCAL, ICAL) for information on the next meeting.

Overview Video

Working towards all the Geophysics, but Backwards - SciPy 2016

Citing SimPEG

There is a paper about SimPEG!

Cockett, R., Kang, S., Heagy, L. J., Pidlisecky, A., & Oldenburg, D. W. (2015). SimPEG: An open source framework for simulation and gradient based parameter estimation in geophysical applications. Computers & Geosciences.

BibTex:

@article{cockett2015simpeg,
  title={SimPEG: An open source framework for simulation and gradient based parameter estimation in geophysical applications},
  author={Cockett, Rowan and Kang, Seogi and Heagy, Lindsey J and Pidlisecky, Adam and Oldenburg, Douglas W},
  journal={Computers \& Geosciences},
  year={2015},
  publisher={Elsevier}
}

Electromagnetics

If you are using the electromagnetics module of SimPEG, please cite:

Lindsey J. Heagy, Rowan Cockett, Seogi Kang, Gudni K. Rosenkjaer, Douglas W. Oldenburg, A framework for simulation and inversion in electromagnetics, Computers & Geosciences, Volume 107, 2017, Pages 1-19, ISSN 0098-3004, http://dx.doi.org/10.1016/j.cageo.2017.06.018.

BibTex:

@article{heagy2017,
    title= "A framework for simulation and inversion in electromagnetics",
    author= "Lindsey J. Heagy and Rowan Cockett and Seogi Kang and Gudni K. Rosenkjaer and Douglas W. Oldenburg",
    journal= "Computers & Geosciences",
    volume = "107",
    pages = "1 - 19",
    year = "2017",
    note = "",
    issn = "0098-3004",
    doi = "http://dx.doi.org/10.1016/j.cageo.2017.06.018"
}

Questions

If you have a question regarding a specific use of SimPEG, the fastest way to get a response is by posting on our Discourse discussion forum: https://simpeg.discourse.group/. Alternatively, if you prefer real-time chat, you can join our Mattermost Team at https://mattermost.softwareunderground.org/simpeg. Please do not create an issue to ask a question.

Meetings

SimPEG hosts weekly meetings for users to interact with each other, for developers to discuss upcoming changes to the code base, and for discussing topics related to geophysics in general. Currently our meetings are held every Wednesday, alternating between a mornings (10:30 am pacific time) and afternoons (3:00 pm pacific time) on even numbered Wednesdays. Find more info on our Mattermost.

Links

Website: https://simpeg.xyz

Forums: https://simpeg.discourse.group/

Mattermost (real time chat): https://mattermost.softwareunderground.org/simpeg

Documentation: https://docs.simpeg.xyz

Code: https://github.com/simpeg/simpeg

Tests: https://dev.azure.com/simpeg/simpeg/_build

Bugs & Issues: https://github.com/simpeg/simpeg/issues

Contributing

We always welcome contributions towards SimPEG whether they are adding new code, suggesting improvements to existing codes, identifying bugs, providing examples, or anything that will improve SimPEG. Please checkout the contributing guide for more information on how to contribute.

pymatsolver's Issues

MKL bug for reordering.

https://software.intel.com/en-us/forums/intel-math-kernel-library/topic/657934

`is_summetric` with `SolverLU`

Copying over from Slack:

Dieter:

Alfredo Sanchez brought a problem to my attention, and I think it is a bug: It is about the following example: https://docs.simpeg.xyz/content/tutorials/05-dcr/plot_inv_2_dcr2d.html Running it with PARDISO works fine. However, running it with SolverLU results in:

Lindsey:

hmm, yes, I think there is an inconsistency between the kwargs that splu / pardiso expect (the is_symmetric flag in Pardiso changes how the matrix is factored, whereas it is not an option in splu https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.splu.html)
This is something we should fix in the pymatsolver layer

Sparse RHS

Mumps should take a sparse vector for the RHS.

https://github.com/JuliaSparse/MUMPS.jl/blob/master/src/MUMPS_4.10.0/SCILAB/examples/sparseRHS_example.sce

instantiate solvers without setting a system matrix?

Motivation

In developing a Simulation class in SimPEG (simpeg/simpeg#672), I want to be able to readily serialize and deserialize an instance of a Simulation, so we are injecting properties throughout, which gives us serialization, deserialization and validation (thanks @fwkoch!!). Right now, the solvers are attached as the simulation as a class, and the solver_opts are attached to the simulation as a dict. We then instantiate a solver as needed in SimPEG.

State the problem

There are a couple potential snags with this approach

• if your solver_opts are not valid, we won't find out until we try and use the solver
• It is not readily clear how to serialize the solver class (a Class property is one option seequent/properties#163, we could also think about serializing / deserializing the name of the class as a string)

Another approach?

We could instantiate the solver with its options and then call it to create Ainv. This would solve the serialization problem and also allow immediate validation of the solver options on its creation.

import pymatsolver
solver = pymatsolver.Pardiso(**opts)
Ainv = solver(A)

so the base class might look like:

class Base(properties.HasProperties):

    check_accuracy = properties.Bool(
        "check the accuracy of the solve?",
        default = False
    )

    accuracy_tol = properties.Float(
        "tolerance on the accuracy of the solver",
        default=1e-6
    )

    def __init__(self, **kwargs):
        super(Base, self).__init__(**kwargs)

    def __call__(self, A):
        self.A = A.tocsr()

There are a couple impacts of this within SimPEG. The one that comes to mind is:

• right now in the inversion, if a solver is not set for the data misfit, we use the same one as on the simulation, here we would likely want to make a copy of the instance, rather than re-use the same instance (so that factors are not cleared)

tests failing / not running

It looks like we need to go in and update our github workflows for treating as seen in #44

pymatsolver docs

The pymatsolver docs do not appear to be building? (https://readthedocs.org/pymatsolver). @rowanc1: is this a leftover from moving it to the SimPEG org?

move to SimPEG org?

Since we are moving towards having SimPEG rely on pymatsolver for all of the solvers (see simpeg/simpeg#672), would it make sense to bring this repo over to the SimPEG org? We can make sure the docs are linked, pypi deploys automated, etc.

PyLops and spgl1

Would spgl1 be something for pymatsolver? Might be less memory-intensive than PARDISO; we would have to check runtime.

Improve warning when not enough memory (output is all zero)

Issue carried over from discussion on Slack with @lheagy, @grosenkj (also pinging @rowanc1).

Using Pardiso (it might or might not affect other solvers as well, I don't know), my output was all zeros (array of correct size, all with 0+0j). No error message printed, everything run smooth as usual. (You only note it when matplotlib complains about plotting it on log-scale "Data has no positive values, and therefore cannot be log-scaled" 😄)

Turns out that there was not enough memory. Making the model smaller resolved it. However, while running the notebook the memory usage never went up, so it didn't occur to me that it could be a memory-related issue. I assume Pardiso checks how much memory is available with the model size, and doesn't even start if there is not sufficient.

The obvious solution might be to check output of Pardiso and pass it along. However, that might not be that straight forward or the output might be cryptic. In that case one could make a check at the end if the result is all 0's, and print a warning (related with some hints, as mkl-installation or mesh-size/memory-requirements).

PARDISO compatibility issue with MKL 2021

It seems that the current wrapper for PARDISO doesn't work for MKL 2021.

BUG? Pydiso not a requirement

Not sure if this is a bug or by design, but pydiso is not a requirement in the setup and therefore needs a separate install of pydiso to use the Paradiso solver. It would be nice to install pymatsolver and get the Paradiso solver with it.

issues with python setup.py install

Hello,

first, thank you for making this package available. I managed to get the tests running in the pymatsolver directory without problems, using the instructions in the README,rst (not using anaconda). There are, however, some issues, when trying to use the python setup.py install command.

1. This does not build/install the MUMPS wrappers. (I could make and copy them manually, but...).
2. The tests directory is installed into .../python/dist-packages directly (on ubuntu, this is site-packages) and not into .../python/dist-packages/pymatsolver.

Use newer, non-mkl, version of PARDISO

To follow up on the Slack-discussion and to not forget about it.

For more info: https://www.pardiso-project.org

From the website:

Important: Please note that the Intel MKL version of PARDISO is based on our version from 2006 and that a lot of new features and improvements of PARDISO are not available in the Intel MKL library.

In their examples this means a speed-up factor of anything between 1.7 and 26, see the comparison charts in the feature-list.

Error running setup.py

When running setup.py, finish with:
error: file '[...]/pymatsolver/pymatsolver/Mumps/MumpsInterface.so' does not exist

Transfer to the SimPEG organization.

This should allow us to collaborate on things.

@jcapriot @sgkang @bsmithyman