From a5b2e5c8a7dcf3c8fe1109fdd62d3bcbcc5a7bdb Mon Sep 17 00:00:00 2001
From: Robert Izzard <r.izzard@surrey.ac.uk>
Date: Sat, 10 Jul 2021 17:21:28 +0100
Subject: [PATCH] fix errors in he-burning jacobian: works up to 1e9K now
---
src/nucsyn/nucsyn_burn_helium.c | 12 ++++++------
src/nucsyn/nucsyn_hbb.c | 8 ++++----
2 files changed, 10 insertions(+), 10 deletions(-)
diff --git a/src/nucsyn/nucsyn_burn_helium.c b/src/nucsyn/nucsyn_burn_helium.c
index 99466f496..e0ba2a650 100644
--- a/src/nucsyn/nucsyn_burn_helium.c
+++ b/src/nucsyn/nucsyn_burn_helium.c
@@ -109,7 +109,7 @@ static void helium_Jacobian(double *const*const RESTRICT J,
const double * RESTRICT const N,
const double * RESTRICT const sigmav)
{
- const double y2 = 3.0 * Pow2(N[XHe4]) * sigmav[SIGMAV_3He4];
+ const double dummy0 = 3.0 * Pow2(N[XHe4]) * sigmav[SIGMAV_3He4];
const double dummy1 = sigmav[SIGMAV_C12_He4] * N[XC12];
const double dummy2 = sigmav[SIGMAV_O16_He4] * N[XO16];
const double dummy3 = sigmav[SIGMAV_Ne20_He4] * N[XNe20];
@@ -118,7 +118,7 @@ static void helium_Jacobian(double *const*const RESTRICT J,
const double dummy6 = sigmav[SIGMAV_Ne20_He4] * N[XHe4];
J[1][1] =
- -3.0 * y2
+ -3.0 * dummy0
-dummy1
-dummy2
-dummy3;
@@ -128,10 +128,10 @@ static void helium_Jacobian(double *const*const RESTRICT J,
J[1][5] = 0.0;
J[2][1] =
- + y2
+ + dummy0
- dummy1;
J[2][2] = - dummy4;
- J[2][3] = - dummy5;
+ J[2][3] = 0.0;
J[2][4] = 0.0;
J[2][5] = 0.0;
@@ -150,9 +150,9 @@ static void helium_Jacobian(double *const*const RESTRICT J,
- dummy3;
J[4][2] = 0.0;
J[4][3] =
- - dummy5;
+ + dummy5;
J[4][4] =
- + dummy6;
+ - dummy6;
J[4][5] = 0.0;
J[5][1] = dummy3;
diff --git a/src/nucsyn/nucsyn_hbb.c b/src/nucsyn/nucsyn_hbb.c
index f7972f4ed..2a9bf9e55 100644
--- a/src/nucsyn/nucsyn_hbb.c
+++ b/src/nucsyn/nucsyn_hbb.c
@@ -29,7 +29,7 @@ void nucsyn_hbb(const double T, /* (log10) temperature */
const Nuclear_mass * mnuc = stardata->store->mnuc;
const Nuclear_mass *imnuc = stardata->store->imnuc;
-
+//#define __TEST_HE_BURN
#ifdef __TEST_HE_BURN
Abundance * X = New_clear_isotope_array;
Number_density * N = New_isotope_array;
@@ -39,10 +39,10 @@ void nucsyn_hbb(const double T, /* (log10) temperature */
X[XHe4] = 1.0;
X_to_N(imnuc,testrho,N,X,ISOTOPE_ARRAY_SIZE);
- const double timemax = 1000;
+ const double timemax = 100;
double time = 0.0;
- double dtime = timemax/1000.0;
- double logtemp = log10(3e8);
+ double dtime = timemax/100.0;
+ double logtemp = log10(1e9);
while(time < timemax)
{
--
GitLab