Skip to content

Allows complex mu values in rTE kernel #86

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jcapriot merged 3 commits into from
Dec 2, 2025
Merged

Allows complex mu values in rTE kernel #86

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
f32ba92
Allows complex mu values in rTE kernel
jcapriot Dec 2, 2025
c363f17
Explicitly cast mu to complex_t pointer
jcapriot Dec 2, 2025
1bea0df
test kernels with complex conductivity and complex permeability
jcapriot Dec 2, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 8 additions & 7 deletions geoana/kernels/_extensions/_rTE.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,19 +9,20 @@ void funcs::rTE(
double * frequencies,
double * lambdas,
complex_t * sigmas,
double * mus,
complex_t * mus,
double * h,
size_t n_frequency,
size_t n_filter,
size_t n_layers
){
) noexcept(true){
double mu_0 = 4.0E-7*M_PI;
complex_t j(0.0, 1.0);

complex_t k2, u, Yh, Y, tanhuh, Y0;
complex_t *sigi;
double *mui;
double omega, mu_c, l2;
complex_t *mui;
double omega, l2;
complex_t mu_c;

for(size_t i_filt=0, i=0; i_filt<n_filter; ++i_filt){
l2 = lambdas[i_filt] * lambdas[i_filt];
Expand Down Expand Up @@ -55,12 +56,12 @@ void funcs::rTEgrad(
double * frequencies,
double * lambdas,
complex_t * sigmas,
double * mus,
complex_t * mus,
double * h,
size_t n_frequency,
size_t n_filter,
size_t n_layers
){
) noexcept(true){
double mu_0 = 4.0E-7*M_PI;
complex_t j(0.0, 1.0);

Expand All @@ -74,7 +75,7 @@ void funcs::rTEgrad(
vec_complex_t k2s(n_layers), us(n_layers), Yhs(n_layers);
vec_complex_t Ys(n_layers-1), tanhs(n_layers-1);

double *mui;
complex_t *mui;
complex_t *sigi, *TE_dhi, *TE_dmui, *TE_dsigmai;
complex_t gyh0, bot, gy, gtanh, Y0;
complex_t gu, gmu, gk2;
Expand Down
8 changes: 4 additions & 4 deletions geoana/kernels/_extensions/_rTE.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,12 +14,12 @@ namespace funcs {
double *frequencies,
double *lambdas,
complex_t *sigmas,
double *mus,
complex_t *mus,
double *depths,
size_t n_frequency,
size_t n_filter,
size_t n_layers
);
) noexcept(true);

void rTEgrad(
complex_t * TE_dsigma,
Expand All @@ -28,12 +28,12 @@ namespace funcs {
double * frequencies,
double * lambdas,
complex_t * sigmas,
double * mus,
complex_t * mus,
double * h,
size_t n_frequency,
size_t n_filter,
size_t n_layers
);
) noexcept(true);
}

#endif
22 changes: 12 additions & 10 deletions geoana/kernels/_extensions/rTE.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,12 +19,12 @@ cdef extern from "_rTE.h" namespace "funcs":
double *frequencies,
double *lambdas,
complex_t *sigmas,
double *mus,
complex_t *mus,
double *thicks,
SIZE_t n_frequency,
SIZE_t n_filter,
SIZE_t n_layers
) nogil
) noexcept nogil

void rTEgrad(
complex_t * TE_dsigma,
Expand All @@ -33,12 +33,12 @@ cdef extern from "_rTE.h" namespace "funcs":
double * frequencies,
double * lambdas,
complex_t * sigmas,
double * mus,
complex_t * mus,
double * h,
SIZE_t n_frequency,
SIZE_t n_filter,
SIZE_t n_layers
) nogil
) noexcept nogil

def rTE_forward(frequencies, lamb, sigma, mu, thicknesses):
"""Compute reflection coefficients for Transverse Electric (TE) mode.
Expand Down Expand Up @@ -67,7 +67,7 @@ def rTE_forward(frequencies, lamb, sigma, mu, thicknesses):
"""
# Sigma and mu must be fortran contiguous
sigma = np.require(sigma, dtype=np.complex128, requirements="F")
mu = np.require(mu, dtype=np.float64, requirements="F")
mu = np.require(mu, dtype=np.complex128, requirements="F")

# These just must be contiguous
lamb = np.require(lamb, dtype=np.float64, requirements="C")
Expand All @@ -92,7 +92,7 @@ def rTE_forward(frequencies, lamb, sigma, mu, thicknesses):
REAL_t[:] f = frequencies
REAL_t[:] lam = lamb
COMPLEX_t[:, :] sig = sigma
REAL_t[:, :] c_mu = mu
COMPLEX_t[:, :] c_mu = mu
REAL_t[:] hs = thicknesses

