Overriding the inverter in a constant entropy network calculation

In my last post, I showed how to run a network calculation that begins in NSE (nuclear statistical equilibrium) and evolve it at constant entropy. The calculation worked by finding the temperature at each time step consistent with the input entropy. In this calculation, at each time step, a root-finding routine varies the temperature. The abundances are evolved over the time step with the current temperature guess. The evolved abundances are then used to compute the entropy. The temperature guess is then varied until the computed entropy agrees with the input entropy. In this way, the network is called for each root-finding iteration, which means it is called many times for each time step. This makes the calculation rather slow.

In this post, I will override the default entropy inverter. I will compute the temperature from the entropy using the fixed abundances from the previous time step. I will then use the temperature to evolve the abundances over the current time step. The abundances will therefore not be consistent with the temperature as determined from the entropy, but, as I will show, when this treatment works, it still gives accurate results.

I first run a calculation with the "default" inverter in NucNet Tools. I follow the instructions in my previous post except I choose for this post an entropy per nucleon of s/k = 15. I change into my examples/network directory by typing

cd nucnet-tools-code/examples/network

I will run an exponential expansion of the density, so here I do not need to set the NNT_HYDRO_CODE environment variable. I compile the code by typing

make clean

make run_constant_entropy

Once the code successfully compiles, I edit the file ../../data_pub/zone.xml until it reads

 <?xml version="1.0" encoding="UTF-8"?>

   <zone_data>

       <zone>

         <optional_properties>

           <property name="entropy per nucleon">15.</property>
           <property name="t9">10.</property>
          <property name="Ye">0.5</property>
          <property name="tend">100.</property>
          <property name="steps">5</property>
          <property name="time">0.</property>
          <property name="tau">0.1</property>
          <property name="dt">1.e-5</property>

       </optional_properties>

       <mass_fractions/>

     </zone>

   </zone_data>

I run the code by typing

./run_constant_entropy ../../data_pub/my_net.xml ../../data_pub/zone.xml out_15.xml "[z <= 60]"

The calculation takes several hours on my computer.

I next override the inverter (the routine to compute the temperature from the entropy). As before, I work in the my_user and my_examples directories I previously created since I am modifying the NucNet Tools codes. I first change into my_users by typing

cd ../../my_users

I next ensure that I have the my_hydro.cpp and my_hydro.h from my previous trajectory post and
my_thermo.cpp and my_thermo.h files from my previous function post. I now need the code to override the default calculation of the network temperature from the entropy. I create (or download) the file my_evolve.cpp to my_user. I also create (or download) the corresponding header file my_evolve.h to my_users. Finally, I update the include file in my_users until it looks like Makefile.inc.

I next change into my_examples/network by typing

cd ../my_examples/network

I create (or download) the main code run_constant_entropy.cpp and the updated Makefile in my_examples/network. The new version of run_constant_entropy.cpp calls the inverter function my_user::my_network_t9_from_entropy in my_users/my_evolve.cpp. The new inverter function, which overrides the user::network_t9_from_entropy we used before, does not evolve the network during each rooting finding iteration. One might naturally wonder why I put the routine in my_evolve.cpp if I'm not evolving the network. I chose to do this simply to mirror the default structure since user::network_t9_from_entropy is in users/evolve.cpp. I could just as well have put my_user::network_t9_from_entropy in my_thermo.cpp and added the prototype to my_thermo.h.

Since I will use an exponential expansion trajectory, I compile the executable by typing (for the bash shell)

export NNT_HYDRO_CODE=exponential_expansion

make run_constant_entropy

Note that setting the NNT_HYDRO_CODE environment variable is not optional in this case. If the variable is not set, the hydro routine used will be my user-defined trajectory function.

Once the code has successfully compiled, I can execute it. I type

./run_constant_entropy ../../data_pub/my_net.xml ../../data_pub/zone.xml out_15.xml "[z <= 60]"

This calculation takes less than 30 minutes on my computer, so there is a considerable speed up in the calculation if we do not evolve the network during each root finding iteration.

I can plot the entropy per nucleon in each of the components as I did in the previous post. For the overridden inverter calculation in my_examples/network this is the resulting figure:

