QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Beyn wrappers for periodic system. 

Former-commit-id: c4c21de7d02af36d35133ff9ef0c426742eab6ff
This commit is contained in:
Marek Nečada 2020-02-28 21:22:36 +02:00
parent 1fb4e5760a
commit 78c793fb68
2 changed files with 67 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
qpms/scatsystem.c
View File

@ -2128,7 +2128,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);
qpms_ss_ensure_nonperiodic_a(ss, "qpms_scatsys_periodic_find_eigenmodes()");
size_t n; // matrix dimension
if (qpms_iri_is_valid(ss->sym, iri)) {
n = ss->saecv_sizes[iri];
@ -2149,3 +2149,47 @@ beyn_result_t *qpms_scatsys_finite_find_eigenmodes(
return result;
}
struct qpms_scatsys_periodic_eval_Beyn_ImTW_param {
const qpms_scatsys_t *ss;
const double *k; ///< Wavevector in cartesian coordinates.
};
/// Wrapper for Beyn algorithm (periodic system)
static int qpms_scatsys_periodic_eval_Beyn_ImTW(complex double *target,
size_t m, complex double omega, void *params){
const struct qpms_scatsys_periodic_eval_Beyn_ImTW_param *p = params;
qpms_scatsys_at_omega_t *ssw = qpms_scatsys_at_omega(p->ss, omega);
QPMS_ENSURE(ssw != NULL, "qpms_scatsys_at_omega() returned NULL");
qpms_scatsys_at_omega_k_t sswk = {
.ssw = ssw,
.k = {p->k[0], p->k[1], p->k[2]}
};
QPMS_ASSERT(m == p->ss->fecv_size);
QPMS_ENSURE(NULL !=
qpms_scatsyswk_build_modeproblem_matrix_full(target, &sswk),
"qpms_scatsyswk_build_modeproblem_matrix_full() returned NULL");
qpms_scatsys_at_omega_free(ssw);
return QPMS_SUCCESS;
}
beyn_result_t *qpms_scatsys_periodic_find_eigenmodes(
const qpms_scatsys_t * const ss, const double k[3],
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_a(ss, "qpms_scatsys_finite_find_eigenmodes()");
size_t n = ss->fecv_size; // matrix dimension
beyn_contour_t *contour = beyn_contour_ellipse(omega_centre,
omega_rr, omega_ri, contour_npoints);
struct qpms_scatsys_periodic_eval_Beyn_ImTW_param p = {ss, k};
beyn_result_t *result = beyn_solve(n, n,
qpms_scatsys_periodic_eval_Beyn_ImTW, NULL, (void *) &p,
contour, rank_tol, rank_sel_min, res_tol);
QPMS_ENSURE(result != NULL, "beyn_solve() returned NULL");
free(contour);
return result;
}

22
qpms/scatsystem.h
View File

@ -519,6 +519,28 @@ qpms_ss_LU qpms_scatsyswk_build_modeproblem_matrix_full_LU(
const qpms_scatsys_at_omega_k_t *sswk
);
/// Searches for periodic scattering system's eigenmodes using Beyn's algorithm.
/**
* Currently, elliptical contour is used.
*
* TODO In the future, this will probably support irrep decomposition as well,
* but for the case of periodic / small systems,
* the bottleneck is the T-matrix and translation matrix evaluation
* rather than the linear algebra.
*/
struct beyn_result_t *qpms_scatsys_periodic_find_eigenmodes(
const qpms_scatsys_t *ss,
/// Wavevector in cartesian coordinates (must lie in the lattice plane).
const double k[3],
complex double omega_centre, ///< Center of the ellipse inside which the eigenfreqs are searched for.
double omega_rr, ///< Real half-axis of the ellipse inside which the eigenfreqs are searched for.
double omega_ri, ///< Imaginary half-axis of the ellipse inside which the eigenfreqs are searched for.
size_t contour_npoints, ///< Number of elliptic contour discretisation points (preferably even number)
double rank_tol, ///< (default: `1e-4`) TODO DOC.
size_t rank_sel_min, ///< Minimum number of eigenvalue candidates, even if they don't pass \a rank_tol.
double res_tol ///< (default: `0.0`) TODO DOC.
);
// ======================= Periodic system -only related stuff end =========================