Skip to content

Commit

Permalink
clean up linear interpolation (#2664)
Browse files Browse the repository at this point in the history 
The locate() function was returning the index of the zone that is higher than the input coordinate r. The logic in interpolate accounted for this, but things were more complex than they need be. This simplifies locate and the linear interpolation.

It also adds a unit test to help understand what's happening.
  • Loading branch information
@zingale
zingale authored Dec 26, 2023
1 parent 3095255 commit c3509d4
Show file tree
Hide file tree
Showing 10 changed files with 4,322 additions and 82 deletions.
34 changes: 34 additions & 0 deletions .github/workflows/model_parser.yml
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,34 @@
name: model parser

on: [pull_request]
jobs:
model_parser:
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v4
with:
fetch-depth: 0

- name: Get submodules
run: |
git submodule update --init
cd external/Microphysics
git fetch; git checkout development
cd ../amrex
git fetch; git checkout development
cd ../..
- name: Install dependencies
run: |
sudo apt-get update -y -qq
sudo apt-get -qq -y install curl cmake jq clang g++>=9.3.0
- name: Compile model_parser test
run: |
cd Util/model_parser/test
make -j 4
- name: Run model_parser test
run: |
cd Util/model_parser/test
./Castro3d.gnu.ex
4 changes: 2 additions & 2 deletions Exec/science/flame_wave/ci-benchmarks/grid_diag.out
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# COLUMN 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
# TIMESTEP TIME MASS XMOM YMOM ZMOM ANG. MOM. X ANG. MOM. Y ANG. MOM. Z KIN. ENERGY INT. ENERGY GAS ENERGY GRAV. ENERGY TOTAL ENERGY CENTER OF MASS X-LOC CENTER OF MASS Y-LOC CENTER OF MASS Z-LOC CENTER OF MASS X-VEL CENTER OF MASS Y-VEL CENTER OF MASS Z-VEL MAXIMUM TEMPERATURE MAXIMUM DENSITY MAXIMUM T_S / T_E
0 0.0000000000000000 1.4902024346464787e+20 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 2.9583337617665660e+37 2.9583337617665660e+37 0.0000000000000000e+00 2.9583337617665660e+37 2.7306640614148750e+04 6.8067773166795473e+02 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 1.0000000000156899e+08 3.1017131702472486e+07 0.0000000000000000e+00
1 5.9694928185132600e-08 1.4902038799260451e+20 -9.8852813993167483e+18 7.0758681407801181e+23 -3.7077151117729030e+15 -7.8559831399077499e+18 7.9789428545426686e+19 1.9321854546506737e+28 2.7622571688141327e+29 2.9583392229687520e+37 2.9583392505913144e+37 0.0000000000000000e+00 2.9583392505913144e+37 2.7306640614147251e+04 6.8067725884159722e+02 0.0000000000000000e+00 -6.6335093690719216e-02 4.7482550784469004e+03 -2.4880589573803199e-05 1.0000972452748042e+08 3.1017431302318867e+07 1.9852339015703824e-12
2 1.2535934918877847e-07 1.4902056479398973e+20 -2.0759055840781603e+19 1.4862415462814482e+24 -1.6351019859135748e+16 -3.4644556495600996e+19 3.5187152213757382e+20 4.0584338401784393e+28 5.5823502529134965e+29 2.9583452181108320e+37 2.9583452739343134e+37 0.0000000000000000e+00 2.9583452739343134e+37 2.7306640614140473e+04 6.8067697458181988e+02 0.0000000000000000e+00 -1.3930329595435043e-01 9.9733989623249036e+03 -1.0972324445112568e-04 1.0003863648020169e+08 3.1017667549506564e+07 1.9917944188178796e-12
1 5.9694928185132600e-08 1.4902038799260451e+20 -9.8852814000373391e+18 7.0758681407801342e+23 -3.7077151119947915e+15 -7.8559831403383900e+18 7.9789428550115738e+19 1.9321854546506743e+28 2.7622571688115688e+29 2.9583392229687520e+37 2.9583392505913144e+37 0.0000000000000000e+00 2.9583392505913144e+37 2.7306640614147251e+04 6.8067725884159722e+02 0.0000000000000000e+00 -6.6335093695554737e-02 4.7482550784469113e+03 -2.4880589575292177e-05 1.0000972452758998e+08 3.1017431302318867e+07 2.1834522208251888e-12
2 1.2535934918877847e-07 1.4902056479398969e+20 -2.0759055857352421e+19 1.4862415462814437e+24 -1.6351019861504854e+16 -3.4644556500407476e+19 3.5187152219157529e+20 4.0584338401784253e+28 5.5823502529036681e+29 2.9583452181108315e+37 2.9583452739343125e+37 0.0000000000000000e+00 2.9583452739343125e+37 2.7306640614140484e+04 6.8067697458182033e+02 0.0000000000000000e+00 -1.3930329606554864e-01 9.9733989623248763e+03 -1.0972324446702355e-04 1.0003863648118678e+08 3.1017667549506564e+07 1.8217100523183610e-12
2 changes: 1 addition & 1 deletion Exec/science/flame_wave/ci-benchmarks/species_diag.out
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,4 +2,4 @@
# TIMESTEP TIME Mass He4 Mass C12 Mass O16 Mass Ne20 Mass Mg24 Mass Si28 Mass S32 Mass Ar36 Mass Ca40 Mass Ti44 Mass Cr48 Mass Fe52 Mass Ni56
0 0.0000000000000000 3.5899765861250859e-15 7.4944801491560746e-24 7.4944801491560746e-24 7.4944801491560746e-24 7.4944801491560746e-24 7.4944801491560746e-24 7.4944801491560746e-24 7.4944801491560746e-24 7.4944801491560746e-24 7.4944801491560746e-24 7.4944801491560746e-24 7.4944801491560746e-24 7.1354824912929246e-14
1 5.9694928185132600e-08 3.5899765861251656e-15 7.4944874211898590e-24 7.4944874177116286e-24 7.4944874177116286e-24 7.4944874177116286e-24 7.4944874177116286e-24 7.4944874177116286e-24 7.4944874177116286e-24 7.4944874177116286e-24 7.4944874177116286e-24 7.4944874177116286e-24 7.4944874177116286e-24 7.1354897598484802e-14
2 1.2535934918877847e-07 3.5899765861253060e-15 7.4944963166333664e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.1354986514892610e-14
2 1.2535934918877847e-07 3.5899765861253060e-15 7.4944963166333664e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.4944963093524191e-24 7.1354986514892597e-14
46 changes: 23 additions & 23 deletions Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out
View file
Original file line number Diff line number Diff line change
@@ -1,29 +1,29 @@
plotfile = plt00095
time = 1.25
variables minimum value maximum value
density 8.7136176158e-05 13348283.786
xmom -4.4636648292e+14 1.4969028735e+15
ymom -1.8931343553e+15 1.9807937913e+15
density 8.7136170328e-05 13347626.831
xmom -4.4623942798e+14 1.4982725823e+15
ymom -1.8879302767e+15 1.980304858e+15
zmom 0 0
rho_E 7.7947390767e+11 5.6568077669e+24
rho_e 7.4321344551e+11 5.6343409475e+24
Temp 100000 3972527783.9
rho_He4 8.7136176158e-17 1473.9666386
rho_C12 3.4854470463e-05 5030539.1488
rho_O16 5.2281705694e-05 7778301.6377
rho_Ne20 8.7136176158e-17 1023673.1153
rho_Mg24 8.7136176158e-17 1040419.2782
rho_Si28 8.7136176158e-17 4251082.0739
rho_S32 8.7136176158e-17 2179431.2961
rho_Ar36 8.7136176158e-17 497747.48798
rho_Ca40 8.7136176158e-17 382056.037
rho_Ti44 8.7136176158e-17 1576.0930955
rho_Cr48 8.7136176158e-17 1467.9139369
rho_Fe52 8.7136176158e-17 14831.710059
rho_Ni56 8.7136176158e-17 182702.27304
phiGrav -4.6147467267e+17 -2.2055818332e+16
grav_x -461195258.85 -48603.568291
grav_y -444709392.81 392306861.64
rho_E 7.7948360809e+11 5.64675518e+24
rho_e 7.4322366401e+11 5.6244082614e+24
Temp 100000 3969215375.1
rho_He4 8.7136170328e-17 1473.6068478
rho_C12 3.4854468131e-05 5028784.2519
rho_O16 5.2281702196e-05 7774667.7686
rho_Ne20 8.7136170328e-17 448739.14288
rho_Mg24 8.7136170328e-17 1131012.4412
rho_Si28 8.7136170328e-17 4244213.8084
rho_S32 8.7136170328e-17 2178795.6667
rho_Ar36 8.7136170328e-17 496791.83986
rho_Ca40 8.7136170328e-17 381536.89513
rho_Ti44 8.7136170328e-17 1576.1652647
rho_Cr48 8.7136170328e-17 1467.4931808
rho_Fe52 8.7136170328e-17 14841.830382
rho_Ni56 8.7136170328e-17 182905.91385
phiGrav -4.6146743474e+17 -2.2055817929e+16
grav_x -461199998.92 -48603.52543
grav_y -444710966.19 392312565.37
grav_z 0 0
rho_enuc -7.6356851771e+21 5.7259582003e+26
rho_enuc -7.4954583352e+21 6.7150736332e+26

84 changes: 28 additions & 56 deletions Util/model_parser/model_parser.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,13 @@ namespace model_string
}


