CCSN Initilization #135

carlnotsagan · 2022-10-27T14:10:58Z

Reads in 1D stellar profile.
Maps to phoebus grid and fills matter variables
Assumes dr is the same for Monopole solver and input file
Has check for npoints not equal to input file
Converts units to scale free assumed by monopole solver.

Yurlungur · 2022-12-13T06:37:57Z

I didn't make it to this yet. Sorry @carlnotsagan I was sick last week and on travel before that. This is on my priority list.

Yurlungur

This actually all looks reasonable. I don't see the units issue off the top of my head though I do have some nitpicks. Might need to debug with the code + the input progenitor.

Yurlungur · 2022-12-13T17:30:00Z

src/ccsn/ccsn.cpp

+  std::pair<int, int> p = Get1DProfileNumZones(model_filename);
+  const int num_zones = p.first - 1;
+  const int num_comments = p.second;
+  const int num_vars = 11;


Since this is always the same, I suggest moving it to the ccsn.hpp header file and defining it as a global constant

constexpr int NUM_VARS = 11;

inside the CCSN namespace. Unless you think this might change later? In which case, putting it here and making it a param is the right thing.

Yurlungur · 2022-12-13T17:30:46Z

src/ccsn/ccsn.cpp

+  if (num_zones > 0 && npoints != num_zones) {
+    std::stringstream msg;
+    msg << "When using a stellar input profile Monopole GR npoints must equal the number "
+           "of zones in input file:"
+        << num_zones << std::endl;
+    PARTHENON_THROW(msg);
+  }


You're doing interpolation so is this actually true? I think this restriction can be relaxed.

Yurlungur · 2022-12-13T17:31:59Z

src/ccsn/ccsn.cpp

+  auto matter_h = monopolepkg->Param<MonopoleGR::Matter_host_t>("matter_h");
+  auto state_h = params.Get<CCSN::State_host_t>("ccsn_state_interp_h");
+
+  // printf("Current working dir: %s\n", get_current_dir_name());


Suggested change

// printf("Current working dir: %s\n", get_current_dir_name());

Yurlungur · 2022-12-13T17:34:44Z

src/ccsn/ccsn.cpp

+      Real yint = MonopoleGR::Interpolate(xa, state_raw, radius,
+                                          j); //  call to some interp routine


Weird spacing here. Also this is a linear interpolation routine.

Yurlungur · 2022-12-13T17:38:09Z

src/ccsn/ccsn.cpp

+  printf("MESA model interpolated to monopoleGR rgrid. Central density = %.14e (g/cc)\n",
+         ccsn_state_interp_d(CCSN::RHO, 0));
+  printf(
+      "MESA model interpolated to monopoleGR rgrid. Central electron fraction = %.14e\n",
+      ccsn_state_interp_d(CCSN::YE, 0));
+  printf("MESA model interpolated to monopoleGR rgrid. Central pressure = %.14e "
+         "(dyn/cm**2)\n",
+         ccsn_state_interp_d(CCSN::PRES, 0));
+  printf(
+      "MESA model interpolated to monopoleGR rgrid. Central temperature = %.14e (GK)\n",
+      ccsn_state_interp_d(CCSN::TEMP, 0));


These should be behind a check for whether or not you're on the root MPI rank. You need
#include <globals.hpp>
which includes the relevant header file from Parthenon and then

