Skip to content

Commit

Permalink
Add test for C++ network with derived rates (pynucastro#572)
Browse files Browse the repository at this point in the history
  • Loading branch information
yut23 authored Jul 17, 2023
1 parent a352912 commit 4222442
Show file tree
Hide file tree
Showing 14 changed files with 1,588 additions and 0 deletions.
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
CEXE_headers += network_properties.H

ifeq ($(USE_REACT),TRUE)
CEXE_sources += actual_network_data.cpp
CEXE_headers += actual_network.H
CEXE_headers += tfactors.H
CEXE_headers += partition_functions.H
CEXE_headers += actual_rhs.H
CEXE_headers += reaclib_rates.H
CEXE_headers += table_rates.H
CEXE_sources += table_rates_data.cpp
USE_SCREENING = TRUE
USE_NEUTRINOS = TRUE
endif
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
NSCREEN := 2
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
@namespace: network

Original file line number Diff line number Diff line change
@@ -0,0 +1,76 @@
#ifndef actual_network_H
#define actual_network_H

#include <AMReX_REAL.H>
#include <AMReX_Array.H>

#include <fundamental_constants.H>
#include <network_properties.H>

using namespace amrex;

void actual_network_init();

const std::string network_name = "pynucastro-cxx";

namespace network
{
extern AMREX_GPU_MANAGED amrex::Array1D<amrex::Real, 1, NumSpec> bion;
extern AMREX_GPU_MANAGED amrex::Array1D<amrex::Real, 1, NumSpec> mion;
}

namespace Rates
{

enum NetworkRates
{
k_he4_fe52_to_ni56 = 1,
k_p_co55_to_ni56 = 2,
k_he4_fe52_to_p_co55 = 3,
k_ni56_to_he4_fe52_derived = 4,
k_ni56_to_p_co55_derived = 5,
k_p_co55_to_he4_fe52_derived = 6,
NumRates = k_p_co55_to_he4_fe52_derived
};

// number of reaclib rates

const int NrateReaclib = 6;

// number of tabular rates

const int NrateTabular = 0;

// rate names -- note: the rates are 1-based, not zero-based, so we pad
// this vector with rate_names[0] = "" so the indices line up with the
// NetworkRates enum

static const std::vector<std::string> rate_names = {
"", // 0
"he4_fe52_to_ni56", // 1,
"p_co55_to_ni56", // 2,
"he4_fe52_to_p_co55", // 3,
"ni56_to_he4_fe52_derived", // 4,
"ni56_to_p_co55_derived", // 5,
"p_co55_to_he4_fe52_derived" // 6,
};

}

#ifdef NSE_NET
namespace NSE_INDEX
{
constexpr int h1_index = 0;
constexpr int n_index = -1;
constexpr int he4_index = 1;

// Each row corresponds to the rate in NetworkRates enum
// First 3 row indices for reactants, followed by 3 product indices
// last index is the corresponding reverse rate index.

extern AMREX_GPU_MANAGED amrex::Array2D<int, 1, Rates::NumRates, 1, 7, Order::C> rate_indices;
extern AMREX_GPU_MANAGED bool initialized;
}
#endif

#endif
Original file line number Diff line number Diff line change
@@ -0,0 +1,49 @@
#include <actual_network.H>


namespace network
{
AMREX_GPU_MANAGED amrex::Array1D<amrex::Real, 1, NumSpec> bion;
AMREX_GPU_MANAGED amrex::Array1D<amrex::Real, 1, NumSpec> mion;
}

#ifdef NSE_NET
namespace NSE_INDEX
{
AMREX_GPU_MANAGED amrex::Array2D<int, 1, Rates::NumRates, 1, 7, Order::C> rate_indices {
-1, 1, 2, -1, -1, 4, 4,
-1, 0, 3, -1, -1, 4, 5,
-1, 1, 2, -1, 0, 3, 6,
-1, -1, 4, -1, 1, 2, -1,
-1, -1, 4, -1, 0, 3, -1,
-1, 0, 3, -1, 1, 2, -1
};
AMREX_GPU_MANAGED bool initialized = false;
}
#endif

void actual_network_init()
{
using namespace Species;
using namespace network;

// binding energies per nucleon in MeV
amrex::Array1D<amrex::Real, 1, NumSpec> ebind_per_nucleon;

ebind_per_nucleon(H1) = 0.0_rt;
ebind_per_nucleon(He4) = 7.073915_rt;
ebind_per_nucleon(Fe52) = 8.609574_rt;
ebind_per_nucleon(Co55) = 8.669618_rt;
ebind_per_nucleon(Ni56) = 8.642779_rt;

// convert to binding energies per nucleus in MeV
for (int i = 1; i <= NumSpec; ++i) {
bion(i) = ebind_per_nucleon(i) * aion[i-1];
}

// Set the mass -- this will be in grams
for (int i = 1; i <= NumSpec; ++i) {
mion(i) = (aion[i-1] - zion[i-1]) * C::Legacy::m_n + zion[i-1] * (C::Legacy::m_p + C::Legacy::m_e) - bion(i) * C::Legacy::MeV2gr;
}

}
Loading

0 comments on commit 4222442

Please sign in to comment.