Skip to content

Commit

Permalink
fix test for small rho
Browse files Browse the repository at this point in the history
Do a check for rho < sqrt(a**2 + z**2)
  • Loading branch information
jcapriot committed Oct 24, 2024
1 parent 3b1f8d5 commit 4cb4a17
Show file tree
Hide file tree
Showing 3 changed files with 10 additions and 16 deletions.
1 change: 1 addition & 0 deletions .github/workflows/release.yml
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,7 @@ jobs:
with:
path: ./dist/geoana*.whl


distribute:
name: Distribute documentation
runs-on: ubuntu-latest
Expand Down
12 changes: 4 additions & 8 deletions geoana/em/static/wholespace.py
Original file line number Diff line number Diff line change
Expand Up @@ -597,10 +597,8 @@ def vector_potential(self, xyz, coordinates="cartesian"):

A_cyl = np.zeros_like(r_vec)

# when rho is small relative to the radius and z, k_sq -> 0
# this is unstable so use a small argument approximation.
# check k_sq for the relative sizes
small_rho = k_sq < 1E-3
# when rho is small relative to the radius and z
small_rho = rho**2/(a**2 + z**2) < 1E-6

temp = np.sqrt(a**2 + z[small_rho]**2)
A_cyl[small_rho, 1] = np.pi * C * (
Expand Down Expand Up @@ -731,10 +729,8 @@ def magnetic_flux_density(self, xyz, coordinates="cartesian"):

B_cyl = np.zeros_like(r_cyl)

# when rho is small relative to the radius and z, k_sq -> 0
# this is unstable so use a small argument approximation.
# check k_sq for the relative sizes
small_rho = k_sq < 1E-3
# when rho is small relative to the radius and z
small_rho = rho**2/(a**2 + z**2) < 1E-6

temp = np.sqrt(a**2 + z[small_rho]**2)
B_cyl[small_rho, 0] = 3 * C * np.pi * a**2 * z[small_rho] * (
Expand Down
13 changes: 5 additions & 8 deletions tests/em/static/test_current_loop.py
Original file line number Diff line number Diff line change
Expand Up @@ -24,10 +24,10 @@ def circle_loop():

@pytest.fixture()
def xyz():
nx, ny, nz = (11, 13, 9)
x = np.linspace(-100, 100, nx)
y = np.linspace(-90, 90, ny)
z = np.linspace(-80, 80, nz)
nx, ny, nz = (32, 32, 32)
x = np.linspace(-50, 50, nx)
y = np.linspace(-50, 50, ny)
z = np.linspace(-50, 50, nz)
xyz = np.meshgrid(x, y, z)
return xyz

Expand Down Expand Up @@ -87,8 +87,5 @@ def test_circular_loop(circle_loop, orient, method, xyz):
test = getattr(circle_loop, method)(xyz)
other = getattr(loop_approx, method)(xyz)

atol = 1E-16
if method == 'magnetic_field':
atol /= mu_0

npt.assert_allclose(test, other, rtol=1E-3, atol=atol)
npt.assert_allclose(test, other, rtol=1E-3, atol=1E-20)

0 comments on commit 4cb4a17

Please sign in to comment.