diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 8b4e7fa..67c64b6 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -43,6 +43,22 @@ void qpms_scatsystem_set_nthreads(long n) { qpms_scatsystem_nthreads_override = n; } +static inline void qpms_ss_ensure_periodic(qpms_scatsys_t *ss) { + QPMS_ENSURE(ss->lattice_dimension > 0, "This method is applicable only to periodic systems."); +} + +static inline void qpms_ss_ensure_periodic_a(qpms_scatsys_t *ss, const char *s) { + QPMS_ENSURE(ss->lattice_dimension > 0, "This method is applicable only to periodic systems. Use %s instead.", s); +} + +static inline void qpms_ss_ensure_nonperiodic(qpms_scatsys_t *ss) { + QPMS_ENSURE(ss->lattice_dimension == 0, "This method is applicable only to nonperiodic systems."); +} + +static inline void qpms_ss_ensure_nonperiodic_a(qpms_scatsys_t *ss, const char *s) { + QPMS_ENSURE(ss->lattice_dimension == 0, "This method is applicable only to nonperiodic systems. Use %s instead.", s); +} + // ------------ Stupid implementation of qpms_scatsys_apply_symmetry() ------------- #define MIN(x,y) (((x)<(y))?(x):(y)) @@ -151,6 +167,7 @@ static void process_lattice_bases(qpms_scatsys_t *ss, const qpms_tolerance_spec_ normsq = cart3_normsq(ss->per->lattice_basis[0]); ss->per->unitcell_volume = sqrt(normsq); ss->per->reciprocal_basis[0] = cart3_divscale(ss->per->lattice_basis[0], normsq); + ss->per->eta = 2. / M_2_SQRTPI / ss->per->unitcell_volume; } break; case 2: { // (I) Create orthonormal basis of the plane in which the basis vectors lie @@ -169,6 +186,7 @@ static void process_lattice_bases(qpms_scatsys_t *ss, const qpms_tolerance_spec_ // (III) Reduce lattice basis and get reciprocal bases in the 2D plane l2d_reduceBasis(b2d[0], b2d[1], &b2d[0], &b2d[1]); ss->per->unitcell_volume = l2d_unitcell_area(b2d[0], b2d[1]); + ss->per->eta = 2. / M_2_SQRTPI / sqrt(ss->per->unitcell_volume); cart2_t r2d[2]; QPMS_ENSURE_SUCCESS(l2d_reciprocalBasis1(b2d[0], b2d[1], &r2d[0], &r2d[1])); // (IV) Rotate everything back to the original 3D space. @@ -185,6 +203,7 @@ static void process_lattice_bases(qpms_scatsys_t *ss, const qpms_tolerance_spec_ l3d_reduceBasis(ss->per->lattice_basis, ss->per->lattice_basis); ss->per->unitcell_volume = fabs(l3d_reciprocalBasis1(ss->per->lattice_basis, ss->per->reciprocal_basis)); + ss->per->eta = 2. / M_2_SQRTPI / cbrt(ss->per->unitcell_volume); // TODO check unitcell_volume for sanity w.r.t tolerance } break; default: @@ -1134,6 +1153,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_full( qpms_bessel_t J ///< Bessel function type. ) { + qpms_ss_ensure_nonperiodic(ss); const size_t full_len = ss->fecv_size; if(!target) QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double)); @@ -1164,7 +1184,50 @@ complex double *qpms_scatsys_build_translation_matrix_e_full( return target; } +static inline int qpms_ss_ppair_W32xy(qpms_scatsys_t *ss, + qpms_ss_pi_t pdest, qpms_ss_pi_t psrc, complex double wavenumber, + const cart2_t kvector, + complex double *target, const ptrdiff_t deststride, const ptrdiff_t srcstride, + qpms_ewald_part parts) { + const qpms_vswf_set_spec_t *srcspec = qpms_ss_bspec_pi(ss, psrc); + const qpms_vswf_set_spec_t *destspec = qpms_ss_bspec_pi(ss, pdest); + // This might be a bit arbitrary, roughly "copied" from Unitcell constructor. TODO review. + const double maxR = sqrt(ss->per->unitcell_volume) * 32; + const double maxK = 1024 * 2 * M_PI / maxR; + + return qpms_trans_calculator_get_trans_array_e32_e(ss->c, + target, NULL /*err*/, destspec, deststride, srcspec, srcstride, + ss->eta, wavenumber, + cart3xy2cart2(ss->per->lattice_basis[0]), cart3xy2cart2(ss->per->lattice_basis[1]), + kvector, + cart2_substract(cart3xy2cart2(ss->p[pdest].pos), cart3xy2cart2(ss->p[psrc].pos)) + u->maxR, u->maxK, parts); +} + +complex double *qpms_scatsys_periodic_build_translation_matrix_full( + complex double *target, const qpms_scatsys_t *ss, + complex double wavenumber, const cart3_t *wavevector) { + QPMS_UNTESTED; + qpms_ss_ensure_periodic(ss); + const size_t full_len = ss->fecv_size; + if(!target) + QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double)); + const ptrdiff_t deststride = ss->fecv_size, srcstride = 1; + // We have some limitations in the current implementation + if(ss->lattice_dimension == 2 && // Currently, we can only the xy-plane + !ss->per->lattice_basis[0].z && !ss->per_lattice_basis[1].z && + !wavevector->z) { + for (qpms_ss_pi_t pd = 0; pd < ss->p_count; ++pd) + for (qpms_ss_pi_t ps = 0; ps < ss->p_count; ++ps) { + QPMS_ENSURE_SUCCESS(qpms_ss_ppair_W32xy(ss, pd, ps, wavenumber, cart3xy2cart2(*wavevector), + target + deststride * ss->fecv_pstarts[pd] + srcstride * ss->fecv_pstarts[ps], + deststride, srcstride, QPMS_EWALD_FULL)); + } + } else + QPMS_NOT_IMPLEMENTED("Only 2D xy-lattices currently supported"); + return target; +} complex double *qpms_scatsysw_build_modeproblem_matrix_full( /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. @@ -1174,6 +1237,7 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_full( { const complex double k = ssw->wavenumber; const qpms_scatsys_t *ss = ssw->ss; + qpms_ss_ensure_nonperiodic(ss); const size_t full_len = ss->fecv_size; if(!target) QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double)); @@ -1216,6 +1280,24 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_full( return target; } +complex double *qpms_scatsyswk_build_modeproblem_matrix_full( + /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. + complex double *target, + const qpms_scatsys_at_omega_k_t *sswk + ) +{ + const qpms_scatsys_at_omega_t *ssw = sswk->ssw; + const complex double k = ssw->wavenumber; + const qpms_scatsys_t *ss = ssw->ss; + qpms_ss_ensure_periodic(ss); + const size_t full_len = ss->fecv_size; + if(!target) + QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double)); + //////////// TODO //////////////// + QPMS_NOT_IMPLEMENTED("TODO"); + return target; +} + // Serial reference implementation. complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial( /// Target memory with capacity for ss->saecv_sizes[iri]**2 elements. If NULL, new will be allocated. @@ -1225,6 +1307,7 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial( ) { const qpms_scatsys_t *ss = ssw->ss; + qpms_ss_ensure_nonperiodic(ss); const complex double k = ssw->wavenumber; const size_t packedlen = ss->saecv_sizes[iri]; if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps? @@ -1335,6 +1418,7 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_orbitorderR( ) { const qpms_scatsys_t *ss = ssw->ss; + qpms_ss_ensure_nonperiodic(ss); const complex double k = ssw->wavenumber; const size_t packedlen = ss->saecv_sizes[iri]; if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps? @@ -1726,6 +1810,7 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed( qpms_bessel_t J ) { + qpms_ss_ensure_nonperiodic(ss); QPMS_UNTESTED; const size_t packedlen = ss->saecv_sizes[iri]; if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps? @@ -1778,6 +1863,7 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed( const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri ) { + qpms_ss_ensure_nonperiodic(ssw->ss); const size_t packedlen = ssw->ss->saecv_sizes[iri]; if (!packedlen) // THIS IS A BIT PROBLEMATIC, TODO how to deal with empty irreps? return target_packed; @@ -1904,7 +1990,14 @@ void qpms_ss_LU_free(qpms_ss_LU lu) { } qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise(complex double *mpmatrix_full, - int *target_piv, const qpms_scatsys_at_omega_t *ssw) { + int *target_piv, const qpms_scatsys_at_omega_t *ssw, const qpms_scatsys_at_k_omega_t *sswk) { + if (sswk) { + QPMS_ASSERT(sswk->ssw == ssw || !ssw); + ssw = sswk->ssw; + QPMS_ASSERT(ssw->ss->lattice_dimension > 0); + } else { + QPMS_ASSERT(ssw->ss->lattice_dimension == 0); + } const qpms_scatsys_t *ss = ssw->ss; QPMS_ENSURE(mpmatrix_full, "A non-NULL pointer to the pre-calculated mode matrix is required"); if (!target_piv) QPMS_CRASHING_MALLOC(target_piv, ss->fecv_size * sizeof(int)); @@ -1914,6 +2007,7 @@ qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise(complex double *mpmat lu.a = mpmatrix_full; lu.ipiv = target_piv; lu.ssw = ssw; + lu.sswk = sswk; lu.full = true; lu.iri = -1; return lu; @@ -1922,6 +2016,7 @@ qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise(complex double *mpmat qpms_ss_LU qpms_scatsysw_modeproblem_matrix_irrep_packed_factorise(complex double *mpmatrix_packed, int *target_piv, const qpms_scatsys_at_omega_t *ssw, qpms_iri_t iri) { QPMS_ENSURE(mpmatrix_packed, "A non-NULL pointer to the pre-calculated mode matrix is required"); + qpms_scatsys_ensure_nonperiodic(ssw->ss); size_t n = ssw->ss->saecv_sizes[iri]; if (!target_piv) QPMS_CRASHING_MALLOC(target_piv, n * sizeof(int)); QPMS_ENSURE_SUCCESS(LAPACKE_zgetrf(LAPACK_ROW_MAJOR, n, n, @@ -1938,8 +2033,16 @@ qpms_ss_LU qpms_scatsysw_modeproblem_matrix_irrep_packed_factorise(complex doubl qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_full_LU( complex double *target, int *target_piv, const qpms_scatsys_at_omega_t *ssw){ + qpms_scatsys_ensure_nonperiodic_a(ssw->ss, "qpms_scatsyswk_build_modeproblem_matrix_full_LU()"); target = qpms_scatsysw_build_modeproblem_matrix_full(target, ssw); - return qpms_scatsysw_modeproblem_matrix_full_factorise(target, target_piv, ssw); + return qpms_scatsysw_modeproblem_matrix_full_factorise(target, target_piv, ssw, NULL); +} + +qpms_ss_LU qpms_scatsyswk_build_modeproblem_matrix_full_LU( + complex double *target, int *target_piv, + const qpms_scatsys_at_omega_k_t *sswk){ + target = qpms_scatsyswk_build_modeproblem_matrix_full(target, sswk); + return qpms_scatsysw_modeproblem_matrix_full_factorise(target, target_piv, sswk->ssw, sswk); } qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU( @@ -1990,6 +2093,7 @@ beyn_result_t *qpms_scatsys_finite_find_eigenmodes( complex double omega_centre, double omega_rr, double omega_ri, size_t contour_npoints, double rank_tol, size_t rank_sel_min, double res_tol) { + qpms_ss_ensure_nonperiodic(ss); size_t n; // matrix dimension if (qpms_iri_is_valid(ss->sym, iri)) { n = ss->saecv_sizes[iri]; diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 051d9c4..d82ff40 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -146,6 +146,9 @@ typedef struct qpms_scatsys_periodic_info_t { * lattice_dimension == 2, a 2D area. */ double unitcell_volume; + + /// Ewald parameter \f$ \eta \f$. + double eta; } qpms_scatsys_periodic_info_t; @@ -240,6 +243,7 @@ typedef struct qpms_scatsys_at_omega_t { complex double wavenumber; ///< Background medium wavenumber } qpms_scatsys_at_omega_t; + /// Creates a new scatsys by applying a symmetry group onto a "proto-scatsys", copying particles if needed. /** In fact, it copies everything except the vswf set specs and qpms_abstract_tmatrix_t instances, * so keep them alive until scatsys is destroyed. @@ -352,6 +356,16 @@ complex double *qpms_scatsys_build_translation_matrix_full( complex double k ///< Wave number to use in the translation matrix. ); +/// Creates the full \f$ (I - WS) \f$ matrix of the periodic scattering system. NI +/** + * \returns \a target on success, NULL on error. + */ +complex double *qpms_scatsysw_build_modeproblem_matrix_full( + /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. + complex double *target, + const qpms_scatsys_at_omega_t *ssw + ); + /// As qpms_scatsys_build_translation_full() but with choice of Bessel function type. /** Might be useful for evaluation of cross sections and testing. */ @@ -377,15 +391,6 @@ complex double *qpms_scatsys_build_translation_matrix_e_irrep_packed( qpms_bessel_t J ); -/// Creates the full \f$ (I - TS) \f$ matrix of the scattering system. -/** - * \returns \a target on success, NULL on error. - */ -complex double *qpms_scatsysw_build_modeproblem_matrix_full( - /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. - complex double *target, - const qpms_scatsys_at_omega_t *ssw - ); /// Creates the mode problem matrix \f$ (I - TS) \f$ directly in the irrep-packed form. /** * \returns \a target on success, NULL on error. @@ -420,6 +425,7 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_irrep_packed_serial( /// LU factorisation (LAPACKE_zgetrf) result holder. typedef struct qpms_ss_LU { const qpms_scatsys_at_omega_t *ssw; + const qpms_scatsys_at_omega_k_t *sswk; ///< Only for periodic systems, otherwise NULL. bool full; ///< true if full matrix; false if irrep-packed. qpms_iri_t iri; ///< Irrep index if `full == false`. /// LU decomposition array. @@ -429,7 +435,7 @@ typedef struct qpms_ss_LU { } qpms_ss_LU; void qpms_ss_LU_free(qpms_ss_LU); -/// Builds an LU-factorised mode/scattering problem \f$ (I - TS) \f$ matrix from scratch. +/// Builds an LU-factorised mode/scattering problem \f$ (I - TS) \f$ matrix from scratch. Nonperiodic systems only. qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_full_LU( complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated). int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated). @@ -448,7 +454,8 @@ qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_irrep_packed_LU( qpms_ss_LU qpms_scatsysw_modeproblem_matrix_full_factorise( complex double *modeproblem_matrix_full, ///< Pre-calculated mode problem matrix (I-TS). Mandatory. int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated). - const qpms_scatsys_at_omega_t *ssw + const qpms_scatsys_at_omega_t *ssw, ///< Must be filled for non-periodic systems. + const qpms_scatsys_at_omega_k_t *sswk ///< Must be filled for periodic systems, otherwise must be NULL. ); /// Computes LU factorisation of a pre-calculated irrep-packed mode/scattering problem matrix, replacing its contents. @@ -466,6 +473,50 @@ complex double *qpms_scatsys_scatter_solve( qpms_ss_LU ludata ///< Pre-factorised \f$ I - TS \f$ matrix data. ); +// ======================= Periodic system -only related stuff ============================= + +/// Scattering system at a given frequency and k-vector. Used only with periodic systems. +typedef struct qpms_scatsys_at_omega_k_t { + const qpms_scatsys_at_omega_t *ssw; + double k[3]; ///< The k-vector's cartesian coordinates. +} qpms_scatsys_at_omega_k_t; + +/// Creates the full \f$ (I - WS) \f$ matrix of the periodic scattering system. NI +/** + * \returns \a target on success, NULL on error. + */ +complex double *qpms_scatsyswk_build_modeproblem_matrix_full( + /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. + complex double *target, + const qpms_scatsys_at_omega_k_t *sswk + ); + +/// Global translation matrix. +complex double *qpms_scatsys_periodic_build_translation_matrix_full( + /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. + complex double *target, + const qpms_scatsys_t *ss, + complex double wavenumber, ///< Wave number to use in the translation matrix. + const double wavevector[] ///< Wavevector / pseudomomentum in cartesian coordinates. + ); + +/// Global translation matrix. NI +complex double *qpms_scatsyswk_build_translation_matrix_full( + /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. + complex double *target, + const qpms_scatsys_omega_k_t *sswk + ); + + +/// Builds an LU-factorised mode/scattering problem \f$ (I - TS) \f$ matrix from scratch. Periodic systems only. +qpms_ss_LU qpms_scatsysw_build_modeproblem_matrix_full_LU( + complex double *target, ///< Pre-allocated target array. Optional (if NULL, new one is allocated). + int *target_piv, ///< Pre-allocated pivot array. Optional (if NULL, new one is allocated). + const qpms_scatsys_at_omega_k_t *sswk + ); + +// ======================= Periodic system -only related stuff end ========================= + /// NOT IMPLEMENTED Dumps a qpms_scatsys_t structure to a file. qpms_errno_t qpms_scatsys_dump(qpms_scatsys_t *ss, char *path);