cdef:
Expand All @@ -109,12 +109,13 @@ def rTE_forward(frequencies, lamb, sigma, mu, thicknesses):
# Types are same size and both pairs of numbers (r, i), so can cast one to the other
complex_t *out_p = <complex_t *> &out[0, 0]
complex_t *sig_p = <complex_t *> &sig[0, 0]
complex_t *mu_p = <complex_t *> &c_mu[0, 0]
REAL_t *h_p = NULL
if n_layers > 1:
h_p = &hs[0]

with nogil:
rTE(out_p, &f[0], &lam[0], sig_p, &c_mu[0, 0], h_p,
rTE(out_p, &f[0], &lam[0], sig_p, mu_p, h_p,
n_frequency, n_filter, n_layers)

return np.array(out)
Expand Down Expand Up @@ -153,7 +154,7 @@ def rTE_gradient(frequencies, lamb, sigma, mu, thicknesses):
"""
# Require that they are all the same ordering as sigma
sigma = np.require(sigma, dtype=np.complex128, requirements="F")
mu = np.require(mu, dtype=np.float64, requirements="F")
mu = np.require(mu, dtype=np.complex128, requirements="F")
lamb = np.require(lamb, dtype=np.float64, requirements="C")
frequencies = np.require(frequencies, dtype=np.float64, requirements="C")
thicknesses = np.require(thicknesses, dtype=np.float64, requirements="C")
Expand All @@ -176,7 +177,7 @@ def rTE_gradient(frequencies, lamb, sigma, mu, thicknesses):
REAL_t[:] f = frequencies
REAL_t[:] lam = lamb
COMPLEX_t[:, :] sig = sigma
REAL_t[:, :] c_mu = mu
COMPLEX_t[:, :] c_mu = mu
REAL_t[:] hs = thicknesses

cdef:
Expand All @@ -194,6 +195,7 @@ def rTE_gradient(frequencies, lamb, sigma, mu, thicknesses):
cdef:
# Types are same size and both pairs of numbers (r, i), so can cast one to the other
complex_t *sig_p = <complex_t *> &sig[0, 0]
complex_t *mu_p = <complex_t *> &c_mu[0, 0]
complex_t *gsig_p = <complex_t *> &gsig[0, 0, 0]
complex_t *gmu_p = <complex_t *> &gmu[0, 0, 0]
complex_t *gh_p = NULL
Expand All @@ -203,7 +205,7 @@ def rTE_gradient(frequencies, lamb, sigma, mu, thicknesses):
gh_p = <complex_t *> &gh[0, 0, 0]

with nogil:
rTEgrad(gsig_p, gmu_p, gh_p, &f[0], &lam[0], sig_p, &c_mu[0, 0], h_p,
rTEgrad(gsig_p, gmu_p, gh_p, &f[0], &lam[0], sig_p, mu_p, h_p,
n_frequency, n_filter, n_layers)

return np.array(gsig), np.array(gh), np.array(gmu)
10 changes: 5 additions & 5 deletions tests/em/fdem/test_em_fdem_layered.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -294,9 +294,9 @@ def test_rTE_jacobian(self):
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)
sigma = 1E-1*(1 + 1E-4j + 0.1 * np.random.rand(n_layer, n_frequency))
mu = mu_0 * (1 + 1.0E-4j + 0.1 * np.random.rand(n_layer, n_frequency))
dmu = mu_0 * (0.1 + 1.0E-5j) * np.random.rand(n_layer, n_frequency)

def rte_sigma(x):
sigma = x.reshape(n_layer, n_frequency)
Expand Down Expand Up @@ -354,8 +354,8 @@ def setUp(self):
thicknesses = np.ones(n_layer-1)
lamb = np.logspace(-5, -3, n_lambda)
np.random.seed(123)
sigma = 1E-1 * (1 + 1.0/(n_layer*n_frequency) * np.arange(n_layer*n_frequency).reshape(n_layer, n_frequency))
mu = mu_0 * (1 + 1.0/(n_layer*n_frequency) * np.arange(n_layer*n_frequency).reshape(n_layer, n_frequency))
sigma = 1E-1 * (1 + 1E-4j + 1.0/(n_layer*n_frequency) * np.arange(n_layer*n_frequency).reshape(n_layer, n_frequency))
mu = mu_0 * (1 + 1.0E-4j + 1.0/(n_layer*n_frequency) * np.arange(n_layer*n_frequency).reshape(n_layer, n_frequency))

self.rTE1 = rTE_forward(frequencies, lamb, sigma, mu, thicknesses)
self.rTE2 = _rTE_forward(frequencies, lamb, sigma, mu, thicknesses)
Expand Down
Loading