From 468e04ce0b380b04dd7a2ed875d57d3aed9aadd3 Mon Sep 17 00:00:00 2001 From: Ryan Thomas Wollaeger Date: Tue, 19 Nov 2024 21:04:53 -0700 Subject: [PATCH] Modify ASCII parser to use min/max header data. + Make changes to parser name and comments to remove references. + Remove power-law opacity generation script. Note: the new metadata is redundant with the rho-T column data, so either might be removed. --- .../example_ascii/create_mean_opac_ascii.py | 56 ------------------- .../photons/example_ascii/kap_plaw.txt | 12 ++-- .../photons/mean_opacity_photons.hpp | 47 +++++++++++----- 3 files changed, 41 insertions(+), 74 deletions(-) delete mode 100644 singularity-opac/photons/example_ascii/create_mean_opac_ascii.py diff --git a/singularity-opac/photons/example_ascii/create_mean_opac_ascii.py b/singularity-opac/photons/example_ascii/create_mean_opac_ascii.py deleted file mode 100644 index 370bdbd..0000000 --- a/singularity-opac/photons/example_ascii/create_mean_opac_ascii.py +++ /dev/null @@ -1,56 +0,0 @@ -#--------------------------------------------------------------------------------------------------# -# This script defaults to using the temperature and density bounds for opacity data found -# in Zhaohuan Zhu et al (2021), "Global 3D Radiation Hydrodynamic Simulations of Proto-Jupiter’s -# Convective EnvelopeConvective Envelope". Otherwise, it merely generates a power law for Rosseland -# and Planck-averaged opacities, and produces a table in an ASCII format equivalent to that of Zhu. -#--------------------------------------------------------------------------------------------------# -import numpy as np -import argparse - -#-- parse command line arguments -parser = argparse.ArgumentParser(description="Create power-law grey opacity ASCII table (for testing).") -parser.add_argument("fname", type=str, help="base opacity file name to create (.txt is appended).") -parser.add_argument("--Nrho", type=int, default=4, help="Number of density points (assumes log space)") -parser.add_argument("--rho", type=float, nargs=2, default=[1e-14, 0.7943282347241912], - help="min, max densities [g/cc]") -parser.add_argument("--NT", type=int, default=4, help="Number of temperature points (assumes log space)") -parser.add_argument("--T", type=float, nargs=2, default=[1.0, 7943282.347242886], - help="min, max temperatures [K]") -parser.add_argument("--rho_ref", type=float, default=1e-6, help="Reference density [g/cc]") -parser.add_argument("--rho_exp", type=float, default=0.0, help="Exponent of density / ref. density") -parser.add_argument("--T_ref", type=float, default=11604, help="Reference density [K]") -parser.add_argument("--T_exp", type=float, default=0.0, help="Exponent of temperature / ref. temperature") -parser.add_argument("--Pkap_coef", type=float, default=0.01, help="Planck opacity coefficient [cm^2/g]") -parser.add_argument("--Rkap_coef", type=float, default=0.01, help="Rosseland opacity coefficient [cm^2/g]") -args = parser.parse_args() - -#-- implement simple power law -def power_law_kappa(r, t): - '''power_law_kappa (cgs): r=density (1st arg), t=temperature (2nd arg)''' - kap_coef = [args.Rkap_coef, args.Pkap_coef] - kap_pwrs = (r / args.rho_ref)**(args.rho_exp) * (t / args.T_ref)**(args.T_exp) - return [kapc * kap_pwrs for kapc in kap_coef] - -#-- set rho and T arrays -rho = np.logspace(np.log10(args.rho[0]), np.log10(args.rho[1]), args.Nrho) -T = np.logspace(np.log10(args.T[0]), np.log10(args.T[1]), args.NT) - -#-- initialize the opacity array - columns: density, temperature, Rosseland, Planck -Nrow = args.Nrho * args.NT -opac = np.zeros((Nrow, 4)) - -#-- populate opacity array -for i in range(args.Nrho): - for j in range(args.NT): - k = j + args.NT * i - opac[k,0] = rho[i] - opac[k,1] = T[j] - opac[k,2:] = power_law_kappa(rho[i], T[j]) - -#-- save opacity array -np.savetxt(args.fname+'.txt', opac, delimiter=" ", - header="rho "+str(args.Nrho)+" T "+str(args.NT), comments=" ", fmt="%.8e") - -#--------------------------------------------------------------------------------------------------# -# End of create_mean_opac_ascii.py -#--------------------------------------------------------------------------------------------------# diff --git a/singularity-opac/photons/example_ascii/kap_plaw.txt b/singularity-opac/photons/example_ascii/kap_plaw.txt index 28c41b3..123abd7 100644 --- a/singularity-opac/photons/example_ascii/kap_plaw.txt +++ b/singularity-opac/photons/example_ascii/kap_plaw.txt @@ -1,17 +1,19 @@ - rho 4 T 4 + NRho 4 NT 4 + rho_min 1.00000000e-14 rho_max 7.94328235e-01 + T_min 1.00000000e+00 T_max 7.94328235e+06 1.00000000e-14 1.00000000e+00 1.00000000e-03 1.00000000e-01 -1.00000000e-14 1.99526231e+02 1.00000000e-03 1.00000000e-01 +1.00000000e-14 1.99526232e+02 1.00000000e-03 1.00000000e-01 1.00000000e-14 3.98107171e+04 1.00000000e-03 1.00000000e-01 1.00000000e-14 7.94328235e+06 1.00000000e-03 1.00000000e-01 4.29866235e-10 1.00000000e+00 1.00000000e-03 1.00000000e-01 -4.29866235e-10 1.99526231e+02 1.00000000e-03 1.00000000e-01 +4.29866235e-10 1.99526232e+02 1.00000000e-03 1.00000000e-01 4.29866235e-10 3.98107171e+04 1.00000000e-03 1.00000000e-01 4.29866235e-10 7.94328235e+06 1.00000000e-03 1.00000000e-01 1.84784980e-05 1.00000000e+00 1.00000000e-03 1.00000000e-01 -1.84784980e-05 1.99526231e+02 1.00000000e-03 1.00000000e-01 +1.84784980e-05 1.99526232e+02 1.00000000e-03 1.00000000e-01 1.84784980e-05 3.98107171e+04 1.00000000e-03 1.00000000e-01 1.84784980e-05 7.94328235e+06 1.00000000e-03 1.00000000e-01 7.94328235e-01 1.00000000e+00 1.00000000e-03 1.00000000e-01 -7.94328235e-01 1.99526231e+02 1.00000000e-03 1.00000000e-01 +7.94328235e-01 1.99526232e+02 1.00000000e-03 1.00000000e-01 7.94328235e-01 3.98107171e+04 1.00000000e-03 1.00000000e-01 7.94328235e-01 7.94328235e+06 1.00000000e-03 1.00000000e-01 diff --git a/singularity-opac/photons/mean_opacity_photons.hpp b/singularity-opac/photons/mean_opacity_photons.hpp index 5536a09..7c9984c 100644 --- a/singularity-opac/photons/mean_opacity_photons.hpp +++ b/singularity-opac/photons/mean_opacity_photons.hpp @@ -59,7 +59,7 @@ class MeanOpacity { NT, lNuMin, lNuMax, NNu, lambda); } - // construct Planck/Rosseland DataBox from Zhu-formatted ascii file + // construct Planck/Rosseland DataBox from ascii file MeanOpacity(const std::string &filename) : filename_(filename.c_str()) { // get number of density and temperature points @@ -72,8 +72,7 @@ class MeanOpacity { std::string extension = filePath.extension().string(); if (extension == ".txt") { - // assume this is one of the original Zhu et al (2021) ASCII files - loadZhuASCII(ff); + loadASCII(ff); #ifdef SPINER_USE_HDF } else if (extension == ".hdf5" || extension == ".h5" || extension == ".sp5") { herr_t status = H5_SUCCESS; @@ -202,8 +201,8 @@ class MeanOpacity { } } - // if .txt file, assume it is the original Zhu dust opacity file - void loadZhuASCII(std::ifstream &ff) { + // ASCII opacity file reader + void loadASCII(std::ifstream &ff) { int NRho = -1; int NT = -1; @@ -211,7 +210,7 @@ class MeanOpacity { // line read from file std::string fline; - // read 1-line header to get sizes + // read 1st line of header to get sizes std::getline(ff, fline); // tokenize fline @@ -222,19 +221,41 @@ class MeanOpacity { fl_tok = std::strtok(nullptr, " "); NRho = std::stoi(fl_tok); - // move to next token to get number of density points + // move to next token to get number of temperature points fl_tok = std::strtok(nullptr, " "); fl_tok = std::strtok(nullptr, " "); NT = std::stoi(fl_tok); - // reseize the Planck and Rosseland databoxes (number of types of opac=2) + // read 2nd line of header to get min/max density + std::getline(ff, fline); + // tokenize fline + cfline = const_cast(fline.c_str()); + fl_tok = std::strtok(cfline, " "); + fl_tok = std::strtok(nullptr, " "); + const Real RhoMin = std::stod(fl_tok); + fl_tok = std::strtok(nullptr, " "); + fl_tok = std::strtok(nullptr, " "); + const Real RhoMax = std::stod(fl_tok); + + // read 3nd line of header to get min/max temperature + std::getline(ff, fline); + // tokenize fline + cfline = const_cast(fline.c_str()); + fl_tok = std::strtok(cfline, " "); + fl_tok = std::strtok(nullptr, " "); + const Real TMin = std::stod(fl_tok); + fl_tok = std::strtok(nullptr, " "); + fl_tok = std::strtok(nullptr, " "); + const Real TMax = std::stod(fl_tok); + + // reseize the Planck and Rosseland databox lkappa_.resize(NRho, NT, 2); - // set rho-T rankges (Zhu tables are uniform in log-log rho-T space) - const Real lTMin = toLog_(1.0); - const Real lTMax = toLog_(7943282.347242886); - const Real lRhoMin = toLog_(1.0e-14); - const Real lRhoMax = toLog_(0.7943282347241912); + // set rho-T ranges + const Real lTMin = toLog_(TMin); + const Real lTMax = toLog_(TMax); + const Real lRhoMin = toLog_(RhoMin); + const Real lRhoMax = toLog_(RhoMax); lkappa_.setRange(1, lTMin, lTMax, NT); lkappa_.setRange(2, lRhoMin, lRhoMax, NRho);