From b3d15a1bb71eb10a7deffce7be2857b51bcc5877 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Sat, 21 Mar 2020 17:56:44 +0200 Subject: [PATCH] Alternative implementation of qpms_scatsys_scattered_E() For testing purposes. Seems to work OK. Former-commit-id: 897e687d27dbd81b2aaac17fb8f19bc4257dc887 --- qpms/qpms_c.pyx | 44 ++++++++++-------- qpms/qpms_cdefs.pxd | 4 ++ qpms/scatsystem.c | 63 ++++++++++++++++++++++++++ qpms/scatsystem.h | 31 +++++++++++++ tests/test_vswf_transop_consistency.py | 45 ++++++++++++++++++ 5 files changed, 169 insertions(+), 18 deletions(-) create mode 100755 tests/test_vswf_transop_consistency.py diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 26a0ad4..4e6da2c 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -710,25 +710,29 @@ cdef class ScatteringSystem: return retdict - def scattered_E(self, cdouble wavenumber, scatcoeffvector_full, evalpos): - evalpos = np.array(evalpos, 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 + def scattered_E(self, cdouble wavenumber, scatcoeffvector_full, evalpos, bint alt=False): + evalpos = np.array(evalpos, dtype=float, copy=False) if evalpos.shape[-1] != 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 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.y = evalpos_a[i,1] 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,1] = res.y results[i,2] = res.z - return results + return results.reshape(evalpos.shape) 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. @@ -863,25 +867,29 @@ cdef class _ScatteringSystemAtOmega: def translation_matrix_full(self, blochvector = None): return self.ss_pyref.translation_matrix_full(wavenumber=self.wavenumber, blochvector=blochvector) - def scattered_E(self, scatcoeffvector_full, evalpos): - evalpos = np.array(evalpos, 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 + def scattered_E(self, scatcoeffvector_full, evalpos, bint alt=False): + evalpos = np.array(evalpos, dtype=float, copy=False) if evalpos.shape[-1] != 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 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.y = evalpos_a[i,1] 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,1] = res.y results[i,2] = res.z - return results + return results.reshape(evalpos.shape) cdef class ScatteringMatrix: diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 61ea266..c6d806a 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -629,6 +629,10 @@ cdef extern from "scatsystem.h": const cdouble *f_excitation_vector_full, cart3_t where) ccart3_t qpms_scatsysw_scattered_E(const qpms_scatsys_at_omega_t *ssw, 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": struct qpms_csf_result: diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 12cf040..8f0844d 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -2018,6 +2018,69 @@ ccart3_t qpms_scatsysw_scattered_E(const qpms_scatsys_at_omega_t *ssw, 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 ccart3_t qpms_scatsys_scattered_E_irrep(const qpms_scatsys_t *ss, diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 1dbf952..d7260a0 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -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. ); +/// 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 /** Evaluates partial scattered fields (corresponding to a given irrep-reduced excitation vector) * at a given point. diff --git a/tests/test_vswf_transop_consistency.py b/tests/test_vswf_transop_consistency.py new file mode 100755 index 0000000..d6ea7ba --- /dev/null +++ b/tests/test_vswf_transop_consistency.py @@ -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)