Benchmark testing

[christopher.jessop]

Hello,

I am trying to benchmark one of the built-in networks in KROME, specifically react_primordialZ. I am specifically trying to plot figure 5 from the Omukai 2005 paper, which the reactions in the network are taken from, so I ran the collapseZ test which worked fine for its intended purpose of recreating figure 9 from Omukai 2005 for the thermal evolution. The problem comes when I try to plot the abundances of H2O, O, O2 and OH against number density. The output fort.22 seems to have the columns for these abundances, but plotting them gives massively different results than the figures in the paper. I altered the plot.gps script to plot columns 9-12, so the plot for Z=-3:

reset
set logscale
set autoscale
set format y "%L"
set format x "%L"

set ylabel "log y(Abundance H2O, O, O2, OH)"
set xlabel "log(n/cm-3)"
set title "Z = -3"
set yrange [1e-13:1e-6]


set key b r
set grid

plot './fort.22' u 2:($1==-3?$11:1/0) w l t "H2O",\
'' u 2:($1==-3?$9:1/0) w l t "O",\
'' u 2:($1==-3?$10:1/0) w l t "OH",\
'' u 2:($1==-3?$12:1/0) w l t "O2",\

In addition to this I altered the test.f90 script to only run for one metallicity, which produced the same plot I get when I try the above script for the standard test.f90, just to double check. I understand that this was not the intended purpose of the test, and dust cooling is not present so there may be differences, but I'd have thought the plots would have been closer. If anyone has any suggestions as to why this may be, what I am doing incorrectly regarding looking at the data output or the way I'm plotting, it would be greatly appreciated. I will also attach the graph I produce and also the specific figure I am trying to reproduce.

I have tried using python, such as the script plot.py but I've no idea what the module kp is and can not find any information online about it. Trying to load the data using techniques I know in python seems more confusing than gnuplot because of the data output format, so that is why I have tried to stick to using gnuplot for this.

Cheers,

Chris Jessop

[tommaso grassi]

Hi Chris,

thanks for using KROME.

The reason of the discrepancy is that H2 formation on dust is not implemented in the test, so there is too much H that destroys oxygen compounds. This was beyond the aims of the test. You can have an idea by changing rate 25 in krome_subs.f90 as follows

   !H + H + H -> H2 + H

    Td = 3d0

    fa = 1d0 / (1 + exp(7.5d2*(1d0/75. - 1d0/Td)))

    k(25) = small + (6.d-32&

         *Tgas**(-0.25d0)+2.d-31*Tgas**(-0.5d0)) &

         + 6d-17*sqrt(Tgas/3e2) * fa / (1 + &

         4d-2 * sqrt(Tgas + Td) + 2d-3*Tgas + 8d-6*Tgas**2) * 1d-3 / n(idx_H)

(this is a quick trick to reuse 3H->H2+H as 2H->H2). I use 3K, that holds up to about 1e8 cm-3, after that density Td should be computed properly (more details here https://arxiv.org/abs/1606.01229).

Additionally, initial conditions for metals are slightly different wrt Omukai, since he has

     x(krome_idx_C) = 0.927d-4 * 1d1**zs(jz)

     x(krome_idx_O) = 3.568d-4 * 1d1**zs(jz)

Hope this explains what you found.

Cheers,

-tg

[christopher.jessop]

Thanks for the explanation, I'll have a look at those changes.

Cheers,

Chris

[christopher.jessop]

Hello,

Just to follow up, I implemented the changes to krome_subs.f90 and changed the initial conditions in test.f90 to match that of Omukai, and the figures do look much more similar. Of course there are still quite large differences but I expected that to be the case. Is there a way to potentially build in H2 formation on dust into the react_primordialZ network? I notice that the change to krome_subs is the H2 formation on dust reaction from the paper, would there be a way to implement this into the reaction network itself rather than using the trick to reuse 3H->H2+H as 2H->H2, such as just adding the reaction to the network after defining new variables for fa and Td?

Thanks,

Chris

[tommaso grassi]

You can add the reaction in the chemical network as written in Omukai's paper.

The trick with 3H I used in my previous reply was for a quick test.

However, you need to know Tdust, that usually comes from DUST cooling, that is not implemented in the collapseZ test.

Alternatively use dust tables as in collapseCO test, i.e. the option -dustTabs=H2,COOL,MC

[christopher.jessop]

Thank you for your suggestions. I have tried a couple of things, firstly I altered the react_primordialZ network to include:

@format:idx,R,R,R,P,rate
@var:Td = 3d0
@var:fa = 1d0 / (1 + exp(7.5d2*(1d0/75. - 1d0/Td)))
75,H,H,GRAIN0,H2,(6.d-32 * Tgas**(-0.25d0)+2.d-31*Tgas**(-0.5d0)) + 6d-17*sqrt(Tgas/3e2) * fa / (1 + 4d-2 * sqrt(Tgas + Td) + 2d-3*Tgas + 8d-6*Tgas**2) * 1d-3 / n(idx_H)

Using Td = 3d0 just as a placeholder, as you said I need that from dust cooling, but it doesn't look to have the same effect as changing k(25) from the workaround, and it also produces a mass conservation warning. Secondly, I tried to change options.opt to read:

-dustTabs=H2,COOL,MC
-useTabs

but that again does not seem to give nearly the same as the workaround. Could you please point out to me where I may be going incorrect trying to add H2 formation on dust into the reaction network? Apologies if some of these questions seem quite naive/simple.

Thanks,

Chris

[tommaso grassi]

Hi,

you have to remove GRAIN0 from the reactants to conserve the mass, or add GRAIN0 in the products.

Then remove (6.d-32 * Tgas**(-0.25d0)+2.d-31*Tgas**(-0.5d0)), that was the original H + H + H -> H2 + H

I admit that the workaround is confusing.

H2 formation on dust is a complicated issue that depends a lot on the boundary conditions and on the assumptions.

I recommend this paper for a good review on the problem https://arxiv.org/abs/1711.10568

For the dust tables please take a look to the method paper, because there are some assumptions that has to be considered

https://arxiv.org/abs/1606.01229

cheers,

-tg