if (parthenon::Globals::my_rank == 0) { // The print statements }

this will ensure that when you run the code in parallel, the messages are only printed once.

Yurlungur · 2022-12-13T17:43:30Z

src/ccsn/ccsn.cpp

+    // convert from specific eint then to code units
+    Real mass = ccsn_state_interp_h(CCSN::RHO, i) * (4. / 3.) *
+                (pi)*std::pow(ccsn_state_interp_h(CCSN::R, i), 3.);
+
+    ccsn_state_conv_interp_h(CCSN::EPS, i) =
+        ccsn_state_interp_h(CCSN::EPS, i) * mass * unit_conv.GetEnergyCGSToCode();


We will never want total energy in a cell (and what you've computed here is total energy enclosed by a given radius, which is definitely not what we want). We either want specific internal energy or we want internal energy by volume:
$u = \varepsilon \rho$
also the monopole solver actually wants specific energy $\varepsilon$, not energy density $u$. So I would leave this as $\varepsilon$ with appropriate units of unit_conv.GetEnergyCGSToCode()/unit_conv.GetMassCGSToCode() and only convert to energy density $u$ in the problem generator when you set the 3D grid variables.

Yurlungur · 2022-12-13T17:48:51Z

src/ccsn/ccsn.cpp

+    Real pres = ccsn_state_conv_interp_h(CCSN::PRES, i); // pressure
+    Real rho = ccsn_state_conv_interp_h(CCSN::RHO, i);   // mass density
+    Real eps = ccsn_state_conv_interp_h(
+        CCSN::EPS, i); // specific internal energy - need to multiply by mass?


These aren't used here so no reason to pull them out.

Yurlungur · 2022-12-13T17:49:18Z

src/ccsn/ccsn.cpp

+    // floor based on pressure
+    if (pres <= 1.1 * Pmin) {
+      pres = rho = eps = 0;
+    } else {
+      // PolytropeThermoFromP(press, K, Gamma, rho, eps);
+    }


since you're passing in the monopoleGR quantities directly into the monopole solver, you probably don't need this.

Yurlungur · 2022-12-13T17:49:50Z

src/ccsn/ccsn.cpp

+  // Pressure floor
+  Real pmin = pin->GetOrAddReal("ccsn", "Pmin", 1e-9);
+  params.Add("pmin", pmin);
+
+  // Mass density floor
+  Real rhomin = pin->GetOrAddReal("ccsn", "rhomin", 1e-12);
+  params.Add("rhomin", rhomin);
+
+  // Internal energy floor
+  Real epsmin = pin->GetOrAddReal("ccsn", "epsmin", 1e-12);
+  params.Add("epsmin", epsmin);


... which probably means you don't need these floors.

Yurlungur · 2022-12-13T23:41:55Z

src/ccsn/ccsn.hpp

+KOKKOS_INLINE_FUNCTION
+std::pair<int, int> Get1DProfileNumZones(const std::string model_filename) {
+
+  // open file
+  std::ifstream inputfile(model_filename);
+  const std::string whitespace(" \t\n\r");
+
+  // error check
+  if (!inputfile.is_open()) {
+    std::cout << model_filename << " not found :( \n.";
+    exit(-1);
+  }
+
+  int nz = 0;
+  int nc = 0;
+  std::string line;
+  std::string comment1("#");
+  std::string comment2("//");
+
+  // get number of zones from 1d file
+  while (!inputfile.eof()) {
+
+    getline(inputfile, line);
+
+    std::size_t first_nonws = line.find_first_not_of(whitespace);
+
+    // skip empty lines
+    if (first_nonws == std::string::npos) {
+      continue;
+    }
+
+    // skip comments
+    if (line.find(comment1) == first_nonws || line.find(comment2) == first_nonws) {
+      nc++;
+      continue;
+    }
+
+    nz++;
+  }
+
+  inputfile.close();
+
+  return std::make_pair(nz + 1, nc);
+}
+
+template <typename ParArray2D>
+KOKKOS_INLINE_FUNCTION void Get1DProfileData(const std::string model_filename,
+                                             const int num_zones, const int num_comments,
+                                             const int num_vars,
+                                             ParArray2D &ccsn_state_raw_h) {
+
+  std::ifstream inputfile(model_filename);
+
+  Real val = 0;
+  std::string line;
+  int line_num = 0;
+
+  // skip comment lines
+  while (line_num < num_comments) {
+    getline(inputfile, line);
+    line_num++;
+  }
+
+  // read file into model_1d
+  for (int i = 0; i < num_zones; i++) // number of zones
+  {
+    for (int j = 0; j < num_vars; j++) //  number of vars
+    {
+      inputfile >> val;
+      ccsn_state_raw_h(j, i) = val;
+    }
+  }
+
+  std::cout << "Read 1D profile into state array.\n";
+
+  inputfile.close();
+}


I think there's probably a way to read the data without looping through the file twice. But that's a problem for future us. This is fine for now.

carlnotsagan · 2022-12-14T18:56:32Z

Thanks for these comments. I have been working on chicoma so may need to wait until that is back up. Sounds like there is still a units issue and these are mostly design based. i'll go through them when i can and update once done!

Yurlungur · 2022-12-14T19:00:43Z

I thought I might try to download and poke at it myself to see if I could find the issue. Might ping you for the progenitor file.

update ccsn files

a3f87df

carlnotsagan added the enhancement New feature or request label Oct 27, 2022

carlnotsagan self-assigned this Oct 27, 2022

carlnotsagan added 25 commits November 1, 2022 11:03

ccsn edits to driver

54986a4

ccsn input file

9b2aef8

ccsn cmakelist

d901e9b

ccsn pgen

7747c29

clang format

e43078a

try to use monopole interp

ee7ceee

clang format

bf5fb09

update vars to match new input format

98d1fd0

fill matter arrays with appropriate data from input files

51adf88

update var names

f02188a

update txt dump routine

0d588b5

update pgen

b1d87bf

update call to interp

b1762f3

clang format

d18c681

update input deck

139ffe4

Merge branch 'main' into ccsn_init

f32f5e3

update pin

1ca6d60

add ccsn to pgen

62875b9

actually initialize CCSN pin

fed5e39

update pin defaults, pass raw radius instead to interpoltor

dda402e

formatting

7c47686

fix comment support, add throw error if npoints dont match

da25da7

start converting from cgs to code units

e4d0d84

set fluid mass and length scales

0760072

convert all units

210bb80

carlnotsagan changed the title ~~WIP: CCSN Initilization~~ CCSN Initilization Nov 21, 2022

carlnotsagan requested a review from jdolence November 21, 2022 22:08

carlnotsagan requested a review from Yurlungur November 21, 2022 22:08

use converted data in pgen

72b21dc

carlnotsagan requested a review from AstroBarker December 1, 2022 20:06

Yurlungur reviewed Dec 13, 2022

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

CCSN Initilization #135

CCSN Initilization #135

carlnotsagan commented Oct 27, 2022 •

edited

Loading

Yurlungur commented Dec 13, 2022

Yurlungur left a comment

Yurlungur Dec 13, 2022

Yurlungur Dec 13, 2022

Yurlungur Dec 13, 2022

Yurlungur Dec 13, 2022

Yurlungur Dec 13, 2022

Yurlungur Dec 13, 2022

Yurlungur Dec 13, 2022

Yurlungur Dec 13, 2022

Yurlungur Dec 13, 2022

Yurlungur Dec 13, 2022

carlnotsagan commented Dec 14, 2022

Yurlungur commented Dec 14, 2022

		Real yint = MonopoleGR::Interpolate(xa, state_raw, radius,
		j); // call to some interp routine

CCSN Initilization #135

Are you sure you want to change the base?

CCSN Initilization #135

Conversation

carlnotsagan commented Oct 27, 2022 • edited Loading

Yurlungur commented Dec 13, 2022

Yurlungur left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

carlnotsagan commented Dec 14, 2022

Yurlungur commented Dec 14, 2022

carlnotsagan commented Oct 27, 2022 •

edited

Loading