From 1158e116d2dbf62a6266f393e807b64ee8531b0b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 24 Jul 2020 09:43:19 +0300 Subject: [PATCH] "Basis fields" for finite systems. There seem to be race conditions in the prange'd cython parts. --- qpms/qpms_c.py | 0 qpms/qpms_c.pyx | 39 +++++++++++++++++++-------- qpms/qpms_cdefs.pxd | 4 +++ qpms/scatsystem.c | 48 +++++++++++++++++++++++---------- qpms/scatsystem.h | 65 +++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 131 insertions(+), 25 deletions(-) create mode 100644 qpms/qpms_c.py diff --git a/qpms/qpms_c.py b/qpms/qpms_c.py new file mode 100644 index 0000000..e69de29 diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 10ed7e0..4733667 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -1286,8 +1286,9 @@ cdef class _ScatteringSystemAtOmega: return results.reshape(evalpos.shape) @boundscheck(False) - def scattered_field_basis(self, evalpos, blochvector=None, btyp=QPMS_HANKEL_PLUS): + def scattered_field_basis(self, evalpos, blochvector=None, particle=None, btyp=QPMS_HANKEL_PLUS): # TODO examples + # FIXME periodic case not implemented """Evaluate scattered field "basis" This function enables the evaluation of "scattered" fields @@ -1301,25 +1302,31 @@ cdef class _ScatteringSystemAtOmega: Evaluation points in cartesian coordinates. blochvector: array_like or None Bloch vector, must be supplied (non-None) for periodic systems, else None. + particle_index: int or None (default), optional + A valid particle index; if specified, only the contribution of the given particle + is evaluated. Not applicable for periodic arrays. btyp: BesselType, optional Kind of the waves. Defaults to BesselType.HANKEL_PLUS. Returns ------- - ndarray of complex, with the shape `evalpos.shape[:-1] + (self.fecv_size, 3)` + ndarray of complex, with the shape `evalpos.shape[:-1] + (n, 3)` "Basis" fields at the positions given in `evalpos`, in cartesian coordinates. + `n` here stays either for `self.fecv_size` (if `particle==None`) + or `len(self.bspec_pi(particle))`. """ - if(btyp != QPMS_HANKEL_PLUS): - raise NotImplementedError("Only first kind Bessel function-based fields are supported") + #if(btyp != QPMS_HANKEL_PLUS): # #TODO, IIRC not supported only for periodic systems + # raise NotImplementedError("Only first kind Bessel function-based fields are supported") cdef qpms_bessel_t btyp_c = BesselType(btyp) - cdef Py_ssize_t fecv_size = self.fecv_size + cdef Py_ssize_t basissize = self.fecv_size if particle is None else len(self.bspec_pi(particle)) + cdef qpms_ss_pi_t pi = particle 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[double,ndim=2] evalpos_a = evalpos.reshape(-1,3) - cdef np.ndarray[complex, ndim=3] results = np.empty((evalpos_a.shape[0], fecv_size, 3), dtype=complex) + cdef np.ndarray[complex, ndim=3] results = np.empty((evalpos_a.shape[0], basissize, 3), dtype=complex) cdef ccart3_t *res - res = malloc(fecv_size*sizeof(ccart3_t)) + res = malloc(basissize*sizeof(ccart3_t)) cdef cart3_t pos cdef Py_ssize_t i, j with nogil, wraparound(False), parallel(): @@ -1327,14 +1334,24 @@ cdef class _ScatteringSystemAtOmega: pos.x = evalpos_a[i,0] pos.y = evalpos_a[i,1] pos.z = evalpos_a[i,2] - qpms_scatsysw_scattered_field_basis(res, self.ssw, btyp_c, pos) - for j in range(fecv_size): + if particle is None: + qpms_scatsysw_scattered_field_basis(res, self.ssw, btyp_c, pos) + else: + qpms_scatsysw_scattered_field_basis_pi(res, self.ssw, pi, btyp_c, pos) + for j in range(basissize): results[i,j,0] = res[j].x results[i,j,1] = res[j].y results[i,j,2] = res[j].z free(res) - return results.reshape(evalpos.shape[:-1] + (self.fecv_size, 3)) - + return results.reshape(evalpos.shape[:-1] + (basissize, 3)) + + def bspec_pi(self, pi): + return self.ss_pyref.bspec_pi(pi) + + property bspecs: + def __get__(self): + return self.ss_pyref.bspecs + cdef class ScatteringMatrix: ''' diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 945def5..3b4d67d 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -709,6 +709,10 @@ cdef extern from "scatsystem.h": qpms_bessel_t btyp, cdouble wavenumber, cart3_t where) nogil qpms_errno_t qpms_scatsysw_scattered_field_basis(ccart3_t *target, const qpms_scatsys_at_omega_t *ssw, qpms_bessel_t btyp, cart3_t where) nogil + qpms_errno_t qpms_scatsys_scattered_field_basis_pi(ccart3_t *target, const qpms_scatsys_t *ss, + qpms_ss_pi_t pi, qpms_bessel_t btyp, cdouble wavenumber, cart3_t where) nogil + qpms_errno_t qpms_scatsysw_scattered_field_basis_pi(ccart3_t *target, const qpms_scatsys_at_omega_t *ssw, + qpms_ss_pi_t pi, qpms_bessel_t btyp, cart3_t where) nogil qpms_errno_t qpms_scatsyswk_scattered_field_basis(ccart3_t *target, const qpms_scatsys_at_omega_k_t *sswk, qpms_bessel_t btyp, cart3_t where) nogil double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, cdouble wavenumber, const double *wavevector) nogil diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index fc5c14a..c3e8ad0 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -2043,6 +2043,37 @@ ccart3_t qpms_scatsysw_scattered_E(const qpms_scatsys_at_omega_t *ssw, cvf, where); } +qpms_errno_t qpms_scatsys_scattered_field_basis_pi( + ccart3_t *target, + const qpms_scatsys_t *ss, + const qpms_ss_pi_t pi, + const qpms_bessel_t btyp, + const complex double k, + const cart3_t where + ) { + qpms_ss_ensure_nonperiodic_a(ss, "qpms_scatsyswk_scattered_field_basis()"); + QPMS_UNTESTED; + const qpms_vswf_set_spec_t *bspec = qpms_ss_bspec_pi(ss, pi); + csphvec_t *vswfs_sph; //Single particle contributions in spherical coordinates + QPMS_CRASHING_MALLOC(vswfs_sph, bspec->n * sizeof(*vswfs_sph)); + const cart3_t particle_pos = ss->p[pi].pos; + const csph_t kr = sph_cscale(k, cart2sph( + cart3_substract(where, particle_pos))); + QPMS_ENSURE_SUCCESS(qpms_uvswf_fill(vswfs_sph, bspec, kr, btyp)); + for(size_t i = 0; i < bspec->n; ++i) + target[i] = csphvec2ccart_csph(vswfs_sph[i], kr); + + free(vswfs_sph); + return QPMS_SUCCESS; +} + +qpms_errno_t qpms_scatsysw_scattered_field_basis_pi( + ccart3_t *target, const qpms_scatsys_at_omega_t *ssw, const qpms_ss_pi_t pi, + const qpms_bessel_t btyp, const cart3_t where) { + return qpms_scatsys_scattered_field_basis_pi(target, ssw->ss, pi, btyp, + ssw->wavenumber, where); +} + qpms_errno_t qpms_scatsys_scattered_field_basis( ccart3_t *target, const qpms_scatsys_t *ss, @@ -2052,20 +2083,9 @@ qpms_errno_t qpms_scatsys_scattered_field_basis( ) { qpms_ss_ensure_nonperiodic_a(ss, "qpms_scatsyswk_scattered_field_basis()"); QPMS_UNTESTED; - csphvec_t *vswfs_sph; //Single particle contributions in spherical coordinates - QPMS_CRASHING_CALLOC(vswfs_sph, ss->max_bspecn, sizeof(*vswfs_sph)); - - 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 csph_t kr = sph_cscale(k, cart2sph( - cart3_substract(where, particle_pos))); - QPMS_ENSURE_SUCCESS(qpms_uvswf_fill(vswfs_sph, bspec, kr, btyp)); - for(size_t i = 0; i < bspec->n; ++i) - target[ss->fecv_pstarts[pi] + i] = csphvec2ccart_csph(vswfs_sph[i], kr); - } - - free(vswfs_sph); + for (qpms_ss_pi_t pi = 0; pi < ss->p_count; ++pi) + QPMS_ENSURE_SUCCESS(qpms_scatsys_scattered_field_basis_pi( + target + ss->fecv_pstarts[pi], ss, pi, btyp, k, where)); return QPMS_SUCCESS; } diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 31596fa..3704ac9 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -757,12 +757,77 @@ qpms_errno_t qpms_scatsys_scattered_field_basis( cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated. ); +/// Evaluates a "basis" for electric field of a given particle at a given point. +/** + * This function evaluates all the included VSWFs from a particle in the system, evaluated + * at a given point. + * + * Finite systems only. + * + * \see qpms_scatsys_scattered_field_basis() for the whole-array equivalent. + * + * \see qpms_scatsysw_scattered_field_basis_pi() + */ +qpms_errno_t qpms_scatsys_scattered_field_basis_pi( + ccart3_t *target, ///< Target array of length at least `qpms_ss_bspec_pi(ssw->ss, pi)->n` + const qpms_scatsys_t *ss, + qpms_ss_pi_t pi, ///< Particle index + qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS). + complex double wavenumber, ///< Wavenumber of the background medium + cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated. + ); + /// Evaluates a "basis" for electric field at a given point. /** * This function evaluates all the included VSWFs from the particles in the system, evaluated * at a given point. Taking a linear combination of these with the coefficients \a scattcoeff_full[] * would be equivalent to the result of qpms_scatsysw_scattered_E(). * + * Note that this might require a relatively large amount of memory. Depending + * on the application, qpms_scatsysw_scattered_field_basis_pi() might be a better + * alternative. + * + * \see qpms_scatsysw_scattered_field_basis_pi() for the single-particle equivalent. + * + * \see qpms_scatsys_scattered_field_basis() + * + * \see qpms_scatsyswk_scattered_field_basis() for periodic systems. + * + */ +qpms_errno_t qpms_scatsysw_scattered_field_basis( + ccart3_t *target, ///< Target array of length \a ss->fecv_size + const qpms_scatsys_at_omega_t *ssw, + qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS). + cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated. + ); + +/// Evaluates a "basis" for electric field of a given particle at a given point. +/** + * This function evaluates all the included VSWFs from a particle in the system, evaluated + * at a given point. + * + * Finite systems only. + * + * \see qpms_scatsysw_scattered_field_basis() for the whole-array equivalent. + * + * \see qpms_scatsys_scattered_field_basis_pi() + */ +qpms_errno_t qpms_scatsysw_scattered_field_basis_pi( + ccart3_t *target, ///< Target array of length at least `qpms_ss_bspec_pi(ssw->ss, pi)->n` + const qpms_scatsys_at_omega_t *ssw, + qpms_ss_pi_t pi, ///< Particle index + qpms_bessel_t typ, ///< Bessel function kind to use (for scattered fields, use QPMS_HANKEL_PLUS). + cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated. + ); + +/// Evaluates a "basis" for electric field at a given point. +/** + * This function evaluates all the included VSWFs from the particles in the system, evaluated + * at a given point. Taking a linear combination of these with the coefficients \a scattcoeff_full[] + * would be equivalent to the result of qpms_scatsysw_scattered_E(). + * + * \see qpms_scatsysw_scattered_field_basis_pi() + * * \see qpms_scatsys_scattered_field_basis() * * \see qpms_scatsyswk_scattered_field_basis() for periodic systems.