From 6c3c841a401ede0d4566bc71db5a779bbc3ffe61 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dieter=20Werthm=C3=BCller?= Date: Wed, 9 Feb 2022 15:19:29 +0100 Subject: [PATCH] v1.4.0 (#270) Allow vector's to be smaller than domain --- CHANGELOG.rst | 79 ++++++++++++++++++++++++-------------------- emg3d/meshes.py | 71 ++++++++++++++++++++++----------------- tests/test_meshes.py | 46 +++++++++++++++----------- 3 files changed, 111 insertions(+), 85 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 908e2615..eeea3447 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -6,28 +6,37 @@ Changelog """""""""" -*latest* will become v1.4.0 ---------------------------- +v1.4.0 : Meshing: improve vector +-------------------------------- + +**2022-02-09** + +- Meshes: Non-backwards compatible changes in ``construct_mesh`` + (``origin_and_widths``; ``estimate_gridding_options``) when providing + ``vector``'s (implemented non-backwards compatible as the old rules were not + intuitive nor logic; previous meshes can still be obtained, mostly, by + setting the parameters carefully). -- Meshes: Non-backwards compatible change in ``construct_mesh`` - (``origin_and_widths``; ``estimate_gridding_options``). It is implemented - non-backwards compatible as the old rules were not intuitive nor logic. The - previous meshes can still be obtained by setting the parameters carefully. - Priority-order changed to ``domain > distance > vector`` (before it was ``domain > vector > distance``). - - A ``vector`` is new cut to the corresponding domain, if ``domain`` or - ``distance`` was defined as well (cut at the first point where - ``vector <= domain[0]``, ``vector >= domain[1]``). + - A provided ``vector`` is new trimmed to the corresponding domain if it is + larger than a also provided domain (from ``domain`` or ``distance``); + trimmed at the first point where + ``vector <= domain[0]``, ``vector >= domain[1]``. + - A ``vector`` can new also be smaller than the defined domain, and the + domain is then filled according to the normal rules; the last cell of + ``vector`` in each direction is taken as starting width for the expansion. -- Removed functions and modules that were deprecated in v1.2.1. +- Bugfixes and maintenance: -- Bugfix when adding ``add_noise`` explicitly to ``Simulation.compute()``. + - Removed functions and modules that were deprecated in v1.2.1. + - Fixed kwargs-error when adding ``add_noise`` explicitly to + ``Simulation.compute()``. + - Python 3.10 added to tests; Python 3.7 tests reduced to minimum. -- Maintenance: Python 3.10 added to tests; Python 3.7 reduced to minimum. - -v1.3.2: Bugfix CLI-select -------------------------- +v1.3.2 : Bugfix CLI-select +-------------------------- **2021-12-01** @@ -35,8 +44,8 @@ CLI: Add ``remove_empty`` to parameter file; set to ``False`` by default (pre-v1.3.1 behaviour, and therefore backwards compatible). -v1.3.1: Select: remove empty pairs ----------------------------------- +v1.3.1 : Select: remove empty pairs +----------------------------------- **2021-11-20** @@ -46,8 +55,8 @@ v1.3.1: Select: remove empty pairs - Maintenance: Added a cron to GHA; 20th of every month at 14:14. -v1.3.0: File-based computations -------------------------------- +v1.3.0 : File-based computations +-------------------------------- **2021-10-27** @@ -81,8 +90,8 @@ v1.3.0: File-based computations -v1.2.1: Remove optimize & bug fix ---------------------------------- +v1.2.1 : Remove optimize & bug fix +---------------------------------- **2021-08-22** @@ -104,8 +113,8 @@ v1.2.1: Remove optimize & bug fix ``Simulation.{misfit;gradient}``. -v1.2.0: White noise -------------------- +v1.2.0 : White noise +-------------------- **2021-07-27** @@ -151,8 +160,8 @@ v1.2.0: White noise - Simplified (unified) ``_edge_curl_factor`` (private fct). -v1.1.0: Adjoint-fix for electric receivers ------------------------------------------- +v1.1.0 : Adjoint-fix for electric receivers +------------------------------------------- **2021-06-30** @@ -210,8 +219,8 @@ The changes in more detail: set in the parameter file, which overwrites the described default behaviour. -v1.0.0: Stable API ------------------- +v1.0.0 : Stable API +------------------- **2021-05-28** @@ -511,8 +520,8 @@ Detailed changes by module """""""""" -v0.17.0: Magnetics in Simulation --------------------------------- +v0.17.0 : Magnetics in Simulation +--------------------------------- **2021-03-03** @@ -544,8 +553,8 @@ v0.17.0: Magnetics in Simulation stuff. Currently only stores the data selection to output data. -v0.16.1: Verbosity & Logging ----------------------------- +v0.16.1 : Verbosity & Logging +----------------------------- **2021-02-09** @@ -577,8 +586,8 @@ v0.16.1: Verbosity & Logging above). -v0.16.0: Arbitrarily shaped sources ------------------------------------ +v0.16.0 : Arbitrarily shaped sources +------------------------------------ **2021-01-13** @@ -628,8 +637,8 @@ v0.16.0: Arbitrarily shaped sources ``map_coordinates``. -v0.15.3: Move to EMSiG ----------------------- +v0.15.3 : Move to EMSiG +----------------------- **2020-12-09** diff --git a/emg3d/meshes.py b/emg3d/meshes.py index edae4294..fabd94ea 100644 --- a/emg3d/meshes.py +++ b/emg3d/meshes.py @@ -371,11 +371,19 @@ def construct_mesh(frequency, properties, center, domain=None, vector=None, where ``center=(cx, cy, cz)``. vector : {tuple, ndarray, dict, None}, optional - Contains vectors of mesh-edges that should be used. If provided, the - vector *must* at least include all of the survey domain, to which - extent it will be cut. Format: ``(xvector, yvector, zvector)`` or + Contains vectors of mesh-edges that should be used. Format: + ``(xvector, yvector, zvector)`` or ``{'x': xvector, 'y': yvector, 'z': zvector}``. + If also ``domain`` or ``distance`` is defined, the following happens: + - There must be at least two cells within the domain, otherwise + ``vector`` is set to None. + - Where the ``vector`` is outside the domain, it will be trimmed to it. + - Where the ``vector`` is smaller than the domain, it will be extended + following the normal rules. The last cell-width in each direction + will be taken as starting cell width, together with the domain + stretching factor. + It can be None, or individual ndarrays can be None (e.g., ``(xvector, yvector, None)``), in which case you have to provide a ``domain`` or ``distance``. If only one ndarray is provided it is applied to all @@ -627,8 +635,7 @@ def origin_and_widths(frequency, properties, center, domain=None, vector=None, else: raise ValueError( - "At least one of `domain`, `distance`, " - "and `vector` must be provided." + "At least one of `domain`/`distance`/`vector` must be provided." ) # Check vector if provided @@ -643,12 +650,9 @@ def origin_and_widths(frequency, properties, center, domain=None, vector=None, if vmax.size > 1: vector = vector[:vmax[1]] - # Ensure vector spans whole domain. - if vmin.size == 0 or vmax.size == 0: - raise ValueError( - "Provided vector MUST at least include " - "all of the survey domain." - ) + # If vector is outside domain, set to None. + if len(vector) < 3: + vector = None # Seasurface related checks. if seasurface is not None: @@ -700,32 +704,37 @@ def origin_and_widths(frequency, properties, center, domain=None, vector=None, for sa in np.linspace(1.0, stretching[0], nsa): if vector is None: + cpart = [center, center] + dlr = [dmin, dmin] + nx0 = 0 + else: + cpart = [vector[0], vector[-1]] + dlr = [vector[1]-vector[0], vector[-1]-vector[-2]] + nx0 = 1 - # Get current stretched grid cell sizes. - thxl = dmin*sa**np.arange(nx) # Left of origin. - thxr = dmin*sa**np.arange(nx) # Right of origin. + # Get current stretched grid cell sizes. + thxl = dlr[0]*sa**np.arange(nx0, nx) # Left of origin. + thxr = dlr[1]*sa**np.arange(nx0, nx) # Right of origin. - # Adjust stretching for seasurface if required. - if seasurface is not None and seasurface > center: - t_nx = np.r_[center, center+np.cumsum(thxr)] - ii = np.argmin(abs(t_nx-seasurface)) - thxr[:ii] *= abs(seasurface-center)/np.sum(thxr[:ii]) + # Adjust stretching for seasurface if required. + if seasurface is not None and seasurface > cpart[1]: + t_nx = np.r_[cpart[1], cpart[1]+np.cumsum(thxr)] + ii = np.argmin(abs(t_nx-seasurface)) + thxr[:ii] *= abs(seasurface-cpart[1])/np.sum(thxr[:ii]) - # Fill from center to left and right domain. - nl = np.sum((center-np.cumsum(thxl)) > domain[0]) + 1 - nr = np.sum((center+np.cumsum(thxr)) < domain[1]) + 1 + # Fill from center to left and right domain. + nl = np.sum((cpart[0]-np.cumsum(thxl)) >= domain[0]) + (1-nx0) + nr = np.sum((cpart[1]+np.cumsum(thxr)) <= domain[1]) + (1-nx0) - # Create the current hx-array. + # Create the current hx-array. + if vector is None: hx = np.r_[thxl[:nl][::-1], thxr[:nr]] - - # Get actual domain: - asurv_domain = [center - np.sum(thxl[:nl]), - center + np.sum(thxr[:nr])] - else: - # Store actual domain, current hx-array, and number of cells. - asurv_domain = [vector[0], vector[-1]] - hx = np.diff(vector) + hx = np.r_[thxl[:nl][::-1], np.diff(vector), thxr[:nr]] + + # Get actual domain: + asurv_domain = [cpart[0] - np.sum(thxl[:nl]), + cpart[1] + np.sum(thxr[:nr])] # Expand for seasurface if necessary. if seasurface is not None and seasurface > asurv_domain[-1]: diff --git a/tests/test_meshes.py b/tests/test_meshes.py index d2b27738..aaa5d322 100644 --- a/tests/test_meshes.py +++ b/tests/test_meshes.py @@ -251,12 +251,9 @@ def test_errors(self, capsys): with pytest.raises(TypeError, match='Unexpected '): meshes.origin_and_widths(1, 1, 0, [-1, 1], unknown=True) - with pytest.raises(ValueError, match="At least one of `domain`, `d"): + with pytest.raises(ValueError, match="At least one of `domain`/`d"): meshes.origin_and_widths(1, 1, 0) - with pytest.raises(ValueError, match="Provided vector MUST at least"): - meshes.origin_and_widths(1, 1, 0, [-1, 1], np.array([0, 1, 2])) - with pytest.raises(ValueError, match="The `seasurface` must be bigge"): meshes.origin_and_widths(1, 1, 0, [-1, 1], seasurface=-2) @@ -312,24 +309,35 @@ def test_basics(self, capsys): def test_domain_vector(self): inp = {'frequency': 1/np.pi, 'properties': 9*mu_0, 'center': 0.0, - 'stretching': [1, 1]} + 'stretching': [1, 1.5]} + x01, hx1 = meshes.origin_and_widths(domain=[-1, 1], **inp) - x02, hx2 = meshes.origin_and_widths(vector=np.array([-1, 0, 1]), **inp) + + x01b, hx1b = meshes.origin_and_widths( # vector will be set to None + domain=[-1, 1], vector=[-2, -1, 0], **inp) + assert_allclose(x01, x01b) + assert_allclose(hx1, hx1b) + + vector2 = np.array([-1, 0, 1]) + x02, hx2 = meshes.origin_and_widths(vector=vector2, **inp) + assert np.in1d(vector2, x02 + np.cumsum(hx2)).all() + + vector3 = np.array([-2, -1, 0, 1, 2]) x03, hx3 = meshes.origin_and_widths( # vector will be cut - domain=[-1, 1], vector=np.array([-2, -1, 0, 1, 2]), **inp) - x04, hx4 = meshes.origin_and_widths( # vector will be cut - distance=[1, 1], vector=np.array([-2, -1, 0, 1, 2]), **inp) - assert_allclose(x01, x02) - assert_allclose(x01, x03) - assert_allclose(x01, x04) - assert_allclose(hx1, hx2) - assert_allclose(hx1, hx3) - assert_allclose(hx1, hx4) + domain=[-1, 1], vector=vector3, **inp) + assert np.in1d(vector3[1:-1], x03 + np.cumsum(hx3)).all() + assert not np.in1d(vector3[0], x03 + np.cumsum(hx3)) + assert not np.in1d(vector3[-1], x03 + np.cumsum(hx3)) - x03, hx3 = meshes.origin_and_widths( - 1/np.pi, 9*mu_0, 0.0, distance=[1, 1], stretching=[1, 1]) - assert_allclose(x01, x03) - assert_allclose(hx1, hx3) + x04, hx4 = meshes.origin_and_widths( # vector will be cut + distance=[1, 1], vector=vector3, **inp) + assert_allclose(x03, x04) + assert_allclose(hx3, hx4) + + x05, hx5 = meshes.origin_and_widths( # vector will be expanded + distance=[2.0, 2.0], vector=vector2, **inp) + print(x05 + np.cumsum(hx5)) + assert np.in1d(np.array([-2, -1, 0, 1, 2]), x05 + np.cumsum(hx5)).all() def test_seasurface(self): x01, hx1 = meshes.origin_and_widths(