From 1821a0d8f58a2234f54c95ed29c0f4b901596cd1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 12 Nov 2019 07:18:20 +0200 Subject: [PATCH] Try to recover from failed quadrature. Former-commit-id: 746709b224d825baf4a41ddfcd040a0c80fc45c8 --- qpms/cycommon.pyx | 2 ++ qpms/qpms_cdefs.pxd | 1 + qpms/qpms_error.h | 3 ++- qpms/tmatrices.c | 46 +++++++++++++++++++++++++++++++++++++-------- 4 files changed, 43 insertions(+), 9 deletions(-) diff --git a/qpms/cycommon.pyx b/qpms/cycommon.pyx index 118c047..88c557d 100644 --- a/qpms/cycommon.pyx +++ b/qpms/cycommon.pyx @@ -47,11 +47,13 @@ try: class DebugFlags(enum.IntFlag): # Should be IntFlag if python version >= 3.6 MISC = QPMS_DBGMSG_MISC THREADS = QPMS_DBGMSG_THREADS + INTEGRATION = QPMS_DBGMSG_INTEGRATION has_IntFlag = True except AttributeError: # For old versions of enum, use IntEnum instead class DebugFlags(enum.IntEnum): MISC = QPMS_DBGMSG_MISC THREADS = QPMS_DBGMSG_THREADS + INTEGRATION = QPMS_DBGMSG_INTEGRATION has_IntFlag = False def dbgmsg_enable(qpms_dbgmsg_flags types): diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 10ce555..ee02028 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -146,6 +146,7 @@ cdef extern from "qpms_error.h": ctypedef enum qpms_dbgmsg_flags: QPMS_DBGMSG_MISC QPMS_DBGMSG_THREADS + QPMS_DBGMSG_INTEGRATION qpms_dbgmsg_flags qpms_dbgmsg_enable(qpms_dbgmsg_flags types) qpms_dbgmsg_flags qpms_dbgmsg_disable(qpms_dbgmsg_flags types) diff --git a/qpms/qpms_error.h b/qpms/qpms_error.h index 9919296..10fb7b8 100644 --- a/qpms/qpms_error.h +++ b/qpms/qpms_error.h @@ -29,7 +29,8 @@ void qpms_pr_debug_at_flf(const char *filename, unsigned int linenum, typedef enum { QPMS_DBGMSG_MISC = 1, - QPMS_DBGMSG_THREADS = 2 // Multithreading-related debug messages. + QPMS_DBGMSG_THREADS = 2, // Multithreading-related debug messages. + QPMS_DBGMSG_INTEGRATION = 4 // Quadrature-related debug messages. } qpms_dbgmsg_flags; void qpms_debug_at_flf(const char *filename, unsigned int linenum, diff --git a/qpms/tmatrices.c b/qpms/tmatrices.c index 78a1eea..0ee828a 100644 --- a/qpms/tmatrices.c +++ b/qpms/tmatrices.c @@ -25,6 +25,8 @@ #include #include #include +#include +#include "optim.h" // These are quite arbitrarily chosen constants for the quadrature in qpms_tmatrix_axialsym_fill() #define TMATRIX_AXIALSYM_INTEGRAL_EPSREL (1e-5) @@ -712,6 +714,38 @@ struct qpms_tmatrix_axialsym_fill_integration_thread_arg { pthread_mutex_t *i2start_mutex; }; +// This wrapper retries to calculate the T-matrix integral and if it fails, reduces torelance. +static inline int axint_wrapper(const gsl_function *f, double epsabs, double epsrel, + size_t intlimit, gsl_integration_workspace *w, double *result, double *abserr) { + static int tol_exceeded = 0; + double errfac; + int retval; + for (errfac = 1; epsabs*errfac < 0.01 && epsrel*errfac < 0.01; errfac *= 8) { + retval = gsl_integration_qag(f, 0, M_PI, epsabs*errfac, epsrel*errfac, + intlimit, 2, w, result, abserr); + if(QPMS_LIKELY(!retval)) { + if(QPMS_LIKELY(errfac == 1)) + return 0; + else { + QPMS_DEBUG(QPMS_DBGMSG_INTEGRATION, + "T-matrix quadrature succeeded only after %g-fold" + "tolerance reduction."); // TODO more details here? + return GSL_EROUND; + } + } + else if(retval == GSL_EROUND) { + if (!tol_exceeded) { + QPMS_WARN("T-matrix quadrature failed; retrying with worse tolerances." + " (this warning will be shown only once)."); + } + tol_exceeded += 1; + } + else QPMS_ENSURE_SUCCESS(retval); + } + QPMS_PR_ERROR("T-matrix quadrature failed even after harsh tolerance reduction."); +} + + static void * qpms_tmatrix_axialsym_fill_integration_thread(void *arg) { const struct qpms_tmatrix_axialsym_fill_integration_thread_arg *a = arg; struct tmatrix_axialsym_integral_param_t p = *a->p_template; @@ -749,27 +783,23 @@ static void * qpms_tmatrix_axialsym_fill_integration_thread(void *arg) { // Re(R') p.btype = QPMS_BESSEL_REGULAR; p.realpart = true; - QPMS_ENSURE_SUCCESS(gsl_integration_qag(&f, 0, M_PI, a->epsabs, a->epsrel, - a->intlimit, 2, w, &result, &abserr)); + axint_wrapper(&f, a->epsabs, a->epsrel, a->intlimit, w, &result, &abserr); a->R[iQR] = result; // Im(R') p.realpart = false; - QPMS_ENSURE_SUCCESS(gsl_integration_qag(&f, 0, M_PI, a->epsabs, a->epsrel, - a->intlimit, 2, w, &result, &abserr)); + axint_wrapper(&f, a->epsabs, a->epsrel, a->intlimit, w, &result, &abserr); a->R[iQR] += I*result; // Re(Q') p.btype = QPMS_HANKEL_PLUS; p.realpart = true; - QPMS_ENSURE_SUCCESS(gsl_integration_qag(&f, 0, M_PI, a->epsabs, a->epsrel, - a->intlimit, 2, w, &result, &abserr)); + axint_wrapper(&f, a->epsabs, a->epsrel, a->intlimit, w, &result, &abserr); a->Q[iQR] = result; // Im(Q') p.realpart = false; - QPMS_ENSURE_SUCCESS(gsl_integration_qag(&f, 0, M_PI, a->epsabs, a->epsrel, - a->intlimit, 2, w, &result, &abserr)); + axint_wrapper(&f, a->epsabs, a->epsrel, a->intlimit, w, &result, &abserr); a->Q[iQR] += I*result; } }