IndexError when calling emg3d.construct_mesh
jingxu94 opened this issue · 6 comments
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
Thanks for reporting @jingxu94 . The entire construct_mesh()-functionality is very new, so unfortunately it is possible that there are bugs, even though I tried to avoid them ;-)
As I explained in our earlier discussion, I can currently not reproduce the problem, so I cannot do debugging myself. I will therefore ask some questions, it would be great if you could help me with that. Hopefully we will get to the bottom of it.
Could you try to run the following and post the output?
import emg3d
z0, hz = emg3d.meshes.get_origin_widths(
frequency=1.0,
min_width_limits=100.0,
properties=[0.3, 2, 100],
center=-1000,
seasurface=0.0,
domain=[-2500, 0],
verb=2,
)
Just to ensure if the issue is in construct_mesh or in get_origin_widths.
Similar error as follows:
Traceback (most recent call last):
File "c:/Users/JingXu/Desktop/code/start_emg3d/model_property_mapping/test_issue.py", line 9, in <module>
verb=2,
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
That is really strange, I get the following without any issues:
Skin depth [m] : 276 / 712 / 5033 [corresponding to `properties`]
Survey domain DS [m] : -2500 - 0
Comp. domain DC [m] : -6974 - 31623
Final extent [m] : -8729 - 32628
Cell widths [m] : 100 / 100 / 9085 [min(DS) / max(DS) / max(DC)]
Number of cells : 48 (25 / 23 / 0) [Total (DS/DC/remain)]
Max stretching : 1.000 (1.000) / 1.380 [DS (seasurface) / DC]
I even recreated an environment fixing all the package versions to match yours, but I do not encounter the problem. I also tried it with pip and with conda, both work fine. The only difference that remains between your system and mine is Windows/Linux.
I will try and get hold of a Windows machine to test it.
OK, I found the problem, it is a Windows issue. Int's default to 64 on Linux/Mac, but on Windows they default to 32. This caused an overflow problem in emg3d.meshes.good_mg_cell_nr.
If you change line 780 (emg3d v0.14.1) from
primes = np.array([2, 3, 5, 7, 11, 13, 17, 19])
to
primes = np.array([2, 3, 5, 7, 11, 13, 17, 19], dtype=np.int64)
the problem will be fixed.
I will push an update shortly.
Thanks for reporting!
Thinks, it works. That's a helpful tip for programming!