Unrealistically large abundances

[Bradley S. Meyer]

Moved post:

Dear Dr brad:
I ran the code just like you metioned in this post and I got a correct result.But,when I change the mass excess in my_net.xml and renamed it as my_net_WS.xml,I got somethig like this:
time(s) = 1000000 t9 = 1.00E-04 rho(g/cc) = 1.50E-09
0 0.00E+00
1 -4.85E+98
2 7.76E+111
3 2.46E+118
4 -4.13E+89
5 0.00E+00
6 0.00E+00
7 4.36E+117
8 -7.58E+158
9 2.65E+205
10 -1.86E+126
11 1.79E+140
12 1.98E+144
13 1.52E+142
14 2.43E+147
15 2.84E+140
16 -3.88E+152
17 1.21E+202
18 -6.24E+144
19 3.14E+127
20 -8.21E+129
21 3.27E+197

[Bradley S. Meyer]

It is likely that changing a mass excess leads to some rate blowing up somewhere during the calculation. You can check that by using the libnucnet example code check_rates. Do

cd nucnet-tools-code/build

make -f Makefile.libnucnet all_libnucnet

cd ../libnucnet

./check_rates ../data_pub/my_net_WS.xml 1.e-5 10 0. 1.e30

That will compute the rates at a variety of temperatures between T₉ = 1.e-5 and 10 and check that any are less than zero or greater than 1.e30. You can play with the temperature range and rate range to identify which, if any, rates blow up.

It may be that your mass excess change makes a previously exothermic rate endothermic. This will cause the reverse rate to be exothermic and blow up at low termperature. Let me know what happens. Best wishes.

Brad Meyer

[Anonymous]

Dear Brad:
Thank you for your reply.
After I ran
cd nucnet-tools-code/build
make -f Makefile.libnucnet all_libnucnet
cd ../libnucnet
./check_rates ../data_pub/my_net_WS.xml 1.e-5 10 0. 1.e30
I got some information like the following:
n + nb136 -> nb137 + gamma: upper bound is 1e+30 and reverse rate at t9 = 0.178989 is 9.99079e+34
n + cf310 -> cf311 + gamma: upper bound is 1e+30 and reverse rate at t9 = 0.0122474 is 7.96014e+35
n + fm322 -> fm323 + gamma: upper bound is 1e+30 and reverse rate at t9 = 0.0468206 is 4.15357e+31
n + ac281 -> ac282 + gamma: upper bound is 1e+30 and reverse rate at t9 = 0.0468206 is 5.49025e+35
n + co90 -> co91 + gamma: upper bound is 1e+30 and reverse rate at t9 = 0.0468206 is 4.87642e+32
n + xe174 -> xe175 + gamma: upper bound is 1e+30 and reverse rate at t9 = 0.0468206 is 1.54452e+30
n + cf316 -> cf317 + gamma: upper bound is 1e+30 and reverse rate at t9 = 0.0468206 is 1.05439e+38
n + ge100 -> ge101 + gamma: upper bound is 1e+30 and reverse rate at t9 = 0.0468206 is 7.25023e+35
n + th282 -> th283 + gamma: upper bound is 1e+30 and reverse rate at t9 = 0.0468206 is 4.14908e+42
n + ne37 -> ne38 + gamma: upper bound is 1e+30 and reverse rate at t9 = 0.0122474 is 3.08414e+49
......
What can I do to make it right?
Thanks again!Best wishes!

[Bradley S. Meyer]

Could you try increasing the rate cutoff until you find the largest rates? Try rerunning check_rates but with 1.e40 instead of 1.e30. If there are still large rates try 1.e50 instead of 1.e40, etc. What are the largest rates you find?

[Anonymous]

When I ran check_rates with cutoff 1.e45,I still find large rate,like:
n + zr131 -> zr132 + gamma: upper bound is 1e+45 and reverse rate at t9 = 0.0398107 is 1.57629e+45
n + th282 -> th283 + gamma: upper bound is 1e+45 and reverse rate at t9 = 0.0398107 is 4.47886e+48
n + th286 -> th287 + gamma: upper bound is 1e+45 and reverse rate at t9 = 0.00251189 is 6.07903e+47
n + lu232 -> lu233 + gamma: upper bound is 1e+45 and reverse rate at t9 = 0.01 is 1.07371e+46
n + mn84 -> mn85 + gamma: upper bound is 1e+45 and reverse rate at t9 = 0.158489 is 3.29673e+46
n + se114 -> se115 + gamma: upper bound is 1e+45 and reverse rate at t9 = 0.158489 is 6.58889e+46
n + ce189 -> ce190 + gamma: upper bound is 1e+45 and reverse rate at t9 = 0.00251189 is 1.5158e+49
But,with 1.e50,it is no more exsit.

