Skip to content

Commit

Permalink
Merge branch 'development'
Browse files Browse the repository at this point in the history
  • Loading branch information
@zingale
zingale committed May 1, 2020
2 parents 6acbe6f + 55ce324 commit bbfabaa
Show file tree
Hide file tree
Showing 156 changed files with 32,537 additions and 15,087 deletions.
3 changes: 3 additions & 0 deletions CHANGES.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,9 @@

* species names are available as an enum in network_properties.H (#304)

* The screening on O16+O16 in iso7 was fixed (#302)

* The VODE integrator is now available in C++ (#299)

# 20.04

Expand Down
2 changes: 1 addition & 1 deletion EOS/breakout/actual_eos.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ AMREX_GPU_HOST_DEVICE inline
void actual_eos (eos_input_t input, eos_t& state)
{

const Real R = k_B * n_A;
const Real R = C::k_B * C::n_A;

Real poverrho;

Expand Down
2 changes: 1 addition & 1 deletion EOS/gamma_law/actual_eos.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,7 @@ AMREX_GPU_HOST_DEVICE inline
void actual_eos (eos_input_t input, eos_t& state)
{

const Real R = k_B * n_A;
const Real R = C::k_B * C::n_A;

// Calculate mu.

Expand Down
30 changes: 15 additions & 15 deletions EOS/gamma_law_general/actual_eos.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,8 @@ AMREX_GPU_HOST_DEVICE inline
void actual_eos(const eos_input_t input, eos_t& state) {

// Get the mass of a nucleon from Avogadro's number.
const Real m_nucleon = 1.0 / n_A;
const Real fac = 1.0 / std::pow(2.0 * M_PI * hbar * hbar, 1.5);
const Real m_nucleon = 1.0 / C::n_A;
const Real fac = 1.0 / std::pow(2.0 * M_PI * C::hbar * C::hbar, 1.5);

if (eos_assume_neutral) {
state.mu = state.abar;
Expand Down Expand Up @@ -79,7 +79,7 @@ void actual_eos(const eos_input_t input, eos_t& state) {
// Solve for the temperature:
// h = e + p/rho = (p/rho)*[1 + 1/(gamma-1)] = (p/rho)*gamma/(gamma-1)

state.T = (state.h * state.mu * m_nucleon / k_B)*(eos_gamma - 1.0)/eos_gamma;
state.T = (state.h * state.mu * m_nucleon / C::k_B)*(eos_gamma - 1.0)/eos_gamma;

break;

Expand All @@ -90,7 +90,7 @@ void actual_eos(const eos_input_t input, eos_t& state) {
// Solve for the density:
// p = rho k T / (mu m_nucleon)

state.rho = state.p * state.mu * m_nucleon / (k_B * state.T);
state.rho = state.p * state.mu * m_nucleon / (C::k_B * state.T);

break;

Expand All @@ -101,7 +101,7 @@ void actual_eos(const eos_input_t input, eos_t& state) {
// Solve for the temperature:
// p = rho k T / (mu m_nucleon)

state.T = state.p * state.mu * m_nucleon / (k_B * state.rho);
state.T = state.p * state.mu * m_nucleon / (C::k_B * state.rho);

break;

Expand All @@ -112,7 +112,7 @@ void actual_eos(const eos_input_t input, eos_t& state) {
// Solve for the temperature
// e = k T / [(mu m_nucleon)*(gamma-1)]

state.T = state.e * state.mu * m_nucleon * (eos_gamma - 1.0) / k_B;
state.T = state.e * state.mu * m_nucleon * (eos_gamma - 1.0) / C::k_B;

break;

Expand All @@ -125,13 +125,13 @@ void actual_eos(const eos_input_t input, eos_t& state) {
// rho = p mu m_nucleon / (k T)

state.T = std::pow(state.p, 2.0/5.0) *
std::pow(2.0*M_PI*hbar*hbar/(state.mu*m_nucleon), 3.0/5.0) *
exp(2.0*state.mu*m_nucleon*state.s/(5.0*k_B) - 1.0) / k_B;
std::pow(2.0*M_PI*C::hbar*C::hbar/(state.mu*m_nucleon), 3.0/5.0) *
exp(2.0*state.mu*m_nucleon*state.s/(5.0*C::k_B) - 1.0) / C::k_B;

// Solve for the density
// rho = p mu m_nucleon / (k T)

state.rho = state.p * state.mu * m_nucleon / (k_B * state.T);
state.rho = state.p * state.mu * m_nucleon / (C::k_B * state.T);

break;

Expand All @@ -142,7 +142,7 @@ void actual_eos(const eos_input_t input, eos_t& state) {
// Solve for temperature and density

state.rho = state.p / state.h * eos_gamma / (eos_gamma - 1.0);
state.T = state.p * state.mu * m_nucleon / (k_B * state.rho);
state.T = state.p * state.mu * m_nucleon / (C::k_B * state.rho);

break;

Expand Down Expand Up @@ -176,25 +176,25 @@ void actual_eos(const eos_input_t input, eos_t& state) {

// Compute the pressure simply from the ideal gas law, and the
// specific internal energy using the gamma-law EOS relation.
state.p = state.rho * state.T * (k_B / (state.mu*m_nucleon));
state.p = state.rho * state.T * (C::k_B / (state.mu*m_nucleon));
state.e = state.p/(eos_gamma - 1.0) * rhoinv;

// enthalpy is h = e + p/rho
state.h = state.e + state.p * rhoinv;

// entropy (per gram) of an ideal monoatomic gas (the Sackur-Tetrode equation)
// NOTE: this expression is only valid for gamma = 5/3.
state.s = (k_B/(state.mu*m_nucleon))*
state.s = (C::k_B/(state.mu*m_nucleon))*
(2.5 + log((std::pow(state.mu*m_nucleon, 2.5) * rhoinv ) *
std::pow(k_B * state.T, 1.5) * fac ) );
std::pow(C::k_B * state.T, 1.5) * fac ) );

// Compute the thermodynamic derivatives and specific heats
state.dpdT = state.p * Tinv;
state.dpdr = state.p * rhoinv;
state.dedT = state.e * Tinv;
state.dedr = 0.0;
state.dsdT = 1.5 * (k_B / (state.mu * m_nucleon)) * Tinv;
state.dsdr = - (k_B / (state.mu * m_nucleon)) * rhoinv;
state.dsdT = 1.5 * (C::k_B / (state.mu * m_nucleon)) * Tinv;
state.dsdr = - (C::k_B / (state.mu * m_nucleon)) * rhoinv;
state.dhdT = state.dedT + state.dpdT * rhoinv;
state.dhdr = 0.0;

Expand Down
22 changes: 11 additions & 11 deletions EOS/multigamma/actual_eos.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@ AMREX_GPU_HOST_DEVICE inline
void actual_eos (eos_input_t input, eos_t& state)
{
// Get the mass of a nucleon from Avogadro's number.
const Real m_nucleon = 1.0_rt / n_A;
const Real m_nucleon = 1.0_rt / C::n_A;

// Special gamma factors
Real sumY_gm1 = 0.0_rt;
Expand Down Expand Up @@ -101,7 +101,7 @@ void actual_eos (eos_input_t input, eos_t& state)
// Solve for the temperature:
// h = e + p/rho = (p/rho)*[1 + 1/(gamma-1)] = (p/rho)*gamma/(gamma-1)
dens = state.rho;
temp = state.h * (m_nucleon / k_B) / sumYg_gm1;
temp = state.h * (m_nucleon / C::k_B) / sumYg_gm1;

break;

Expand All @@ -111,7 +111,7 @@ void actual_eos (eos_input_t input, eos_t& state)

// Solve for the density:
// p = rho k T / (abar m_nucleon)
dens = state.p * state.abar * (m_nucleon / k_B) / state.T;
dens = state.p * state.abar * (m_nucleon / C::k_B) / state.T;
temp = state.T;

break;
Expand All @@ -123,7 +123,7 @@ void actual_eos (eos_input_t input, eos_t& state)
// Solve for the temperature:
// p = rho k T / (mu m_nucleon)
dens = state.rho;
temp = state.p * state.abar * (m_nucleon / k_B) / state.rho;
temp = state.p * state.abar * (m_nucleon / C::k_B) / state.rho;

break;

Expand All @@ -134,7 +134,7 @@ void actual_eos (eos_input_t input, eos_t& state)
// Solve for the temperature
// e = k T / [(mu m_nucleon)*(gamma-1)]
dens = state.rho;
temp = state.e * (m_nucleon / k_B) / sumY_gm1;
temp = state.e * (m_nucleon / C::k_B) / sumY_gm1;

break;

Expand All @@ -153,7 +153,7 @@ void actual_eos (eos_input_t input, eos_t& state)

// Solve for temperature and density
dens = state.p * state.abar / state.h * sumYg_gm1;
temp = state.p * state.abar * (m_nucleon / k_B) / dens;
temp = state.p * state.abar * (m_nucleon / C::k_B) / dens;

break;

Expand Down Expand Up @@ -186,18 +186,18 @@ void actual_eos (eos_input_t input, eos_t& state)

// Compute the pressure simply from the ideal gas law, and the
// specific internal energy using the gamma-law EOS relation.
state.p = dens * (k_B / m_nucleon) * temp / state.abar;
state.e = (k_B / m_nucleon) * temp * sumY_gm1;
state.p = dens * (C::k_B / m_nucleon) * temp / state.abar;
state.e = (C::k_B / m_nucleon) * temp * sumY_gm1;

// Enthalpy is h = e + p/rho
state.h = state.e + state.p / dens;

// entropy (per gram) -- this is wrong. Not sure what the expression
// is for a multigamma gas
state.s = ((k_B / m_nucleon) / state.abar) *
state.s = ((C::k_B / m_nucleon) / state.abar) *
(2.5_rt + std::log((std::pow(state.abar * m_nucleon, 2.5_rt) / dens) *
std::pow(k_B * temp, 1.5_rt) /
std::pow(2.0_rt * M_PI * hbar * hbar, 1.5_rt)));
std::pow(C::k_B * temp, 1.5_rt) /
std::pow(2.0_rt * M_PI * C::hbar * C::hbar, 1.5_rt)));

// Compute the thermodynamic derivatives and specific heats
state.dpdT = state.p / temp;
Expand Down
5 changes: 5 additions & 0 deletions Make.Microphysics
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,11 @@ ifeq ($(USE_REACT), TRUE)
DEFINES += -DREACTIONS
endif

USE_NETWORK_SOLVER ?= FALSE
ifeq ($(USE_NETWORK_SOLVER), TRUE)
DEFINES += -DNETWORK_SOLVER
endif

ifeq ($(USE_REACT_SPARSE_JACOBIAN), TRUE)
DEFINES += -DREACT_SPARSE_JACOBIAN

Expand Down
6 changes: 3 additions & 3 deletions conductivity/stellar/actual_conductivity.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ actual_conductivity(eos_t& state) {

const Real zbound = 0.1e0_rt;
const Real t7peek = 1.0e20_rt;
const Real k2c = 4.0_rt/3.0_rt*a_rad*c_light;
const Real k2c = 4.0_rt/3.0_rt*C::a_rad*C::c_light;
const Real meff = 1.194648642401440e-10_rt;
const Real weid = 6.884326138694269e-5_rt;
const Real iec = 1.754582332329132e16_rt;
Expand Down Expand Up @@ -217,7 +217,7 @@ actual_conductivity(eos_t& state) {

// from iben apj 196 525 1975 for non-degenerate regimes
if (dlog10 < drelim) {
Real zdel = state.xne/(n_A*t6*std::sqrt(t6));
Real zdel = state.xne/(C::n_A*t6*std::sqrt(t6));
Real zdell10 = std::log10(zdel);
Real eta0 = std::exp(-1.20322_rt + twoth * std::log(zdel));
Real eta02 = eta0*eta0;
Expand Down Expand Up @@ -261,7 +261,7 @@ actual_conductivity(eos_t& state) {
dnefac = 1.5_rt/eta0 * (1.0_rt - 0.8225_rt/eta02);
}
Real wpar2 = 9.24735e-3_rt * zdel *
(state.rho*n_A*(w[3]+w[4]+w[5])/state.xne + dnefac)/(std::sqrt(t6)*pefac);
(state.rho*C::n_A*(w[3]+w[4]+w[5])/state.xne + dnefac)/(std::sqrt(t6)*pefac);
Real walf = 0.5_rt * std::log(wpar2);
Real walf10 = 0.5_rt * std::log10(wpar2);

Expand Down
104 changes: 56 additions & 48 deletions constants/fundamental_constants.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,72 +3,80 @@
#ifndef _fundamental_constants_H_
#define _fundamental_constants_H_

namespace C
{
// speed of light in vacuum
constexpr amrex::Real c_light = 2.99792458e10; // cm/s

// newton's gravitational constant
constexpr amrex::Real Gconst = 6.67428e-8; // cm^3/g/s^2
// newton's gravitational constant
constexpr amrex::Real Gconst = 6.67428e-8; // cm^3/g/s^2

// new value; if uncommented initial models will need to be re-HSE'ed
// constexpr amrex::Real Gconst = 6.67384e-8; // cm^3/g/s^2
// new value; if uncommented initial models will need to be re-HSE'ed
// constexpr amrex::Real Gconst = 6.67384e-8; // cm^3/g/s^2

// boltzmann's constant
constexpr amrex::Real k_B = 1.3806488e-16; // erg/K
// boltzmann's constant
constexpr amrex::Real k_B = 1.3806488e-16; // erg/K

// planck's constant over 2pi
constexpr amrex::Real hbar = 1.054571726e-27; // erg s
// planck's constant over 2pi
constexpr amrex::Real hbar = 1.054571726e-27; // erg s

// planck's constant
constexpr amrex::Real hplanck = 6.62606957e-27; // erg s
// planck's constant
constexpr amrex::Real hplanck = 6.62606957e-27; // erg s

// avogradro's Number
constexpr amrex::Real n_A = 6.02214129e23; // mol^-1
// avogradro's Number
constexpr amrex::Real n_A = 6.02214129e23; // mol^-1

// convert eV to erg
constexpr amrex::Real ev2erg = 1.602176487e-12;
// convert eV to erg
constexpr amrex::Real ev2erg = 1.602176487e-12;

// convert MeV to eV
constexpr amrex::Real MeV2eV = 1.0e6;
// convert MeV to eV
constexpr amrex::Real MeV2eV = 1.0e6;

// mass of proton
constexpr amrex::Real m_p = 1.672621777e-24; // g
// convert MeV to grams
constexpr amrex::Real MeV2gr = (MeV2eV * ev2erg) / (c_light * c_light);

// mass of neutron
constexpr amrex::Real m_n = 1.674927351e-24; // g
// conversion factor for nuclear energy generation rate
constexpr amrex::Real enuc_conv2 = -n_A * c_light * c_light;

// mass of electron
constexpr amrex::Real m_e = 9.10938291e-28; // g
// mass of proton
constexpr amrex::Real m_p = 1.672621777e-24; // g

// speed of light in vacuum
constexpr amrex::Real c_light = 2.99792458e10; // cm/s
// mass of neutron
constexpr amrex::Real m_n = 1.674927351e-24; // g

// electron charge
// NIST: q_e = 1.602176565e-19 C
//
// C is the SI unit Coulomb; in cgs we have the definition:
// 1 C = 0.1 * |c_light| * 1 statC
// where statC is the cgs unit statCoulomb; 1 statC = 1 erg^1/2 cm^1/2
// and |c_light| is the speed of light in cgs (but without units)
constexpr amrex::Real q_e = 4.80320451e-10; // erg^1/2 cm^1/2
// mass of electron
constexpr amrex::Real m_e = 9.10938291e-28; // g

// stefan-boltzmann constant
constexpr amrex::Real sigma_SB = 5.670373e-5; // erg/s/cm^2/K^4
// electron charge
// NIST: q_e = 1.602176565e-19 C
//
// C is the SI unit Coulomb; in cgs we have the definition:
// 1 C = 0.1 * |c_light| * 1 statC
// where statC is the cgs unit statCoulomb; 1 statC = 1 erg^1/2 cm^1/2
// and |c_light| is the speed of light in cgs (but without units)
constexpr amrex::Real q_e = 4.80320451e-10; // erg^1/2 cm^1/2

// radiation constant
constexpr amrex::Real a_rad = 4.0*sigma_SB/c_light;
// stefan-boltzmann constant
constexpr amrex::Real sigma_SB = 5.670373e-5; // erg/s/cm^2/K^4

// Number of centimeters in a parsec and an AU.
// Note that since the length of an AU is defined exactly
// by IAU convention, the length of a parsec is also
// defined exactly as (6.48e5 / pi) AU.
constexpr amrex::Real AU = 1.49597871e13; // cm
constexpr amrex::Real parsec = 3.085677587679311e18; // cm
// radiation constant
constexpr amrex::Real a_rad = 4.0*sigma_SB/c_light;

// Hubble constant (in s^{-1}, converted from 100 (km/s)/Mpc by dividing by 3.08568025e19km/Mpc)
constexpr amrex::Real Hubble_const = 32.407764868e-19;
// Number of centimeters in a parsec and an AU.
// Note that since the length of an AU is defined exactly
// by IAU convention, the length of a parsec is also
// defined exactly as (6.48e5 / pi) AU.
constexpr amrex::Real AU = 1.49597871e13; // cm
constexpr amrex::Real parsec = 3.085677587679311e18; // cm

// solar mass (from http://asa.usno.navy.mil/SecK/Constants.html)
constexpr amrex::Real M_solar = 1.9884e33;
// Hubble constant (in s^{-1}, converted from 100 (km/s)/Mpc by dividing by 3.08568025e19km/Mpc)
constexpr amrex::Real Hubble_const = 32.407764868e-19;

// solar radius
constexpr amrex::Real R_solar = 6.957e10;
// solar mass (from http://asa.usno.navy.mil/SecK/Constants.html)
constexpr amrex::Real M_solar = 1.9884e33;

// solar radius
constexpr amrex::Real R_solar = 6.957e10;
};

#endif
2 changes: 1 addition & 1 deletion integration/Make.package
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,11 +87,11 @@ ifeq ($(USE_SIMPLIFIED_SDC), TRUE)
F90EXE_sources += integrator_sdc.F90
else
F90EXE_sources += integrator.F90
CEXE_HEADERS += integrator.H
endif
F90EXE_sources += integrator_scaling.F90
f90EXE_sources += integration_data.f90


INCLUDE_LOCATIONS += $(MICROPHYSICS_HOME)/integration/utils
VPATH_LOCATIONS += $(MICROPHYSICS_HOME)/integration/utils
EXTERN_CORE += $(MICROPHYSICS_HOME)/integration/utils
Expand Down
Loading

0 comments on commit bbfabaa

Please sign in to comment.