diff --git a/.github/environment_ci.yml b/.github/environment_ci.yml index 9e03b20d..cc85be11 100644 --- a/.github/environment_ci.yml +++ b/.github/environment_ci.yml @@ -15,7 +15,6 @@ dependencies: - pydata-sphinx-theme==0.15.4 - sphinx-gallery>=0.1.13 - numpydoc>=1.5 - - discretize - jupyter - graphviz - pillow diff --git a/.github/workflows/test_with_conda.yml b/.github/workflows/test_with_conda.yml index 3de54d21..784f0339 100644 --- a/.github/workflows/test_with_conda.yml +++ b/.github/workflows/test_with_conda.yml @@ -18,10 +18,12 @@ jobs: fail-fast: false matrix: os: ["ubuntu-latest"] - python-version: ["3.10", "3.11", "3.12"] + python-version: ["3.10", "3.11", "3.12", "3.13"] steps: - uses: actions/checkout@v4 + with: + fetch-depth: 0 - name: Setup Conda uses: conda-incubator/setup-miniconda@v3 @@ -31,6 +33,25 @@ jobs: activate-environment: geoana-test python-version: ${{ matrix.python-version }} + - name: Install numba + if: matrix.python-version != '3.13' + # Numba doesn't work on python 3.13 just yet: + run: | + conda install --yes -c conda-forge numba + + # Install discretize from it's repo until wheels/conda-forge built against numpy2.0 available: + - name: Pull Discretize + uses: actions/checkout@v4 + with: + fetch-depth: 0 + repository: 'simpeg/discretize' + ref: 'main' + path: 'discretize' + + - name: Install discretize + run: | + pip install ./discretize + - name: Conda information run: | conda info @@ -39,11 +60,11 @@ jobs: - name: Install Our Package run: | - pip install . --config-settings=setup-args="-Dwith_extensions=true" + pip install --no-build-isolation --editable . --config-settings=setup-args="-Dwith_extensions=true" - name: Run Tests run: | - pytest --cov-config=.coveragerc --cov=geoana --cov-report=xml -s -v -W ignore::DeprecationWarning + pytest tests --cov-config=.coveragerc --cov=geoana --cov-report=xml -s -v -W ignore::DeprecationWarning - name: "Upload coverage to Codecov" if: matrix.python-version == '3.11' diff --git a/geoana/gravity.py b/geoana/gravity.py index 8d11e2e8..71dd2d6e 100644 --- a/geoana/gravity.py +++ b/geoana/gravity.py @@ -555,7 +555,7 @@ def gravitational_gradient(self, xyz): ind1 = r == self.radius g_tens[ind0] = super().gravitational_gradient(xyz[ind0]) g_tens[~ind0] = -G * 4 / 3 * np.pi * self.rho * np.eye(3) - g_tens[ind1] = np.NaN + g_tens[ind1] = np.nan return g_tens diff --git a/tests/test_em_fdem_electric_dipole.py b/tests/test_em_fdem_electric_dipole.py index 17cef136..237736e7 100644 --- a/tests/test_em_fdem_electric_dipole.py +++ b/tests/test_em_fdem_electric_dipole.py @@ -4,7 +4,6 @@ from scipy.constants import mu_0, epsilon_0 from geoana.em import fdem -import discretize # from SimPEG.EM import FDEM # from SimPEG import Maps @@ -31,7 +30,7 @@ def E_from_EDWS( epsilon = epsilon_0*epsr sig_hat = sig + 1j*fdem.omega(f)*epsilon - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) dx = XYZ[:, 0] - srcLoc[0] dy = XYZ[:, 1] - srcLoc[1] @@ -133,7 +132,7 @@ def electric_dipole_e(self, orientation): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) extest, eytest, eztest = E_from_EDWS( xyz, edws.location, edws.sigma, edws.frequency, @@ -170,7 +169,7 @@ def test_electric_dipole_tilted_e(self): y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) extest0, eytest0, eztest0 = E_from_EDWS( xyz, edws.location, edws.sigma, edws.frequency, diff --git a/tests/test_em_fdem_layered.py b/tests/test_em_fdem_layered.py index e9409c16..1f6bcf77 100644 --- a/tests/test_em_fdem_layered.py +++ b/tests/test_em_fdem_layered.py @@ -1,4 +1,6 @@ import unittest +from tabnanny import check + import pytest import numpy as np from numpy.testing import assert_allclose @@ -8,11 +10,13 @@ from geoana.kernels.tranverse_electric_reflections import ( rTE_forward, rTE_gradient, _rTE_forward, _rTE_gradient ) -import discretize from scipy.special import iv, kv from geoana.em.fdem.base import sigma_hat -from discretize.tests import check_derivative +try: + from discretize.tests import check_derivative +except ImportError: + check_derivative = None class TestHalfSpace: @@ -180,7 +184,7 @@ def test_magnetic_field_errors(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) mag_layer.magnetic_field(xyz) def test_magnetic_field(self): @@ -279,64 +283,65 @@ def test_magnetic_dipole(self): assert_allclose(em_layer_sec, em_half_sec, rtol=1e-05, atol=1e-08) -class TestrTEGradient(unittest.TestCase): - def test_rTE_jacobian(self): - """Test to make sure numpy and compiled give same results""" - n_layer = 11 - n_frequency = 5 - n_lambda = 8 - frequencies = np.logspace(1, 4, 5) - thicknesses = np.ones(n_layer-1) - lamb = np.logspace(0, 3, n_lambda) - np.random.seed(0) - sigma = 1E-1*(1 + 0.1 * np.random.rand(n_layer, n_frequency)) - mu = mu_0 * (1 + 0.1 * np.random.rand(n_layer, n_frequency)) - dmu = mu_0 * 0.1 * np.random.rand(n_layer, n_frequency) - - def rte_sigma(x): - sigma = x.reshape(n_layer, n_frequency) - rTE = rTE_forward(frequencies, lamb, sigma, mu, thicknesses) - - J_sigma, _, _ = rTE_gradient(frequencies, lamb, sigma, mu, thicknesses) - - def J(y): - y = y.reshape(n_layer, n_frequency) - # do summation over layers, broadcast over frequencies only - # equivalent to: - # out = np.empty_like(rTE) - # for i in range(n_frequency): - # out[i] = J_sigma[:, i, :]).T@y[:, i] - return np.einsum('i...k,i...', J_sigma, y) - return rTE, J - - def rte_h(x): - thicknesses = x - rTE = rTE_forward(frequencies, lamb, sigma, mu, thicknesses) - _, J_h, _ = rTE_gradient(frequencies, lamb, sigma, mu, thicknesses) - - def J(y): - #(J_h.T@y).T - out = np.einsum('i...,i', J_h, y) - # do sum over n_layer, broadcast over all others - return out - - return rTE, J - - def rte_mu(x): - mu = x.reshape(n_layer, n_frequency) - rTE = rTE_forward(frequencies, lamb, sigma, mu, thicknesses) - - _, _, J_mu = rTE_gradient(frequencies, lamb, sigma, mu, thicknesses) - - def J(y): - y = y.reshape(n_layer, n_frequency) - # do summation over layers, broadcast over frequencies only - return np.einsum('i...k,i...', J_mu, y) - return rTE, J - - self.assertTrue(check_derivative(rte_sigma, sigma.reshape(-1), num=4, plotIt=False)) - self.assertTrue(check_derivative(rte_h, thicknesses, num=4, plotIt=False)) - self.assertTrue(check_derivative(rte_mu, mu.reshape(-1), dx=dmu.reshape(-1), num=4, plotIt=False)) +if check_derivative is not None: + class TestrTEGradient(unittest.TestCase): + def test_rTE_jacobian(self): + """Test to make sure numpy and compiled give same results""" + n_layer = 11 + n_frequency = 5 + n_lambda = 8 + frequencies = np.logspace(1, 4, 5) + thicknesses = np.ones(n_layer-1) + lamb = np.logspace(0, 3, n_lambda) + np.random.seed(0) + sigma = 1E-1*(1 + 0.1 * np.random.rand(n_layer, n_frequency)) + mu = mu_0 * (1 + 0.1 * np.random.rand(n_layer, n_frequency)) + dmu = mu_0 * 0.1 * np.random.rand(n_layer, n_frequency) + + def rte_sigma(x): + sigma = x.reshape(n_layer, n_frequency) + rTE = rTE_forward(frequencies, lamb, sigma, mu, thicknesses) + + J_sigma, _, _ = rTE_gradient(frequencies, lamb, sigma, mu, thicknesses) + + def J(y): + y = y.reshape(n_layer, n_frequency) + # do summation over layers, broadcast over frequencies only + # equivalent to: + # out = np.empty_like(rTE) + # for i in range(n_frequency): + # out[i] = J_sigma[:, i, :]).T@y[:, i] + return np.einsum('i...k,i...', J_sigma, y) + return rTE, J + + def rte_h(x): + thicknesses = x + rTE = rTE_forward(frequencies, lamb, sigma, mu, thicknesses) + _, J_h, _ = rTE_gradient(frequencies, lamb, sigma, mu, thicknesses) + + def J(y): + #(J_h.T@y).T + out = np.einsum('i...,i', J_h, y) + # do sum over n_layer, broadcast over all others + return out + + return rTE, J + + def rte_mu(x): + mu = x.reshape(n_layer, n_frequency) + rTE = rTE_forward(frequencies, lamb, sigma, mu, thicknesses) + + _, _, J_mu = rTE_gradient(frequencies, lamb, sigma, mu, thicknesses) + + def J(y): + y = y.reshape(n_layer, n_frequency) + # do summation over layers, broadcast over frequencies only + return np.einsum('i...k,i...', J_mu, y) + return rTE, J + + self.assertTrue(check_derivative(rte_sigma, sigma.reshape(-1), num=4, plotIt=False)) + self.assertTrue(check_derivative(rte_h, thicknesses, num=4, plotIt=False)) + self.assertTrue(check_derivative(rte_mu, mu.reshape(-1), dx=dmu.reshape(-1), num=4, plotIt=False)) class TestCompiledVsNumpy(unittest.TestCase): diff --git a/tests/test_em_fdem_magnetic_dipole.py b/tests/test_em_fdem_magnetic_dipole.py index 5db21d5a..48e8757e 100644 --- a/tests/test_em_fdem_magnetic_dipole.py +++ b/tests/test_em_fdem_magnetic_dipole.py @@ -1,7 +1,5 @@ import unittest import numpy as np -import discretize - from scipy.constants import mu_0, epsilon_0 from geoana.em import fdem @@ -26,7 +24,7 @@ def H_from_MagneticDipoleWholeSpace( omega = lambda f: 2*np.pi*f - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) # Check dx = XYZ[:, 0]-srcLoc[0] @@ -88,7 +86,7 @@ def E_from_MagneticDipoleWholeSpace( omega = lambda f: 2 * np.pi * f - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) dx = XYZ[:, 0]-srcLoc[0] dy = XYZ[:, 1]-srcLoc[1] @@ -188,7 +186,7 @@ def magnetic_dipole_b(self, orientation): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) # srcLoc, obsLoc, component, orientation='Z', moment=1., mu=mu_0 @@ -221,7 +219,7 @@ def magnetic_dipole_e(self, orientation): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) extest, eytest, eztest = E_from_MagneticDipoleWholeSpace( xyz, mdws.location, mdws.sigma, mdws.frequency, @@ -258,7 +256,7 @@ def test_magnetic_dipole_tilted_b(self): y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) bxtest0, bytest0, bztest0 = B_from_MagneticDipoleWholeSpace( xyz, mdws.location, mdws.sigma, mdws.frequency, @@ -312,7 +310,7 @@ def test_magnetic_dipole_tilted_e(self): y = np.linspace(-30., 30., 10) z = np.linspace(-40., 40., 10) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) extest0, eytest0, eztest0 = E_from_MagneticDipoleWholeSpace( xyz, mdws.location, mdws.sigma, mdws.frequency, diff --git a/tests/test_em_fdem_planewave.py b/tests/test_em_fdem_planewave.py index fc1386a8..03ad0fb3 100644 --- a/tests/test_em_fdem_planewave.py +++ b/tests/test_em_fdem_planewave.py @@ -1,6 +1,5 @@ import pytest import numpy as np -import discretize from geoana.em import fdem from scipy.constants import mu_0, epsilon_0 @@ -48,7 +47,7 @@ def test_electric_field(): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) z = xyz[:, 2] kz = np.outer(k, z) @@ -90,7 +89,7 @@ def test_current_density(): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) z = xyz[:, 2] kz = np.outer(k, z) @@ -130,7 +129,7 @@ def test_magnetic_field(): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) z = xyz[:, 2] kz = z * k @@ -171,7 +170,7 @@ def test_magnetic_flux_density(): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) z = xyz[:, 2] kz = z * k diff --git a/tests/test_em_static.py b/tests/test_em_static.py index 8caac7c0..96a0533f 100644 --- a/tests/test_em_static.py +++ b/tests/test_em_static.py @@ -2,7 +2,10 @@ import unittest import numpy as np from scipy.constants import mu_0, epsilon_0 -import discretize +try: + import discretize +except ImportError: + discretize = None from scipy.special import ellipk, ellipe from geoana.em import static, fdem @@ -72,7 +75,7 @@ def test_errors(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) self.clws.magnetic_flux_density(xyz, coordinates='square') with pytest.raises(ValueError): self.lcfs.nodes = [0, 1, 2, 3, 4, 5] @@ -85,26 +88,21 @@ def test_errors(self): def test_vector_potential(self): n = 50 - mesh = discretize.TensorMesh( - [np.ones(n), np.ones(n), np.ones(n)], x0="CCC" - ) - - # radius = 1.0 - # self.clws.radius = radius - # self.clws.current = 1./(np.pi * radius**2) + x = y = z = np.linspace(-25, 25, n) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) for orientation in ["x", "y", "z"]: self.clws.orientation = orientation self.mdws.orientation = orientation inds = ( - (np.absolute(mesh.cell_centers[:, 0]) > 5) & - (np.absolute(mesh.cell_centers[:, 1]) > 5) & - (np.absolute(mesh.cell_centers[:, 2]) > 5) + (np.absolute(xyz[:, 0]) > 5) & + (np.absolute(xyz[:, 1]) > 5) & + (np.absolute(xyz[:, 2]) > 5) ) - a_clws = self.clws.vector_potential(mesh.cell_centers)[inds] - a_mdws = self.mdws.vector_potential(mesh.cell_centers)[inds] + a_clws = self.clws.vector_potential(xyz)[inds] + a_mdws = self.mdws.vector_potential(xyz)[inds] self.assertTrue(isinstance(a_clws, np.ndarray)) self.assertTrue(isinstance(a_mdws, np.ndarray)) @@ -114,128 +112,130 @@ def test_vector_potential(self): 0.5 * TOL * (np.linalg.norm(a_clws) + np.linalg.norm(a_mdws)) ) - def test_magnetic_field_tensor(self): - print("\n === Testing Tensor Mesh === \n") - n = 30 - h = 2. - mesh = discretize.TensorMesh( - [h*np.ones(n), h*np.ones(n), h*np.ones(n)], x0="CCC" - ) + if discretize is not None: + def test_magnetic_field_tensor(self): + print("\n === Testing Tensor Mesh === \n") + n = 30 + h = 2. + mesh = discretize.TensorMesh( + [np.full(n, h), np.full(n, h), np.full(n, h)], + origin='CCC' + ) - for radius in [0.5, 1, 1.5]: - self.clws.radius = radius - self.clws.current = 1./(np.pi * radius**2) - - fdem_dipole = fdem.MagneticDipoleWholeSpace(frequency=0) - - for location in [ - np.r_[0, 0, 0], np.r_[10, 0, 0], np.r_[0, -10, 0], - np.r_[0, 0, 10], np.r_[10, 10, 10] - ]: - self.clws.location = location - self.mdws.location = location - fdem_dipole.location = location - - for orientation in ["x", "y", "z"]: - self.clws.orientation = orientation - self.mdws.orientation = orientation - fdem_dipole.orientation = orientation - - a_clws = np.hstack([ - self.clws.vector_potential(mesh.gridEx)[:, 0], - self.clws.vector_potential(mesh.gridEy)[:, 1], - self.clws.vector_potential(mesh.gridEz)[:, 2] - ]) - - a_mdws = np.hstack([ - self.mdws.vector_potential(mesh.gridEx)[:, 0], - self.mdws.vector_potential(mesh.gridEy)[:, 1], - self.mdws.vector_potential(mesh.gridEz)[:, 2] - ]) - - b_clws = mesh.edge_curl * a_clws - b_mdws = mesh.edge_curl * a_mdws - - b_clws_true = np.hstack([ - self.clws.magnetic_flux_density(mesh.gridFx)[:, 0], - self.clws.magnetic_flux_density(mesh.gridFy)[:, 1], - self.clws.magnetic_flux_density(mesh.gridFz)[:, 2] - ]) - b_mdws_true = np.hstack([ - self.mdws.magnetic_flux_density(mesh.gridFx)[:, 0], - self.mdws.magnetic_flux_density(mesh.gridFy)[:, 1], - self.mdws.magnetic_flux_density(mesh.gridFz)[:, 2] - ]) - b_fdem = np.hstack([ - fdem_dipole.magnetic_flux_density(mesh.gridFx)[:, 0], - fdem_dipole.magnetic_flux_density(mesh.gridFy)[:, 1], - fdem_dipole.magnetic_flux_density(mesh.gridFz)[:, 2] - ]) - - self.assertTrue(isinstance(b_fdem, np.ndarray)) - - inds = (np.hstack([ - (np.absolute(mesh.gridFx[:, 0]) > h*2 + location[0]) & - (np.absolute(mesh.gridFx[:, 1]) > h*2 + location[1]) & - (np.absolute(mesh.gridFx[:, 2]) > h*2 + location[2]), - (np.absolute(mesh.gridFy[:, 0]) > h*2 + location[0]) & - (np.absolute(mesh.gridFy[:, 1]) > h*2 + location[1]) & - (np.absolute(mesh.gridFy[:, 2]) > h*2 + location[2]), - (np.absolute(mesh.gridFz[:, 0]) > h*2 + location[0]) & - (np.absolute(mesh.gridFz[:, 1]) > h*2 + location[1]) & - (np.absolute(mesh.gridFz[:, 2]) > h*2 + location[2]) - ])) - - loop_passed_1 = ( - np.linalg.norm(b_fdem[inds] - b_clws[inds]) < - 0.5 * TOL * ( - np.linalg.norm(b_fdem[inds]) + - np.linalg.norm(b_clws[inds]) + for radius in [0.5, 1, 1.5]: + self.clws.radius = radius + self.clws.current = 1./(np.pi * radius**2) + + fdem_dipole = fdem.MagneticDipoleWholeSpace(frequency=0) + + for location in [ + np.r_[0, 0, 0], np.r_[10, 0, 0], np.r_[0, -10, 0], + np.r_[0, 0, 10], np.r_[10, 10, 10] + ]: + self.clws.location = location + self.mdws.location = location + fdem_dipole.location = location + + for orientation in ["x", "y", "z"]: + self.clws.orientation = orientation + self.mdws.orientation = orientation + fdem_dipole.orientation = orientation + + a_clws = np.hstack([ + self.clws.vector_potential(mesh.gridEx)[:, 0], + self.clws.vector_potential(mesh.gridEy)[:, 1], + self.clws.vector_potential(mesh.gridEz)[:, 2] + ]) + + a_mdws = np.hstack([ + self.mdws.vector_potential(mesh.gridEx)[:, 0], + self.mdws.vector_potential(mesh.gridEy)[:, 1], + self.mdws.vector_potential(mesh.gridEz)[:, 2] + ]) + + b_clws = mesh.edge_curl * a_clws + b_mdws = mesh.edge_curl * a_mdws + + b_clws_true = np.hstack([ + self.clws.magnetic_flux_density(mesh.gridFx)[:, 0], + self.clws.magnetic_flux_density(mesh.gridFy)[:, 1], + self.clws.magnetic_flux_density(mesh.gridFz)[:, 2] + ]) + b_mdws_true = np.hstack([ + self.mdws.magnetic_flux_density(mesh.gridFx)[:, 0], + self.mdws.magnetic_flux_density(mesh.gridFy)[:, 1], + self.mdws.magnetic_flux_density(mesh.gridFz)[:, 2] + ]) + b_fdem = np.hstack([ + fdem_dipole.magnetic_flux_density(mesh.gridFx)[:, 0], + fdem_dipole.magnetic_flux_density(mesh.gridFy)[:, 1], + fdem_dipole.magnetic_flux_density(mesh.gridFz)[:, 2] + ]) + + self.assertTrue(isinstance(b_fdem, np.ndarray)) + + inds = (np.hstack([ + (np.absolute(mesh.gridFx[:, 0]) > h*2 + location[0]) & + (np.absolute(mesh.gridFx[:, 1]) > h*2 + location[1]) & + (np.absolute(mesh.gridFx[:, 2]) > h*2 + location[2]), + (np.absolute(mesh.gridFy[:, 0]) > h*2 + location[0]) & + (np.absolute(mesh.gridFy[:, 1]) > h*2 + location[1]) & + (np.absolute(mesh.gridFy[:, 2]) > h*2 + location[2]), + (np.absolute(mesh.gridFz[:, 0]) > h*2 + location[0]) & + (np.absolute(mesh.gridFz[:, 1]) > h*2 + location[1]) & + (np.absolute(mesh.gridFz[:, 2]) > h*2 + location[2]) + ])) + + loop_passed_1 = ( + np.linalg.norm(b_fdem[inds] - b_clws[inds]) < + 0.5 * TOL * ( + np.linalg.norm(b_fdem[inds]) + + np.linalg.norm(b_clws[inds]) + ) ) - ) - loop_passed_2 = ( - np.linalg.norm(b_fdem[inds] - b_clws_true[inds]) < - 0.5 * TOL * ( - np.linalg.norm(b_fdem[inds]) + - np.linalg.norm(b_clws_true[inds]) + loop_passed_2 = ( + np.linalg.norm(b_fdem[inds] - b_clws_true[inds]) < + 0.5 * TOL * ( + np.linalg.norm(b_fdem[inds]) + + np.linalg.norm(b_clws_true[inds]) + ) ) - ) - dipole_passed = ( - np.linalg.norm(b_fdem[inds] - b_mdws[inds]) < - 0.5 * TOL * ( - np.linalg.norm(b_fdem[inds]) + - np.linalg.norm(b_mdws[inds]) + dipole_passed = ( + np.linalg.norm(b_fdem[inds] - b_mdws[inds]) < + 0.5 * TOL * ( + np.linalg.norm(b_fdem[inds]) + + np.linalg.norm(b_mdws[inds]) + ) ) - ) - print( - "Testing r = {}, loc = {}, orientation = {}".format( - radius, location, orientation + print( + "Testing r = {}, loc = {}, orientation = {}".format( + radius, location, orientation + ) ) - ) - print( - " fdem: {:1.4e}, loop: {:1.4e}, dipole: {:1.4e}" - " Passed? loop: {}, dipole: {} \n".format( - np.linalg.norm(b_fdem[inds]), - np.linalg.norm(b_clws[inds]), - np.linalg.norm(b_mdws[inds]), - loop_passed_1, - loop_passed_2, - dipole_passed + print( + " fdem: {:1.4e}, loop: {:1.4e}, dipole: {:1.4e}" + " Passed? loop: {}, dipole: {} \n".format( + np.linalg.norm(b_fdem[inds]), + np.linalg.norm(b_clws[inds]), + np.linalg.norm(b_mdws[inds]), + loop_passed_1, + loop_passed_2, + dipole_passed + ) ) - ) - self.assertTrue(loop_passed_1) - self.assertTrue(loop_passed_2) - self.assertTrue(dipole_passed) + self.assertTrue(loop_passed_1) + self.assertTrue(loop_passed_2) + self.assertTrue(dipole_passed) def test_dipole_magnetic_flux_density(self): mdws = static.MagneticDipoleWholeSpace() x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) n_obs = xyz.shape[0] xyz_ = spatial.cylindrical_2_cartesian(xyz) @@ -256,7 +256,7 @@ def test_pole_magnetic_flux_density(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) xyz_ = spatial.cylindrical_2_cartesian(xyz) r = mpws.vector_distance(xyz_) @@ -273,7 +273,7 @@ def test_pole_magnetic_field(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) xyz_ = spatial.cylindrical_2_cartesian(xyz) r = mpws.vector_distance(xyz_) @@ -291,7 +291,7 @@ def test_loop_magnetic_flux_density(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) xyz_ = np.atleast_2d(xyz) xyz_ = spatial.cylindrical_2_cartesian(xyz_) @@ -320,132 +320,133 @@ def test_loop_magnetic_flux_density(self): b_test = clws.magnetic_flux_density(xyz, coordinates='cylindrical') np.testing.assert_equal(b_test, b) - def test_magnetic_field_3Dcyl(self): - print("\n === Testing 3D Cyl Mesh === \n") - n = 50 - ny = 10 - h = 2. - mesh = discretize.CylindricalMesh( - [h*np.ones(n), np.ones(ny) * 2 * np.pi / ny, h*np.ones(n)], x0="00C" - ) + if discretize is not None: + def test_magnetic_field_3Dcyl(self): + print("\n === Testing 3D Cyl Mesh === \n") + n = 50 + ny = 10 + h = 2. + mesh = discretize.CylindricalMesh( + [h*np.ones(n), np.ones(ny) * 2 * np.pi / ny, h*np.ones(n)], x0="00C" + ) - for radius in [0.5, 1, 1.5]: - self.clws.radius = radius - self.clws.current = 1./(np.pi * radius**2) # So that all have moment of 1.0 - - fdem_dipole = fdem.MagneticDipoleWholeSpace(frequency=0, quasistatic=True) - - for location in [ - np.r_[0, 0, 0], np.r_[0, 4, 0], np.r_[4, 4., 0], - np.r_[4, 4, 4], np.r_[4, 0, -4] - ]: - self.clws.location = location - self.mdws.location = location - fdem_dipole.location = location - - for orientation in ["z"]: - self.clws.orientation = orientation - self.mdws.orientation = orientation - fdem_dipole.orientation = orientation - - a_clws = np.hstack([ - self.clws.vector_potential( - mesh.gridEx, coordinates="cylindrical" - )[:, 0], - self.clws.vector_potential( - mesh.gridEy, coordinates="cylindrical" - )[:, 1], - self.clws.vector_potential( - mesh.gridEz, coordinates="cylindrical" - )[:, 2] - ]) - - a_mdws = np.hstack([ - self.mdws.vector_potential( - mesh.gridEx, coordinates="cylindrical" - )[:, 0], - self.mdws.vector_potential( - mesh.gridEy, coordinates="cylindrical" - )[:, 1], - self.mdws.vector_potential( - mesh.gridEz, coordinates="cylindrical" - )[:, 2] - ]) - - b_clws = mesh.edge_curl * a_clws - b_mdws = mesh.edge_curl * a_mdws - - b_fdem = np.hstack([ - spatial.cartesian_2_cylindrical( - spatial.cylindrical_2_cartesian(mesh.gridFx), - fdem_dipole.magnetic_flux_density( - spatial.cylindrical_2_cartesian(mesh.gridFx) - ) - )[:, 0], - spatial.cartesian_2_cylindrical( - spatial.cylindrical_2_cartesian(mesh.gridFy), - fdem_dipole.magnetic_flux_density( - spatial.cylindrical_2_cartesian(mesh.gridFy) + for radius in [0.5, 1, 1.5]: + self.clws.radius = radius + self.clws.current = 1./(np.pi * radius**2) # So that all have moment of 1.0 + + fdem_dipole = fdem.MagneticDipoleWholeSpace(frequency=0, quasistatic=True) + + for location in [ + np.r_[0, 0, 0], np.r_[0, 4, 0], np.r_[4, 4., 0], + np.r_[4, 4, 4], np.r_[4, 0, -4] + ]: + self.clws.location = location + self.mdws.location = location + fdem_dipole.location = location + + for orientation in ["z"]: + self.clws.orientation = orientation + self.mdws.orientation = orientation + fdem_dipole.orientation = orientation + + a_clws = np.hstack([ + self.clws.vector_potential( + mesh.gridEx, coordinates="cylindrical" + )[:, 0], + self.clws.vector_potential( + mesh.gridEy, coordinates="cylindrical" + )[:, 1], + self.clws.vector_potential( + mesh.gridEz, coordinates="cylindrical" + )[:, 2] + ]) + + a_mdws = np.hstack([ + self.mdws.vector_potential( + mesh.gridEx, coordinates="cylindrical" + )[:, 0], + self.mdws.vector_potential( + mesh.gridEy, coordinates="cylindrical" + )[:, 1], + self.mdws.vector_potential( + mesh.gridEz, coordinates="cylindrical" + )[:, 2] + ]) + + b_clws = mesh.edge_curl * a_clws + b_mdws = mesh.edge_curl * a_mdws + + b_fdem = np.hstack([ + spatial.cartesian_2_cylindrical( + spatial.cylindrical_2_cartesian(mesh.gridFx), + fdem_dipole.magnetic_flux_density( + spatial.cylindrical_2_cartesian(mesh.gridFx) + ) + )[:, 0], + spatial.cartesian_2_cylindrical( + spatial.cylindrical_2_cartesian(mesh.gridFy), + fdem_dipole.magnetic_flux_density( + spatial.cylindrical_2_cartesian(mesh.gridFy) + ) + )[:, 1], + spatial.cartesian_2_cylindrical( + spatial.cylindrical_2_cartesian(mesh.gridFz), + fdem_dipole.magnetic_flux_density( + spatial.cylindrical_2_cartesian(mesh.gridFz) + ) + )[:, 2] + ]) + + + inds = (np.hstack([ + (np.absolute(mesh.gridFx[:, 0]) > h*4 + location[0]) & + (np.absolute(mesh.gridFx[:, 2]) > h*4 + location[2]), + (np.absolute(mesh.gridFy[:, 0]) > h*4 + location[0]) & + (np.absolute(mesh.gridFy[:, 2]) > h*4 + location[2]), + (np.absolute(mesh.gridFz[:, 0]) > h*4 + location[0]) & + (np.absolute(mesh.gridFz[:, 2]) > h*4 + location[2]) + ])) + + loop_passed = ( + np.linalg.norm(b_fdem[inds] - b_clws[inds]) < + 0.5 * TOL * ( + np.linalg.norm(b_fdem[inds]) + + np.linalg.norm(b_clws[inds]) ) - )[:, 1], - spatial.cartesian_2_cylindrical( - spatial.cylindrical_2_cartesian(mesh.gridFz), - fdem_dipole.magnetic_flux_density( - spatial.cylindrical_2_cartesian(mesh.gridFz) - ) - )[:, 2] - ]) - - - inds = (np.hstack([ - (np.absolute(mesh.gridFx[:, 0]) > h*4 + location[0]) & - (np.absolute(mesh.gridFx[:, 2]) > h*4 + location[2]), - (np.absolute(mesh.gridFy[:, 0]) > h*4 + location[0]) & - (np.absolute(mesh.gridFy[:, 2]) > h*4 + location[2]), - (np.absolute(mesh.gridFz[:, 0]) > h*4 + location[0]) & - (np.absolute(mesh.gridFz[:, 2]) > h*4 + location[2]) - ])) - - loop_passed = ( - np.linalg.norm(b_fdem[inds] - b_clws[inds]) < - 0.5 * TOL * ( - np.linalg.norm(b_fdem[inds]) + - np.linalg.norm(b_clws[inds]) ) - ) - dipole_passed = ( - np.linalg.norm(b_fdem[inds] - b_mdws[inds]) < - 0.5 * TOL * ( - np.linalg.norm(b_fdem[inds]) + - np.linalg.norm(b_mdws[inds]) + dipole_passed = ( + np.linalg.norm(b_fdem[inds] - b_mdws[inds]) < + 0.5 * TOL * ( + np.linalg.norm(b_fdem[inds]) + + np.linalg.norm(b_mdws[inds]) + ) ) - ) - print( - "Testing r = {}, loc = {}, orientation = {}".format( - radius, location, orientation + print( + "Testing r = {}, loc = {}, orientation = {}".format( + radius, location, orientation + ) ) - ) - print( - " fdem: {:1.4e}, loop: {:1.4e}, dipole: {:1.4e}" - " Passed? loop: {}, dipole: {} \n".format( - np.linalg.norm(b_fdem[inds]), - np.linalg.norm(b_clws[inds]), - np.linalg.norm(b_mdws[inds]), - loop_passed, - dipole_passed + print( + " fdem: {:1.4e}, loop: {:1.4e}, dipole: {:1.4e}" + " Passed? loop: {}, dipole: {} \n".format( + np.linalg.norm(b_fdem[inds]), + np.linalg.norm(b_clws[inds]), + np.linalg.norm(b_mdws[inds]), + loop_passed, + dipole_passed + ) ) - ) - self.assertTrue(loop_passed) - self.assertTrue(dipole_passed) + self.assertTrue(loop_passed) + self.assertTrue(dipole_passed) def Vt_from_ESphere( XYZ, loc, sig_s, sig_b, radius, amp ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) sig_cur = (sig_s - sig_b) / (sig_s + 2 * sig_b) r_vec = XYZ - loc @@ -462,7 +463,7 @@ def Vt_from_ESphere( def Vp_from_ESphere( XYZ, loc, sig_s, sig_b, radius, amp ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) r_vec = XYZ - loc x = r_vec[:, 0] @@ -483,7 +484,7 @@ def Et_from_ESphere( XYZ, loc, sig_s, sig_b, radius, amp ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) sig_cur = (sig_s - sig_b) / (sig_s + 2 * sig_b) r_vec = XYZ - loc @@ -505,7 +506,7 @@ def Et_from_ESphere( def Ep_from_ESphere( XYZ, loc, sig_s, sig_b, radius, amp ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) r_vec = XYZ - loc r = np.linalg.norm(r_vec, axis=-1) @@ -527,7 +528,7 @@ def Jt_from_ESphere( XYZ, loc, sig_s, sig_b, radius, amp ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) r_vec = XYZ - loc r = np.linalg.norm(r_vec, axis=-1) @@ -557,7 +558,7 @@ def Js_from_ESphere( def rho_from_ESphere( XYZ, dr, loc, sig_s, sig_b, radius, amp ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) sig_cur = (sig_s - sig_b) / (sig_s + 2 * sig_b) r_vec = XYZ - loc @@ -625,7 +626,7 @@ def testV(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) vttest = Vt_from_ESphere( xyz, ess.location, ess.sigma_sphere, ess.sigma_background, ess.radius, ess.primary_field @@ -668,7 +669,7 @@ def testE(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) ettest = Et_from_ESphere( xyz, ess.location, ess.sigma_sphere, ess.sigma_background, ess.radius, ess.primary_field @@ -711,7 +712,7 @@ def testJ(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) jttest = Jt_from_ESphere( xyz, ess.location, ess.sigma_sphere, ess.sigma_background, ess.radius, ess.primary_field @@ -754,7 +755,7 @@ def test_rho(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) rho_test = rho_from_ESphere( xyz, None, ess.location, ess.sigma_sphere, ess.sigma_background, ess.radius, ess.primary_field @@ -768,7 +769,7 @@ def Vt_from_Sphere( XYZ, loc, mu_s, mu_b, radius, amp ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) mu_cur = (mu_s - mu_b) / (mu_s + 2 * mu_b) r_vec = XYZ - loc @@ -785,7 +786,7 @@ def Vt_from_Sphere( def Vp_from_Sphere( XYZ, loc, mu_s, mu_b, radius, amp ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) r_vec = XYZ - loc x = r_vec[:, 0] @@ -806,7 +807,7 @@ def Ht_from_Sphere( XYZ, loc, mu_s, mu_b, radius, amp ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) mu_cur = (mu_s - mu_b) / (mu_s + 2 * mu_b) r_vec = XYZ - loc @@ -828,7 +829,7 @@ def Ht_from_Sphere( def Hp_from_Sphere( XYZ, loc, mu_s, mu_b, radius, amp ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) r_vec = XYZ - loc x = r_vec[:, 0] @@ -850,7 +851,7 @@ def Bt_from_Sphere( XYZ, loc, mu_s, mu_b, radius, amp ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) r_vec = XYZ - loc r = np.linalg.norm(r_vec, axis=-1) @@ -928,7 +929,7 @@ def testV(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) vttest = Vt_from_Sphere( xyz, mss.location, mss.mu_sphere, mss.mu_background, mss.radius, mss.primary_field @@ -971,7 +972,7 @@ def testH(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) httest = Ht_from_Sphere( xyz, mss.location, mss.mu_sphere, mss.mu_background, mss.radius, mss.primary_field @@ -1014,7 +1015,7 @@ def testB(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) btest = Bt_from_Sphere( xyz, mss.location, mss.mu_sphere, mss.mu_background, mss.radius, mss.primary_field @@ -1046,7 +1047,7 @@ def V_from_PointCurrentW( XYZ, loc, rho, cur ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) r_vec = XYZ - loc r = np.linalg.norm(r_vec, axis=-1) @@ -1067,7 +1068,7 @@ def E_from_PointCurrentW( XYZ, loc, rho, cur ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) r_vec = XYZ - loc r = np.linalg.norm(r_vec, axis=-1) @@ -1113,7 +1114,7 @@ def test_potential(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) vtest = V_from_PointCurrentW( xyz, pcws.location, pcws.rho, pcws.current @@ -1137,7 +1138,7 @@ def test_current_density(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) jtest = J_from_PointCurrentW( xyz, pcws.location, pcws.rho, pcws.current @@ -1161,7 +1162,7 @@ def test_electric_field(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) etest = E_from_PointCurrentW( xyz, pcws.location, pcws.rho, pcws.current @@ -1178,7 +1179,7 @@ def V_from_PointCurrentH( XYZ, loc, rho, cur ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) r_vec = XYZ - loc r = np.linalg.norm(r_vec, axis=-1) @@ -1191,7 +1192,7 @@ def E_from_PointCurrentH( XYZ, loc, rho, cur ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) r_vec = XYZ - loc r = np.linalg.norm(r_vec, axis=-1) @@ -1268,7 +1269,7 @@ def test_potential(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 0., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) vtest = V_from_PointCurrentH( xyz, pchs.location, pchs.rho, pchs.current @@ -1283,7 +1284,7 @@ def test_potential(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) with pytest.raises(ValueError): pchs.potential(xyz) @@ -1300,7 +1301,7 @@ def test_electric_field(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 0., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) etest = E_from_PointCurrentH( xyz, pchs.location, pchs.rho, pchs.current @@ -1315,7 +1316,7 @@ def test_electric_field(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) with pytest.raises(ValueError): pchs.electric_field(xyz) @@ -1332,7 +1333,7 @@ def test_current_density(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 0., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) jtest = J_from_PointCurrentH( xyz, pchs.location, pchs.rho, pchs.current @@ -1347,7 +1348,7 @@ def test_current_density(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) with pytest.raises(ValueError): pchs.current_density(xyz) @@ -1356,7 +1357,7 @@ def test_current_density(self): def V_from_Dipole1( XYZ_M, XYZ_N, rho, cur, loc_a, loc_b ): - XYZ_M = discretize.utils.as_array_n_by_dim(XYZ_M, 3) + XYZ_M = np.atleast_2d(XYZ_M) r_vec1 = XYZ_M - loc_a r_vec2 = XYZ_M - loc_b @@ -1370,8 +1371,8 @@ def V_from_Dipole1( def V_from_Dipole2( XYZ_M, XYZ_N, rho, cur, loc_a, loc_b ): - XYZ_M = discretize.utils.as_array_n_by_dim(XYZ_M, 3) - XYZ_N = discretize.utils.as_array_n_by_dim(XYZ_N, 3) + XYZ_M = np.atleast_2d(XYZ_M) + XYZ_N = np.atleast_2d(XYZ_N) r_vec1 = XYZ_M - loc_a r_vec2 = XYZ_M - loc_b @@ -1392,7 +1393,7 @@ def V_from_Dipole2( def E_from_Dipole1( XYZ_M, XYZ_N, rho, cur, loc_a, loc_b ): - XYZ_M = discretize.utils.as_array_n_by_dim(XYZ_M, 3) + XYZ_M = np.atleast_2d(XYZ_M) r_vec1 = XYZ_M - loc_a r_vec2 = XYZ_M - loc_b @@ -1406,8 +1407,8 @@ def E_from_Dipole1( def E_from_Dipole2( XYZ_M, XYZ_N, rho, cur, loc_a, loc_b ): - XYZ_M = discretize.utils.as_array_n_by_dim(XYZ_M, 3) - XYZ_N = discretize.utils.as_array_n_by_dim(XYZ_N, 3) + XYZ_M = np.atleast_2d(XYZ_M) + XYZ_N = np.atleast_2d(XYZ_N) r_vec1 = XYZ_M - loc_a r_vec2 = XYZ_M - loc_b @@ -1498,12 +1499,12 @@ def test_potential(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 0., 50) - xyz1 = discretize.utils.ndgrid([x, y, z]) + xyz1 = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) x = np.linspace(-30., 20., 50) y = np.linspace(-20., 30., 50) z = np.linspace(-30., 0., 50) - xyz2 = discretize.utils.ndgrid([x, y, z]) + xyz2 = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) vtest1 = V_from_Dipole1( xyz1, None, dhs.rho, dhs.current, dhs.location_a, dhs.location_b @@ -1525,8 +1526,8 @@ def test_potential(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz3 = discretize.utils.ndgrid([x, y, z]) - xyz4 = discretize.utils.ndgrid([x, y, z]) + xyz3 = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) + xyz4 = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) with pytest.raises(ValueError): dhs.potential(xyz3) @@ -1538,12 +1539,12 @@ def test_electric_field(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 0., 50) - xyz1 = discretize.utils.ndgrid([x, y, z]) + xyz1 = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) x = np.linspace(-30., 20., 50) y = np.linspace(-20., 30., 50) z = np.linspace(-30., 0., 50) - xyz2 = discretize.utils.ndgrid([x, y, z]) + xyz2 = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) etest1 = E_from_Dipole1( xyz1, None, dhs.rho, dhs.current, dhs.location_a, dhs.location_b @@ -1565,8 +1566,8 @@ def test_electric_field(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz3 = discretize.utils.ndgrid([x, y, z]) - xyz4 = discretize.utils.ndgrid([x, y, z]) + xyz3 = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) + xyz4 = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) with pytest.raises(ValueError): dhs.electric_field(xyz3) @@ -1578,12 +1579,12 @@ def test_current_density(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 0., 50) - xyz1 = discretize.utils.ndgrid([x, y, z]) + xyz1 = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) x = np.linspace(-30., 20., 50) y = np.linspace(-20., 30., 50) z = np.linspace(-30., 0., 50) - xyz2 = discretize.utils.ndgrid([x, y, z]) + xyz2 = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) jtest1 = J_from_Dipole1( xyz1, None, dhs.rho, dhs.current, dhs.location_a, dhs.location_b @@ -1605,8 +1606,8 @@ def test_current_density(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz3 = discretize.utils.ndgrid([x, y, z]) - xyz4 = discretize.utils.ndgrid([x, y, z]) + xyz3 = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) + xyz4 = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) with pytest.raises(ValueError): dhs.current_density(xyz3) diff --git a/tests/test_em_tdem_halfspace.py b/tests/test_em_tdem_halfspace.py index 60953219..4ee3516b 100644 --- a/tests/test_em_tdem_halfspace.py +++ b/tests/test_em_tdem_halfspace.py @@ -1,7 +1,6 @@ import numpy as np from geoana.em.tdem import VerticalMagneticDipoleHalfSpace -import discretize from geoana.em.tdem.base import theta from scipy.special import erf, ive @@ -10,7 +9,7 @@ def H_from_Vertical( XYZ, loc, time, sigma, mu, moment ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) r_vec = XYZ - loc r = np.linalg.norm(r_vec[:, :2], axis=-1) @@ -38,8 +37,7 @@ def H_from_Vertical( def dH_from_Vertical( XYZ, loc, time, sigma, mu, moment ): - - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) r_vec = XYZ - loc r = np.linalg.norm(r_vec[:, :2], axis=-1) @@ -85,7 +83,7 @@ def test_magnetic_field(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) htest = H_from_Vertical( xyz, vmdhs.location, vmdhs.time, vmdhs.sigma, vmdhs.mu, vmdhs.moment @@ -103,7 +101,7 @@ def test_magnetic_flux(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) btest = B_from_Vertical( xyz, vmdhs.location, vmdhs.time, vmdhs.sigma, vmdhs.mu, vmdhs.moment @@ -121,7 +119,7 @@ def test_magnetic_field_dt(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) dh_test = dH_from_Vertical( xyz, vmdhs.location, vmdhs.time, vmdhs.sigma, vmdhs.mu, vmdhs.moment @@ -139,7 +137,7 @@ def test_magnetic_flux_dt(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) db_test = dB_from_Vertical( xyz, vmdhs.location, vmdhs.time, vmdhs.sigma, vmdhs.mu, vmdhs.moment diff --git a/tests/test_em_tdem_planewave.py b/tests/test_em_tdem_planewave.py index ab208ac4..d9807b2d 100644 --- a/tests/test_em_tdem_planewave.py +++ b/tests/test_em_tdem_planewave.py @@ -1,6 +1,5 @@ import pytest import numpy as np -import discretize from geoana.em import tdem from scipy.constants import mu_0, epsilon_0 @@ -40,7 +39,7 @@ def test_electric_field(): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) z = xyz[:, 2] bunja = -mu_0 ** 0.5 * z * np.exp(-(mu_0 * z ** 2) / 4) @@ -77,7 +76,7 @@ def test_current_density(): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) z = xyz[:, 2] bunja = -mu_0 ** 0.5 * z * np.exp(-(mu_0 * z ** 2) / 4) @@ -115,7 +114,7 @@ def test_magnetic_field(): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) z = xyz[:, 2] hy = -(np.sqrt(1 / (np.pi * mu_0)) * np.exp(-(mu_0 * z ** 2) / 4)) @@ -150,7 +149,7 @@ def test_magnetic_flux_density(): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) z = xyz[:, 2] by = -mu_0 * (np.sqrt(1 / (np.pi * mu_0)) * np.exp(-(mu_0 * z ** 2) / 4)) diff --git a/tests/test_gravity.py b/tests/test_gravity.py index 70be7178..72d6bf94 100644 --- a/tests/test_gravity.py +++ b/tests/test_gravity.py @@ -1,5 +1,4 @@ import pytest -import discretize import numpy as np from scipy.constants import G @@ -10,7 +9,7 @@ def U_from_PointMass( XYZ, loc, m ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) r_vec = XYZ - loc r = np.linalg.norm(r_vec, axis=-1) @@ -23,7 +22,7 @@ def g_from_PointMass( XYZ, loc, m ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) r_vec = XYZ - loc r = np.linalg.norm(r_vec, axis=-1) @@ -36,7 +35,7 @@ def gtens_from_PointMass( XYZ, loc, m ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) r_vec = XYZ - loc r = np.linalg.norm(r_vec, axis=-1) @@ -72,7 +71,7 @@ def test_gravitational_potential(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) utest = U_from_PointMass( xyz, pm.location, pm.mass @@ -93,7 +92,7 @@ def test_gravitational_field(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) gtest = g_from_PointMass( xyz, pm.location, pm.mass @@ -114,7 +113,7 @@ def test_gravitational_gradient(self): x = np.linspace(-20., 20., 5) y = np.linspace(-30., 30., 5) z = np.linspace(-40., 40., 5) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) g_tenstest = gtens_from_PointMass( xyz, pm.location, pm.mass @@ -132,7 +131,7 @@ def U_from_Sphere( XYZ, loc, m, rho, radius ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) r_vec = XYZ - loc r = np.linalg.norm(r_vec, axis=-1) @@ -148,7 +147,7 @@ def g_from_Sphere( XYZ, loc, m, rho, radius ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) r_vec = XYZ - loc r = np.linalg.norm(r_vec, axis=-1) @@ -164,7 +163,7 @@ def gtens_from_Sphere( XYZ, loc, m, rho, radius ): - XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3) + XYZ = np.atleast_2d(XYZ) r_vec = XYZ - loc r = np.linalg.norm(r_vec, axis=-1) @@ -175,7 +174,7 @@ def gtens_from_Sphere( g_tens[ind0] = -G * m * (np.eye(3) / r[ind0, None, None] ** 3 - 3 * r_vec[ind0, None] * r_vec[ind0, None, :] / r[ind0, None, None] ** 5) g_tens[~ind0] = -G * 4 / 3 * np.pi * rho * np.eye(3) - g_tens[ind1] = np.NaN + g_tens[ind1] = np.nan return g_tens @@ -215,7 +214,7 @@ def test_gravitational_potential(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) utest = U_from_Sphere( xyz, s.location, s.mass, s.rho, s.radius @@ -239,7 +238,7 @@ def test_gravitational_field(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) gtest = g_from_Sphere( xyz, s.location, s.mass, s.rho, s.radius @@ -263,7 +262,7 @@ def test_gravitational_gradient(self): x = np.linspace(-20., 20., 50) y = np.linspace(-30., 30., 50) z = np.linspace(-40., 40., 50) - xyz = discretize.utils.ndgrid([x, y, z]) + xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3) g_tens_test = gtens_from_Sphere( xyz, s.location, s.mass, s.rho, s.radius diff --git a/tests/test_potential_prism.py b/tests/test_potential_prism.py index 87773975..01097725 100644 --- a/tests/test_potential_prism.py +++ b/tests/test_potential_prism.py @@ -1,6 +1,5 @@ import numpy as np from numpy.testing import assert_allclose -from discretize.tests import check_derivative import pytest import geoana.kernels.potential_field_prism as pf @@ -11,6 +10,11 @@ except ImportError: njit = None +try: + from discretize.tests import check_derivative +except ImportError: + check_derivative = None + class TestCompiledVsNumpy(): xyz = np.mgrid[-50:50:51j, -50:50:51j, -50:50:51j] @@ -73,6 +77,7 @@ def test_fxyz(self): assert_allclose(v0, v1) +@pytest.mark.skipif(check_derivative is None, reason="discretize not installed.") class TestGravityPrismDerivatives(): xyz = np.mgrid[-100:100:26j, -100:100:26j, -100:100:26j] x = xyz[0].ravel() @@ -144,6 +149,7 @@ def J(v): assert check_derivative(grav_field, self.z, num=3, plotIt=False) +@pytest.mark.skipif(check_derivative is None, reason="discretize not installed.") class TestMagPrismDerivatives(): xyz = np.mgrid[-100:100:26j, -100:100:26j, -100:100:26j] x = xyz[0].ravel() @@ -288,34 +294,34 @@ def test_grav_init_and_errors(): prism.rho = 'abc' -if njit is not None: - @pytest.mark.parametrize('function', [ - pf.prism_f, - pf.prism_fz, - pf.prism_fzx, - pf.prism_fzy, - pf.prism_fzz, - pf.prism_fzzz, - pf.prism_fxxy, - pf.prism_fxxz, - pf.prism_fxyz - ]) - def test_numba_jitting_nopython(function): - - # create a vectorized jit function of it: - - x = np.random.rand(10) - y = np.random.rand(10) - z = np.random.rand(10) - @njit - def jitted_func(x, y, z): - n = len(x) - out = np.empty_like(x) - for i in range(n): - out[i] = function(x[i], y[i], z[i]) - return out - - v1 = jitted_func(x, y, z) - v2 = function(x, y, z) - - assert_allclose(v1, v2) +@pytest.mark.skipif(njit is None, reason="numba is not installed.") +@pytest.mark.parametrize('function', [ + pf.prism_f, + pf.prism_fz, + pf.prism_fzx, + pf.prism_fzy, + pf.prism_fzz, + pf.prism_fzzz, + pf.prism_fxxy, + pf.prism_fxxz, + pf.prism_fxyz +]) +def test_numba_jitting_nopython(function): + + # create a vectorized jit function of it: + + x = np.random.rand(10) + y = np.random.rand(10) + z = np.random.rand(10) + @njit + def jitted_func(x, y, z): + n = len(x) + out = np.empty_like(x) + for i in range(n): + out[i] = function(x[i], y[i], z[i]) + return out + + v1 = jitted_func(x, y, z) + v2 = function(x, y, z) + + assert_allclose(v1, v2) diff --git a/tests/test_utils.py b/tests/test_utils.py index bba5c81b..79b2a149 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -3,6 +3,11 @@ from geoana.utils import check_xyz_dim, mkvc, ndgrid import pytest +try: + from numpy.exceptions import ComplexWarning +except ImportError: + from numpy import ComplexWarning + def test_config_info(): info = geoana.show_config() @@ -95,7 +100,7 @@ def test_dtype_cast_good(): def test_bad_dtype_cast(): xyz = np.random.rand(10, 3) + 1j* np.random.rand(10, 3) - with pytest.warns(np.ComplexWarning): + with pytest.warns(ComplexWarning): xyz2 = check_xyz_dim(xyz) xyz = [['0', 'O', 'o']] with pytest.raises(ValueError):