[Bradley S. Meyer]

I don't see anything horribly out of the ordinary there. Let's try zeroing in on where things go awry. Try:

cd nucnet-tools-code/examples/analysis

make print_zone_abundances

./print_zone_abundances ../../my_examples/network/my_output.xml "" > out.txt

Add out.txt as an attachment to your reply.

[Anonymous]

I have ran these code and got out.txt,but I realized that my last caculation didin't totally finish because of time.Is the out.txt still helpful?Or I need to rerun my caculation.
Thank you for your reply and patience.Best wishes!

[Bradley S. Meyer]

Go ahead and run your calculation first. In the meantime, could you attach your zone.xml file, or is it the same as in the blog post?

[Anonymous]

Dear Brad:
I have ran the whole caculation and got the out.txt(in attachments),I checked zone.xml file else,it is the same as in the blog post.

[Bradley S. Meyer]

I have another idea. You could try turning off the network limiter below, say T9 = 0.01. You would replace the lines

 user::limit_evolution_network( zone );

 i_step++;

in run_single_zone.cpp with

if( zone.getProperty<double>( nnt::s_T9 ) > 0.01 )
{
    user::limit_evolution_network( zone );
}

i_step++;

Now recompile and rerun. The calculation will take longer, but it may allow the network to adjust to large rates at low temperature. If that works, there are some diagnostics you can do to find out what was happening in the previous calculation.

[Bradley S. Meyer]

Your calculation is good down to time dump 188 (t9 = 0.00487). After that, the abundances blow up. That's probably because some rate goes bad. If you are willing to share my_net_WS.xml, I'd be happy to take a look at this to see if I can figure out what is going on. The way to do this would be to do

gzip my_net_WS.xml

and attach my_net_WS.xml.gz to your reply to this post. After I've downloaded my_net_WS.xml from your reply, I will delete it before I approve your post so the file won't be publicly available.

[Anonymous]

Dear Brad:
It is OK to share my_net_WS33.xml with you and I add it as a attachment.I will also try the way that you have metioned above.
Thank you for your help!Best wishes!

[Bradley S. Meyer]

I've gotten your my_net_WS33.xml (which I've deleted from above for privacy). The changes in the mass excess changed the sign of a number of reaction Q values. To see this, I wrote and committed a new example code. I first made sure that my_net.xml and my_net_WS33.xml are in the directory nucnet-tools-code/data_pub. I then typed

cd nucnet-tools-code/examples/misc

svn up

make Qvalue_sign_change

I then ran the code by typing

./Qvalue_sign_change ../../data_pub/my_net.xml ../../data_pub/my_net_WS33.xml ../../data_pub/my_net.xml

The reactions in my_net.xml and my_net_WS33.xml are the same. The output from the code shows the reactions that change sign of the Q value when the nuclear data changes from that in my_net.xml to that in my_net_WS33.xml. An example reaction that changes sign of the Q value is n + cu98 -> cu99 + gamma:

 n + cu98 -> cu99 + gamma                                    0.521  -0.506

The Q value for this reaction for the nuclear data in my_net.xml is 0.521 MeV while that for the nuclear data in my_net_WS33.xml is -0.506 MeV. Since the reverse reaction rate is computed from detailed balance, and the appropriate expression includes the term exp(-Q/kT), a negative Q value can cause a reverse rate to blow up at low temperature. Which reaction is causing the problem remains to be seen. I'll be looking at that later today.

You can use options with the Qvalue_sign_change.cpp code. To get a usage statement, type

./Qvalue_sign_change

or

./Qvalue_sign_change --help

To list all reactions except weak decay, you can type

./Qvalue_sign_change ../../data_pub/my_net.xml ../../data_pub/my_net_WS33.xml ../../data_pub/my_net.xml --reac_xpath "[not(product[contains(.,'neutrino')])]"