///
/// return the index into the model coordinate, loc, such that
/// model::profile(model_index).r(loc) < r < model::profile(model_index).r(loc+1)
/// if r < model::profile(model_index).r(0) then we return 0, likewise
/// if r > model::profile(model_index).r(model::npts-2) then we return model::npt-2,
/// since this will give us the interval [npts-2, npts-1] to interpolate in
///
AMREX_INLINE AMREX_GPU_HOST_DEVICE
int
locate(const Real r, const int model_index) {
Expand All @@ -74,7 +81,7 @@ locate(const Real r, const int model_index) {
loc = 0;

} else if (r > model::profile(model_index).r(model::npts-2)) {
loc = model::npts-1;
loc = model::npts-2;

} else {

Expand All @@ -91,7 +98,7 @@ locate(const Real r, const int model_index) {
}
}

loc = ihi;
loc = ilo;
}

return loc;
Expand All @@ -102,73 +109,38 @@ AMREX_INLINE AMREX_GPU_HOST_DEVICE
Real
interpolate(const Real r, const int var_index, const int model_index=0) {

// find the value of model_state component var_index at point r
// using linear interpolation. Eventually, we can do something
// fancier here.
// this gives us an index such that profile.r(id) < r < profile.r(id+1)

int id = locate(r, model_index);

Real slope;
Real interp;

if (id == 0) {

slope = (model::profile(model_index).state(id+1, var_index) -
model::profile(model_index).state(id, var_index)) /
(model::profile(model_index).r(id+1) - model::profile(model_index).r(id));
interp = slope * (r - model::profile(model_index).r(id)) + model::profile(model_index).state(id, var_index);

// safety check to make sure interp lies within the bounding points
Real minvar = amrex::min(model::profile(model_index).state(id+1, var_index),
model::profile(model_index).state(id, var_index));
Real maxvar = amrex::max(model::profile(model_index).state(id+1, var_index),
model::profile(model_index).state(id, var_index));
interp = amrex::max(interp, minvar);
interp = amrex::min(interp, maxvar);

} else if (id == model::npts-1) {

slope = (model::profile(model_index).state(id, var_index) -
model::profile(model_index).state(id-1, var_index)) /
(model::profile(model_index).r(id) - model::profile(model_index).r(id-1));
interp = slope * (r - model::profile(model_index).r(id)) + model::profile(model_index).state(id, var_index);


// safety check to make sure interp lies within the bounding points
Real minvar = amrex::min(model::profile(model_index).state(id-1, var_index),
model::profile(model_index).state(id, var_index));
Real maxvar = amrex::max(model::profile(model_index).state(id-1, var_index),
model::profile(model_index).state(id, var_index));
interp = amrex::max(interp, minvar);
interp = amrex::min(interp, maxvar);

} else {

if (r >= model::profile(model_index).r(id)) {

slope = (model::profile(model_index).state(id+1, var_index) -
model::profile(model_index).state(id, var_index)) /
(model::profile(model_index).r(id+1) - model::profile(model_index).r(id));
interp = slope * (r - model::profile(model_index).r(id)) + model::profile(model_index).state(id, var_index);

} else {

slope = (model::profile(model_index).state(id, var_index) -
model::profile(model_index).state(id-1, var_index)) /
(model::profile(model_index).r(id) - model::profile(model_index).r(id-1));
interp = slope * (r - model::profile(model_index).r(id)) + model::profile(model_index).state(id, var_index);

}

slope = (model::profile(model_index).state(id+1, var_index) -
model::profile(model_index).state(id, var_index)) /
(model::profile(model_index).r(id+1) - model::profile(model_index).r(id));
interp = slope * (r - model::profile(model_index).r(id)) +
model::profile(model_index).state(id, var_index);

// safety check to make sure interp lies within the bounding points. We don't
// do this at the lower boundary, which usually corresponds to the center of the star.
if (r >= model::profile(model_index).r(0)) {
Real minvar = std::min(model::profile(model_index).state(id+1, var_index),
model::profile(model_index).state(id, var_index));
Real maxvar = std::max(model::profile(model_index).state(id+1, var_index),
model::profile(model_index).state(id, var_index));
interp = std::clamp(interp, minvar, maxvar);
}

return interp;

}

// Subsample the interpolation to get an averaged profile. For this we need to know the
// 3D coordinate (relative to the model center) and cell size.

///
/// Subsample the interpolation to get an averaged profile. For this we need to know the
/// 3D coordinate (relative to the model center) and cell size.
///
AMREX_GPU_HOST_DEVICE AMREX_INLINE
Real interpolate_3d (const Real* loc, const Real* dx, int var_index, int nsub = 1, int model_index = 0)
{
Expand Down
33 changes: 33 additions & 0 deletions Util/model_parser/test/GNUmakefile
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,33 @@
PRECISION = DOUBLE
PROFILE = FALSE

DEBUG = FALSE

DIM = 3

COMP = gnu

USE_MPI = FALSE
USE_OMP = FALSE

USE_ALL_CASTRO = FALSE
USE_AMR_CORE = FALSE

# define the location of the CASTRO top directory
CASTRO_HOME ?= ../../..

USE_MODEL_PARSER = TRUE
NUM_MODELS := 2

# This sets the EOS directory in Castro/EOS
EOS_DIR := helmholtz

# This sets the network directory in Castro/Networks
NETWORK_DIR := subch_simple

EXTERN_SEARCH += .

Bpack := ./Make.package
Blocs := . .. $(CASTRO_HOME)/Source/problems

include $(CASTRO_HOME)/Exec/Make.Castro
4 changes: 4 additions & 0 deletions Util/model_parser/test/Make.package
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
CEXE_sources += main.cpp
CEXE_sources += extern_parameters.cpp
CEXE_headers += extern_parameters.H

4 changes: 4 additions & 0 deletions Util/model_parser/test/README.md
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
# Unit test for model_parser

This simply reads in an initial model and does some checks to make
sure the locate and interpolation routines work as expected.
69 changes: 69 additions & 0 deletions Util/model_parser/test/main.cpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,69 @@
#include <iostream>

#include <model_parser.H>


int main(int argc, char *argv[]) {

amrex::Initialize(argc, argv);

// initialize the external runtime parameters in C++

init_extern_parameters();

// now initialize the C++ Microphysics
#ifdef REACTIONS
network_init();
#endif

Real small_temp = 1.e-200;
Real small_dens = 1.e-200;
eos_init(small_temp, small_dens);

std::string model = "sub_chandra.M_WD-1.10.M_He-0.050.hse.CO.N14.N.10.00km";
read_model_file(model);

// test locate and ensure that the index we get is such that
// profile.r(idx) < r < profile.r(idx+1)

Real r{3.89e7};

std::cout << "testing locate" << std::endl;

auto idx = locate(r, 0);
AMREX_ALWAYS_ASSERT(r >= model::profile(0).r(idx) &&
r <= model::profile(0).r(idx+1));

std::cout << "testing interpolate" << std::endl;

// density is monotonically decreasing
// test to make sure that we see that

auto dens = interpolate(r, model::idens);
AMREX_ALWAYS_ASSERT(dens <= model::profile(0).state(idx, model::idens) &&
dens >= model::profile(0).state(idx+1, model::idens));

// test the bounds. Our model spans r = [5.e5, 4.0955e9]

std::cout << "testing boundaries of the model" << std::endl;

auto idx_lo_bnd = locate(-1.0, 0);
AMREX_ALWAYS_ASSERT(idx_lo_bnd == 0);
std::cout << "value a r < 0 = " << interpolate(-1.0, model::idens) << std::endl;

auto idx_hi_bnd = locate(4.1e9, 0);
AMREX_ALWAYS_ASSERT(idx_hi_bnd == model::npts-2);
std::cout << "value a r > r_max = " << interpolate(4.1e9, model::idens) << std::endl;

// test if we interpolate to exactly a point in the model that we
// recover the data there to roundoff

std::cout << "testing interpolating exactly on model point" << std::endl;

int idx_test{100};
Real r_test = model::profile(0).r(idx_test);
auto dens_test = interpolate(r_test, model::idens);

AMREX_ALWAYS_ASSERT(std::abs(dens_test - model::profile(0).state(idx_test, model::idens)) < 1.e-15_rt);

}
Loading

0 comments on commit c3509d4

Please sign in to comment.