tboudreaux/libcomposition
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

feat(Standard-Compositions): Imported standard metal fractions from MESA

Browse Source 
Imported standard metal fractions from MESA to build compositions from one of the standard schemes
This commit is contained in:
poojanagrawal
2026-06-02 15:57:33 +02:00
parent 3741768893
commit 663bdcea03
7 changed files with 1947 additions and 2 deletions
Download Patch File Download Diff File Expand all files Collapse all files

346
src/composition/lib/io/standard_compositions.cpp Normal file
View File

@@ -0,0 +1,346 @@
#include "fourdst/composition/io/standard_compositions.h"
#include "fourdst/composition/io/StandardAbundancesBinary.h"
#include "fourdst/composition/composition.h"
#include "fourdst/atomic/atomicSpecies.h"
#include "fourdst/atomic/species.h"
#include "fourdst/composition/utils.h"
#include <string>
#include <vector>
#include <algorithm>
#include <fstream>
#include <numeric>
#include <stdexcept>
#include <unordered_map>
#include <print>
#include <ranges>
#include <cctype>
namespace fourdst:: composition::io {
namespace {
inline void ltrim(std::string &s) {
s.erase(
s.begin(),
std::ranges::find_if(s,
[](const unsigned char ch) {
return !std::isspace(ch);
})
);
}
inline void rtrim(std::string &s) {
s.erase(
std::find_if(
s.rbegin(),
s.rend(),
[](const unsigned char ch) {
return !std::isspace(ch);
}).base(),
s.end()
);
}
inline void trim(std::string &s) {
ltrim(s);
rtrim(s);
}
}
bool to_bool(std::string s) {
std::transform(s.begin(), s.end(), s.begin(),
[](unsigned char c){ return std::tolower(c); });
return s == "true";
}
CompositionData ChemicalFileParser::parse_compositon_data(const std::vector<char>& data,const std::string& scheme) const {
// get file and metal_fraction_scheme
// Load the file
// find the metal_fraction_scheme
// return abundances
// LOG_TRACE_L1(m_logger, "Parsing chemical abundance for: {}", scheme);
bool debug = false;
if (debug){
std::println("Parsing chemical abundance for: {}", scheme);
}
std::istringstream stream(std::string(data.begin(), data.end()));
// add error message if something goes wrong
std::string line;
int start_line = 0;
int i = 0;
CompositionData comp;
while (std::getline(stream, line)) {
// find where the end of the scheme block is
auto end_pos = std::ranges::search(line,std::format("END {}", scheme));
// exit if have reached the end of block
if (!end_pos.empty()) {
break;
}
if (start_line>0){
const size_t colon_pos = line.find(':');
line = line.substr(colon_pos+1);
line.erase(std::remove_if(line.begin(), line.end(),
[](char c){ return c == '[' || c == ']'; }),
line.end());
trim(line);
std::string item;
std::stringstream ss(line);
double val;
switch(i-start_line){
case 1:
comp.comment_str = line;
break;
case 2:
comp.he_abundance = std::pow(10.0,std::stod(line));
break;
case 3:
comp.requires_atomic_weight = to_bool(line);
break;
case 4:
while(std::getline(ss, item, ',')) {
comp.elements.push_back(item);
}
break;
case 5:
while(std::getline(ss, item, ',')) {
val = std::pow(10.0,std::stod(item));
comp.abundances.push_back(val);
}
break;
}
}
// find where the start of the scheme block is
auto start_pos = std::ranges::search(line, std::format("BEGIN {}", scheme));
if (!start_pos.empty()) {
start_line = i;
}
i+=1;
}
// if (start_pos==0):
// raise error ("Scheme {} not found", scheme)
if (debug){
std::println("he_abundance: {}", comp.he_abundance);
std::println("requires_atomic_weight: {}", comp.requires_atomic_weight);
std::println("elements: {}",comp.elements);
std::println("abundances: {}", comp.abundances);
}
return comp;
}
IsotopicPercentage ChemicalFileParser::parse_isotopic_percentage(const std::vector<char>& data,const std::string& scheme) const {
// get file and iso_scheme
// Load the file
// find the iso_scheme
// get iso_comp data
// IsotopicPercentage object
bool debug = false;
if (debug){
std::println("Parsing Isotopic Percentage for: {}", scheme);
}
std::istringstream stream(std::string(data.begin(), data.end()));
// add error message if something goes wrong
ParsedChemicalData parsed;
std::string line;
int start_line = 0;
int i = 0;
IsotopicPercentage iso;
while (std::getline(stream, line)) {
// find where the end of the scheme block is
auto end_pos = std::ranges::search(line,std::format("END {}", scheme));
// exit if have reached the end of block
if (!end_pos.empty()) {
break;
}
if (start_line>0){
const size_t colon_pos = line.find(':');
line = line.substr(colon_pos+1);
line.erase(std::remove_if(line.begin(), line.end(),
[](char c){ return c == '[' || c == ']'; }),
line.end());
trim(line);
std::string item;
std::stringstream ss(line);
parsed.push_back(line);
switch(i-start_line){
case 1:
iso.comment_str = line;
break;
case 2:
while(std::getline(ss, item, ',')) {
iso.atomic_numbers.push_back(std::stoi(item));
}
break;
case 3:
while(std::getline(ss, item, ',')) {
iso.elements.push_back(item);
}
break;
case 4:
while(std::getline(ss, item, ',')) {
iso.mass_numbers.push_back(std::stoi(item));
}
break;
case 5:
while(std::getline(ss, item, ',')) {
iso.percentages.push_back(std::stod(item));
}
break;
}
}
// find where the start of the scheme block is
auto start_pos = std::ranges::search(line, std::format("BEGIN {}", scheme));
if (!start_pos.empty()) {
start_line = i;
}
i+=1;
}
if (debug){
std::println("atomic_numbers: {}", iso.atomic_numbers);
std::println("elements: {}",iso.elements);
std::println("mass_numbers: {}", iso.mass_numbers);
std::println("percentages: {}", iso.percentages);
}
return iso;
}
}
namespace fourdst::composition {
Composition get_composition_record(const std::string& metal_fraction_scheme,
const std::string& isotopic_percentage_scheme,
double initial_z, double initial_y) {
std::vector<char> data;
io::ChemicalFileParser parser;
io::CompositionData compositions;
io::IsotopicPercentage isotopes;
data = std::ranges::to<std::vector<char>>(StandardAbundances);
compositions = parser.parse_compositon_data(data,metal_fraction_scheme);
isotopes = parser.parse_isotopic_percentage(data,isotopic_percentage_scheme);
std::string name;
std::vector<fourdst::atomic::Species> species;
// construct name of the isotopes
for (const auto [E,A] : std::ranges::views::zip(isotopes.elements, isotopes.mass_numbers)){
name = std::format("{}-{}",E,A);
// std::println("{}", name);
auto SpeciesObject = fourdst::atomic::species.at(name);
species.push_back(SpeciesObject);
// std::println("Species: {} has mass: {}", SpeciesObject.name(), SpeciesObject.mass());
}
std::vector<double> massFracs, metal_fractions;
// hydrogen and helium are treated separately since they are not metals
// H1
massFracs.push_back(std::max(0.0, std::min(1.0, 1.0 - (initial_z + initial_y))));
// H2
massFracs.push_back(0.0);
// He3
// anders & grevesse 1989 solar mass fractions
double xsol_he3,xsol_he4;
xsol_he3=2.9291e-05;
xsol_he4=2.7521e-01;
massFracs.push_back(initial_y*xsol_he3/(xsol_he3 + xsol_he4));
// He4
massFracs.push_back(initial_y*xsol_he4/(xsol_he3 + xsol_he4));
// Metals
// multiply by atomic weight if needed
if (compositions.requires_atomic_weight){
// get isotope with max abundance for each metal
std::vector<double> element_atomic_weight;
for (size_t i = 0; i < isotopes.atomic_numbers.size();) {
size_t Z = isotopes.atomic_numbers[i];
if (Z<2) {
continue;
}
double max_iso = isotopes.percentages[i];
size_t loc = i;
size_t j = i ;
while (j<isotopes.atomic_numbers.size() && isotopes.atomic_numbers[j] == Z) {
if (isotopes.percentages[j]>max_iso) {
loc = j;
}
++j;
}
element_atomic_weight.push_back(species[loc].mass());
i=j;
}
for(size_t i=0;i < compositions.abundances.size();++i){
metal_fractions.push_back(element_atomic_weight[i]*compositions.abundances[i]);
}
} else {
metal_fractions = compositions.abundances;
}
double sum = std::accumulate(metal_fractions.begin(),metal_fractions.end(),0.0);
for(size_t i=0;i < metal_fractions.size();++i){
std::println("testing: {} , {}", i,compositions.elements[i]);
metal_fractions[i] = metal_fractions[i]/sum;
}
// Z = max(0d0, min(1d0, 1d0 - (h1 + h2 + he3 + he4)))
for (size_t i = 0; i < isotopes.atomic_numbers.size();++i) {
if (isotopes.atomic_numbers[i]<2) {
continue;
}
double frac = 1e-2*isotopes.percentages[i]*species[i].mass();
}
//
Composition comp = buildCompositionFromMassFractions(species, massFracs);
return comp;
}
}