WIP 2D in 3D -periodic scatsystem

Former-commit-id: e53b6b1f8361efc6c66260e776688b4e941660ad

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];

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);