[Anonymous]

Dear Brad:
After I changed run_single_zone.cpp as you written above and rerun my caculation,I got the ordinary data,finally.How could I do to find out what is the problem with my caculation.Sorry for too many questions and thank you for your help!Best wishes!

[Bradley S. Meyer]

Thanks for your post. The issue is with those reactions that, due to the new data, significantly changed their Q values.

The first thing you can try is to go back to

user::limit_evolution_network( zone );

i_step++;

in run_single_zone.cpp. Then look for

  //============================================================================
// Evolve abundances.
//============================================================================

 user::evolve( zone );

and change it to

//============================================================================
// Evolve abundances.
//============================================================================

user::safe_evolve( zone );

Now recompile and rerun with your my_net_WS33.xml. You might need to update your nucnet-tools-code by doing

cd nucnet-tools-code

svn up

This will check the abundances in each time step and then cut the time step and recompute if the abundances change too much (because of the changed Q value). That should work.

The second thing you can do is probably even better. That is to go to the NucNet Project jina_to_webnucleo. Follow the steps to create the reaction_data.xml file. You can get that file as described in the wiki, or you can download from the JINA Reaclib Database snapshots page (you can include Chaps. 9, 10, and 11). If you do the latter, be sure to choose the Reaclib2 format. Run the reaction_converter code on the data you downloaded. That resulting file reaction_data.xml has forward and reverse reactions. Now move that to the appropriate data_pub directory (where you have your my_net_WS33.xml file), and create a new xml file my_net_combined.xml:

<?xml version="1.0" encoding="ISO-8859-1"?>

<nuclear_network
  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
  xmlns:xi="http://www.w3.org/2001/XInclude"
>

   <nuclear_data>
     <xi:include href="./my_net_WS33.xml" xpointer="xpointer(//nuclide)" />
   </nuclear_data>

   <reaction_data>
     <xi:include href="./reaction_data.xml" xpointer="xpointer(//reaction)" />
   </reaction_data>

</nuclear_network>

You can run your calculation (with user::evolve() or user::safe_evolve()) with this file as the net input. The calculation will select the nuclear data from my_net_WS33.xml and reaction data from reaction_data.xml. For example, in my_examples/network, you can type:

./run_single_zone ../../data_pub/my_net_combined.xml ../../data_pub/zone_rproc.xml my_output2.xml "[z <= 90]"

That calculation will begin by looking through all reactions are removing the endothermic version. It then won't have a problem at low temperatures. Note that my standard my_net.xml has the endothermic reactions already removed--any reverse reaction rates are computed from detailed balance.

Of course all this hits on the important point that it is ultimately not realistic to change nuclear masses without changing the reaction rates. Changing the relative masses of reactants or products in a reaction certainly would change the reaction rate among them. Nevertheless, in exploring the importance of nuclear masses for, say, the r process, one often wants to change masses but leave rates unchanged. This second approach (allowing the code to pick out the exothermic reactions and compute the reverse rates from detailed balance) seems to me the most reasonable way to do that. Let me know if you have any trouble with this. Best wishes.

[Anonymous]

Dear professor,

I got some problems that maybe have the same reason with the discussion above. When I do the calculation with the mass excesses from WS4 and the (n,gamma) reaction rate calculate from talys using the WS4 mass table, I get unrealistically large abundance and even negative abundance.

I follow this discussion, type
./checkrates ../datapub/my_net_WS4.xml 1.e-5 10 0. 1.e30

and get (I select several lines)

...
Removing te163 + gamma -> n + te162 (Data source: rath)
Removing pr190 + gamma -> h1 + ce189 (Data source: rath)
Removing ir169 + gamma -> h1 + os168 (Data source: rath)
Removing ne32 + gamma -> h1 + f31 (Data source: rpsm)
...
n + nb123 -> nb124 + gamma: upper bound is 1e+30 and rate at t9 = 0.0398107 is 8.59003e+212
n + nb123 -> nb124 + gamma: upper bound is 1e+30 and reverse rate at t9 = 0.0398107 is 2.90451e+225
n + pt240 -> pt241 + gamma: upper bound is 1e+30 and rate at t9 = 0.0398107 is 1.15422e+201
n + pt240 -> pt241 + gamma: upper bound is 1e+30 and reverse rate at t9 = 0.0398107 is 3.49003e+215
n + rh144 -> rh145 + gamma: upper bound is 1e+30 and rate at t9 = 0.0398107 is 5.81618e+212
n + rh144 -> rh145 + gamma: upper bound is 1e+30 and reverse rate at t9 = 0.0398107 is 3.53392e+225
n + ar53 -> ar54 + gamma: upper bound is 1e+30 and reverse rate at t9 = 0.158489 is 2.69691e+36
n + ho211 -> ho212 + gamma: upper bound is 1e+30 and rate at t9 = 0.0398107 is 1.34631e+216
n + ho211 -> ho212 + gamma: upper bound is 1e+30 and reverse rate at t9 = 0.0398107 is 2.58738e+228
...
Reaction: n + hg263 -> hg264 + gamma
Out of bounds for t9 = 2.511887e-03
ERROR: Invalid rate in ../vendor/libnucnet/0.32/src/LibnucnetReac.c on line 2018

It seems that some rate are blow up. I focus on n + nb123 -> nb124 + gamma, both the forward and reverse rate of which are blow up. I type
./printratesbystring ../datapub/mynetWS4.xml "n + nb123 -> nb124 + gamma"
and get
n + nb123 -> nb124 + gamma
T9(K) Rate
===== ============
0.1 1.458660e-07
0.2 1.781160e-05
0.3 2.040530e-04
0.4 9.775290e-04
0.5 2.735710e-03
0.6 5.590050e-03
0.7 9.541830e-03
0.8 1.470660e-02
0.9 2.132680e-02
1.0 2.965150e-02
1.5 9.727840e-02
2.0 1.867370e-01
2.5 2.678350e-01
3.0 3.270910e-01
3.5 3.616630e-01
4.0 3.719020e-01
4.5 3.585852e-01
5.0 3.222640e-01
6.0 2.053640e-01
7.0 9.660460e-02
8.0 3.777110e-02
9.0 1.375940e-02
10.0 4.933990e-03
It seems normal in my point of view, I also print out the rate in the default date in the reaclib, I get
n + nb123 -> nb124 + gamma
T9(K) Rate
===== ============
0.1 1.713904e-01
0.2 2.012593e-01
0.3 3.101040e-01
0.4 4.762934e-01
0.5 7.020363e-01
0.6 9.910647e-01
0.7 1.346407e+00
0.8 1.770013e+00
0.9 2.262741e+00
1.0 2.824434e+00
1.5 6.603994e+00
2.0 1.167789e+01
2.5 1.744627e+01
3.0 2.325660e+01
3.5 2.851864e+01
4.0 3.277220e+01
4.5 3.571959e+01
5.0 3.723069e+01
6.0 3.615252e+01
7.0 3.094792e+01
8.0 2.378403e+01
9.0 1.660176e+01
10.0 1.060907e+01
Which seems even larger than the rate from WS4. So I don't know why the rate in WS4 would blow up.

And I don't know how the rate 8.59003e+212 of n + nb123 -> nb124 + gamma: upper bound is 1e+30 and rate at t9 = 0.0398107 is 8.59003e+212 obtained from the reaction rate above.

[Bradley S. Meyer]

Thanks for your post. It looks like the rates blow up for t9 < 0.1. Could you try editing the nucnet-tools-code/vendor/libnucnet/0.32/examples/print_rates_by_string.c code to change the t9's to compute for t9 < 0.1? For example, you might try:

  double a_t9[23] = { 0.01, 0.02, 0.03, 0.04 ...

Just make sure you have 23 entries, or modify the code to use a different number from 23 (in two places). Then go back to nucnet-tools-code/bulid and remake the libnucnet examples and rerun the print_rates_by_string on your file. It would be interesting to see what the rates do below t9 = 0.1. Let me know how it goes. Best wishes.

[Anonymous]

Dear professor,

If I just type
make print_rates_by_string
it seems fine to make it without download the files from sourceforge, and in the ....../0.32/examples/ (Not in the /nucnet-tools-code/libnucnet/) I get theprint_rates_by_string I run it by type
./print_rates_by_string ../../../../data_pub/my_net.xml "n + nb123 -> nb124 + gamma"
and get
n + nb123 -> nb124 + gamma
T9(K) Rate
===== ============
0.0 9.807410e-30
Reaction: n + nb123 -> nb124 + gamma
Out of bounds for t9 = 2.000000e-02
ERROR: Invalid rate in ../src/LibnucnetReac.c on line 2018

So you are right that the rates blows up for t9<0.1 . Does this due to the function fitting of the rate? My input rate.txt for n + nb123 -> nb124 + gamma is

rate_table
ws4_ng_table
2
n
nb123
2
nb124
gamma
30
0.0001 0.00000E+00 1.00
0.0005 0.00000E+00 1.00
0.0010 0.00000E+00 1.00
0.0050 0.00000E+00 1.00
0.0100 9.80741E-30 1.00
0.0500 9.67296E-11 1.00
0.1000 1.45866E-07 1.00
0.1500 2.91057E-06 1.00
0.2000 1.78116E-05 1.00
0.2500 6.91165E-05 1.00
0.3000 2.04053E-04 1.00
0.4000 9.77529E-04 1.00
0.5000 2.73571E-03 1.00
0.6000 5.59005E-03 1.00
0.7000 9.54183E-03 1.00
0.8000 1.47066E-02 1.00
0.9000 2.13268E-02 1.00
1.0000 2.96515E-02 1.00
1.5000 9.72784E-02 1.00
2.0000 1.86737E-01 1.00
2.5000 2.67835E-01 1.00
3.0000 3.27091E-01 1.00
3.5000 3.61663E-01 1.00
4.0000 3.71902E-01 1.00
5.0000 3.22264E-01 1.00
6.0000 2.05364E-01 1.00
7.0000 9.66046E-02 1.00
8.0000 3.77711E-02 1.00
9.0000 1.37594E-02 1.00
10.0000 4.93399E-03 1.00

In this files, each rate seems right. How could I avoid this problem?

Thank you very much.

[Bradley S. Meyer]

Thanks for this post. The problem is that the spline interpolation routine has a problem for the low temperatures and zero rate. The following worked for me. Change your data not to include the zero rates:

rate_table
ws4_ng_table
2
n
nb123
2
nb124
gamma
26
0.0100 9.80741E-30 1.00
0.0500 9.67296E-11 1.00
0.1000 1.45866E-07 1.00
0.1500 2.91057E-06 1.00
0.2000 1.78116E-05 1.00
0.2500 6.91165E-05 1.00
0.3000 2.04053E-04 1.00
0.4000 9.77529E-04 1.00
0.5000 2.73571E-03 1.00
0.6000 5.59005E-03 1.00
0.7000 9.54183E-03 1.00
0.8000 1.47066E-02 1.00
0.9000 2.13268E-02 1.00
1.0000 2.96515E-02 1.00
1.5000 9.72784E-02 1.00
2.0000 1.86737E-01 1.00
2.5000 2.67835E-01 1.00
3.0000 3.27091E-01 1.00
3.5000 3.61663E-01 1.00
4.0000 3.71902E-01 1.00
5.0000 3.22264E-01 1.00
6.0000 2.05364E-01 1.00
7.0000 9.66046E-02 1.00
8.0000 3.77711E-02 1.00
9.0000 1.37594E-02 1.00
10.0000 4.93399E-03 1.00

Then it works. This is, of course, undesirable behavior on the part of libnucnet, but it is due to the fact that the code interpolates the log10 of the rate. I will see if I can fix that in the next release, but for now, perhaps you can just make your text files not include zero rate. Sorry for the trouble. Best wishes.

[Anonymous]

Thank you very much.

Yesterday, I try it with more temperature points in the range 0.01 and 0.05, and it seems fine in the fitting.
With your reply, I think I now understand why this happens.

[Bradley S. Meyer]

Good. I'll still try to fix the zero rates issue. Best wishes.

[Anonymous]

Dear professor,

When I type
make -f Makefile.libnucnet alllibnucnet
After I modify theprint_rates_by_string.c.
It seems to download the print_rates_by_string.c from sourceforge.net and cover the one I modified.

How to make it without download the new files from sourceforge?

[Bradley S. Meyer]

Sorry about that. As you discovered in your next post, remaking only the specifically changed code works.