ent_full

The total entropy per nucleon is nearly constant throughout the calculation; however, there is some deviation at high temperature. If I plot the total entropy per nucleon versus T9 on a more restricted ordinate scale, here is the resulting figure:

ent

It is clear at high temperature there is some deviation away from a constant entropy per nucleon of s/k = 15. When we do not evolve the abundances during a root finding iteration, the temperature found will not be quite consistent with the value that would be found from the evolved abundances (the value calculated with the compute_thermo_quantity code). Near T9 = 10, the entropy per nucleon in baryons is decreasing as the material expands. This means that lighter nuclei are assembling into heavier ones, thereby decreasing the number of nuclei per nucleon. When we do not allow the abundances to evolve during the root-finding iterations, the number of nuclei per nucleon will be a little larger than if we do evolve the abundances. This means that the entropy in baryons in the former case will be a little larger than in the latter case for the same temperature. Thus, during this high-temperature regime, by not evolving the abundances during the root-finding iterations, the code computes a temperature slightly lower than it should be. Then when the compute_thermo_quantity code calculates the entropy per nucleon from the too low temperature and the evolved abundances, the total entropy per nucleon will be less than the input entropy.

These considerations become evident if I plot the temperature versus time. I get the relevant data for the exact calculation (the one that evolves the abundances during the root-finding iterations) by typing

cd ../../examples/analysis

./examples_make

cd ../../my_examples/network

../../examples/analysis/print_properties ../../examples/network/out_15.xml time t9 > props_exact_15.txt

I then do the same for the approximate (approx) calculation (the one that does not evolve the abundances during the root-finding iterations) by typing

../../examples/analysis/print_properties out_15.xml time t9 > props_approx_15.txt

If I plot column 3 versus column 2 in props_exact_15.txt and props_approx_15.txt, I get this figure:

t9

Here I can see that the temperature in the exact calculation lies slightly above that in approx calculation, at least for a short period of time. Nevertheless, the approximate calculation temperature agrees extremely well over most of the calculation, and where the deviation ends, the temperature is still sufficiently high that the abundances are still in NSE, which means any difference is temperature history is forgotten. We thus expect the final abundances to be essentially the same for the two calculations. I can confirm this by plotting the final elemental abundances. I type

../../examples/analysis/print_abundances_vs_nucleon_number ../../examples/network/out_15.xml z "[last()]" > yz_exact_15.txt

../../examples/analysis/print_abundances_vs_nucleon_number out_15.xml z "[last()]" > yz_approx_15.txt

If I plot column 2 versus column 1 for yz_exact_15.txt and yz_approx_15.txt, I get this figure:

yz

The two curves lie on top of each other; thus, the final abundances are the same. We clearly make little error in neglecting the evolution of the abundances during the iterations to find the temperature from the entropy at each time step.

For higher entropy per nucleon than s/k = 15, the approximate calculation is an even better match to the exact one since the entropy in baryons is an even smaller component of the total. When the code that neglects the evolution of the abundances in the iterations to find the temperature works, it gives a good approximation to the exact result.

For lower entropy, on the other hand, the approximate calculation may not work. For example, if I try using the run_constant_entropy in my_examples/network with an entropy per nucleon s/k = 10, the time step frequently crashes and prevents the code from evolving in a stable fashion. Since baryons make up a dominant fraction of the total entropy, the temperature computed by the root finder can be considerably different from the temperature in the previous time step. When the network evolves the abundances with the new temperature, certain species change their abundance dramatically. In an effort to reduce the large changes in the next time step, the code decreases the time step value considerably. With the time step value experiencing large variations, the network evolution is not stable.

When we use the exact entropy inverter, the network evolution is stable, even at low entropy, because the temperature determined by the inverter is consistent with the evolved abundances. There are possible strategies to alleviate these difficulties while using the approximate entropy inverter (e.g., using equilibrium abundances), but one might simply choose to use the "exact" entropy inverter when evolving abundances at constant low entropy.

Posted by Bradley S. Meyer 2016-02-26


Anonymous

Cancel  Add attachments





Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:





No, thanks