diff --git a/data/material-model/entropy-table/pyrtable/material_table_temperature_pressure_small.txt b/data/material-model/entropy-table/pyrtable/material_table_temperature_pressure_small.txt index 6c440a8238f..94b5752f916 100644 --- a/data/material-model/entropy-table/pyrtable/material_table_temperature_pressure_small.txt +++ b/data/material-model/entropy-table/pyrtable/material_table_temperature_pressure_small.txt @@ -4,5 +4,5 @@ T(K) P(bar) s,J/K/kg rho,kg/m3 alpha,1/K cp,J/K/kg vp,km/s vs,km/s h,J/kg maph,1 250.0 0.0 509.74740059944844 3332.8967443100287 1.9781287327429287e-05 748.5636509347647 8.116376883289359 4.7495273882440685 -13445454.390915 0.0 4000.0 0.0 3804.799798825367 2879.31713365359 0.00012376961201237974 1865.6003296108431 5.89690417370827 3.068046128481904 -8401781.805825338 0.0 -250.0 1500000.0 198.7646619931198 5836.877132241759 3.4121952803538903e-06 376.4469624433934 14.555814776451808 7.855984247196738 16942560.667925145 1.0 -4000.0 1500000.0 3001.438264292299 5612.196624661615 1.174389488674247e-05 1244.420157015243 14.133650137979814 7.318949741082264 21273542.53213353 1.0 +250.0 400781.25 326.1989261514726 4738.970221094591 9.15389351114028e-06 565.4531513648064 11.974450477467947 6.779618381631337 -3722250.1404551915 1.0 +4000.0 400781.25 3394.028743328079 4373.768154238075 2.6303932404721963e-05 1362.8625515701538 10.986106548987486 5.8231315598976465 873079.6281181354 1.0 diff --git a/source/material_model/utilities.cc b/source/material_model/utilities.cc index 651f86f848f..e728e74bce7 100644 --- a/source/material_model/utilities.cc +++ b/source/material_model/utilities.cc @@ -1207,13 +1207,16 @@ namespace aspect Assert (in.pressure/1e5 >= minimum_pressure[i] && in.pressure/1e5 < maximum_pressure[i], ExcInternalError()); // const unsigned maph = unsigned(material_lookup[0]->get_data(temperature_pressure, 7)); - const std::vector& maph_points = material_lookup[0]->get_interpolation_point_coordinates(7); + const std::vector& temperature_points = material_lookup[i]->get_interpolation_point_coordinates(0); + const std::vector& pressure_points = material_lookup[i]->get_interpolation_point_coordinates(1); const unsigned i_T = static_cast(std::round((in.temperature - minimum_temperature[i]) / interval_temperature[i])); const unsigned i_p = static_cast(std::round((in.pressure/1e5 - minimum_pressure[i]) / interval_pressure[i])); - const unsigned i_point = i_p * material_lookup[0]->get_table_points(0) + i_T; + // const unsigned i_point = i_p * material_lookup[0]->get_table_points(0) + i_T; - const unsigned maph = unsigned(maph_points[i_point]); + Point<2> temperature_pressure(temperature_points[i_T], pressure_points[i_p]); + + const unsigned maph = unsigned(material_lookup[i]->get_data(temperature_pressure, 7)); // determine the value if (in.phase_index < maph) @@ -1221,6 +1224,10 @@ namespace aspect else function_value = 0.0; + if (in.pressure > 2e10){ + std::cout << "foo" << std::endl; + } + return function_value; }