[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master d6447d2b: Library (cosmology.h): slower/robust
From: |
Mohammad Akhlaghi |
Subject: |
[gnuastro-commits] master d6447d2b: Library (cosmology.h): slower/robust gsl integration when fast fails |
Date: |
Thu, 1 Feb 2024 19:07:10 -0500 (EST) |
branch: master
commit d6447d2b26eec1577faf91d5150a36e77db29f9b
Author: Thorsten Alteholz <thorsten@alteholz.dev>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Library (cosmology.h): slower/robust gsl integration when fast fails
Until now, when CosmicCalculator was given an input redshift that is higher
than 100, GSL would produce a code-dump. This was due to the choice of
'gsl_integration_qng' for the numerical intergration (because it is fast).
With this commit, after 'gsl_integration_qng' is finished we check whether
it was successful or not. If it failed, the slower, but more robust
'gsl_integration_qag' is used. In order to use the 'error' function to
report the failure, the variable that was previously named 'error' was
renamed to 'gslerr'.
---
lib/cosmology.c | 61 ++++++++++++++++++++++++++++++++++++++++++++++++++-------
1 file changed, 54 insertions(+), 7 deletions(-)
diff --git a/lib/cosmology.c b/lib/cosmology.c
index a3735062..62bf775f 100644
--- a/lib/cosmology.c
+++ b/lib/cosmology.c
@@ -28,6 +28,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <stdio.h>
#include <stdlib.h>
+#include <gsl/gsl_errno.h>
#include <gsl/gsl_const_mksa.h>
#include <gsl/gsl_integration.h>
@@ -249,9 +250,12 @@ gal_cosmology_proper_distance(double z, double H0, double
o_lambda_0,
double o_matter_0, double o_radiation_0)
{
cosmology_density_check(o_lambda_0, o_matter_0, o_radiation_0);
+ int status;
size_t neval;
gsl_function F;
- double result, error, c=GSL_CONST_MKSA_SPEED_OF_LIGHT;
+ gsl_integration_workspace *w;
+ gsl_error_handler_t *gslerrhandler;
+ double result, gslerr, c=GSL_CONST_MKSA_SPEED_OF_LIGHT;
double o_curv_0 = 1.0 - ( o_lambda_0 + o_matter_0 + o_radiation_0 );
double H0s=H0/1000/GSL_CONST_MKSA_PARSEC; /* H0 in units of seconds. */
struct cosmology_integrand_t p={o_lambda_0, o_curv_0, o_matter_0,
@@ -261,9 +265,29 @@ gal_cosmology_proper_distance(double z, double H0, double
o_lambda_0,
F.params=&p;
F.function=&cosmology_integrand_proper_dist;
- /* Do the integration. */
- gsl_integration_qng(&F, 0.0f, z, GSLIEPSABS, GSLIEPSREL, &result,
- &error, &neval);
+ /* Temporarily switch off error handling */
+ gslerrhandler=gsl_set_error_handler_off();
+
+ /* Do the integration (first with the fast GSL function). */
+ status=gsl_integration_qng(&F, 0.0f, z, GSLIEPSABS, GSLIEPSREL, &result,
+ &gslerr, &neval);
+
+ /* If the first integration failed, try a slower, but more robust one. */
+ if (status!=GSL_SUCCESS)
+ {
+ w=gsl_integration_workspace_alloc(GSLILIMIT);
+ status=gsl_integration_qag(&F, 0.0f, z, GSLIEPSABS, GSLIEPSREL,
+ GSLILIMIT, GSL_INTEG_GAUSS21, w,
+ &result, &gslerr);
+ gsl_integration_workspace_free(w);
+ if (status!=GSL_SUCCESS)
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' to "
+ "fix the problem. The status of the second integration is "
+ "%i.", __func__, PACKAGE_BUGREPORT, status);
+ }
+
+ /* Reactivate "normal" error handling */
+ gslerrhandler=gsl_set_error_handler(gslerrhandler);
/* Return the result. */
return result * c / H0s / (1e6 * GSL_CONST_MKSA_PARSEC);
@@ -279,9 +303,12 @@ gal_cosmology_comoving_volume(double z, double H0, double
o_lambda_0,
double o_matter_0, double o_radiation_0)
{
cosmology_density_check(o_lambda_0, o_matter_0, o_radiation_0);
+ int status;
size_t neval;
gsl_function F;
- double result, error;
+ double result, gslerr;
+ gsl_integration_workspace *w;
+ gsl_error_handler_t *gslerrhandler;
double c=GSL_CONST_MKSA_SPEED_OF_LIGHT;
double H0s=H0/1000/GSL_CONST_MKSA_PARSEC; /* H0 in units of seconds. */
double cH = c / H0s / (1e6 * GSL_CONST_MKSA_PARSEC);
@@ -293,9 +320,29 @@ gal_cosmology_comoving_volume(double z, double H0, double
o_lambda_0,
F.params=&p;
F.function=&cosmology_integrand_comoving_volume;
+ /* temporarily switch off error handling */
+ gslerrhandler=gsl_set_error_handler_off();
+
/* Do the integration. */
- gsl_integration_qng(&F, 0.0f, z, GSLIEPSABS, GSLIEPSREL,
- &result, &error, &neval);
+ status=gsl_integration_qng(&F, 0.0f, z, GSLIEPSABS, GSLIEPSREL,
+ &result, &gslerr, &neval);
+
+ /* If the first integration failed, try a slower, but more robust one. */
+ if (status!=GSL_SUCCESS)
+ {
+ w=gsl_integration_workspace_alloc(GSLILIMIT);
+ status=gsl_integration_qag(&F, 0.0f, z, GSLIEPSABS, GSLIEPSREL,
+ GSLILIMIT, GSL_INTEG_GAUSS21, w,
+ &result, &gslerr);
+ gsl_integration_workspace_free(w);
+ if (status!=GSL_SUCCESS)
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' to "
+ "fix the problem. The status of the second integration is "
+ "%i.", __func__, PACKAGE_BUGREPORT, status);
+ }
+
+ /* Reactivate "normal" error handling */
+ gslerrhandler=gsl_set_error_handler(gslerrhandler);
/* Return the result. */
return result * 4 * M_PI * cH*cH*cH;
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master d6447d2b: Library (cosmology.h): slower/robust gsl integration when fast fails,
Mohammad Akhlaghi <=