A multigrid solver for 3D electromagnetic diffusion
Home Page: https://emg3d.emsig.xyz
License: Apache License 2.0
A multigrid solver for 3D electromagnetic diffusion in Python.
Available through conda and pip; consult the documentation for detailed installation instructions.
I could use some help setting up a basic problem on a unit cube (or square). I've been able to run the following but trying a forcing function of say all ones has been difficult.
h0 = 0.1 * np.ones(10)
grid = emg3d.utils.TensorMesh((h0,h0,h0), np.zeros(3))
model = emg3d.utils.Model(grid)
sfield = emg3d.utils.Field(grid)
efield = emg3d.solver.solver(grid, model, sfield)
A couple questions:
Create a tutorial in the emg3d-gallery on how to utilize discretize to plot {electric;magnetic} field lines.
Introduce checks for links in the docs.
EDIT
This issue has two steps:
discretize v6.0.0 => #148 ; this is a bug, as it will currently fail;discretize v1.0.0 => currently I count that as an enhancement.The second part requires to implement all the name changes.
There need to be some fixes to make emg3d compatible with discretize>=0.6.0. Most things should work, but there are some tests that fail. Currently it is fixed in:
'discretize<0.6.0'discretize<0.6.0'discretize<0.6.0'=> Fix this and undo the above three lines.
The new examples in section 6 of https://github.com/empymod/emg3d-examples introduce transient CSEM using emg3d.
In these examples, empymod is used for the transform. The required steps include:
=> Add some helper routines to emg3d.utils to simplify these steps.
At the moment, tests are run normally and once with the env DISABLE_JIT set to 1, in order to check the jitted functions for coverage for coveralls.
The better way to do that is with
@pytest.mark.parametrize("njit", [True, False])
def test_my_jitted_function(njit):
if njit:
test_fct = my_jitted_function
else:
test_fct = my_jitted_function.py_func
This way, the tests run once with the jitted functions and once the pure python functions. The advantage of that approach is that we can define which testes are run in pure Python mode, we don't have to run the entire tests in pure python mode.
It is currently implemented in the dev-branch. However, More tests are required to cover 100% of njitted:
gauss_seidel (the test is rather bad at the moment, actually)restrict (only sc_dir=0 is tested, not sc_dir=1, 2, 3, 4, 5, 6)_volume_avg_weightsnjitted.gauss_seidel_{x;y;z}.Issues to tackle:
allowed_failure or similar (actions/runner#2347).The first benchmark always has a lower memory profile, so currently the latest commit is never added to the benchmark results.
The reason might be that the data is loaded outside the setup_cache-fct. Move np.load into get_model in the file emg3d-asv/benchmarks/memory.py and re-run the benchmarks.
There is a memory issue in the supposedly stripped-down, very simplified utils.TensorMesh:
import discretize
import numpy as np
import emg3d
%load_ext memory_profiler
nx, ny, nz = 200, 200, 100
h = [np.ones(nx)*20., np.ones(ny)*20., np.ones(nz)*20.]
%memit discretize.TensorMesh(h, x0='00N')
%memit emg3d.utils.TensorMesh(h, x0=np.array([0, 0, -4200]))
Output:
peak memory: 220.55 MiB, increment: 0.33 MiB
peak memory: 506.71 MiB, increment: 286.16 MiB
I tried Primary/secondary field example shown following url.
https://empymod.github.io/emg3d-gallery/gallery/tutorials/total_vs_ps_field.html#sphx-glr-gallery-tutorials-total-vs-ps-field-py
But, I encountered this error. This error will occur other examples too.
This needs to update to using grid.vectorCCx, grid.vectorCCy, grid.vectorCCz methods.
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
<ipython-input-1-acfab86d284a> in <module>
118 # Layered_background
119 res_x = np.ones(grid.nC)*res[0] # Air resistivity
--> 120 res_x[grid.gridCC[:, 2] < 0] = res[1] # Water resistivity
121 res_x[grid.gridCC[:, 2] < -1000] = res[2] # Background resistivity
122
AttributeError: 'TensorMesh' object has no attribute 'gridCC'
Mac OS Catalina 10.15.5
Python 3.7.7
emg3d 0.12.0 (installed via conda)
scipy 1.5.2
numpy 1.19.1
In b4de17e versioning was handed over to setuptools_scm. Unfortunately, setuptools_scm handles what is included in a build, it includes by default everything under version control, which is definitely not, what we want. MANIFEST can be used to exclude stuff, so we have to update that.
#7 showed that forming gridE* uses significant memory. They are ONLY used in prolongation:
# Calculate required points of finer grid.
yz_points = grid.gridEx[::grid.nCx, 1:]
xz_points = grid.gridEy[:, ::2].reshape(grid.nNz, -1, 2)
xz_points = xz_points[:, :grid.nNx, :].reshape(-1, 2)
xy_points = grid.gridEz[:grid.nNx*grid.nNy, :2]
Is there a better way to get *_points etc without gridE*?
The functions
io.save, andio.loaduse a emg3d-specific file format (using h5, npz, or json). These functions are also used by
surveys.Survey.to_file; surveys.Survey.from_file;simulations.Simulation.to_file; simulations.Simulation.from_file;However, the file format is not well documented.
I'm seeing a test failure most of the time (like 9 out of 10). It looks like successive times is around 1e-5 for me, which is triggering the assert. It's not clear why it's so large.
(emg3d1) ± pytest .
================================ test session starts =================================
platform darwin -- Python 3.7.3, pytest-4.6.3, py-1.8.0, pluggy-0.12.0
plugins: cov-2.7.1
collected 28 items
tests/test_njitted.py ...... [ 21%]
tests/test_solver.py ........ [ 50%]
tests/test_utils.py ........F..... [100%]
====================================== FAILURES ======================================
_____________________________________ test_Time ______________________________________
def test_Time():
t0 = utils.timeit() # Create almost at the same time a
time = utils.Time() # t0-stamp and a Time-instance.
> assert_allclose(t0, time.t0)
E AssertionError:
E Not equal to tolerance rtol=1e-07, atol=0
E
E Mismatch: 100%
E Max absolute difference: 1.493e-06
E Max relative difference: 3.3439148e-07
E x: array(4.464825)
E y: array(4.464827)
tests/test_utils.py:383: AssertionError
======================== 1 failed, 27 passed in 5.43 seconds =========================
In [18]: t1=emg3d.utils.timeit();t2=emg3d.utils.timeit();print(t1,t2,t2-t1)
81.447158438 81.447194455 3.6016999999333166e-05
This should return
Field-instance, if input was a Field (including freq)Model-instance, if input was a ModelIn the latter case, it has to do the interpolation recursively for res_x, res_y, res_z, and mu_r.
After this step, tagging is the only thing required to release a version.
There should be a possibility to retrieve some information from the solver, such as number of iterations or final norm (some of the information which is also printed if verbose).
The benchmarks (https://empymod.github.io/emg3d-asv) fail since d8e98c0 because of the introduction of VolumeModel. Update https://github.com/empymod/emg3d-asv to include this change.
The function get_source_field returns the field for a dipole source. Source fields can be summed up to yield arbitrary sources, e.g. a square loop (hence providing a 'magnetic' source). An example of this is given in the example magnetic_source_el_loop in the gallery:
# Initiate a zero-valued source field.
sfield = emg3d.fields.SourceField(pgrid, freq=freq)
# Define the four dipole segments.
srcloop = [
np.r_[src[0]-0.5, src[0]+0.5, src[1]-0.5, src[1]-0.5, src[2], src[2]],
np.r_[src[0]+0.5, src[0]+0.5, src[1]-0.5, src[1]+0.5, src[2], src[2]],
np.r_[src[0]+0.5, src[0]-0.5, src[1]+0.5, src[1]+0.5, src[2], src[2]],
np.r_[src[0]-0.5, src[0]-0.5, src[1]+0.5, src[1]-0.5, src[2], src[2]],
]
# Add the source fields up.
for srcl in srcloop:
sfield += emg3d.fields.get_source_field(pgrid, srcl, freq, strength)
The function get_source_field should be extended to do this in one call, get_source_field(grid, src, freq, strength), where each parameter in the tuple src=(x0, x1, y0, y1, z0, z1) will be an array instead of a float.
emg3d can be used to model electric and magnetic sources with electric or magnetic receivers.
How to get the magnetic response from the electric field is shown here:
How to calculate the electric and magnetic fields due to a magnetic source is shown in:
However, nicely implemented is the default, the electric field to an electric source.
Implement the other cases properly, with msrc and mrec parameters just like in empymod. Steps to take care of in doing so (not all are mandatory; just take magnetic loop sources and magnetic dipole receivers into account, should be fine):
mu_r from isotropic to tri-axial (otherwise duality will only work for isotropic resistivities). => https://github.com/empymod/emg3d/tree/triaxial-mu-r\eta, \mu_r to definitions \eta=V\sigma=V/\rho, \zeta=V/\mu_r.msrc (el. dipole or loop) and mrec (el. or magn. dipole) into emg3d.solver, hence there are four possible fields that can be returned; use duality (swap +/-; \eta/\zeta; src-dir).get_e_field.exp(-iwt) to exp(iwt) (not necessary, just to make it more alike empymod and avoid occasional confusion).This, together with emsig/empymod#47, will make empymod and emg3d very similar in behaviour and possibilities.
MAINTENANCE.rst
MAINTENANCE.rst in documentation (just like README.rst, CHANGELOG.rst)Create an example for MT measurements.
The benchmarks (https://empymod.github.io/emg3d-asv/) show two peculiarities which should be investigated at some point.
vol is not precomputed any longer (not stored in Model either, but available to compute in TensorMesh), so it has to be done within the solver routine, before it was done outside. If that is the case all good. The other possibility is the change from np.outer to broadcasting. If that is the case then bad.Because empymod is now a dependency, but it is not in the main channel of conda, my build recipe fails. This means it is not available as conda install -c prisae emg3d.
=> Time to move to conda-forge!
Current work-around: pip install emg3d.
Include a verb=-1 for a one-liner verbosity; good for inversion logs or if creating multi-source, multi-frequency data or similar.
:: emg3d :: 6.5e-03; 5(3); 0:00:10; CONVERGED
Version 0.9.0 actually applied the diffusive approximation, setting \epsilon_r=0.
Now, the results are controversial. Slightly less RAM is used, and it looks like it is faster on some examples (example-notebooks 1X, for instance), but slower on others (2a). This needs more investigation.
Todo:
Make functions similar to this empymod.utils.get_minimum() and empymod.utils.set_minimum(), but for setting the defaults for emg3d.solve().
E.g.,
emg3d.utils.set_defaults(
sslsolver=True,
linerelaxation=False,
verb=-1,
)
Having this global option saves a lot of typing later on in large examples.
Create an example using the primary/secondary-field approach.
Replace np.int, np.float, and np.complex with either the pure Python equivalents int, float, complex, or with the numpy versions np.int_, np.float_, and np.complex_.
They are likely going to be deprecated and removed in the future: numpy/numpy#14882
Create a notebook explaining the automatic gridding parameters.
Currently, discretize is not a dependency, and all could be done with utils.TensorMesh. However, this leads sometimes to confusion, and the gain is little. The main reasons why discretize has not been made a dependency are removed:
matplotlib is no longer a hard dependencypymatsolver is no longer a dependency (simpeg/discretize#193)condanumpySo I think it is time to go all-in and make it a hard dependency (adjust docs and website) as soon as all boxes above are ticked and a release happened (pin the version then); utils.TensorMesh should then only be used internally.
I think the best would be to do in emg3d.meshes
from discretize import TensorMesh
and then rename the meshes.TensorMesh to meshes._TensorMesh, for internal use of the solver only.
Adjust accordingly:
Hi all,
I just start to learn emg3d following the example gallery.
I got an index error when I run the script in the tutorial.
Here's the code:
import emg3d
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('ggplot')
print(emg3d.Report())
###############################################################################
# Survey
# ------
src = [0, 0, -950, 0, 0] # x-dir. source at the origin, 50 m above seafloor
off = np.arange(5, 81)*100 # Offsets
rec = [off, off*0, -1000] # In-line receivers on the seafloor
res = [1e10, 0.3, 1] # 1D resistivities (Ohm.m): [air, water, backgr.]
freq = 1.0 # Frequency (Hz)
###############################################################################
# Mesh
# ----
#
# We create quite a coarse grid (100 m minimum cell width), to have reasonable
# fast computation times.
grid = emg3d.construct_mesh(
frequency=freq,
min_width_limits=100.0,
properties=[res[1], 100., 2, 100],
center=(src[0], src[1], -1000),
seasurface=0.0,
domain=([-100, 8100], [-500, 500], [-2500, 0]),
verb=2,
)
print(grid)The output information:
--------------------------------------------------------------------------------
Date: Wed Nov 18 14:57:21 2020 China Standard Time
OS : Windows
CPU(s) : 8
Machine : AMD64
Architecture : 64bit
RAM : 63.9 GB
Environment : Python
Python 3.7.9 (default, Aug 31 2020, 17:10:11) [MSC v.1916 64 bit (AMD64)]
numpy : 1.19.2
scipy : 1.5.4
numba : 0.51.2
emg3d : 0.14.1
empymod : 2.0.3
xarray : 0.16.1
discretize : 0.5.1
h5py : 2.10.0
matplotlib : 3.3.2
IPython : 7.13.0
Intel(R) Math Kernel Library Version 2019.0.0 Product Build 20180829 for
Intel(R) 64 architecture applications
--------------------------------------------------------------------------------
== GRIDDING IN X ==
Skin depth [m] : 276 / 5033 / 5033 [corresponding to `properties`]
Survey domain DS [m] : -100 - 8100
Comp. domain DC [m] : -31723 - 39723
Final extent [m] : -33701 - 41701
Cell widths [m] : 100 / 100 / 5465 [min(DS) / max(DS) / max(DC)]
Number of cells : 128 (82 / 46 / 0) [Total (DS/DC/remain)]
Max stretching : 1.000 (1.000) / 1.190 [DS (seasurface) / DC]
== GRIDDING IN Y ==
Skin depth [m] : 276 / 5033 / 5033 [corresponding to `properties`]
Survey domain DS [m] : -500 - 500
Comp. domain DC [m] : -32123 - 32123
Final extent [m] : -34890 - 34890
Cell widths [m] : 100 / 100 / 9016 [min(DS) / max(DS) / max(DC)]
Number of cells : 40 (10 / 30 / 0) [Total (DS/DC/remain)]
Max stretching : 1.000 (1.000) / 1.350 [DS (seasurface) / DC]
== GRIDDING IN Z ==
Traceback (most recent call last):
File "c:/Users/JingXu/Desktop/code/start_emg3d/model_property_mapping/mapping.py", line 30, in <module>
verb=2,
File "D:\Anaconda3\envs\tdem\lib\site-packages\emg3d\meshes.py", line 462, in construct_mesh
z0, hz = get_origin_widths(**kwargs, **zparams)
File "D:\Anaconda3\envs\tdem\lib\site-packages\emg3d\meshes.py", line 626, in get_origin_widths
thxr = hx[-1]*sa**np.arange(nx)
IndexError: index -1 is out of bounds for axis 0 with size 0
The formulation for emg3d uses electric sources and receivers. The routine get_h_field uses Faraday's law to calculate the magnetic field at receiver locations.
It would be nice to also have the possibility to calculate the fields due to a magnetic source.
Dear Dieter,
I installed emg3d by pip and after that I wanted to import it but I faced with this error:
(base) mohammad@ubuntu:~$ pip3 install emg3d
Collecting emg3d
Cache entry deserialization failed, entry ignored
Downloading https://files.pythonhosted.org/packages/c1/d4/7d2a18dc68cb4e65d1b58a7cd501bfd26cbd7e9180736fbe3d02a326d547/emg3d-0.9.1-py3-none-any.whl (65kB)
100% |████████████████████████████████| 71kB 49kB/s
Collecting numba>=0.40.0 (from emg3d)
Downloading https://files.pythonhosted.org/packages/53/34/22b6c2ded558976b5367be01b851ae679a0d1ba994de002d54afe84187b5/numba-0.46.0-cp36-cp36m-manylinux1_x86_64.whl (3.6MB)
100% |████████████████████████████████| 3.6MB 158kB/s
Collecting numpy>=1.15.0 (from emg3d)
Downloading https://files.pythonhosted.org/packages/d2/ab/43e678759326f728de861edbef34b8e2ad1b1490505f20e0d1f0716c3bf4/numpy-1.17.4-cp36-cp36m-manylinux1_x86_64.whl (20.0MB)
100% |████████████████████████████████| 20.0MB 42kB/s
Collecting scipy>=1.1.0 (from emg3d)
Cache entry deserialization failed, entry ignored
Downloading https://files.pythonhosted.org/packages/54/18/d7c101d5e93b6c78dc206fcdf7bd04c1f8138a7b1a93578158fa3b132b08/scipy-1.3.3-cp36-cp36m-manylinux1_x86_64.whl (25.2MB)
100% |████████████████████████████████| 25.2MB 33kB/s
Collecting empymod (from emg3d)
Downloading https://files.pythonhosted.org/packages/ad/01/b257911dbc987a440be8972384798a368b41be870526959eebc48142f49c/empymod-1.10.3-py3-none-any.whl (217kB)
100% |████████████████████████████████| 225kB 304kB/s
Collecting llvmlite>=0.30.0dev0 (from numba>=0.40.0->emg3d)
Downloading https://files.pythonhosted.org/packages/b3/93/924788871a889ead1c115ba2d43e67932a8865d38fdb179505c99e4bb575/llvmlite-0.30.0-cp36-cp36m-manylinux1_x86_64.whl (20.2MB)
100% |████████████████████████████████| 20.2MB 39kB/s
Installing collected packages: numpy, llvmlite, numba, scipy, empymod, emg3d
Successfully installed emg3d-0.9.1 empymod-1.10.3 llvmlite-0.30.0 numba-0.46.0 numpy-1.17.4 scipy-1.3.3
(base) mohammad@ubuntu:~$ python3
Python 3.7.3 (default, Mar 27 2019, 22:11:17)
[GCC 7.3.0] :: Anaconda, Inc. on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import emg3d
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
ModuleNotFoundError: No module named 'emg3d'
>>> exit()
Regards,
Mohammad
The license copyright statement should be updated from the default:
https://github.com/empymod/emg3d/blob/master/LICENSE#L189
Make it possible to add/subtract two models.
Docs of discretize moved from rtfd to simpeg.xyz. Adjust the link.
Since v0.12.0 there is the optimize module. However, this is still considered experimental, and many things have to be improved. This is an issue to keep track of that:
Also, optimize.gradient should be made stand-alone, so it can be called from other packages.
Hi Dieter, congrats for the new release, seems to be a big one.
The example in gempy broke. I tried to just refactor the functions in utils but in the end I faced a reshaping error in the solution.
File "C:/Users/legui/PycharmProjects/gempy/examples/integrations/gempy_discretize_emg3d.py", line 217, in <module>
cgrid.plot_3d_slicer(
File "C:\Users\legui\miniconda3\envs\gp-dev\lib\site-packages\discretize\View.py", line 830, in plot_3d_slicer
tracker = Slicer(
File "C:\Users\legui\miniconda3\envs\gp-dev\lib\site-packages\discretize\View.py", line 1421, in __init__
self.v = mesh.r(v.reshape((mesh.nC, -1), order='F'), 'CC', 'CC', 'M')
ValueError: cannot reshape array of size 245440 into shape (237568,newaxis)
So far I am fixing the version of emg3d to 0.11 for gempy so no rush
Currently, new Field instances are initiated in C-order (Python default). Change this to F and benchmark it.
Check if the same is the case in the Model initiation.
The data_write and data_read needs some improvement or rethinking in general.
Field: implement (de-)serialization; add meta data.Model: implement (de-)serialization; add meta data.TensorMesh: implement (de-)serialization; add meta data.xarray or not? Internally with .npzdata_write and data_read then still required. Potentially. But it calls the functions own serialization routines and collects them into one file. => NO; deprecatedemg3d.io?)emg3d.io?)Meta data: __version__, date (scooby-dump for all dependencies?)
This is a bug which has always been an issue, in one or another way. The Krylov methods in its current implementation have a problem with initial guesses, either they do not converge, or the field is not updated in place. I have currently no idea whatsoever how to fix this.
No, the bug was fixed in #32 by simply resetting initial guesses to zero if an sslsolver is used.
However, I leave this issue open as an 'Enhancement'-request to re-enable that feature again.
Original initial message:
There is a bug in solver.krylov for very rare cases. I thought I fixed it with
d255e2d
but that was wrong (overwriting a potentially provided efield). So that commit was undone again in
f34dcb6
When using
simulation.compute(observed=True)
obtained following message
TypeError: cannot pickle '_thread.RLock' object
Repo: https://github.com/empymod/emg3d-gallery
Deployed to: https://empymod.github.io/emg3d-gallery
Strings are handled differently since h5py version 3 (https://docs.h5py.org/en/latest/whatsnew/3.0.html, released 30/10/2020), which breaks its implementation in emg3d.
Needs fixing (when fixing adjust .travis.yml again).
Current work-around is to use h5py<3.0.
emg3d.utils.get_receiver takes a field (x-, y-, or z-directed, electric or magnetic), and returns the values at receiver location [x, y, z].
Create a thin wrapper for that, which takes an entire Field instance, and either [x, y, z, azimuth, dip] or [x1, x2, y1, y2, z1, z2] and returns the corresponding response.
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
Django
The Web framework for perfectionists with deadlines.
Laravel
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
Microsoft
Open source projects and samples from Microsoft.
Google
Google ❤️ Open Source for everyone.
Alibaba
Alibaba Open Source for everyone
D3
Data-Driven Documents codes.
Tencent
China tencent open source team.