Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix flux eos #387

Merged
merged 3 commits into from
Nov 19, 2021
Merged

Fix flux eos #387

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
e658f6b
Fix asinh
G-071 Nov 17, 2021
dc286a8
Fixed deg calculation order
G-071 Nov 17, 2021
4ca5886
Added tests for eos=wd
G-071 Nov 17, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions octotiger/unitiger/hydro_impl/flux_kernel_interface.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,8 +63,8 @@ CUDA_GLOBAL_METHOD inline T pow_wrapper(const T& tmp1, const double& tmp2) {
return std::pow(tmp1, tmp2);
}
template <typename T>
CUDA_GLOBAL_METHOD inline T asin_wrapper(const T& tmp1) {
return std::asin(tmp1);
CUDA_GLOBAL_METHOD inline T asinh_wrapper(const T& tmp1) {
return std::asinh(tmp1);
}
template <typename T>
CUDA_GLOBAL_METHOD inline bool skippable(const T& tmp1) {
Expand Down Expand Up @@ -121,7 +121,7 @@ CUDA_GLOBAL_METHOD inline double_t inner_flux_loop(const double omega, const siz

const double_t edeg_tmp1 = rho * hdeg - pdeg;
const double_t edeg_tmp2 = 2.4 * A_ * x_pow_5;
const double_t pdeg_tmp1 = A_ * (x * (2 * x_sqr - 3) * x_sqr_sqrt + 3 * asin_wrapper(x));
const double_t pdeg_tmp1 = A_ * (x * (2 * x_sqr - 3) * x_sqr_sqrt + 3 * asinh_wrapper(x));
const double_t pdeg_tmp2 = 1.6 * A_ * x_pow_5;
select_wrapper(edeg, (x > 0.001), edeg_tmp1, edeg_tmp2);
select_wrapper(pdeg, (x > 0.001), pdeg_tmp1, pdeg_tmp2);
Expand Down Expand Up @@ -166,7 +166,7 @@ CUDA_GLOBAL_METHOD inline double_t inner_flux_loop(const double omega, const siz
const double_t hdeg = 8.0 * A_ * Binv * (x_sqr_sqrt - 1.0);
const double_t edeg_tmp1 = rho * hdeg - pdeg;
const double_t edeg_tmp2 = 2.4 * A_ * x_pow_5;
const double_t pdeg_tmp1 = A_ * (x * (2 * x_sqr - 3) * x_sqr_sqrt + 3 * asin_wrapper(x));
const double_t pdeg_tmp1 = A_ * (x * (2 * x_sqr - 3) * x_sqr_sqrt + 3 * asinh_wrapper(x));
const double_t pdeg_tmp2 = 1.6 * A_ * x_pow_5;
select_wrapper(edeg, (x > 0.001), edeg_tmp1, edeg_tmp2);
select_wrapper(pdeg, (x > 0.001), pdeg_tmp1, pdeg_tmp2);
Expand Down
23 changes: 15 additions & 8 deletions octotiger/unitiger/hydro_impl/flux_kernel_templates.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,8 +65,8 @@ CUDA_GLOBAL_METHOD inline double pow_wrapper<double>(const double& tmp1, const d
return std::pow(tmp1, tmp2);
}
template <>
CUDA_GLOBAL_METHOD inline double asin_wrapper<double>(const double& tmp1) {
return std::asin(tmp1);
CUDA_GLOBAL_METHOD inline double asinh_wrapper<double>(const double& tmp1) {
return std::asinh(tmp1);
}
template <>
CUDA_GLOBAL_METHOD inline bool skippable<bool>(const bool& tmp1) {
Expand Down Expand Up @@ -111,12 +111,14 @@ CUDA_GLOBAL_METHOD inline double_t cell_inner_flux_loop(const double omega, cons
const auto x_pow_5 = x_sqr * x_sqr * x;
const double_t hdeg = 8.0 * A_ * Binv * (x_sqr_sqrt - 1.0);


const double_t pdeg_tmp1 = A_ * (x * (2 * x_sqr - 3) * x_sqr_sqrt + 3 * asinh_wrapper(x));
const double_t pdeg_tmp2 = 1.6 * A_ * x_pow_5;
select_wrapper(pdeg, (x < 0.001), pdeg_tmp2, pdeg_tmp1);

const double_t edeg_tmp1 = rho * hdeg - pdeg;
const double_t edeg_tmp2 = 2.4 * A_ * x_pow_5;
const double_t pdeg_tmp1 = A_ * (x * (2 * x_sqr - 3) * x_sqr_sqrt + 3 * asin_wrapper(x));
const double_t pdeg_tmp2 = 1.6 * A_ * x_pow_5;
select_wrapper(edeg, (x > 0.001), edeg_tmp1, edeg_tmp2);
select_wrapper(pdeg, (x > 0.001), pdeg_tmp1, pdeg_tmp2);

dpdeg_drho = 8.0 / 3.0 * A_ * Binv * x_sqr / x_sqr_sqrt;
}
Expand Down Expand Up @@ -158,19 +160,24 @@ CUDA_GLOBAL_METHOD inline double_t cell_inner_flux_loop(const double omega, cons
dpdeg_drho = static_cast<double_t>(0.0);

// all workitems choose the same path
// from to_prim
if (A_ != 0.0) {
const auto Binv = 1.0 / B_;
const auto x = pow_wrapper(rho * Binv, 1.0 / 3.0);
const auto x_sqr = x * x;
const auto x_sqr_sqrt = sqrt_wrapper(x_sqr + 1.0);
const auto x_pow_5 = x_sqr * x_sqr * x;
const double_t hdeg = 8.0 * A_ * Binv * (x_sqr_sqrt - 1.0);


const double_t pdeg_tmp1 = A_ * (x * (2 * x_sqr - 3) * x_sqr_sqrt + 3 * asinh_wrapper(x));
const double_t pdeg_tmp2 = 1.6 * A_ * x_pow_5;
select_wrapper(pdeg, (x < 0.001), pdeg_tmp2, pdeg_tmp1);

const double_t edeg_tmp1 = rho * hdeg - pdeg;
const double_t edeg_tmp2 = 2.4 * A_ * x_pow_5;
const double_t pdeg_tmp1 = A_ * (x * (2 * x_sqr - 3) * x_sqr_sqrt + 3 * asin_wrapper(x));
const double_t pdeg_tmp2 = 1.6 * A_ * x_pow_5;
select_wrapper(edeg, (x > 0.001), edeg_tmp1, edeg_tmp2);
select_wrapper(pdeg, (x > 0.001), pdeg_tmp1, pdeg_tmp2);

dpdeg_drho = 8.0 / 3.0 * A_ * Binv * x_sqr / x_sqr_sqrt;
}
ek = 0.0;
Expand Down
4 changes: 2 additions & 2 deletions octotiger/util/vec_base_wrapper.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,8 +25,8 @@ CUDA_GLOBAL_METHOD inline T pow_wrapper(const T& tmp1, const double& tmp2) {
return pow(tmp1, tmp2);
}
template <typename T>
CUDA_GLOBAL_METHOD inline T asin_wrapper(const T& tmp1) {
return asin(tmp1);
CUDA_GLOBAL_METHOD inline T asinh_wrapper(const T& tmp1) {
return asinh(tmp1);
}
template <typename T>
CUDA_GLOBAL_METHOD inline T copysign_wrapper(const T& tmp1, const T& tmp2) {
Expand Down
4 changes: 2 additions & 2 deletions octotiger/util/vec_scalar_host_wrapper.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,8 @@ inline double limiter_wrapper<double>(const double& a, const double& b) {
return minmod_theta_wrapper<double>(a, b, 64. / 37.);
}
template <>
inline double asin_wrapper<double>(const double& tmp1) {
return std::asin(tmp1);
inline double asinh_wrapper<double>(const double& tmp1) {
return std::asinh(tmp1);
}
template <>
inline bool skippable<double>(const double& tmp1) {
Expand Down
9 changes: 7 additions & 2 deletions octotiger/util/vec_vc_wrapper.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,8 +41,13 @@ inline vc_type sqrt_wrapper<vc_type>(const vc_type& tmp1) {
return Vc::sqrt(tmp1);
}
template <>
inline vc_type asin_wrapper<vc_type>(const vc_type& tmp1) {
return Vc::asin(tmp1);
inline vc_type asinh_wrapper<vc_type>(const vc_type& tmp1) {
// return Vc::asinh(tmp1);
vc_type ret = 0.0;
for (auto vec_i = 0; vec_i < vc_type::size(); vec_i++) {
ret[vec_i] = std::asinh(tmp1[vec_i]);
}
return ret;
}
template <>
inline vc_type copysign_wrapper<vc_type>(const vc_type& tmp1, const vc_type& tmp2) {
Expand Down
11 changes: 9 additions & 2 deletions src/unitiger/hydro_impl/flux_cpu_kernel.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,8 +78,15 @@ inline vc_type pow_wrapper<vc_type>(const vc_type& tmp1, const double& tmp2) {
// return ret;
}
template <>
inline vc_type asin_wrapper<vc_type>(const vc_type& tmp1) {
return Vc::asin(tmp1);
inline vc_type asinh_wrapper<vc_type>(const vc_type& tmp1) {
// not implemented
//return Vc::asinh(tmp1);

vc_type ret = 0.0;
for (auto vec_i = 0; vec_i < vc_type::size(); vec_i++) {
ret[vec_i] = std::asinh(tmp1[vec_i]);
}
return ret;
}
template <>
inline bool skippable<mask_type>(const mask_type& tmp1) {
Expand Down
4 changes: 2 additions & 2 deletions src/unitiger/hydro_impl/flux_kernel_interface.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,8 +23,8 @@ inline double pow_wrapper<double>(const double& tmp1, const double& tmp2) {
return std::pow(tmp1, tmp2);
}
template <>
inline double asin_wrapper<double>(const double& tmp1) {
return std::asin(tmp1);
inline double asinh_wrapper<double>(const double& tmp1) {
return std::asinh(tmp1);
}
template <>
inline bool skippable<double>(const double& tmp1) {
Expand Down
46 changes: 46 additions & 0 deletions test_problems/rotating_star/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -232,3 +232,49 @@ set_tests_properties(test_problems.rotating_star.init.fixture_cleanup PROPERTIES
FIXTURES_CLEANUP test_problems.rotating_star.init
)

if (OCTOTIGER_WITH_GRIDDIM EQUAL 8)
set(rho_regex "rho 1.395105e-03 1.034051e-02")
set(egas_regex "egas 1.282717e-03 1.179200e-02")
set(tau_regex "tau 6.104387e-04 4.524548e-03")
set(pot_regex "pot 4.917696e-03 3.234157e-02")
set(sx_regex "sx 2.162120e-03 1.434366e-02")
set(sy_regex "sy 2.162120e-03 1.434366e-02")
set(sz_regex "sz 2.531845e-03 1.679645e-02")
set(zx_regex "zx 2.003861e-04 1.158221e-03")
set(zy_regex "zy 2.003861e-04 1.158221e-03")
set(zz_regex "zz 5.833372e-04 2.900349e-03")
set(spc1_regex "spc_1 1.395105e-03 1.034051e-02")
set(spc2_regex "spc_2 8.687667e-13 3.239063e-12")
set(spc3_regex "spc_3 0.000000e.00 0.000000e.00")
set(spc4_regex "spc_4 0.000000e.00 0.000000e.00")
set(spc5_regex "spc_5 0.000000e.00 0.000000e.00")
set(silo_scenario_filename "none")

# Rotating Star - CPU
test_rotating_star_scenario(test_problems.cpu.eos_wd.rotating_star_vc rotating_star_eos_wd_log.txt ${silo_scenario_filename}
" --correct_am_hydro=0 --eos=WD --monopole_host_kernel_type=VC --multipole_host_kernel_type=VC --monopole_device_kernel_type=OFF --multipole_device_kernel_type=OFF --hydro_device_kernel_type=OFF --hydro_host_kernel_type=LEGACY --amr_boundary_kernel_type=AMR_LEGACY")
# Rotating Star - OLD CPU
test_rotating_star_scenario(test_problems.cpu.eos_wd.rotating_star_legacy rotating_star_eos_wd_old_log.txt ${silo_scenario_filename}
" --correct_am_hydro=0 --eos=WD --monopole_host_kernel_type=LEGACY --multipole_host_kernel_type=LEGACY --monopole_device_kernel_type=OFF --multipole_device_kernel_type=OFF --hydro_device_kernel_type=OFF --hydro_host_kernel_type=LEGACY")
if(OCTOTIGER_WITH_CUDA)
test_rotating_star_scenario(test_problems.gpu.eos_wd.rotating_star_cuda rotating_star_eos_wd_cuda_log.txt ${silo_scenario_filename}
" --correct_am_hydro=0 --eos=WD --cuda_number_gpus=1 --cuda_streams_per_gpu=32 --cuda_buffer_capacity=1024 --monopole_host_kernel_type=VC --multipole_host_kernel_type=VC --monopole_device_kernel_type=CUDA --multipole_device_kernel_type=CUDA --hydro_device_kernel_type=CUDA --hydro_host_kernel_type=LEGACY")
endif()

if(OCTOTIGER_WITH_KOKKOS)
test_rotating_star_scenario(test_problems.cpu.eos_wd.rotating_star_kokkos rotating_star_eos_wd_kokkos_log.txt ${silo_scenario_filename}
" --correct_am_hydro=0 --eos=WD --monopole_host_kernel_type=KOKKOS --multipole_host_kernel_type=KOKKOS --monopole_device_kernel_type=OFF --multipole_device_kernel_type=OFF --hydro_device_kernel_type=OFF --hydro_host_kernel_type=KOKKOS")
if(OCTOTIGER_WITH_CUDA)
test_rotating_star_scenario(test_problems.gpu.eos_wd.rotating_star_kokkos rotating_star_eos_wd_kokkos_cuda_log.txt ${silo_scenario_filename}
" --correct_am_hydro=0 --eos=WD --cuda_number_gpus=1 --cuda_streams_per_gpu=32 --cuda_buffer_capacity=1024 --monopole_host_kernel_type=KOKKOS --multipole_host_kernel_type=KOKKOS --monopole_device_kernel_type=KOKKOS_CUDA --multipole_device_kernel_type=KOKKOS_CUDA --hydro_device_kernel_type=KOKKOS_CUDA --hydro_host_kernel_type=KOKKOS")
endif()
endif()

endif()


add_test(test_problems.rotating_star.init.fixture_cleanup ${CMAKE_COMMAND} -E remove ${PROJECT_BINARY_DIR}/test_problems/rotating_star/rotating_star.bin)
set_tests_properties(test_problems.rotating_star.init.fixture_cleanup PROPERTIES
FIXTURES_CLEANUP test_problems.rotating_star.init
)