diff --git a/docs/_static/style.css b/docs/_static/style.css new file mode 100644 index 00000000..e3b0f28e --- /dev/null +++ b/docs/_static/style.css @@ -0,0 +1,33 @@ +/* "emsig"-colours + blue : #447fc1 + blue-lighter : #61a2d8 + blue-darker : #285985 + red : #ef413d +*/ + +/* Button colours */ + +.btn-info { + font-weight: bold; + background: #447fc1; + border-color: #447fc1; +} + +.btn-info:hover { + background: #ef413d; + border-color: #ef413d; +} + +/* Card headers */ +.card-header { + font-size: var(--pst-font-size-h5); + font-weight: bold; + padding: 1rem 2rem; + text-align: center; +} + +/* fa stuff */ +.fa { + margin: 0 1rem 0 0; + color: #ef413d; +} diff --git a/docs/api/index.rst b/docs/api/index.rst index 604d4c2e..23bcad2e 100644 --- a/docs/api/index.rst +++ b/docs/api/index.rst @@ -12,7 +12,7 @@ API reference .. module:: emg3d .. toctree:: - :maxdepth: 1 + :maxdepth: 2 :hidden: core @@ -29,11 +29,34 @@ API reference time utils -Direct links to the main user-facing functions and classes: -- Grid: :class:`emg3d.meshes.TensorMesh` -- Model: :class:`emg3d.models.Model` -- Survey: :class:`emg3d.surveys.Survey` -- Simulation: :class:`emg3d.simulations.Simulation` -- Solver: :func:`emg3d.solver.solve` -- All sources and receivers are in :mod:`emg3d.electrodes` +.. panels:: + :container: container-lg pb-1 + :column: col-lg-12 p-2 + + Grid: :class:`emg3d.meshes.TensorMesh` + + --- + :column: col-lg-12 p-2 + + Model: :class:`emg3d.models.Model` + + --- + :column: col-lg-12 p-2 + + Survey: :class:`emg3d.surveys.Survey` + + --- + :column: col-lg-12 p-2 + + Simulation: :class:`emg3d.simulations.Simulation` + + --- + :column: col-lg-12 p-2 + + Solver: :func:`emg3d.solver.solve` + + --- + :column: col-lg-12 p-2 + + All sources and receivers are in :mod:`emg3d.electrodes` diff --git a/docs/conf.py b/docs/conf.py index 612e6932..5b237fa0 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -7,6 +7,7 @@ extensions = [ # 'sphinx.ext.autodoc', 'numpydoc', + 'sphinx_panels', 'sphinx.ext.intersphinx', # 'sphinx.ext.autosummary', 'sphinx.ext.mathjax', @@ -18,6 +19,7 @@ 'IPython.sphinxext.ipython_console_highlighting', 'IPython.sphinxext.ipython_directive', ] +panels_add_bootstrap_css = False autosummary_generate = True add_module_names = True add_function_parentheses = False @@ -55,7 +57,7 @@ # General information about the project. project = 'emg3d' -author = 'The EMSiG community' +author = 'The emsig community' copyright = f'2018-{time.strftime("%Y")}, {author}' # |version| and |today| tags (|release|-tag is not used). @@ -78,8 +80,8 @@ html_theme_options = { "github_url": "https://github.com/emsig/emg3d", "external_links": [ - {"name": "Gallery", "url": "https://dev1.emsig.xyz/gallery/gallery"}, - {"name": "EMSiG", "url": "https://emsig.xyz"}, + {"name": "Gallery", "url": "https://emsig.xyz/emg3d-gallery/gallery"}, + {"name": "emsig", "url": "https://emsig.xyz"}, ], # "use_edit_page_button": True, } @@ -94,6 +96,11 @@ html_use_modindex = True html_file_suffix = '.html' htmlhelp_basename = 'emg3d' +html_css_files = [ + "style.css", + "https://cdnjs.cloudflare.com/ajax/libs/font-awesome/4.7.0/" + + "css/font-awesome.min.css" +] # ==== 4. linkcheck ==== diff --git a/docs/dev/contributing.rst b/docs/dev/contributing.rst index 12cb8695..85dbe643 100644 --- a/docs/dev/contributing.rst +++ b/docs/dev/contributing.rst @@ -1,28 +1,28 @@ +.. _dev-contributing: + Contributing ============ -.. todo:: +emg3d is a community code, and contributions of any kind are welcomed: From +typos in the documentation to additional examples, reporting or fixing bugs in +the code, all the way to new suggestion or implementing new features. + +Good places to get started is to browse the existing issues, check out the +roadmap, or have a look at any open PR: - Rework this for version 1. +- `Issues `_; +- `Roadmap-project `_; +- `PR's `_. -New contributions, bug reports, or any kind of feedback is always welcomed! -Have a look at the `Roadmap-project -`_ to get an idea of things that -could be implemented. The GitHub `issues -`_ and -`PR's `_ are also a good starting -point. The best way for interaction is at https://github.com/emsig or by -joining the `Slack channel `_ «em-x-d» of SimPEG. If -you prefer to get in touch outside of GitHub/Slack use the contact form on -https://werthmuller.org. +There are various different ways to get in touch, see +`emsig.xyz#contributing-contact `_. -To install emg3d from source, you can download the latest version from GitHub -and install it in your python distribution via: +If you think about making changes to the code and contribute code have a look +at :doc:`tests` -.. code-block:: console - make install +.. note:: -Please make sure your code follows the pep8-guidelines by using, for instance, -the python module ``flake8``, and also that your code is covered with -appropriate tests. Just get in touch if you have any doubts. + The community is tiny so far, so there are no former protocols (such as + CoC) in place (yet). Please get in touch if you think it is about time to + implement those. diff --git a/docs/dev/index.rst b/docs/dev/index.rst index 009c9541..cbf82d67 100644 --- a/docs/dev/index.rst +++ b/docs/dev/index.rst @@ -11,9 +11,46 @@ Development .. toctree:: :maxdepth: 2 + :hidden: contributing - tests_benchmarks + tests maintenance - making_a_release + release +.. panels:: + :container: container-lg pb-2 text-center + + :fa:`code,fa-4x` + + .. link-button:: contributing + :type: ref + :text: Contributing + :classes: btn-link stretched-link font-weight-bold + + --- + + :fa:`fas fa-hourglass-half,fa-4x` + + .. link-button:: tests + :type: ref + :text: Tests & Benchmarks + :classes: btn-link stretched-link font-weight-bold + + --- + + :fa:`cogs,fa-4x` + + .. link-button:: maintenance + :type: ref + :text: Maintenance + :classes: btn-link stretched-link font-weight-bold + + --- + + :fa:`rocket,fa-4x` + + .. link-button:: release + :type: ref + :text: Making a release + :classes: btn-link stretched-link font-weight-bold diff --git a/docs/dev/maintenance.rst b/docs/dev/maintenance.rst index 1c1c0f84..a24ad9e8 100644 --- a/docs/dev/maintenance.rst +++ b/docs/dev/maintenance.rst @@ -1,13 +1,15 @@ Maintenance =========== -.. todo:: - - Rework this for version 1. +Status reports and other tools to have the checks all in one place, for quick +QC. Quick overview / QC ------------------- +All possible badges of the CI chain. Definitely check this before making a +release. + - .. image:: https://github.com/emsig/emg3d/actions/workflows/linux.yml/badge.svg :target: https://github.com/emsig/emg3d/actions/workflows/linux.yml :alt: GitHub Actions linux @@ -57,10 +59,12 @@ Quick overview / QC Info from ReadTheDocs --------------------- +To check the environment in which the documentation was built. + .. ipython:: In [1]: import emg3d ...: emg3d.Report( ...: ['sphinx', 'numpydoc', 'ipykernel', 'sphinx_numfig', - ...: 'sphinx_automodapi', 'pydata_sphinx_theme'] + ...: 'sphinx_panels', 'sphinx_automodapi', 'pydata_sphinx_theme'] ...: ) diff --git a/docs/dev/making_a_release.rst b/docs/dev/release.rst similarity index 56% rename from docs/dev/making_a_release.rst rename to docs/dev/release.rst index 49f83493..5dd2df7d 100644 --- a/docs/dev/making_a_release.rst +++ b/docs/dev/release.rst @@ -1,9 +1,9 @@ Making a release ================ -.. todo:: - - Rework this for version 1. +Making a release is by now straight forward. Creating a «New Release» on GitHub +will trigger all the required things. However, there are still a few things to +do beforehand and a few things to check afterwards. 1. Update ``CHANGELOG.rst``. @@ -11,9 +11,12 @@ Making a release 3. Tagging it on GitHub will automatically deploy it to PyPi, which in turn will create a PR for the conda-forge `feedstock - `_. Merge that PR. + `_. The PR will be + automatically merged. (Note: if the dependencies change or the minimum + Python version or other installation things then the feedstock has to be + updated manually!) -4. Check that: +4. After releasing it, check that: - `PyPi `_ deployed; - `conda-forge `_ deployed; @@ -21,55 +24,54 @@ Making a release - `emg3d.emsig.xyz `_ created a tagged version. -Useful things -------------- - -- If there were changes to README, check it with:: - - python setup.py --long-description | rst2html.py --no-raw > index.html - -- If unsure, test it first on testpypi (requires ~/.pypirc):: - - ~/anaconda3/bin/twine upload dist/* -r testpypi - -- If unsure, test the test-pypi for conda if the skeleton builds:: - - conda skeleton pypi --pypi-url https://test.pypi.io/pypi/ emg3d - -- If it fails, you might have to install ``python3-setuptools``:: - - sudo apt install python3-setuptools - - -CI --- +CI automatic and manual bits +---------------------------- -Automatic bits -`````````````` +Automatic +````````` -- Testing on Github Actions includes: +- Testing on `Github Actions `_ + includes: - - Tests using ``pytest`` + - Tests using ``pytest`` (Linux, MacOS, Windows) - Linting / code style with ``pytest-flake8`` - - Ensure all http(s)-links work (``sphinx linkcheck``) + - Ensure all http(s)-links work (``sphinx -b linkcheck``) - Line-coverage with ``pytest-cov`` on `Coveralls `_ - Code-quality on `Codacy - `_ + `_ - Manual on `ReadTheDocs `_ - DOI minting on `Zenodo `_ -Manual things -````````````` +Manual +`````` - Benchmarks with `Airspeed Velocity `_ (``asv``) - Gallery in `emg3d-gallery `_ (``sphinx-gallery``) -Automatically deploys if tagged -``````````````````````````````` -- `PyPi `_ -- `conda -c conda-forge `_ +Useful things +------------- + +The following info was from the time when we still manually deployed. Now +every push to main triggers a test to Test-PyPI, so things can be verified +there. However, these hints my still be useful at some point. + +- If there were changes to README, check it with:: + + python setup.py --long-description | rst2html.py --no-raw > index.html + +- If unsure about something, test it first on testpypi (requires ~/.pypirc):: + + ~/anaconda3/bin/twine upload dist/* -r testpypi + +- If unsure, test the test-pypi for conda if the skeleton builds:: + + conda skeleton pypi --pypi-url https://test.pypi.io/pypi/ emg3d + +- If it fails, you might have to install ``python3-setuptools``:: + + sudo apt install python3-setuptools diff --git a/docs/dev/tests.rst b/docs/dev/tests.rst new file mode 100644 index 00000000..914411b4 --- /dev/null +++ b/docs/dev/tests.rst @@ -0,0 +1,497 @@ +Tests and Benchmarks +==================== + +If you ended up here you probably think about fixing some bugs or contributing +some code. **Awesome!** Just open a PR, and we will guide you through the +process. The following section contains some more detailed information of the +continues integration (CI) procedure we follow. In the end, each commit has to +pass them before it can be merged into the main branch on GitHub. + + +The first step to develop code is to clone the GitHub repo locally: + +.. code-block:: console + + git clone git@github.com:emsig/emg3d.git + +All requirements for the dev-toolchain are collected in the +``requirements-dev.txt`` file, so you can install them all by running + +.. code-block:: console + + pip install -r requirements_dev.txt + +With this you have all the basic tools to run the tests, lint your code, build +the documentation, and so on. + +Continuous Integration +---------------------- + +The CI elements are: + +1. Linting: ``flake8`` +2. Tests: ``pytest`` +3. Code coverage: ``coveralls`` +4. Link checks: ``sphinx`` +5. Code quality: ``codacy`` +6. Documentation: ``sphinx`` +7. Benchmarks: ``asv`` + + +(1) to (6) are run automatically through GitHub actions when committing changes +to GitHub. Any code change should pass these tests. Additionally, it is crucial +that new code comes with the appropriate tests and documentation, and if +applicable also with the appropriate benchmarks. However, you do not need any +of that to start a PR - everything can go step-by-step! + +Many of the tests are set up in the Makefile (only tested on Linux): + +- To install the current branch in editable mode: + + .. code-block:: console + + make install + +- To check linting: + + .. code-block:: console + + make flake8 + +- To run pytest: + + .. code-block:: console + + make pytest + +- To build the documentation: + + .. code-block:: console + + make doc + +- Or to list all the possibilities, simply run: + + .. code-block:: console + + make + +There is also a benchmark suite using *airspeed velocity*, located in the +`emsig/emg3d-asv `_-repository. The results +of my machine can be found in the `emsig/emg3d-bench +`_, its rendered version at +`emsig.xyz/emg3d-asv `_. They ensure that we do +not slow than the computation by introducing regressions, particularly when we +make changes to :mod:`emg3d.core` or :mod:`emg3d.solver`. + + +.. _improve-cpu-ram: + +CPU & RAM +--------- + +Here some information if someone is interested in tackling the very core of +emg3d, trying to make it faster or reduce the memory consumption. The multigrid +method is attractive because it shows optimal scaling for both runtime and +memory consumption. Below some insights about what has been tried and what +still could be tried in order to improve the current code. + + +Runtime +~~~~~~~ + +The costliest functions (for big models) are: + + - >90 %: :func:`emg3d.solver.smoothing` (:func:`emg3d.core.gauss_seidel`) + - <5 % each, in decreasing importance: + + - :func:`emg3d.solver.prolongation` + (:class:`emg3d.solver.RegularGridProlongator`) + - :func:`emg3d.solver.residual` (:func:`emg3d.core.amat_x`) + - :func:`emg3d.solver.restriction` + +Example with 262,144 / 2,097,152 cells (``nu_{i,1,c,2}=0,2,1,2``; +``sslsolver=False``; ``semicoarsening=True``; ``linerelaxation=True``): + + - 93.7 / 95.8 % ``smoothing`` + - 3.6 / 2.0 % ``prolongation`` + - 1.9 / 1.9 % ``residual`` + - 0.6 / 0.4 % ``restriction`` + +The rest can be ignored. For small models, the percentage of ``smoothing`` goes +down and of ``prolongation`` and ``restriction`` go up. But then the modeller +is fast anyway. + +:func:`emg3d.core.gauss_seidel` and :func:`emg3d.core.amat_x` are written +in ``numba``; jitting :class:`emg3d.solver.RegularGridProlongator` turned out +to not improve things, and many functions used in the restriction are jitted +too. The costliest functions (RAM- and CPU-wise) are therefore already written +in ``numba``. + +**Any serious attempt to improve the speed will have to tackle the smoothing +itself.** + + +**Things which could be tried** + +- Not much has been tested with the ``numba``-options ``parallel``; ``prange``; + and ``nogil``. +- There might be an additional gain by making :class:`emg3d.meshes.TensorMesh`, + :class:`emg3d.models.Model`, and :class:`emg3d.fields.Field` instances jitted + classes. + +**Things which have been tried** + +- One important aspect of the smoothing part is the memory layout. + :func:`emg3d.core.gauss_seidel` and :func:`emg3d.core.gauss_seidel_x` + are ideal for F-arrays (loop z-y-x, hence slowest to fastest axis). + :func:`emg3d.core.gauss_seidel_y` and + :func:`emg3d.core.gauss_seidel_z`, however, would be optimal for C-arrays. + But copying the arrays to C-order and afterwards back is costlier in most + cases for both CPU and RAM. The one possible and therefore implemented + solution was to swap the loop-order in :func:`emg3d.core.gauss_seidel_y`. +- Restriction and prolongation information could be saved in a dictionary + instead of recomputing it every time. Turns out to be not worth the + trouble. +- Rewrite :class:`emg3d.solver.RegularGridProlongator` as jitted function, but + the iterator approach seems to be better for large grids. + + +Memory +~~~~~~ + +Most of the memory requirement comes from storing the data itself, mainly the +fields (source field, electric field, and residual field) and the model +parameters (resistivity, eta, mu). For a big model, they some up; e.g., almost +3 GB for an isotropic model with 256 x 256 x 256 cells. Anyhow, memory +consumption is pretty low already, and there is probably not much to gain, at +least in the solver part (:mod:`emg3d.core` and :mod:`emg3d.solver`). That +looks different for some of the interpolation and plotting routines, which +could be improved . + + +Benchmark scripts for status quo +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +To test CPU and RAM on your machine, you can use and adjust the following +script. The old notebooks which were used to generate the above figures in the +manual can be found at + +- RAM: `4a_RAM-requirements.ipynb + `_, +- CPU: `4b_Runtime.ipynb + `_. + +.. ipython:: + :verbatim: + + In [1]: import emg3d + ...: import numpy as np + ...: import matplotlib.pyplot as plt + ...: from memory_profiler import memory_usage + + In [2]: def compute(nx): + ...: """Simple computation routine. + ...: + ...: This is the actual model it runs. Adjust this to your needs. + ...: + ...: - Model size is nx * nx * nx, centered around the origin. + ...: - Source is at the origin, x-directed. + ...: - Frequency is 1 Hz. + ...: - Homogenous space of 1 Ohm.m. + ...: + ...: """ + ...: + ...: # Grid + ...: hx = np.ones(nx)*50 + ...: x0 = -nx//2*50 + ...: grid = emg3d.TensorMesh([hx, hx, hx], x0=(x0, x0, x0)) + ...: + ...: # Model and source field + ...: model = emg3d.Model(grid, property_x=1.0) + ...: sfield = emg3d.get_source_field( + ...: grid, source=[0, 0, 0, 0, 0], frequency=1.0) + ...: + ...: # Compute the field + ...: _, inf = emg3d.solve( + ...: model, sfield, verb=0, plain=True, return_info=True) + ...: + ...: return inf['time'] + + In [3]: # Loop over model sizes (adjust to your needs). + ...: nsizes = np.array([32, 48, 64, 96, 128, 192, 256, 384]) + ...: memory = np.zeros(nsizes.shape) + ...: runtime = np.zeros(nsizes.shape) + ...: + ...: # Loop over nx + ...: for i, nx in enumerate(nsizes): + ...: print(f" => {nx}^3 = {nx**3:12,d} cells") + ...: mem, time = memory_usage((compute, (nx, ), {}), retval=True) + ...: memory[i] = max(mem) + ...: runtime[i] = time + ...: + + In [4]: # Plot CPU + ...: plt.figure() + ...: plt.title('Runtime') + ...: plt.loglog(nsizes**3/1e6, runtime, '.-') + ...: plt.xlabel('Number of cells (in millions)') + ...: plt.ylabel('CPU (s)') + ...: plt.axis('equal') + ...: plt.show() + + In [5]: # Plot RAM + ...: plt.figure() + ...: plt.title('Memory') + ...: plt.loglog(nsizes**3/1e6, memory/1e3, '-', zorder=10) + ...: plt.xlabel('Number of cells (in millions)') + ...: plt.ylabel('RAM (GB)') + ...: plt.axis('equal') + ...: plt.show() + + +Scripts for solver investigations +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The non-standard Cholesky solver, :func:`emg3d.core.solve`, does almost all the +work, in the end. Improving the speed of that part only slightly would have a +huge effect overall. Here some notes from some dabbling. + +**Benchmark Tests for Cholesky Solve** + +- `numba`, `numpy`, `scipy`, `lapack` +- Benchmarks: + - small and big + - real and complex valued + +"Givens": + +- Diagonal != 0 +- Diagonal values are large (no pivoting) +- Only diagonal values would be complex + +.. ipython:: + :verbatim: + + In [1]: import numba as nb + ...: import numpy as np + ...: import scipy as sp + ...: from numpy.testing import assert_allclose + ...: from scipy.linalg.lapack import get_lapack_funcs + ...: + ...: _numba_setting = {'nogil': True, 'fastmath': True, 'cache': True} + +Status quo + +.. ipython:: + :verbatim: + + In [1]: @nb.njit(**_numba_setting) + ...: def _emg3d_solve(amat, bvec): + ...: n = len(bvec) + ...: h = np.zeros(1, dtype=amat.dtype)[0] # Pre-allocate + ...: d = 1./amat[0] + ...: + ...: for i in range(1, min(n, 6)): + ...: amat[i] *= d + ...: + ...: for j in range(1, n): + ...: h *= 0. # Reset h + ...: for k in range(max(0, j-5), j): + ...: h += amat[j+5*k]*amat[j+5*k]*amat[6*k] + ...: amat[6*j] -= h + ...: d = 1./amat[6*j] + ...: for i in range(j+1, min(n, j+6)): + ...: h *= 0. # Reset h + ...: for k in range(max(0, i-5), j): + ...: h += amat[i+5*k]*amat[j+5*k]*amat[6*k] + ...: amat[i+5*j] -= h + ...: amat[i+5*j] *= d + ...: + ...: amat[6*(n-1)] = d + ...: for j in range(n-2, -1, -1): + ...: amat[6*j] = 1./amat[6*j] + ...: + ...: for j in range(1, n): + ...: h *= 0. # Reset h + ...: for k in range(max(0, j-5), j): + ...: h += amat[j+5*k]*bvec[k] + ...: bvec[j] -= h + ...: + ...: for j in range(n): + ...: bvec[j] *= amat[6*j] + ...: + ...: for j in range(n-2, -1, -1): + ...: h *= 0. # Reset h + ...: for k in range(j+1, min(n, j+6)): + ...: h += amat[k+5*j]*bvec[k] + ...: bvec[j] -= h + ...: + ...: def emg3d_solve(amat, bvec): + ...: out = bvec.copy() + ...: _emg3d_solve(amat.copy(), out) + ...: return out + + +An alternative + +.. ipython:: + :verbatim: + + In [1]: @nb.njit(**_numba_setting) + ...: def _emg3d_solve2(amat, bvec): + ...: n = len(bvec) + ...: h = np.zeros(1, dtype=amat.dtype)[0] # Pre-allocate + ...: d = 1./amat[0] + ...: + ...: for i in range(1, min(n, 6)): + ...: amat[i] *= d + ...: + ...: for j in range(1, n): + ...: h *= 0. # Reset h + ...: for k in range(max(0, j-5), j): + ...: h += amat[j+5*k]*amat[j+5*k]*amat[6*k] + ...: amat[6*j] -= h + ...: d = 1./amat[6*j] + ...: for i in range(j+1, min(n, j+6)): + ...: h *= 0. # Reset h + ...: for k in range(max(0, i-5), j): + ...: h += amat[i+5*k]*amat[j+5*k]*amat[6*k] + ...: + ...: amat[i+5*j] = d*(amat[i+5*j] - h) + ...: + ...: amat[6*(n-1)] = d + ...: for j in range(n-2, -1, -1): + ...: amat[6*j] = 1./amat[6*j] + ...: + ...: for j in range(1, n): + ...: h *= 0. # Reset h + ...: for k in range(max(0, j-5), j): + ...: h += amat[j+5*k]*bvec[k] + ...: bvec[j] -= h + ...: + ...: for j in range(n): + ...: bvec[j] *= amat[6*j] + ...: + ...: for j in range(n-2, -1, -1): + ...: h *= 0. # Reset h + ...: for k in range(j+1, min(n, j+6)): + ...: h += amat[k+5*j]*bvec[k] + ...: bvec[j] -= h + ...: + ...: + ...: def emg3d_solve2(amat, bvec): + ...: out = bvec.copy() + ...: _emg3d_solve2(amat.copy(), out) + ...: return out + +SciPy and NumPy solvers + +.. ipython:: + :verbatim: + + In [1]: def np_linalg_solve(A, b): + ...: return np.linalg.solve(A, b) + ...: + ...: def sp_linalg_solve(A, b): + ...: out = b.copy() + ...: sp.linalg.solve(A.copy(), out, overwrite_a=True, + ...: overwrite_b=True, check_finite=False) + ...: return out + ...: + ...: def sp_linalg_lu_solve(A, b): + ...: out = b.copy() + ...: lu_and_piv = sp.linalg.lu_factor(A.copy(), overwrite_a=True, + ...: check_finite=False) + ...: xlu = sp.linalg.lu_solve(lu_and_piv, out, overwrite_b=True, + ...: check_finite=False) + ...: return out + ...: + ...: def sp_linalg_cho_solve(A, b): + ...: amat = A.copy() + ...: clow = sp.linalg.cho_factor(amat, lower=True, overwrite_a=True, + ...: check_finite=False) + ...: out = b.copy() + ...: sp.linalg.cho_solve(clow, out, overwrite_b=True, + ...: check_finite=False) + ...: return out + ...: + ...: def sp_linalg_cho_banded(A, b): + ...: amat = A.copy() + ...: c = sp.linalg.cholesky_banded(amat, overwrite_ab=True, + ...: lower=True, check_finite=False) + ...: out = b.copy() + ...: sp.linalg.cho_solve_banded((c, True), out, overwrite_b=True, + ...: check_finite=False) + ...: return out + + +Measuring them. You can get the data at +https://github.com/emsig/data/raw/main/emg3d/benchmarks/CholeskySolveBenchmark.npz + +.. ipython:: + :verbatim: + + In [1]: data = np.load('CholeskySolveBenchmark.npz') + ...: + ...: for cr in ['real', 'cplx']: + ...: for bs in ['small', 'big']: + ...: + ...: print(f"dtype={cr}; size={bs}") + ...: + ...: # Get test data. + ...: amat = data[cr+'_'+bs+'_'+'amat'] + ...: bvec = data[cr+'_'+bs+'_'+'bvec'] + ...: out = data[cr+'_'+bs+'_'+'out'] + ...: + ...: # Re-arrange to full (symmetric) or banded matrix + ...: # for some solvers. + ...: n = bvec.size + ...: A = np.zeros((n, n), dtype=amat.dtype) + ...: Ab = np.zeros((6, n), dtype=amat.dtype) + ...: for i in range(n): + ...: A[i, i] = amat[i*6] + ...: Ab[0, i] = amat[i*6] + ...: for j in range(1, 6): + ...: if i+j < n: + ...: A[i, i+j] = amat[i*6+j] + ...: A[i+j, i] = amat[i*6+j] + ...: Ab[j, i] = amat[i*6+j] + ...: + ...: # Assert result is correct + ...: assert_allclose(emg3d_solve(amat, bvec), out, rtol=1e-6) + ...: assert_allclose(emg3d_solve2(amat, bvec), out, rtol=1e-6) + ...: assert_allclose(np_linalg_solve(A, bvec), out, rtol=1e-6) + ...: assert_allclose(sp_linalg_solve(A, bvec), out, rtol=1e-6) + ...: assert_allclose(sp_linalg_lu_solve(A, bvec), out, rtol=1e-6) + ...: if cr == 'real': + ...: assert_allclose(sp_linalg_cho_solve(A, bvec), + ...: out, rtol=1e-6) + ...: assert_allclose(sp_linalg_cho_banded(Ab, bvec), + ...: out, rtol=1e-6) + ...: + ...: # Test speed + ...: print(' np.linalg.solve : ', end='') + ...: %timeit np_linalg_solve(A, bvec) + ...: + ...: print(' sp.linalg.solve : ', end='') + ...: %timeit sp_linalg_solve(A, bvec) + ...: + ...: print(' sp.linalg.lu_solve : ', end='') + ...: %timeit sp_linalg_lu_solve(A, bvec) + ...: + ...: if cr == 'real': + ...: + ...: print(' sp.linalg.cho_solve : ', end='') + ...: %timeit sp_linalg_cho_solve(A, bvec) + ...: + ...: print(' sp.linalg.cho_banded : ', end='') + ...: %timeit sp_linalg_cho_banded(Ab, bvec) + ...: + ...: print(' emg3d.solve : ', end='') + ...: %timeit emg3d_solve(amat, bvec) + ...: + ...: print(' emg3d.solve2 : ', end='') + ...: %timeit emg3d_solve2(amat, bvec) + ...: + ...: print(80*'-') diff --git a/docs/dev/tests_benchmarks.rst b/docs/dev/tests_benchmarks.rst deleted file mode 100644 index 7f60affa..00000000 --- a/docs/dev/tests_benchmarks.rst +++ /dev/null @@ -1,25 +0,0 @@ -Tests and Benchmarks -==================== - -.. todo:: - - Rework this for version 1. - -The modeller comes with a test suite using ``pytest``. If you want to run the -tests, just install ``pytest`` and run it within the ``emg3d``-top-directory. - -.. code-block:: console - - > make pytest - -It should run all tests successfully. Please let us know if not! - -Note that installations of ``em3gd`` via conda or pip do not have the -test-suite included. To run the test-suite you must download or clone ``emg3d`` -from GitHub. - -There is also a benchmark suite using *airspeed velocity*, located in the -`emsig/emg3d-asv `_-repository. The results -of my machine can be found in the `emsig/emg3d-bench -`_, its rendered version at -`emsig.xyz/emg3d-asv `_. diff --git a/docs/index.rst b/docs/index.rst index 0f156e1b..91f89953 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -7,7 +7,6 @@ emg3d Documentation :Release: |version| :Date: |today| :Source: `github.com/emsig/emg3d `_ -:EMSiG: `emsig.xyz `_ ---- @@ -18,9 +17,114 @@ emg3d Documentation api/index dev/index -- :ref:`manual` -- `Gallery `_ -- :ref:`api` -- :ref:`development` +.. panels:: + + --- + + :fa:`book,fa-4x` + + User Guide + ^^^^^^^^^^ + + The manual contains installation instructions, information on how to get + started, tips and tricks, background theory, how to cite emg3d, important + references, license information, and much more! + + +++ + + .. link-button:: manual + :type: ref + :text: To the user guide + :classes: btn-block btn-info stretched-link + + --- + + :fa:`fas fa-picture-o,fa-4x` + + Gallery + ^^^^^^^ + + The gallery contains examples and tutorials on the usage of emg3d, and is + generally the best way to get started. Download them and modify them to + your needs! It is also a good place to see what is possible with emg3d. + + +++ + + .. link-button:: https://emsig.xyz/emg3d-gallery + :type: url + :text: To the gallery + :classes: btn-block btn-info stretched-link + + --- + + :fa:`code,fa-4x` + + API reference + ^^^^^^^^^^^^^ + + The API reference of emg3d is extensive, and includes almost every function + and class. A lot of the underlying theory is described in the docstring of + the functions, particularly in the ``core`` module. + + +++ + + .. link-button:: api + :type: ref + :text: To the API reference + :classes: btn-block btn-info stretched-link + + --- + + :fa:`terminal,fa-4x` + + Developer guide + ^^^^^^^^^^^^^^^ + + emg3d is a community code, please help to shape its future! From typos to + bugs to new developments, every contribution is welcomed. This section is + the best way to get started. + + +++ + + .. link-button:: development + :type: ref + :text: To the development guide + :classes: btn-block btn-info stretched-link + + --- + :column: col-lg-12 p-2 + + :fa:`comments,fa-4x` + + Contact + ^^^^^^^ + + .. link-button:: https://emsig.xyz/#contributing-contact + :type: url + :text: Go to «Contributing & Contact» @ emsig.xyz + :classes: btn-block btn-info stretched-link + --- + :column: col-lg-12 p-2 + + .. dropdown:: About the name and logo of **emg3d** + + The name **emg3d** is a portmanteau of *electromagnetic* (em), + *multigrid* (mg), and *three-dimensional* (3D): emg3d is a + three-dimensional multigrid solver for electromagnetic diffusion. + + The symbol stands for the full ('F') multigrid cycling scheme through + the different grid sizes. This is given, for a maximum of three + coarsening steps, by:: + + h_ + 2h_ \ / + 4h_ \ /\ / + 8h_ \/\/ \/ + + + See, for example, the theory section around :numref:`Figure %s + `. + + .. image:: ./_static/emg3d-logo.svg diff --git a/docs/manual/about.rst b/docs/manual/about.rst index 17d46b27..88cd63d2 100644 --- a/docs/manual/about.rst +++ b/docs/manual/about.rst @@ -1,21 +1,20 @@ -.. _about: - About ===== The code ``emg3d`` [WeMS19]_ is a three-dimensional modeller for -electromagnetic (EM) diffusion as used, for instance, in controlled-source EM -(CSEM) surveys frequently applied in the search for resources such as -groundwater, geothermal energy, hydrocarbons, and minerals. +electromagnetic (EM) diffusion as used for instance in geophysical +controlled-source EM (CSEM) surveys. This includes use cases in the search for +resources such as groundwater, geothermal energy, hydrocarbons, and minerals, +or civil engineering and environmental applications. The core of the code is primarily based on [Muld06]_, [Muld07]_, and [Muld08]_. -You can read more about the background of the code in the chapter -:doc:`credits`. An introduction to the underlying theory of multigrid methods -is given in the chapter :doc:`theory`, and further literature is provided in -the :doc:`references`. The code is currently restricted to regular, stretched -grids. As a matrix-free multigrid solver it scales linearly with the number of -cells for both CPU and RAM. This makes it possible to use emg3d for models with -several millions of cells on a regular laptop. +You can read more about the background of the code in the :doc:`credits`. An +introduction to the underlying theory of multigrid methods is given in the +:doc:`theory`, and further literature is provided in the :doc:`references`. The +code is currently restricted to regular, stretched grids. As a matrix-free +multigrid solver it scales linearly with the number of cells for both CPU and +RAM. This makes it possible to use emg3d for models with several millions of +cells on a regular laptop. @@ -48,11 +47,12 @@ What is it _not_? ----------------- - The code is meant to be used in Python or in a terminal. There is **no** GUI. -- Some knowledge of EM fields is definitely helpful, as GIGO applies («garbage - in, garbage out»). For example, placing your receivers very close to the - computational boundary *will* result in bad or wrong responses. +- Some knowledge of EM fields in particular and numerical modelling in general + is definitely helpful, as GIGO applies («garbage in, garbage out»). For + example, placing your receivers very close to the computational boundary + *will* result in bad or wrong responses. - It is not a model builder; there are other tools that can be used to generate - complex geological models, e.g., `GemPy `_. + complex geological models, for instance `GemPy `_. Related ecosystem @@ -65,28 +65,9 @@ parameters and resulting fields. Furthermore, it can serve as a link to use `PyVista `_ to create nice 3D plots even within a notebook. -`EMSiG `_ with its codes `empymod -`_ and `emg3d `_ is part -of a bigger, fast growing, open-source **EM & Potential Geo-Exploration Python -Ecosystem**: - -.. raw:: html - -

- - - - - - - - - - - - - - - -

+There are some first successful attempts of using emsig as a forward modeller +in both `SimPEG `_ and `pyGIMLi `_ +inversions. Get in touch if you are interested in these developments. +See also the note about the `EM & Potential Geo-Exploration Python Ecosystem +`_ on `emsig.xyz `_. diff --git a/docs/manual/citation.rst b/docs/manual/citation.rst index 089d386e..e8a35878 100644 --- a/docs/manual/citation.rst +++ b/docs/manual/citation.rst @@ -1,15 +1,17 @@ -.. _citation: - Citation ======== If you publish results for which you used `emg3d`, please give credit by citing `Werthmüller et al. (2019) `_: - Werthmüller, D., W. A. Mulder, and E. C. Slob, 2019, - emg3d: A multigrid solver for 3D electromagnetic diffusion: - Journal of Open Source Software, 4(39), 1463; - DOI: `10.21105/joss.01463 `_. +.. panels:: + :container: container-lg pb-1 + :column: col-lg-12 p-2 + + Werthmüller, D., W. A. Mulder, and E. C. Slob, 2019, + emg3d: A multigrid solver for 3D electromagnetic diffusion: + Journal of Open Source Software, 4(39), 1463; + DOI: `10.21105/joss.01463 `_. All releases have a Zenodo-DOI, which can be found on `10.5281/zenodo.3229006 `_. diff --git a/docs/manual/cli.rst b/docs/manual/cli.rst index a4762173..519de36c 100644 --- a/docs/manual/cli.rst +++ b/docs/manual/cli.rst @@ -1,31 +1,30 @@ CLI interface ============= -.. todo:: - - The CLI section needs rework. - -Command-line interface for certain specific tasks, such as forward modelling -and gradient computation of the misfit function. The command is ``emg3d``, -consult the inbuilt help to get started: +The command-line interface can be used for certain specific tasks, such as +forward modelling and the computation of the gradient. The command is +``emg3d``, consult the inbuilt help in your terminal to get started: .. code-block:: console emg3d --help The CLI is driven by command-line parameters and a configuration file. The -default configuration file is ``emg3d.cfg``, but another name can be provided -as first positional argument to ``emg3d``. Note that arguments provided in the -command line overwrite the settings in the configuration file. +default configuration filename is ``emg3d.cfg``, but another name can be +provided as first positional argument to ``emg3d``. Note that arguments +provided in the command line overwrite the settings in the configuration file. + +For an example see +https://emsig.xyz/emg3d-gallery/gallery/tutorials/cli.html>`_ in the gallery. Format of the config file ------------------------- -The shown values are the defaults. All values are commented out in this -example; remove the comment signs to use them. +The shown values are examples. All values are commented out in this example; +remove the comment signs to use them. -:: +``emg3d.cfg``:: # Files # ----- @@ -45,8 +44,9 @@ example; remove the comment signs to use them. # for `compute()`. [simulation] # max_workers = 4 # Also via `-n` or `--nproc`. - # gridding = single # One grid for all sources and frequencies. + # gridding = single # min_offset = 0.0 # Only relevant if `observed=True` (r`. It -is a right-handed system (RHS) with x pointing East, y pointing North, and z -pointing upwards. The azimuth is defined as the anticlockwise rotation from -Easting towards Northing, and elevation is defined as the anticlockwise -rotation from the horizontal plane up. +The used coordinate system is shown in :numref:`Figure %s `. +It is a right-handed system (RHS) with x pointing East, y pointing North, and z +pointing upwards. Azimuth :math:`\theta` is defined as the anticlockwise +rotation from Easting towards Northing, and elevation :math:`\varphi` is +defined as the anticlockwise rotation from the horizontal plane up. .. figure:: ../_static/coordinate_system.svg :align: center @@ -28,41 +24,53 @@ rotation from the horizontal plane up. Grid dimension -------------- -The function :func:`emg3d.solver.solve` is the main entry point, and it takes -care whether multigrid is used as a solver or as a preconditioner (or not at -all), while the actual multigrid solver is :func:`emg3d.solver.multigrid`. Most -input parameters for :func:`emg3d.solver.solve` are sufficiently described in -its docstring. Here a few additional information. - -- You can input any three-dimensional tensor mesh into `emg3d`. However, the - implemented multigrid technique works with the existing nodes, meaning there - are no new nodes created as coarsening is done by combining adjacent - cells. The more times the grid dimension can be divided by two the better it - is suited for MG. Ideally, the number should be dividable by two a few times - and the dimension of the coarsest grid should be 2 or a small, odd number - :math:`p`, for which good sizes can then be computed with :math:`p 2^n`. Good - grid sizes (in each direction) up to 1024 are - - - :math:`2·2^{3, 4, ..., 9}`: 16, 32, 64, 128, 256, 512, 1024, - - :math:`3·2^{3, 4, ..., 8}`: 24, 48, 96, 192, 384, 768, - - :math:`5·2^{3, 4, ..., 7}`: 40, 80, 160, 320, 640, - - :math:`7·2^{3, 4, ..., 7}`: 56, 112, 224, 448, 896, - - and preference decreases from top to bottom row (stick to the first two or - three rows if possible). Good grid sizes in sequential order, excluding p=7: - 16, 24, 32, 40, 48, 64, 80, 96, 128, 160, 192, 256, 320, 384, 512, 640, 768, - 1024. You can get this list via :func:`emg3d.meshes.good_mg_cell_nr()`. - -- The multigrid method can be used as a solver or as a preconditioner, for - instance for BiCGSTAB. Using multigrid as a preconditioner for BiCGSTAB - together with semicoarsening and line relaxation is the most stable version, - but expensive, and therefore only recommended on highly stretched grids. - Which combination of solver is best (fastest) depends to a large extent on - the grid stretching, but also on anisotropy and general model complexity. - See `«Parameter tests» - `_ - in the gallery for an example how to run some tests on your particular - problem. +The multigrid method, as implemented, puts certain restrictions to the grid +dimension. + +.. note:: + + If you use emg3d through the :ref:`high-level-usage` + (:class:`emg3d.simulations.Simulation`) with automatic gridding, or use + :func:`emg3d.meshes.construct_mesh`, then the following is taken care of by + emg3d itself. However, if you define the computational grids yourself, the + following section is important. + +You can provide any three-dimensional regular (stretched) grid into emg3d. +However, the implemented multigrid technique works with the existing nodes, +meaning there are no new nodes created as coarsening is done by combining +adjacent cells. The more times the grid dimension can be divided by two the +better it is therefore suited for MG. Ideally, the number should be dividable +by two a few times and the dimension of the coarsest grid should be 2 or a +small, odd number :math:`p`, for which good sizes can then be computed with +:math:`p\,2^n`. Good grid sizes (in each direction) up to 1024 are + +- :math:`2·2^{3, 4, ..., 9}`: 16, 32, 64, 128, 256, 512, 1024, +- :math:`3·2^{3, 4, ..., 8}`: 24, 48, 96, 192, 384, 768, +- :math:`5·2^{3, 4, ..., 7}`: 40, 80, 160, 320, 640, + +and preference decreases from top to bottom row. In sequential order: 16, 24, +32, 40, 48, 64, 80, 96, 128, 160, 192, 256, 320, 384, 512, 640, 768, 1024. You +can obtain the good cell number via :func:`emg3d.meshes.good_mg_cell_nr()`. + + +Solver or Preconditioner +------------------------- + +The multigrid method can be used as a solver on its own or as a preconditioner +for a Krylov subspace solver such as BiCGSTAB. This can be controlled through +the parameter ``sslsolver`` in :func:`emg3d.solver.solve`. *SSL* stands here +for the module :mod:`scipy.sparse.linalg`, and other solvers provided by this +module can be used as well. + +Using multigrid as a preconditioner for BiCGSTAB together with semicoarsening +and line relaxation is the **most stable combination**, for which it is the +default setting. However, it is also the **most expensive** setting, and often +you can obtain faster results by adjusting the combination of solver, +semicoarsening, and line relaxation. Which combination is best (fastest) +depends to a large extent on the grid stretching, but also on anisotropy and +general model complexity. See `«Parameter tests» +`_ in +the gallery for an example how to run some tests on your particular problem. CPU & RAM @@ -70,14 +78,19 @@ CPU & RAM The multigrid method is attractive because it shows optimal scaling for both runtime and memory consumption. In the following are a few notes regarding -memory and runtime requirements. It also contains information about what has -been tried and what still could be tried in order to improve the current code. +memory and runtime requirements. If you are interested in helping to improve +either have a look at :ref:`improve-cpu-ram`. Runtime ``````` -An example of a runtime test is shown in :numref:`Figure %s `. +An example of a runtime test is shown in :numref:`Figure %s `. The +example shows the runtime to solve for a source of 1.0 Hz at the origin of a +homogeneous space of 1.0 Ohm.m, where the grid starts at 32 x 32 x 32 (32,768) +to 384 x 384 x 384 (56,623,104). (You can find the script in +:ref:`improve-cpu-ram`.) + .. figure:: ../_static/CPU.png :scale: 80 % @@ -88,70 +101,19 @@ An example of a runtime test is shown in :numref:`Figure %s `. Runtime as a function of cell size, which shows nicely the linear scaling of multigrid solvers (using a single thread). -The costliest functions (for big models) are: - - - >90 %: :func:`emg3d.solver.smoothing` (:func:`emg3d.core.gauss_seidel`) - - <5 % each, in decreasing importance: - - - :func:`emg3d.solver.prolongation` - (:class:`emg3d.solver.RegularGridProlongator`) - - :func:`emg3d.solver.residual` (:func:`emg3d.core.amat_x`) - - :func:`emg3d.solver.restriction` -Example with 262,144 / 2,097,152 cells (``nu_{i,1,c,2}=0,2,1,2``; -``sslsolver=False``; ``semicoarsening=True``; ``linerelaxation=True``): - - - 93.7 / 95.8 % ``smoothing`` - - 3.6 / 2.0 % ``prolongation`` - - 1.9 / 1.9 % ``residual`` - - 0.6 / 0.4 % ``restriction`` - -The rest can be ignored. For small models, the percentage of ``smoothing`` goes -down and of ``prolongation`` and ``restriction`` go up. But then the modeller -is fast anyway. - -:func:`emg3d.core.gauss_seidel` and :func:`emg3d.core.amat_x` are written -in ``numba``; jitting :class:`emg3d.solver.RegularGridProlongator` turned out -to not improve things, and many functions used in the restriction are jitted -too. The costliest functions (RAM- and CPU-wise) are therefore already written -in ``numba``. - -**Any serious attempt to improve the speed will have to tackle the smoothing -itself.** - - -**Things which could be tried** - -- Not much has been tested with the ``numba``-options ``parallel``; ``prange``; - and ``nogil``. -- There might be an additional gain by making :class:`emg3d.meshes.TensorMesh`, - :class:`emg3d.models.Model`, and :class:`emg3d.fields.Field` instances jitted - classes. - -**Things which have been tried** - -- One important aspect of the smoothing part is the memory layout. - :func:`emg3d.core.gauss_seidel` and :func:`emg3d.core.gauss_seidel_x` - are ideal for F-arrays (loop z-y-x, hence slowest to fastest axis). - :func:`emg3d.core.gauss_seidel_y` and - :func:`emg3d.core.gauss_seidel_z`, however, would be optimal for C-arrays. - But copying the arrays to C-order and afterwards back is costlier in most - cases for both CPU and RAM. The one possible and therefore implemented - solution was to swap the loop-order in :func:`emg3d.core.gauss_seidel_y`. -- Restriction and prolongation information could be saved in a dictionary - instead of recomputing it every time. Turns out to be not worth the - trouble. -- Rewrite :class:`emg3d.solver.RegularGridProlongator` as jitted function, but - the iterator approach seems to be better for large grids. +The result shows the linear scaling: if you double the number of cells, you +double the runtime. Memory `````` -Most of the memory requirement comes from storing the data itself, mainly the -fields (source field, electric field, and residual field) and the model -parameters (resistivity, eta, mu). For a big model, they some up; e.g., almost -3 GB for an isotropic model with 256x256x256 cells. +Most of the memory requirement in emg3d comes from storing the data itself, +mainly the fields (source field, electric field, and residual field) and the +model parameters (resistivity, eta, mu). For a big model, they some up; e.g., +almost 3 GB for an isotropic model with 256 x 256 x 256 cells. The overhead +from the computation is small in comparison. An example of a memory test is shown in :numref:`Figure %s `. @@ -167,99 +129,6 @@ An example of a memory test is shown in :numref:`Figure %s `. for solving one multigrid F-Cycle. -The theory of multigrid says that in an ideal scenario, multigrid requires -8/7 (a bit over 1.14) the memory requirement of carrying out one Gauss-Seidel -step on the finest grid. As can be seen in the figure, for models up to 2 -million cells that holds pretty much, afterwards it becomes a bit worse. - -However, for this estimation one has to run the model first. Another way to -estimate the requirement is by starting from the RAM used to store the fields -and parameters. As can be seen in the figure, for big models one is on the -save side estimating the required RAM as 1.35 times the storage required for -the fields and model parameters. - -The figure also shows nicely the linear behaviour of multigrid; for twice the -number of cells twice the memory is required (from a certain size onwards). - -**Attempts at improving memory usage should focus on the difference between the -red line (actual usage) and the dashed black line (1.14 x base usage).** - -Scripts -``````` - -To test CPU and RAM on your machine, you can use and adjust the following -script. The old notebooks which were used to generate the above figures can be -found at - -- RAM: `4a_RAM-requirements.ipynb - `_, -- CPU: `4b_Runtime.ipynb - `_. - -.. ipython:: - :verbatim: - - In [1]: import emg3d - ...: import numpy as np - ...: import matplotlib.pyplot as plt - ...: from memory_profiler import memory_usage - - In [2]: def compute(nx): - ...: """Simple computation routine. - ...: - ...: This is the actual model it runs. Adjust this to your needs. - ...: - ...: - Model size is nx * nx * nx, centered around the origin. - ...: - Source is at the origin, x-directed. - ...: - Frequency is 1 Hz. - ...: - Homogenous space of 1 Ohm.m. - ...: - ...: """ - ...: - ...: # Grid - ...: hx = np.ones(nx)*50 - ...: x0 = -nx//2*50 - ...: grid = emg3d.TensorMesh([hx, hx, hx], x0=(x0, x0, x0)) - ...: - ...: # Model and source field - ...: model = emg3d.Model(grid, property_x=1.0) - ...: sfield = emg3d.get_source_field( - ...: grid, source=[0, 0, 0, 0, 0], frequency=1.0) - ...: - ...: # Compute the field - ...: _, inf = emg3d.solve( - ...: model, sfield, verb=0, plain=True, return_info=True) - ...: - ...: return inf['time'] - - In [3]: # Loop over model sizes (adjust to your needs). - ...: nsizes = np.array([32, 48, 64, 96, 128, 192, 256, - ...: 384, 512, 768, 1024]) - ...: memory = np.zeros(nsizes.shape) - ...: runtime = np.zeros(nsizes.shape) - ...: - ...: # Loop over nx - ...: for i, nx in enumerate(nsizes): - ...: print(f" => {nx}^3 = {nx**3:12,d} cells") - ...: mem, time = memory_usage((compute, (nx, ), {}), retval=True) - ...: memory[i] = max(mem) - ...: runtime[i] = time - ...: - - In [4]: # Plot CPU - ...: plt.figure() - ...: plt.title('Runtime') - ...: plt.loglog(nsizes**3/1e6, runtime, '.-') - ...: plt.xlabel('Number of cells (in millions)') - ...: plt.ylabel('CPU (s)') - ...: plt.axis('equal') - ...: plt.show() - - In [5]: # Plot RAM - ...: plt.figure() - ...: plt.title('Memory') - ...: plt.loglog(nsizes**3/1e6, memory/1e3, '-', zorder=10) - ...: plt.xlabel('Number of cells (in millions)') - ...: plt.ylabel('RAM (GB)') - ...: plt.axis('equal') - ...: plt.show() +The results show again nicely the linear behaviour of multigrid; for twice the +number of cells twice the memory is required (from a certain size onwards, for +small models there is an non-negligible overhead). diff --git a/docs/manual/installation.rst b/docs/manual/installation.rst index fcb04e53..025b0764 100644 --- a/docs/manual/installation.rst +++ b/docs/manual/installation.rst @@ -1,5 +1,3 @@ -.. _installation: - Installation ============ @@ -20,14 +18,17 @@ and ``numba``. Various other packages are recommended or required for some advanced functionalities, namely: - ``xarray``: For the :class:`emg3d.surveys.Survey` and - :class:`emg3d.simulations.Simulation` classes (many sources and receivers at - once). + :class:`emg3d.simulations.Simulation` classes (model many sources and + frequencies at once). - ``discretize``: For advanced meshing tools (fancy mesh-representations and plotting utilities). - ``matplotlib``: To use the plotting utilities within ``discretize``. - ``h5py``: Save and load data in the HDF5 format. -- ``empymod``: Time-domain modelling (:class:`emg3d.utils.Fourier`). +- ``empymod``: Time-domain modelling (:class:`emg3d.time.Fourier`). - ``scooby``: For the version and system report (:class:`emg3d.utils.Report`). +- ``tqdm``: For nice progress bars when computing many sources and frequencies. + +All soft dependencies are also available both on ``conda-forge`` and ``pip``. If you are new to Python we recommend using a Python distribution, which will ensure that all dependencies are met, specifically properly compiled versions diff --git a/docs/manual/io.rst b/docs/manual/io.rst new file mode 100644 index 00000000..13d15617 --- /dev/null +++ b/docs/manual/io.rst @@ -0,0 +1,163 @@ +.. _io-persistence: + +I/O & Persistence +================= + + +Saving and loading +------------------ + +emg3d has functions to store data to disk and to read data from disk, in +different file formats. Currently three file formats are supported, each +with its own advantages and disadvantages: + +- ``.h5``: + Uses h5py to store inputs to a hierarchical, compressed binary HDF5 file. + **Recommended file format.** + + - Advantage: Widely used, compressed file format, which can be read and + written in many programs. + - Disadvantage: You have to install ``h5py``. + + +- ``.npz``: + Uses numpy to store inputs to a flat, compressed binary file. + + - Advantage: No extra installation is required, and the outputs are + compressed. + - Disadvantage: Only useful within the Python ecosystem. + + +- ``.json``: + Uses json to store inputs to a hierarchical, plain text file. + + - Advantage: No extra installation is required, and the output is a plain + text file that can be viewed in any editor; good for developing and + debugging. + - Disadvantage: Not compressed (files can become huge). + + +Example +~~~~~~~ + +You should be able to save and load everything you do in emg3d with these +functions. Please have a look at the API of :func:`emg3d.io.save` and +:func:`emg3d.io.load`. But in a nutshell, the first argument is a string +containing the relative or absolute path, file name, and the appropriate suffix +indicating the file format. Afterwards it is simply a ``name=value`` list, +where the name can be anything, and the value must be an existing variable. +(There are a few more options, see the API.) + +.. ipython:: + :verbatim: + + In [1]: emg3d.save( + ...: '/path/to/filename.ending', + ...: inp_model=model1, + ...: out_model=model2, + ...: survey=survey, + ...: efield=efield, + ...: ) + Out[1]: Data saved to «/path/to/filename.ending» + +When you load such a file it will give you a dictionary containing as keys the +names you have defined: + +.. ipython:: + :verbatim: + + In [1]: data = emg3d.load('/path/to/filename.ending') + Out[1]: Data loaded from «/path/to/filename.ending» + + In [2]: data.keys() + Out[2]: dict_keys(['_date', '_format', '_version', 'efield', 'inp_model', 'out_model', 'survey']) + +In addition to the variables you have defined there are a few other, "private" +(starting with an underscore) variables such as the date, format, and version +of emg3d with which the archive was created. + + +``{to;from}_file`` +~~~~~~~~~~~~~~~~~~ + +The two classes :class:`emg3d.surveys.Survey` and +:class:`emg3d.simulations.Simulation` have ``to_file`` and ``from_file`` +methods, which are basically wrappers around the saving and loading functions. +They can be used in the following way: + +Storing to disk + +.. ipython:: + :verbatim: + + In [1]: my_survey.to_file('mydata.h5') + + +and loading from disk + +.. ipython:: + :verbatim: + + In [1]: my_survey = emg3d.Survey.from_file('mydata.h5') + + + + +Serialization +------------- + +The following are advanced information if you want to read data created with +emg3d outside of Python or if you want to create data outside of Python which +you can read subsequently with emg3d. As a pure end-user of emg3d you can +ignore this section. + +Here a few info with regards to the (de-)serialization used in emg3d. + +- When invoking ``emg3d.save('filename.ending', a=a, b=something, foo=bar)``, + the data is collected in a dict ``{'a': a, 'b': something, 'foo': bar}``. +- Afterwards the dict is serialized. Instances of emg3d + (:class:`emg3d.meshes.TensorMesh`, :class:`emg3d.fields.Field`, + :class:`emg3d.surveys.Survey`, :class:`emg3d.simulations.Simulation`) have + ``to_dict`` and ``from_dict`` methods to (de-)serialize themselves. These are + used when saving and loading them. In principal emg3d can save everything + that is either serialized already or is present in + ``emg3d.utils._KNOWN_CLASSES``. You can define your own classes which have + ``{to;from}_dict`` methods, and add them to the known classes with the + decorator ``@utils._known_class``. + + - Things which are done when serializing and undone when de-serializing: + + - ``None`` is saved as a string ``'NoneType'``. + + - Things done when serializing: + + - Dictionary key names are converted to strings + - Grids generated with discretize are stored as if they were created using + emg3d. + + - Things done when de-serializing: + + - ``np.bool_`` is returned as ``bool``. + + + + +These first two points are always carried out. After this it depends on the +file format, as different file formats have different limitations. + + + +- ``.h5``: + Each nesting level creates a new data set. + + +- ``.npz``: + The serialized dict is converted into a flattened dict, where the keys are + separated with ``'>'``. + +- ``.json``: + + - NumPy-arrays are turned into lists, where ``'__array-'`` plus the ``dtype`` + are added to the key. + - Complex numbers are stacked, real values followed by imaginary values; + ``__complex`` is added to the key. diff --git a/docs/manual/license.rst b/docs/manual/license.rst index 684d38d1..efd02762 100644 --- a/docs/manual/license.rst +++ b/docs/manual/license.rst @@ -18,7 +18,7 @@ License of Text and Website License of Code --------------- -Copyright 2018-2021 The EMSiG community. +Copyright 2018-2021 The emsig community. Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. diff --git a/docs/manual/mgwhat.rst b/docs/manual/mgwhat.rst index 7f3bee4a..9233e99f 100644 --- a/docs/manual/mgwhat.rst +++ b/docs/manual/mgwhat.rst @@ -5,8 +5,8 @@ If you have never heard of the multigrid method before you might ask yourself "*multi-what?*" The following is an intent to describe the multigrid method without the maths; just some keywords and some figures. **It is a heavily simplified intro, using a 2D grid for simplicity.** Have a look at the -:doc:`theory`-section for more details. A good, four-page intro with some maths -is given by [Muld11]_. More in-depth information can be found, e.g., in +:doc:`theory`-section for more details. A good, six-page intro with some maths +is given by [Muld20]_. More in-depth information can be found, e.g., in [BrHM00]_, [Hack85]_, and [Wess91]_. The multigrid method ([Fedo64]_) @@ -27,7 +27,7 @@ To solve the system, it solves for all fields adjacent to one node, moves then to the next node, and so on until it reaches the last node, see :numref:`Figure %s `, where the red lines indicate the fields which are solved simultaneously per step (the fields on the boundaries are never computed, as -they are assumed to be 0). +they are assumed to be zero). .. figure:: ../_static/schematics2.svg :width: 60 % @@ -60,7 +60,7 @@ has two advantages: - Coarser grid size transforms lower frequency error to higher frequency error, relatively to cell size, which means faster convergence. -The implemented multigrid method simply joins two adjacent cells to get from +The implemented multigrid method simply joins two adjacent cells to go from finer to coarser grids, see :numref:`Figure %s ` for an example coarsening starting with a 16 cells by 16 cells grid. diff --git a/docs/manual/references.rst b/docs/manual/references.rst index 62a1b6d1..71eb23cf 100644 --- a/docs/manual/references.rst +++ b/docs/manual/references.rst @@ -40,9 +40,9 @@ References diffusion with multigrid: Computing and Visualization in Science, 11, 29--138; DOI: `10.1007/s00791-007-0064-y `_. -.. [Muld11] Mulder, W. A., 2011, *in* Numerical Methods, Multigrid: Springer - Netherlands, 895--900; DOI `10.1007/978-90-481-8702-7_153 - `_. +.. [Muld20] Mulder, W. A., 2020, *in* Numerical Methods, Multigrid: Springer + International Publishing, 1--6; DOI `10.1007/978-3-030-10475-7_153-1 + `_. .. [MuWS08] Mulder, W. A., M. Wirianto, and E. C. Slob, 2008, Time-domain modeling of electromagnetic diffusion with a frequency-domain code: Geophysics, 73, F1--F8; DOI: `10.1190/1.2799093 diff --git a/docs/manual/usage.rst b/docs/manual/usage.rst index 409e2070..73e98961 100644 --- a/docs/manual/usage.rst +++ b/docs/manual/usage.rst @@ -1,133 +1,252 @@ -.. _usage: - Getting started =============== -.. todo:: +There are different usages of ``emg3d``. As an end user you might be interested +in the :ref:`high-level-usage`: define your survey and your model, and +simulate electromagnetic data for them. As a developer you might be interested +in the :ref:`solver-level-usage` or :ref:`cli-level-usage`, to implement emg3d +for instance as a forward modelling kernel in you inversion code. And then +there is also :ref:`time-level-usage`. - The basic example needs complete rework. Include all source and receiver - types. Include saving and loading data: - - :class:`emg3d.electrodes.TxElectricDipole`; - :class:`emg3d.electrodes.TxMagneticDipole`; - :class:`emg3d.electrodes.TxElectricWire`; - - :class:`emg3d.electrodes.RxElectricPoint`; - :class:`emg3d.electrodes.RxElectricPoint`; - - :func:`emg3d.io.save`; - :func:`emg3d.io.load`. +Here we show first an introductory example of the :ref:`high-level-usage`. +Further down you find an overview and more detailed explanation of the +different :ref:`usages`. Basic Example ------------- -Here we show a *very* basic example. To see some more realistic models have a -look at the `gallery `_. This particular -example is also there, with some further explanations and examples to show how -to plot the model and the data; see `«Minimum working example» -`_. It -also contains an example without using ``discretize``. - -First, we load ``emg3d`` and ``discretize`` (to create a mesh), along with -``numpy``: +The following is a *very* basic example. To see some more realistic models have +a look at the `gallery `_. First, we load +``emg3d``, ``numpy``, and ``matplotlib``. Note that this example requires that +you have installed ``discretize`` and ``xarray`` as well. .. ipython:: In [1]: import emg3d ...: import numpy as np - ...: from matplotlib.colors import LogNorm + ...: import matplotlib as mpl + +One of the elementary ingredients for modelling is a model, in our case a +resistivity or conductivity model. emg3d is *not* a model builder. We just +construct here a very simple dummy model. You can see more realistic models in +the `models section `_ of the +gallery. -First, we define the mesh (see :class:`discretize.TensorMesh` for more info). -In reality, this task requires some careful considerations. E.g., to avoid edge -effects, the mesh should be large enough in order for the fields to dissipate, -yet fine enough around source and receiver to accurately model them. This grid -is too small, but serves as a minimal example. +Grid +~~~~ + +We first build a simple four-quadrant grid, centred around the origin. .. ipython:: - In [2]: grid = emg3d.TensorMesh( - ...: h=[[(25, 10, -1.04), (25, 28), (25, 10, 1.04)], - ...: [(50, 8, -1.03), (50, 16), (50, 8, 1.03)], - ...: [(30, 8, -1.05), (30, 16), (30, 8, 1.05)]], - ...: origin='CCC') - ...: grid - Out[2]: - ...: TensorMesh: 49,152 cells + In [1]: hx = np.ones(2)*5000 + ...: grid = emg3d.TensorMesh( + ...: [[5000, 5000], [10000], [5000, 5000]], origin='CCC') + ...: grid # QC + Out[1]: + ...: + ...: TensorMesh: 4 cells ...: - ...: MESH EXTENT CELL WIDTH FACTOR - ...: dir nC min max min max max - ...: --- --- --------------------------- ------------------ ------ - ...: x 48 -662.16 662.16 25.00 37.01 1.04 - ...: y 32 -857.96 857.96 50.00 63.34 1.03 - ...: z 32 -540.80 540.80 30.00 44.32 1.05 - -Next we define a very simple fullspace model with -:math:`\rho_x=1.5\,\Omega\,\text{m}`, :math:`\rho_y=1.8\,\Omega\,\text{m}`, and -:math:`\rho_z=3.3\,\Omega\,\text{m}`. The source is an x-directed dipole at the -origin, with a 10 Hz signal of 1 A. + ...: MESH EXTENT CELL WIDTH FACTOR + ...: dir nC min max min max max + ...: --- --- --------------------------- ------------------ ------ + ...: x 2 -5,000.00 5,000.00 5,000.00 5,000.00 1.00 + ...: y 1 -5,000.00 5,000.00 10,000.00 10,000.00 1.00 + ...: z 2 -5,000.00 5,000.00 5,000.00 5,000.00 1.00 + +Model +~~~~~ + +We use the constructed grid to create a simple resistivity model: .. ipython:: - In [3]: model = emg3d.Model(grid, property_x=1.5, property_y=1.8, - ...: property_z=3.3, mapping='Resistivity') - ...: model - Out[3]: Model [resistivity]; triaxial; 48 x 32 x 32 (49,152) + In [1]: model = emg3d.Model( + ...: grid=grid, + ...: property_x=[1, 10, 0.3, 0.3], + ...: mapping='Resistivity' + ...: ) + ...: model # QC + Out[1]: Model: resistivity; isotropic; 2 x 2 x 2 (4) - In [4]: sfield = emg3d.get_source_field(grid=grid, source=[0, 0, 0, 0, 0], frequency=10) +We defined an isotropic model in terms of resistivities, but through the +``mapping`` parameter one can also define a model in terms of conductivities or +the logarithms thereof. -Now we can compute the electric field with ``emg3d``: +Quick check how the model looks like: .. ipython:: - In [5]: efield = emg3d.solve(model=model, sfield=sfield, verb=4) - Out[5]: - ...: :: emg3d START :: 13:56:59 :: v0.17.1.dev18+gf20d741.d20210309 + @savefig basic_model.png width=4in + In [1]: fig, ax = plt.subplots() + ...: cf = ax.pcolormesh( + ...: grid.cell_centers_x/1e3, + ...: grid.cell_centers_z/1e3, + ...: model.property_x[:, 0, :].T, + ...: shading='nearest', + ...: norm=mpl.colors.LogNorm(), + ...: ) + ...: fig.colorbar(cf) + ...: ax.set_title(r'Depth slice ($\Omega$ m)'); + ...: ax.set_xlabel('Easting (km)'); + ...: ax.set_ylabel('Depth (km)'); + + +So we have an upper halfspace of 0.3 Ohm.m, a lower-left quadrant of 1 Ohm.m, +and a lower-right quadrant of 10 Ohm.m. + +Survey +~~~~~~ + +Now that we have a model we need to define our survey. Currently there are +three source types implemented, + +- :class:`emg3d.electrodes.TxElectricDipole`; +- :class:`emg3d.electrodes.TxMagneticDipole`; +- :class:`emg3d.electrodes.TxElectricWire`; + +and two receiver types, + +- :class:`emg3d.electrodes.RxElectricPoint`; +- :class:`emg3d.electrodes.RxElectricPoint`. + +We are going to define a simple survey with an electric dipole source and a +line of electric point receivers. + +.. ipython:: + + In [1]: source = emg3d.TxElectricDipole( + ...: coordinates=(-3000, 0, 0, 0, 0) # x, y, z, azimuth, elevation + ...: ) ...: - ...: MG-cycle : 'F' sslsolver : False - ...: semicoarsening : False [0] tol : 1e-06 - ...: linerelaxation : False [0] maxit : 50 - ...: nu_{i,1,c,2} : 0, 2, 1, 2 verb : 4 - ...: Original grid : 48 x 32 x 32 => 49,152 cells - ...: Coarsest grid : 3 x 2 x 2 => 12 cells - ...: Coarsest level : 4 ; 4 ; 4 + ...: offsets = np.linspace(-2000, 3000, 21) + ...: receivers = emg3d.surveys.txrx_coordinates_to_dict( + ...: emg3d.RxElectricPoint, + ...: coordinates=(offsets, 0, 0, 0, 0), # x, y, z, azimuth, elevation + ...: ) ...: - ...: [hh:mm:ss] rel. error [abs. error, last/prev] l s + ...: survey = emg3d.Survey( + ...: sources=source, + ...: receivers=receivers, + ...: frequencies=1.0, # Hz + ...: ) ...: - ...: h_ - ...: 2h_ \ / - ...: 4h_ \ /\ / - ...: 8h_ \ /\ / \ / - ...: 16h_ \/\/ \/ \/ + ...: survey # QC + + +Simulation +~~~~~~~~~~ + +Now that we have a model and a survey we can define our simulation: + +.. ipython:: + + In [1]: sim = emg3d.Simulation( + ...: survey=survey, + ...: model=model, + ...: ) ...: - ...: [13:56:59] 2.623e-02 after 1 F-cycles [1.464e-06, 0.026] 0 0 - ...: [13:57:00] 2.253e-03 after 2 F-cycles [1.258e-07, 0.086] 0 0 - ...: [13:57:00] 3.051e-04 after 3 F-cycles [1.704e-08, 0.135] 0 0 - ...: [13:57:00] 5.500e-05 after 4 F-cycles [3.071e-09, 0.180] 0 0 - ...: [13:57:01] 1.170e-05 after 5 F-cycles [6.531e-10, 0.213] 0 0 - ...: [13:57:01] 2.745e-06 after 6 F-cycles [1.532e-10, 0.235] 0 0 - ...: [13:57:01] 6.873e-07 after 7 F-cycles [3.837e-11, 0.250] 0 0 + ...: sim # QC + +From the output we see that we defined a survey with one source, 21 receivers, +and one frequency. We see that we have a model, defined as resistivity, which +consists of four cells. And we see that the simulation created a computational +grid of 96x48x64 cells. + +A simulation takes many input parameters, and you should really read the API +reference of :class:`emg3d.simulations.Simulation`. Particularly important are +the gridding parameters, and as such it is highly recommended to also read +:func:`emg3d.meshes.construct_mesh`. The simulation will try its best to create +suitable computational grids. However, these are not necessarily the best ones, +and the user has to pay attention to the grid creation. Also, the default grids +will most likely be quite big, to be on the safe side. Providing suitable user +inputs can significantly reduce the grid sizes and therefore computational +time! + +So lets check in more detail the computational grid it created: + +.. ipython:: + + In [1]: comp_model = sim.get_model(source='TxED-1', frequency=1.0) + ...: comp_model.grid + Out[1]: ...: - ...: > CONVERGED - ...: > MG cycles : 7 - ...: > Final rel. error : 6.873e-07 + ...: TensorMesh: 294,912 cells ...: - ...: :: emg3d END :: 13:57:01 :: runtime = 0:00:02 + ...: MESH EXTENT CELL WIDTH FACTOR + ...: dir nC min max min max max + ...: --- --- --------------------------- ------------------ ------ + ...: x 96 -6,939.37 13,615.81 91.89 3,045.59 1.42 + ...: y 48 -11,286.30 11,286.30 91.89 3,045.59 1.42 + ...: z 64 -13,651.53 2,245.90 91.89 1,866.31 1.21 +We can see that the grid extends further in the directions where there are +higher resistivities. For instance, in positive z we have 0.3 Ohm.m, so the +grid does only extend roughly to +2.2 km. In positive x and z as well as in +positive and negative y we have 10 Ohm.m, so the grid goes to over 10 km. -So the computation required seven multigrid F-cycles and took just a bit more -than 2 seconds. It was able to coarsen in each dimension four times, where the -input grid had 49,152 cells, and the coarsest grid had 12 cells. +Computing the fields is now a simple command, .. ipython:: - @savefig basic_example.png width=4in - In [6]: grid.plot_slice(efield.field, normal='Y', v_type='Ex', view='abs', - ...: pcolor_opts={'norm': LogNorm()}); + In [1]: sim.compute() + +Results +~~~~~~~ + +Let's plot the electric field at receiver locations (the responses): +.. ipython:: + @savefig basic_receivers.png width=4in + In [1]: responses = sim.data.synthetic.data.squeeze() + ...: fig, ax = plt.subplots() + ...: ax.semilogy(offsets/1e3, abs(responses.real), 'C0o-', label='|Real|') + ...: ax.semilogy(offsets/1e3, abs(responses.imag), 'C1o-', label='|Imag|') + ...: ax.legend() + ...: ax.set_title('Electric field at receivers') + ...: ax.set_xlabel('Easting (km)') + ...: ax.set_ylabel('E-field (V/m)') -Usages ------- +We can also get the entire electric field in a similar way as we got the +computational model, and plot it for QC: +.. ipython:: + + @savefig basic_efield.png width=4in + In [1]: efield = sim.get_efield(source='TxED-1', frequency=1.0) + ...: + ...: fig, ax = plt.subplots() + ...: cf = ax.pcolormesh( + ...: efield.grid.cell_centers_x/1e3, + ...: efield.grid.nodes_z/1e3, + ...: abs(efield.fx[:, 24, :].T), + ...: shading='gouraud', + ...: norm=mpl.colors.LogNorm(vmin=1e-16, vmax=1e-8), + ...: ) + ...: fig.colorbar(cf) + ...: ax.set_xlim([-3.5, 4]) + ...: ax.set_ylim([-3, 1]) + ...: ax.set_title(r'|$E_x$| Field (V/m)'); + ...: ax.set_xlabel('Easting (km)'); + ...: ax.set_ylabel('Depth (km)'); + + +This is frankly a *very* fast rundown, and many things are only scratched at +the surface or not explained at all. However, you should find much more +information and explanation in the rest of the manual, and many examples in the +gallery. + + + +.. _usages: + +Usage levels +------------ + +.. _high-level-usage: Simulations / High-level usage ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -152,8 +271,10 @@ and frequency dependent gridding. (``tqdm`` and ``discretize`` are recommended). +.. _solver-level-usage: + Solver-level usage -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +~~~~~~~~~~~~~~~~~~ .. figure:: ../_static/levels2.svg :align: center @@ -186,6 +307,8 @@ source field manually, as shown in :numref:`Figure %s `. *Note:* This requires only ``emg3d`` (``discretize`` is recommended). +.. _cli-level-usage: + Command-line interface (CLI-level) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -200,12 +323,14 @@ Command-line interface (CLI-level) The command-line interface is a terminal utility for the high-level (Simulation) usage of emg3d. The model and the survey have to be provided as files (HDF5, npz, or json), various settings can be defined in the config file, -and the output will be written to the output file. +and the output will be written to the output file (see also +:ref:`io-persistence`). *Note:* In addition to ``emg3d`` this requires the soft dependency ``xarray`` (``tqdm`` and ``discretize`` are recommended), and ``h5py`` if the provided files are in the HDF5 format. +.. _time-level-usage: Time-domain modelling ~~~~~~~~~~~~~~~~~~~~~ diff --git a/emg3d/__init__.py b/emg3d/__init__.py index 2bf1846d..8b574957 100644 --- a/emg3d/__init__.py +++ b/emg3d/__init__.py @@ -1,4 +1,4 @@ -# Copyright 2018-2021 The EMSiG community. +# Copyright 2018-2021 The emsig community. # # This file is part of emg3d. # diff --git a/emg3d/__main__.py b/emg3d/__main__.py index a08289b6..e747f111 100644 --- a/emg3d/__main__.py +++ b/emg3d/__main__.py @@ -1,4 +1,4 @@ -# Copyright 2018-2021 The EMSiG community. +# Copyright 2018-2021 The emsig community. # # This file is part of emg3d. # diff --git a/emg3d/cli/__init__.py b/emg3d/cli/__init__.py index 19ed3fc8..d152399d 100644 --- a/emg3d/cli/__init__.py +++ b/emg3d/cli/__init__.py @@ -1,4 +1,4 @@ -# Copyright 2018-2021 The EMSiG community. +# Copyright 2018-2021 The emsig community. # # This file is part of emg3d. # diff --git a/emg3d/cli/main.py b/emg3d/cli/main.py index cc0753d6..30b26205 100644 --- a/emg3d/cli/main.py +++ b/emg3d/cli/main.py @@ -1,7 +1,7 @@ """ Entry point for the command-line interface (CLI). """ -# Copyright 2018-2021 The EMSiG community. +# Copyright 2018-2021 The emsig community. # # This file is part of emg3d. # diff --git a/emg3d/cli/parser.py b/emg3d/cli/parser.py index 921d649d..ae035b51 100644 --- a/emg3d/cli/parser.py +++ b/emg3d/cli/parser.py @@ -1,7 +1,7 @@ """ Parser for the configuration file of the command-line interface. """ -# Copyright 2018-2021 The EMSiG community. +# Copyright 2018-2021 The emsig community. # # This file is part of emg3d. # @@ -121,7 +121,7 @@ def parse_config_file(args_dict): # Get absolute paths. ffile = Path(os.path.join(path, fname)) - # Ensure there is a file ending, if not, fall back to h5p + # Ensure there is a file ending, if not, fall back to h5. if ffile.suffix not in ['.h5', '.json', '.npz']: ffile = ffile.with_suffix('.h5') diff --git a/emg3d/cli/run.py b/emg3d/cli/run.py index 662a39d6..e9ee3cd5 100644 --- a/emg3d/cli/run.py +++ b/emg3d/cli/run.py @@ -1,7 +1,7 @@ """ Functions that actually call emg3d within the CLI interface. """ -# Copyright 2018-2021 The EMSiG community. +# Copyright 2018-2021 The emsig community. # # This file is part of emg3d. # diff --git a/emg3d/core.py b/emg3d/core.py index ef0ef1dc..e6bc7ea5 100644 --- a/emg3d/core.py +++ b/emg3d/core.py @@ -20,7 +20,7 @@ start, as by far the most time is spent in these functions, particularly in :func:`solve`. """ -# Copyright 2018-2021 The EMSiG community. +# Copyright 2018-2021 The emsig community. # # This file is part of emg3d. # diff --git a/emg3d/electrodes.py b/emg3d/electrodes.py index e04a8043..b1dbb7b8 100644 --- a/emg3d/electrodes.py +++ b/emg3d/electrodes.py @@ -1,7 +1,7 @@ """ Electrodes define any type of sources and receivers used in a survey. """ -# Copyright 2018-2021 The EMSiG community. +# Copyright 2018-2021 The emsig community. # # This file is part of emg3d. # diff --git a/emg3d/fields.py b/emg3d/fields.py index fc0e02aa..5f844e44 100644 --- a/emg3d/fields.py +++ b/emg3d/fields.py @@ -3,7 +3,7 @@ magnetic field, generating the source field; obtaining the fields at receiver locations. """ -# Copyright 2018-2021 The EMSiG community. +# Copyright 2018-2021 The emsig community. # # This file is part of emg3d. # diff --git a/emg3d/io.py b/emg3d/io.py index da23ba38..5335009f 100644 --- a/emg3d/io.py +++ b/emg3d/io.py @@ -1,7 +1,7 @@ """ Utility functions for writing and reading data. """ -# Copyright 2018-2021 The EMSiG community. +# Copyright 2018-2021 The emsig community. # # This file is part of emg3d. # diff --git a/emg3d/maps.py b/emg3d/maps.py index 5e5bfa5b..f429bbbc 100644 --- a/emg3d/maps.py +++ b/emg3d/maps.py @@ -5,7 +5,7 @@ Interpolation routines mapping values between different grids. """ -# Copyright 2018-2021 The EMSiG community. +# Copyright 2018-2021 The emsig community. # # This file is part of emg3d. # diff --git a/emg3d/meshes.py b/emg3d/meshes.py index 4d75e580..342a43e3 100644 --- a/emg3d/meshes.py +++ b/emg3d/meshes.py @@ -1,7 +1,7 @@ """ Everything related to meshes appropriate for the multigrid solver. """ -# Copyright 2018-2021 The EMSiG community. +# Copyright 2018-2021 The emsig community. # # This file is part of emg3d. # @@ -681,7 +681,8 @@ def origin_and_widths(frequency, properties, center, domain=None, vector=None, for nx in np.unique(cell_numbers): # Loop over possible alphas for domain. - for sa in np.linspace(1.0, stretching[0], 100): + nsa = max(1, min(100, int((stretching[0] - 1) / 0.001))) + for sa in np.linspace(1.0, stretching[0], nsa): if vector is None: @@ -736,7 +737,8 @@ def origin_and_widths(frequency, properties, center, domain=None, vector=None, hxo = hx # Loop over possible alphas for buffer. - for ca in np.linspace(sa, stretching[1], 100): + nca = max(1, min(100, int((stretching[1] - sa) / 0.001))) + for ca in np.linspace(sa, stretching[1], nca): # Get current stretched grid cell sizes. thxl = hx[0]*ca**np.arange(1, nx_remain+1) # Left of survey. diff --git a/emg3d/models.py b/emg3d/models.py index fd6e0d56..0cf728bf 100644 --- a/emg3d/models.py +++ b/emg3d/models.py @@ -1,7 +1,7 @@ """ Everything to store electromagnetic material properties for the solver. """ -# Copyright 2018-2021 The EMSiG community. +# Copyright 2018-2021 The emsig community. # # This file is part of emg3d. # diff --git a/emg3d/optimize.py b/emg3d/optimize.py index 3d10a0f7..254ef127 100644 --- a/emg3d/optimize.py +++ b/emg3d/optimize.py @@ -2,7 +2,7 @@ Functionalities related to optimization (minimization, inversion), such as the misfit function and its gradient. """ -# Copyright 2018-2021 The EMSiG community. +# Copyright 2018-2021 The emsig community. # # This file is part of emg3d. # diff --git a/emg3d/simulations.py b/emg3d/simulations.py index bf63b095..83f75806 100644 --- a/emg3d/simulations.py +++ b/emg3d/simulations.py @@ -5,7 +5,7 @@ The simulation module combines the different pieces of ``emg3d`` providing a high-level, specialised modelling tool for the end user. """ -# Copyright 2018-2021 The EMSiG community. +# Copyright 2018-2021 The emsig community. # # This file is part of emg3d. # diff --git a/emg3d/solver.py b/emg3d/solver.py index 08464e6a..96d90991 100644 --- a/emg3d/solver.py +++ b/emg3d/solver.py @@ -10,7 +10,7 @@ and its implementation. """ -# Copyright 2018-2021 The EMSiG community. +# Copyright 2018-2021 The emsig community. # # This file is part of emg3d. # diff --git a/emg3d/surveys.py b/emg3d/surveys.py index 7507a6d7..22a04348 100644 --- a/emg3d/surveys.py +++ b/emg3d/surveys.py @@ -2,7 +2,7 @@ A survey stores a set of sources and their frequencies, receivers, and the measured data. """ -# Copyright 2018-2021 The EMSiG community. +# Copyright 2018-2021 The emsig community. # # This file is part of emg3d. # diff --git a/emg3d/time.py b/emg3d/time.py index 56f17c94..09871957 100644 --- a/emg3d/time.py +++ b/emg3d/time.py @@ -1,7 +1,7 @@ """ Functionalities related to time-domain modelling using a frequency-domain code. """ -# Copyright 2018-2021 The EMSiG community. +# Copyright 2018-2021 The emsig community. # # This file is part of emg3d. # diff --git a/emg3d/utils.py b/emg3d/utils.py index 8e9dd449..1e35b866 100644 --- a/emg3d/utils.py +++ b/emg3d/utils.py @@ -1,7 +1,7 @@ """ Utility functions for the multigrid solver. """ -# Copyright 2018-2021 The EMSiG community. +# Copyright 2018-2021 The emsig community. # # This file is part of emg3d. # diff --git a/requirements-dev.txt b/requirements-dev.txt index e980b2ff..c34d5cd4 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -2,6 +2,7 @@ -r requirements.txt # SOFT DEPENDENCIES +tqdm h5py scooby xarray @@ -15,12 +16,14 @@ setuptools_scm # FOR DOCUMENTATION sphinx > 3 numpydoc >= 0.9 +sphinx_panels sphinx_numfig pydata_sphinx_theme sphinx_automodapi ipykernel # FOR TESTING +asv pytest coveralls pytest-cov diff --git a/setup.py b/setup.py index 2a20068e..1fc226db 100644 --- a/setup.py +++ b/setup.py @@ -16,7 +16,7 @@ name="emg3d", description="A multigrid solver for 3D electromagnetic diffusion.", long_description=readme, - author="The EMSiG community", + author="The emsig community", author_email="info@emsig.xyz", url="https://emsig.xyz", license="Apache-2.0", diff --git a/tests/alternatives.py b/tests/alternatives.py index 7741ff99..4e70b19d 100644 --- a/tests/alternatives.py +++ b/tests/alternatives.py @@ -13,7 +13,7 @@ not) be easier to understand. """ -# Copyright 2018-2021 The EMSiG community. +# Copyright 2018-2021 The emsig community. # # This file is part of emg3d. # diff --git a/tests/unused.py b/tests/unused.py index 14b7e737..fb8b4209 100644 --- a/tests/unused.py +++ b/tests/unused.py @@ -7,7 +7,7 @@ even just for testing purposes. """ -# Copyright 2018-2021 The EMSiG community. +# Copyright 2018-2021 The emsig community. # # This file is part of emg3d. #