Alternative implementation of qpms_scatsys_scattered_E()

For testing purposes. Seems to work OK.


Former-commit-id: 897e687d27dbd81b2aaac17fb8f19bc4257dc887
This commit is contained in:
Marek Nečada 2020-03-21 17:56:44 +02:00
parent aa82b3a01a
commit b3d15a1bb7
5 changed files with 169 additions and 18 deletions

View File

@ -710,25 +710,29 @@ cdef class ScatteringSystem:
return retdict return retdict
def scattered_E(self, cdouble wavenumber, scatcoeffvector_full, evalpos): def scattered_E(self, cdouble wavenumber, scatcoeffvector_full, evalpos, bint alt=False):
evalpos = np.array(evalpos, copy=False) evalpos = np.array(evalpos, dtype=float, copy=False)
cdef np.ndarray evalpos_a = evalpos
cdef np.ndarray[dtype=complex] scv = np.array(scatcoeffvector_full, copy=False)
cdef cdouble[::1] scv_view = scv
if evalpos.shape[-1] != 3: if evalpos.shape[-1] != 3:
raise ValueError("Last dimension of evalpos has to be 3") raise ValueError("Last dimension of evalpos has to be 3")
cdef np.ndarray[complex] results = np.empty(evalpos.shape, dtype=complex) cdef np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
cdef np.ndarray[dtype=complex, ndim=1] scv = np.array(scatcoeffvector_full, copy=False)
cdef cdouble[::1] scv_view = scv
cdef np.ndarray[complex, ndim=2] results = np.empty((evalpos_a.shape[0],3), dtype=complex)
cdef ccart3_t res cdef ccart3_t res
cdef cart3_t pos cdef cart3_t pos
for i in np.ndindex(evalpos.shape[:-1]): cdef size_t i
for i in range(evalpos_a.shape[0]):
pos.x = evalpos_a[i,0] pos.x = evalpos_a[i,0]
pos.y = evalpos_a[i,1] pos.y = evalpos_a[i,1]
pos.z = evalpos_a[i,2] pos.z = evalpos_a[i,2]
res = qpms_scatsys_scattered_E(self.s, wavenumber, &scv_view[0], pos) if alt:
res = qpms_scatsys_scattered_E__alt(self.s, wavenumber, &scv_view[0], pos)
else:
res = qpms_scatsys_scattered_E(self.s, wavenumber, &scv_view[0], pos)
results[i,0] = res.x results[i,0] = res.x
results[i,1] = res.y results[i,1] = res.y
results[i,2] = res.z results[i,2] = res.z
return results return results.reshape(evalpos.shape)
def empty_lattice_modes_xy(EpsMu epsmu, reciprocal_basis, wavevector, double maxomega): def empty_lattice_modes_xy(EpsMu epsmu, reciprocal_basis, wavevector, double maxomega):
'''Empty (2D, xy-plane) lattice mode (diffraction order) frequencies of a non-dispersive medium. '''Empty (2D, xy-plane) lattice mode (diffraction order) frequencies of a non-dispersive medium.
@ -863,25 +867,29 @@ cdef class _ScatteringSystemAtOmega:
def translation_matrix_full(self, blochvector = None): def translation_matrix_full(self, blochvector = None):
return self.ss_pyref.translation_matrix_full(wavenumber=self.wavenumber, blochvector=blochvector) return self.ss_pyref.translation_matrix_full(wavenumber=self.wavenumber, blochvector=blochvector)
def scattered_E(self, scatcoeffvector_full, evalpos): def scattered_E(self, scatcoeffvector_full, evalpos, bint alt=False):
evalpos = np.array(evalpos, copy=False) evalpos = np.array(evalpos, dtype=float, copy=False)
cdef np.ndarray evalpos_a = evalpos
cdef np.ndarray[dtype=complex] scv = np.array(scatcoeffvector_full, copy=False)
cdef cdouble[::1] scv_view = scv
if evalpos.shape[-1] != 3: if evalpos.shape[-1] != 3:
raise ValueError("Last dimension of evalpos has to be 3") raise ValueError("Last dimension of evalpos has to be 3")
cdef np.ndarray[complex] results = np.empty(evalpos.shape, dtype=complex) cdef np.ndarray[double,ndim=2] evalpos_a = evalpos.reshape(-1,3)
cdef np.ndarray[dtype=complex, ndim=1] scv = np.array(scatcoeffvector_full, copy=False)
cdef cdouble[::1] scv_view = scv
cdef np.ndarray[complex, ndim=2] results = np.empty((evalpos_a.shape[0],3), dtype=complex)
cdef ccart3_t res cdef ccart3_t res
cdef cart3_t pos cdef cart3_t pos
for i in np.ndindex(evalpos.shape[:-1]): cdef size_t i
for i in range(evalpos_a.shape[0]):
pos.x = evalpos_a[i,0] pos.x = evalpos_a[i,0]
pos.y = evalpos_a[i,1] pos.y = evalpos_a[i,1]
pos.z = evalpos_a[i,2] pos.z = evalpos_a[i,2]
res = qpms_scatsysw_scattered_E(self.ssw, &scv_view[0], pos) if alt:
res = qpms_scatsysw_scattered_E__alt(self.ssw, &scv_view[0], pos)
else:
res = qpms_scatsysw_scattered_E(self.ssw, &scv_view[0], pos)
results[i,0] = res.x results[i,0] = res.x
results[i,1] = res.y results[i,1] = res.y
results[i,2] = res.z results[i,2] = res.z
return results return results.reshape(evalpos.shape)
cdef class ScatteringMatrix: cdef class ScatteringMatrix:

