Skip to content
Snippets Groups Projects
Commit a5b2e5c8 authored by Izzard, Robert Dr (Maths & Physics)'s avatar Izzard, Robert Dr (Maths & Physics)
Browse files

fix errors in he-burning jacobian: works up to 1e9K now

parent 8e383ae0
No related branches found
No related tags found
No related merge requests found
...@@ -109,7 +109,7 @@ static void helium_Jacobian(double *const*const RESTRICT J, ...@@ -109,7 +109,7 @@ static void helium_Jacobian(double *const*const RESTRICT J,
const double * RESTRICT const N, const double * RESTRICT const N,
const double * RESTRICT const sigmav) 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 dummy1 = sigmav[SIGMAV_C12_He4] * N[XC12];
const double dummy2 = sigmav[SIGMAV_O16_He4] * N[XO16]; const double dummy2 = sigmav[SIGMAV_O16_He4] * N[XO16];
const double dummy3 = sigmav[SIGMAV_Ne20_He4] * N[XNe20]; const double dummy3 = sigmav[SIGMAV_Ne20_He4] * N[XNe20];
...@@ -118,7 +118,7 @@ static void helium_Jacobian(double *const*const RESTRICT J, ...@@ -118,7 +118,7 @@ static void helium_Jacobian(double *const*const RESTRICT J,
const double dummy6 = sigmav[SIGMAV_Ne20_He4] * N[XHe4]; const double dummy6 = sigmav[SIGMAV_Ne20_He4] * N[XHe4];
J[1][1] = J[1][1] =
-3.0 * y2 -3.0 * dummy0
-dummy1 -dummy1
-dummy2 -dummy2
-dummy3; -dummy3;
...@@ -128,10 +128,10 @@ static void helium_Jacobian(double *const*const RESTRICT J, ...@@ -128,10 +128,10 @@ static void helium_Jacobian(double *const*const RESTRICT J,
J[1][5] = 0.0; J[1][5] = 0.0;
J[2][1] = J[2][1] =
+ y2 + dummy0
- dummy1; - dummy1;
J[2][2] = - dummy4; J[2][2] = - dummy4;
J[2][3] = - dummy5; J[2][3] = 0.0;
J[2][4] = 0.0; J[2][4] = 0.0;
J[2][5] = 0.0; J[2][5] = 0.0;
...@@ -150,9 +150,9 @@ static void helium_Jacobian(double *const*const RESTRICT J, ...@@ -150,9 +150,9 @@ static void helium_Jacobian(double *const*const RESTRICT J,
- dummy3; - dummy3;
J[4][2] = 0.0; J[4][2] = 0.0;
J[4][3] = J[4][3] =
- dummy5; + dummy5;
J[4][4] = J[4][4] =
+ dummy6; - dummy6;
J[4][5] = 0.0; J[4][5] = 0.0;
J[5][1] = dummy3; J[5][1] = dummy3;
......
...@@ -29,7 +29,7 @@ void nucsyn_hbb(const double T, /* (log10) temperature */ ...@@ -29,7 +29,7 @@ void nucsyn_hbb(const double T, /* (log10) temperature */
const Nuclear_mass * mnuc = stardata->store->mnuc; const Nuclear_mass * mnuc = stardata->store->mnuc;
const Nuclear_mass *imnuc = stardata->store->imnuc; const Nuclear_mass *imnuc = stardata->store->imnuc;
//#define __TEST_HE_BURN
#ifdef __TEST_HE_BURN #ifdef __TEST_HE_BURN
Abundance * X = New_clear_isotope_array; Abundance * X = New_clear_isotope_array;
Number_density * N = New_isotope_array; Number_density * N = New_isotope_array;
...@@ -39,10 +39,10 @@ void nucsyn_hbb(const double T, /* (log10) temperature */ ...@@ -39,10 +39,10 @@ void nucsyn_hbb(const double T, /* (log10) temperature */
X[XHe4] = 1.0; X[XHe4] = 1.0;
X_to_N(imnuc,testrho,N,X,ISOTOPE_ARRAY_SIZE); X_to_N(imnuc,testrho,N,X,ISOTOPE_ARRAY_SIZE);
const double timemax = 1000; const double timemax = 100;
double time = 0.0; double time = 0.0;
double dtime = timemax/1000.0; double dtime = timemax/100.0;
double logtemp = log10(3e8); double logtemp = log10(1e9);
while(time < timemax) while(time < timemax)
{ {
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment