In a previous post, I showed how one can modify reaction rates during a network calculation. This was done with the default rate modifier, which uses views to select out reactions and then multiply them by a constant factor. In this post, I will show how you can modify rates with your own rate modification routine.
I will use the constant entropy calculation to demonstrate how to modify rates. Those calculations showed how one could include screening only or screening and NSE corrections. They did not, however, show how to include NSE corrections only (without screening). This is because the default routines include application of NSE corrections within the application of screening. This makes sense because Coulomb corrections should affect both the rate of interaction of charged particles and the chemical potential of species (and hence their equilibrium abundance). Nevertheless, one might want to include the change in the chemical potential without changing the rate of the forward reactions.
I could include NSE corrections without screening the reaction rates by changing my screening function. In this case, I would supply that does not change the forward and reverse reaction rates but then applies a reverse-ratio correction to the reverse reaction. I will instead apply a rate modification function that modifies reaction rates for the network calculation after they have been computed but before the network equations have been solved.
I do this by downloading the file 2016-08-11.tar.gz and placing in in my nucnet-tools-code directory. In that directory, I then type
tar xvf 2016-08-11.tar
or, more simply,
tar zxvf 2016-08-11.tar
This creates and places codes in the directories my_examples and my_users. I now change into the directory my_user
and look at the new files. my_rate_modifies.cpp/h are the code and header for the routine for the rate modification. In this case, the rate modifier is the function modify_rates, which loops over all reactions in the evolution view of the network calculation and applies the correct reverse ratio correction factor.
I next change into my_examples/network by typing
and compute the new run_entropy.cpp to the default version by typing:
diff -Naur ../../examples/run_entropy.cpp run_entropy.cpp
The lines beginning with a - are only in the old version of the code (examples/network/run_entropy.cpp) while those beginning with a + are only in the new version. As in my previous post, I am including the Coulomb correction to the entropy during the calculation. In this particular calculation, I am also setting the rate modification function if the evolution zone has the property use nse_correction set to yes.
The essential prototype of the rate modification function is
void rate_modification_function( void );
In my example, my modification function is my_user::modify_rates, which is in the my_users namespace. I bind to the function a reference to the evolution zone. These are the lines
zone.updateFunction( nnt::s_RATE_MODIFICATION_FUNCTION, static_cast<boost::function<void()> >( boost::bind( my_user::modify_rates, boost::ref( zone ) ) ) );
I now compile the code by typing
Once that is done, I edit the ../../data_pub/zone.xml file to read
<zone_data> <zone> <optional_properties> <property name="tend">1.</property> <property name="tau">inf</property> <property name="munuekT">-inf</property> <property name="t9_0">2.</property> <property name="steps">5</property> <property name="rho_0">3.e9</property> <property name="x" tag1="0">1</property> <property name="x" tag1="1">0</property> <property name="use nse correction">yes</property> </optional_properties> <mass_fractions> <nuclide name="he4"> <z>2</z> <a>4</a> <x>1</x> </nuclide> </mass_fractions> </zone> </zone_data>
I now run the constant entropy calculation by typing
./run_entropy ../../data_pub/my_net.xml ../../data_pub/zone.xml my_output_3e9_nse_correction_only.xml --nuc_xpath "[z <= 60]" --reac_xpath "[not(product[contains(.,'neutrino')])]" --sdot_reac_xpath none
When this is done, I get the time and temperature by typing
../../examples/analysis/print_properties my_output_3e9_nse_correction_only.xml time t9 > props_nse_correction_only.txt
I plot the temperature T9 versus time (column 3 vs. column 2 of props_nse_correction_only.txt) along with the same quantities from the other calculations to get this figure:
In this figure, my new calculation is shown as the green curve. I see that the temperature rise is later than in the calculation that had NSE corrections and screening. This is expected since without screening, the reaction rates are not as fast as with screening. On the other hand, this new calculation goes to the same final temperature since it attains the same final abundance distribution.
The steps presented here are those one could use for any rate modification (including screening, NSE corrections, temperature- or time-dependent modifications, etc.). One would create a .cpp and .h file with any name that includes a routine to loop over the evolution reactions in a zone and apply a modification. One would then include the file in a Makefile.inc. Next one would update the zone rate modification function in the main program, include the appropriate header, and include the file in the make process. Extra data to the rate modifier can be bound to the function when it is registered with the zone.