View File

@ -629,6 +629,10 @@ cdef extern from "scatsystem.h":
const cdouble *f_excitation_vector_full, cart3_t where) const cdouble *f_excitation_vector_full, cart3_t where)
ccart3_t qpms_scatsysw_scattered_E(const qpms_scatsys_at_omega_t *ssw, ccart3_t qpms_scatsysw_scattered_E(const qpms_scatsys_at_omega_t *ssw,
const cdouble *f_excitation_vector_full, cart3_t where) const cdouble *f_excitation_vector_full, cart3_t where)
ccart3_t qpms_scatsys_scattered_E__alt(const qpms_scatsys_t *ss, cdouble wavenumber,
const cdouble *f_excitation_vector_full, cart3_t where)
ccart3_t qpms_scatsysw_scattered_E__alt(const qpms_scatsys_at_omega_t *ssw,
const cdouble *f_excitation_vector_full, cart3_t where)
cdef extern from "ewald.h": cdef extern from "ewald.h":
struct qpms_csf_result: struct qpms_csf_result:

View File

@ -2018,6 +2018,69 @@ ccart3_t qpms_scatsysw_scattered_E(const qpms_scatsys_at_omega_t *ssw,
cvf, where); cvf, where);
} }
// Alternative implementation, using translation operator and regular dipole waves at zero
ccart3_t qpms_scatsys_scattered_E__alt(const qpms_scatsys_t *ss,
const complex double k,
const complex double *cvf,
const cart3_t where
) {
QPMS_UNTESTED;
ccart3_t res = {0,0,0};
ccart3_t res_kc = {0,0,0}; // kahan sum compensation
// bspec containing only electric dipoles
const qpms_vswf_set_spec_t dipspec = {
.n = 3,
.ilist = (qpms_uvswfi_t[]){
qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, -1, 1),
qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, 0, 1),
qpms_tmn2uvswfi(QPMS_VSWF_ELECTRIC, +1, 1),
},
.lMax=1, .lMax_M=0, .lMax_N=1, .lMax_L=-1,
.capacity=0,
.norm = ss->c->normalisation,
};
ccart3_t regdipoles_0[3]; {
const sph_t origin_sph = {.r = 0, .theta = M_PI_2, .phi=0}; // Should work with any theta/phi (TESTWORTHY)
csphvec_t regdipoles_0_sph[3];
QPMS_ENSURE_SUCCESS(qpms_uvswf_fill(regdipoles_0_sph, &dipspec,
sph2csph(origin_sph), QPMS_BESSEL_REGULAR));
for(int i = 0; i < 3; ++i)
regdipoles_0[i] = csphvec2ccart(regdipoles_0_sph[i], origin_sph);
}
complex double *s; // Translation matrix
QPMS_CRASHING_MALLOC(s, ss->max_bspecn * sizeof(*s) * dipspec.n);
for (qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) {
const qpms_vswf_set_spec_t *bspec = qpms_ss_bspec_pi(ss, pi);
const cart3_t particle_pos = ss->p[pi].pos;
const complex double *particle_cv = cvf + ss->fecv_pstarts[pi];
const cart3_t origin_cart = {0, 0, 0};
QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_lc3p(
ss->c, s, &dipspec, 1, bspec, 3, k, particle_pos, where,
QPMS_HANKEL_PLUS));
for(size_t i = 0; i < bspec->n; ++i)
for(size_t j = 0; j < 3; ++j){
ccart3_t summand = ccart3_scale(particle_cv[i] * s[3*i+j], regdipoles_0[j]);
ckahanadd(&(res.x), &(res_kc.x), summand.x);
ckahanadd(&(res.y), &(res_kc.y), summand.y);
ckahanadd(&(res.z), &(res_kc.z), summand.z);
}
}
free(s);
return res;
}
ccart3_t qpms_scatsysw_scattered_E__alt(const qpms_scatsys_at_omega_t *ssw,
const complex double *cvf, const cart3_t where) {
return qpms_scatsys_scattered_E__alt(ssw->ss, ssw->wavenumber,
cvf, where);
}
#if 0 #if 0
ccart3_t qpms_scatsys_scattered_E_irrep(const qpms_scatsys_t *ss, ccart3_t qpms_scatsys_scattered_E_irrep(const qpms_scatsys_t *ss,

View File

@ -727,6 +727,37 @@ ccart3_t qpms_scatsysw_scattered_E(
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated. cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
); );
/// Evaluates scattered electric field at a point, given a full vector of scattered field coefficients.
/**
* This is an alternative implementation of qpms_scatsys_scattered_E(), and should give the same results
* up to rounding errors.
*
* \return Complex electric field at the point defined by \a evalpoint.
*
* \see qpms_scatsys_scattered_E()
*/
ccart3_t qpms_scatsys_scattered_E__alt(
const qpms_scatsys_t *ss,
complex double wavenumber, ///< Wavenumber of the background medium.
const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
/// Evaluates scattered electric field at a point, given a full vector of scattered field coefficients.
/**
* This is an alternative implementation of qpms_scatsys_scattered_E(), and should give the same results
* up to rounding errors.
*
* \return Complex electric field at the point defined by \a evalpoint.
*
* \see qpms_scatsysw_scattered_E()
*/
ccart3_t qpms_scatsysw_scattered_E__alt(
const qpms_scatsys_at_omega_t *ssw,
const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
);
#if 0 #if 0
/** Evaluates partial scattered fields (corresponding to a given irrep-reduced excitation vector) /** Evaluates partial scattered fields (corresponding to a given irrep-reduced excitation vector)
* at a given point. * at a given point.

View File

@ -0,0 +1,45 @@
#!/usr/bin/env python3
"""
Tests whether direct evaluation of VSWFs gives the same result as their
representation in terms of translation operators and regular electric dipole
waves at origin
"""
from qpms import Particle, CTMatrix, lorentz_drude, EpsMuGenerator, TMatrixGenerator, BaseSpec, FinitePointGroup, ScatteringSystem, TMatrixInterpolator, EpsMu, dbgmsg_enable, dbgmsg_disable, dbgmsg_active, BesselType,eV, hbar
from qpms.symmetries import point_group_info
import numpy as np
eh = eV/hbar
np.random.seed(666)
dbgmsg_enable(2)
part_radius = 80e-9
p = 1580e-9
sym = FinitePointGroup(point_group_info['D4h'])
bspec1 = BaseSpec(lMax=3)
medium=EpsMuGenerator(EpsMu(1.52**2))
t1 = TMatrixGenerator.sphere(medium, EpsMuGenerator(lorentz_drude['Au']), r=part_radius)
p1 = Particle((0,0,0),t1,bspec=bspec1)
ss, ssw = ScatteringSystem.create([p1], EpsMuGenerator(EpsMu(1.52**2)), 1.4*eh, sym)
points = np.random.random((100,3)) * p
points = points[np.linalg.norm(points, axis=-1) > part_radius]
t,l,m = bspec1.tlm()
fails=0
for i in range(ss.fecv_size):
fvc = np.zeros((ss.fecv_size,), dtype=complex)
fvc[i] = 1
E = ssw.scattered_E(fvc, points)
E_alt = ssw.scattered_E(fvc, points,alt=True)
diff = abs(E-E_alt)
reldiffavg = np.average(diff/(abs(E)+abs(E_alt)))
fail = reldiffavg > 1e-3
fails += fail
print('E' if t[i] == 2 else 'M', l[i], m[i], np.amax(diff), reldiffavg, 'FAIL!' if fail else 'OK')
exit(fails)