diff --git a/regression/quokka-tests.ini b/regression/quokka-tests.ini index cdee09f6e..47534ffce 100644 --- a/regression/quokka-tests.ini +++ b/regression/quokka-tests.ini @@ -123,3 +123,19 @@ doVis = 1 visVar = temperature testSrcTree = . ignore_return_code = 0 + +[MetalAdvectionProblem-GPU] +buildDir = . +target = test_sne +inputFile = tests/metal_problem_regression.in +link1File = extern/grackle_data_files/input/CloudyData_UVB=HM2012.h5 +dim = 3 +cmakeSetupOpts = -DAMReX_GPU_BACKEND=CUDA +restartTest = 0 +useMPI = 1 +numprocs = 1 +compileTest = 0 +doVis = 1 +visVar = temperature +testSrcTree = . +ignore_return_code = 0 \ No newline at end of file diff --git a/src/problems/CMakeLists.txt b/src/problems/CMakeLists.txt index d710eca96..b8b5df246 100644 --- a/src/problems/CMakeLists.txt +++ b/src/problems/CMakeLists.txt @@ -58,3 +58,5 @@ add_subdirectory(PopIII) add_subdirectory(ShockCloud) add_subdirectory(StarCluster) add_subdirectory(SphericalCollapse) + +add_subdirectory(MetalAdvectionProblem) diff --git a/src/problems/MetalAdvectionProblem/CMakeLists.txt b/src/problems/MetalAdvectionProblem/CMakeLists.txt new file mode 100644 index 000000000..bc4dd63bc --- /dev/null +++ b/src/problems/MetalAdvectionProblem/CMakeLists.txt @@ -0,0 +1,7 @@ +if (AMReX_SPACEDIM EQUAL 3) + add_executable(test_sne test_sne.cpp ${QuokkaObjSources}) + if(AMReX_GPU_BACKEND MATCHES "CUDA") + setup_target_for_cuda_compilation(test_sne) + endif() + +endif() diff --git a/src/problems/MetalAdvectionProblem/data_sets.dat b/src/problems/MetalAdvectionProblem/data_sets.dat new file mode 100644 index 000000000..9e187cc57 --- /dev/null +++ b/src/problems/MetalAdvectionProblem/data_sets.dat @@ -0,0 +1,206 @@ +##----Values for R8 Model-------################ + +AMREX_GPU_MANAGED amrex::GpuArray logphi_data{5.23749982, 5.83925514, 6.19098487, 6.44028658, 6.63341552, + 6.79097415, 6.92395454, 7.03892608, 7.1401333 , 7.23047697, + 7.31202697, 7.38631194, 7.45449324, 7.5174744 , 7.57597231, + 7.63056519, 7.68172614, 7.72984716, 7.77525663, 7.81823237, + 7.85901159, 7.89779849, 7.93477008, 7.9700808 , 8.00386622, + 8.03624594, 8.06732603, 8.09720099, 8.12595536, 8.1536651 , + 8.18039866, 8.206218 , 8.23117933, 8.25533386, 8.2787285 , + 8.30140621, 8.32340645, 8.34476546, 8.36551667, 8.38569095, + 8.40531688, 8.42442095, 8.44302778, 8.46116029, 8.47883984, + 8.49608638, 8.51291857, 8.52935389, 8.54540873, 8.5610985 , + 8.57643767, 8.59143987, 8.60611797, 8.62048409, 8.6345497 , + 8.64832568, 8.66182248, 8.6750501 , 8.68801804, 8.70073526, + 8.71321029, 8.72545119, 8.73746565, 8.74926096, 8.76084405, + 8.77222155, 8.78339974, 8.79438465, 8.80518201, 8.81579732, + 8.8262358 , 8.83650248, 8.84660217, 8.85653946, 8.86631876, + 8.87594432, 8.88542019, 8.89475028, 8.90393833, 8.91298796, + 8.92190263, 8.93068568, 8.93934034, 8.94786969, 8.95627673, + 8.96456434, 8.97273531, 8.98079244, 8.98873851, 8.99657621, + 9.00430812, 9.01193675, 9.01946449, 9.02689367, 9.03422652, + 9.0414652 , 9.04861178, 9.0556683 , 9.06263668, 9.06951882}; +AMREX_GPU_MANAGED amrex::GpuArray logg_data{-9.85457856, -9.39107618, -9.19504351, -9.08451673, -9.01789453, + -8.9773358 , -8.95292764, -8.93860797, -8.93051793, -8.92611324, + -8.92379746, -8.92261559, -8.92202532, -8.92173815, -8.92160118, + -8.92153673, -8.92150669, -8.9214927 , -8.92148604, -8.92148266, + -8.92148075, -8.92147948, -8.92147848, -8.92147761, -8.9214768 , + -8.92147602, -8.92147525, -8.9214745 , -8.92147377, -8.92147304, + -8.92147233, -8.92147163, -8.92147094, -8.92147026, -8.92146959, + -8.92146894, -8.92146829, -8.92146765, -8.92146703, -8.92146641, + -8.92146581, -8.92146521, -8.92146462, -8.92146405, -8.92146348, + -8.92146293, -8.92146238, -8.92146185, -8.92146132, -8.9214608 , + -8.92146029, -8.92145979, -8.9214593 , -8.92145882, -8.92145835, + -8.92145789, -8.92145744, -8.92145699, -8.92145655, -8.92145612, + -8.92145569, -8.92145528, -8.92145487, -8.92145447, -8.92145407, + -8.92145369, -8.92145331, -8.92145294, -8.92145257, -8.92145221, + -8.92145186, -8.92145152, -8.92145118, -8.92145085, -8.92145052, + -8.9214502 , -8.92144989, -8.92144959, -8.92144929, -8.921449 , + -8.92144871, -8.92144844, -8.92144816, -8.9214479 , -8.92144764, + -8.92144738, -8.92144713, -8.92144689, -8.92144665, -8.92144642, + -8.92144619, -8.92144596, -8.92144574, -8.92144553, -8.92144531, + -8.92144511, -8.9214449 , -8.9214447 , -8.92144451, -8.92144431 + }; +AMREX_GPU_MANAGED amrex::GpuArray z_data{6.08467742e+19, 1.82540323e+20, 3.04233871e+20, 4.25927419e+20, + 5.47620968e+20, 6.69314516e+20, 7.91008065e+20, 9.12701613e+20, + 1.03439516e+21, 1.15608871e+21, 1.27778226e+21, 1.39947581e+21, + 1.52116935e+21, 1.64286290e+21, 1.76455645e+21, 1.88625000e+21, + 2.00794355e+21, 2.12963710e+21, 2.25133065e+21, 2.37302419e+21, + 2.49471774e+21, 2.61641129e+21, 2.73810484e+21, 2.85979839e+21, + 2.98149194e+21, 3.10318548e+21, 3.22487903e+21, 3.34657258e+21, + 3.46826613e+21, 3.58995968e+21, 3.71165323e+21, 3.83334677e+21, + 3.95504032e+21, 4.07673387e+21, 4.19842742e+21, 4.32012097e+21, + 4.44181452e+21, 4.56350806e+21, 4.68520161e+21, 4.80689516e+21, + 4.92858871e+21, 5.05028226e+21, 5.17197581e+21, 5.29366935e+21, + 5.41536290e+21, 5.53705645e+21, 5.65875000e+21, 5.78044355e+21, + 5.90213710e+21, 6.02383065e+21, 6.14552419e+21, 6.26721774e+21, + 6.38891129e+21, 6.51060484e+21, 6.63229839e+21, 6.75399194e+21, + 6.87568548e+21, 6.99737903e+21, 7.11907258e+21, 7.24076613e+21, + 7.36245968e+21, 7.48415323e+21, 7.60584677e+21, 7.72754032e+21, + 7.84923387e+21, 7.97092742e+21, 8.09262097e+21, 8.21431452e+21, + 8.33600806e+21, 8.45770161e+21, 8.57939516e+21, 8.70108871e+21, + 8.82278226e+21, 8.94447581e+21, 9.06616935e+21, 9.18786290e+21, + 9.30955645e+21, 9.43125000e+21, 9.55294355e+21, 9.67463710e+21, + 9.79633065e+21, 9.91802419e+21, 1.00397177e+22, 1.01614113e+22, + 1.02831048e+22, 1.04047984e+22, 1.05264919e+22, 1.06481855e+22, + 1.07698790e+22, 1.08915726e+22, 1.10132661e+22, 1.11349597e+22, + 1.12566532e+22, 1.13783468e+22, 1.15000403e+22, 1.16217339e+22, + 1.17434274e+22, 1.18651210e+22, 1.19868145e+22, 1.21085081e+22 + }; + + +AMREX_GPU_MANAGED Real z_star = 245.0 * pc; +AMREX_GPU_MANAGED Real Sigma_star = 42.0 * Msun/pc/pc; +AMREX_GPU_MANAGED Real rho_dm = 0.0064 * Msun/pc/pc/pc; +AMREX_GPU_MANAGED Real R0 = 8.e3 * pc; +AMREX_GPU_MANAGED Real ks_sigma_sfr = 2.088579882548443e-55; +AMREX_GPU_MANAGED Real hscale= 150. * pc; +AMREX_GPU_MANAGED Real sigma1 = 700000.0; +AMREX_GPU_MANAGED Real sigma2 = 7000000.0; +AMREX_GPU_MANAGED Real rho01 = 2.78556e-24; +AMREX_GPU_MANAGED Real rho02 = 2.7855600000000006e-29;; + + +################----------------------################################## + +################-----------Values for R16-----------################################## + +AMREX_GPU_MANAGED amrex::GpuArray logphi_data{0.83638381, 4.50705067, 5.10271383, 5.45268878, 5.70140736, + 5.89447928, 6.05229308, 6.1857468 , 6.30135334, 6.40331919, + 6.49451684, 6.57699822, 6.65227803, 6.72150811, 6.78558362, + 6.84521466, 6.90097381, 6.95332936, 7.00266958, 7.04931937, + 7.09355406, 7.1356083 , 7.17568429, 7.21395706, 7.25057937, + 7.28568514, 7.31939257, 7.35180662, 7.38302073, 7.41311872, + 7.44217572, 7.47025967, 7.49743193, 7.52374829, 7.54925951, + 7.57401195, 7.5980481 , 7.62140688, 7.64412421, 7.66623309, + 7.68776413, 7.70874559, 7.72920373, 7.74916295, 7.76864596, + 7.78767394, 7.80626666, 7.82444266, 7.84221924, 7.85961268, + 7.87663823, 7.89331028, 7.90964238, 7.92564731, 7.94133718, + 7.95672339, 7.9718168 , 7.98662762, 8.00116561, 8.01543998, + 8.02945951, 8.04323254, 8.056767 , 8.07007045}; +AMREX_GPU_MANAGED amrex::GpuArray logg_data{-11.24587667, -10.68390187, -10.45405356, -10.31158848, + -10.21082455, -10.13474898, -10.0750678 , -10.02714739, + -9.98806155, -9.95589857, -9.92931504, -9.90722171, + -9.88888792, -9.87375864, -9.86122421, -9.85088355, + -9.84247276, -9.83565399, -9.83007266, -9.8256095 , + -9.82214665, -9.81939435, -9.81722964, -9.81559821, + -9.81433642, -9.8133616 , -9.81265555, -9.81214602, + -9.81177066, -9.81149643, -9.81129813, -9.81115635, + -9.8110556 , -9.81098525, -9.81093618, -9.81090247, + -9.81087945, -9.81086374, -9.81085319, -9.81084604, + -9.8108412 , -9.81083792, -9.81083566, -9.81083407, + -9.81083293, -9.81083206, -9.81083139, -9.81083083, + -9.81083035, -9.81082991, -9.81082951, -9.81082913, + -9.81082876, -9.8108284 , -9.81082804, -9.81082769, + -9.81082735, -9.810827 , -9.81082666, -9.81082632, + -9.81082599, -9.81082565, -9.81082532, -9.81082499 + }; +AMREX_GPU_MANAGED amrex::GpuArray z_data{1.50900000e+20, 5.55695238e+20, 9.60490476e+20, 1.36528571e+21, + 1.77008095e+21, 2.17487619e+21, 2.57967143e+21, 2.98446667e+21, + 3.38926190e+21, 3.79405714e+21, 4.19885238e+21, 4.60364762e+21, + 5.00844286e+21, 5.41323810e+21, 5.81803333e+21, 6.22282857e+21, + 6.62762381e+21, 7.03241905e+21, 7.43721429e+21, 7.84200952e+21, + 8.24680476e+21, 8.65160000e+21, 9.05639524e+21, 9.46119048e+21, + 9.86598571e+21, 1.02707810e+22, 1.06755762e+22, 1.10803714e+22, + 1.14851667e+22, 1.18899619e+22, 1.22947571e+22, 1.26995524e+22, + 1.31043476e+22, 1.35091429e+22, 1.39139381e+22, 1.43187333e+22, + 1.47235286e+22, 1.51283238e+22, 1.55331190e+22, 1.59379143e+22, + 1.63427095e+22, 1.67475048e+22, 1.71523000e+22, 1.75570952e+22, + 1.79618905e+22, 1.83666857e+22, 1.87714810e+22, 1.91762762e+22, + 1.95810714e+22, 1.99858667e+22, 2.03906619e+22, 2.07954571e+22, + 2.12002524e+22, 2.16050476e+22, 2.20098429e+22, 2.24146381e+22, + 2.28194333e+22, 2.32242286e+22, 2.36290238e+22, 2.40338190e+22, + 2.44386143e+22, 2.48434095e+22, 2.52482048e+22, 2.56530000e+22 + }; + + +AMREX_GPU_MANAGED Real z_star = 245.0 * pc; +AMREX_GPU_MANAGED Real Sigma_star = 1.71 * Msun/pc/pc; +AMREX_GPU_MANAGED Real rho_dm = 1.4e-3 * Msun/pc/pc/pc; +AMREX_GPU_MANAGED Real R0 = 16.e3 * pc; +AMREX_GPU_MANAGED Real ks_sigma_sfr = 5.499927024044233e-57; +AMREX_GPU_MANAGED Real hscale= 150. * pc; +AMREX_GPU_MANAGED Real sigma1 = 11.* kmps; +AMREX_GPU_MANAGED Real sigma2 = 10. * 11.* kmps; +AMREX_GPU_MANAGED Real rho01 = 0.0268988*Const_mH; +AMREX_GPU_MANAGED Real rho02 = 1.e-5 * 0.0268988*Const_mH;; + + +#############----Values for R4 Model-------################ + +AMREX_GPU_MANAGED amrex::GpuArray logphi_data{3.24975864, 6.27954616, 6.8675804 , 7.2142687 , 7.46059677, + 7.65149496, 7.80718072, 7.93849124, 8.05191629, 8.15165531, + 8.24058664, 8.32076511, 8.3937092 , 8.46057409, 8.52225965, + 8.57947934, 8.63280826, 8.682718 , 8.72959937, 8.77377954, + 8.81553494, 8.85510245, 8.89268623, 8.92846337, 8.96258787, + 8.99519552, 9.02640573, 9.05632431, 9.0850453 , 9.1126533 , + 9.13922402, 9.16482567, 9.18951991, 9.21336297, 9.23640588, + 9.25869529, 9.28027478, 9.30118469, 9.32146184, 9.34113998, + 9.36025027, 9.37882149, 9.39688021, 9.41445096, 9.43155666, + 9.44821863, 9.46445669, 9.48028925, 9.49573374, 9.51080637, + 9.52552235, 9.53989592, 9.55394059, 9.56766904, 9.58109321, + 9.59422433, 9.60707313, 9.61964974, 9.63196405, 9.64402569, + 9.6558438 , 9.66742698, 9.67878337, 9.6899207}; +AMREX_GPU_MANAGED amrex::GpuArray logg_data{-8.84486253, -8.53977369, -8.41945731, -8.36855149, -8.34812648, + -8.34077673, -8.33835551, -8.33764733, -8.33744008, -8.33738188, + -8.3373655 , -8.33736028, -8.33735794, -8.33735633, -8.33735492, + -8.33735358, -8.33735228, -8.33735101, -8.33734977, -8.33734856, + -8.33734737, -8.33734622, -8.33734509, -8.33734399, -8.33734291, + -8.33734186, -8.33734084, -8.33733984, -8.33733887, -8.33733792, + -8.337337 , -8.33733609, -8.33733522, -8.33733436, -8.33733353, + -8.33733272, -8.33733193, -8.33733115, -8.3373304 , -8.33732966, + -8.33732894, -8.33732824, -8.33732756, -8.33732689, -8.33732625, + -8.33732561, -8.337325 , -8.3373244 , -8.33732382, -8.33732326, + -8.33732271, -8.33732217, -8.33732165, -8.33732115, -8.33732066, + -8.33732019, -8.33731973, -8.33731928, -8.33731885, -8.33731843, + -8.33731801, -8.33731761, -8.33731722, -8.33731684 + }; +AMREX_GPU_MANAGED amrex::GpuArray z_data{1.50900000e+20, 3.40123810e+20, 5.29347619e+20, 7.18571429e+20, + 9.07795238e+20, 1.09701905e+21, 1.28624286e+21, 1.47546667e+21, + 1.66469048e+21, 1.85391429e+21, 2.04313810e+21, 2.23236190e+21, + 2.42158571e+21, 2.61080952e+21, 2.80003333e+21, 2.98925714e+21, + 3.17848095e+21, 3.36770476e+21, 3.55692857e+21, 3.74615238e+21, + 3.93537619e+21, 4.12460000e+21, 4.31382381e+21, 4.50304762e+21, + 4.69227143e+21, 4.88149524e+21, 5.07071905e+21, 5.25994286e+21, + 5.44916667e+21, 5.63839048e+21, 5.82761429e+21, 6.01683810e+21, + 6.20606190e+21, 6.39528571e+21, 6.58450952e+21, 6.77373333e+21, + 6.96295714e+21, 7.15218095e+21, 7.34140476e+21, 7.53062857e+21, + 7.71985238e+21, 7.90907619e+21, 8.09830000e+21, 8.28752381e+21, + 8.47674762e+21, 8.66597143e+21, 8.85519524e+21, 9.04441905e+21, + 9.23364286e+21, 9.42286667e+21, 9.61209048e+21, 9.80131429e+21, + 9.99053810e+21, 1.01797619e+22, 1.03689857e+22, 1.05582095e+22, + 1.07474333e+22, 1.09366571e+22, 1.11258810e+22, 1.13151048e+22, + 1.15043286e+22, 1.16935524e+22, 1.18827762e+22, 1.20720000e+22 + }; + + +AMREX_GPU_MANAGED Real z_star = 245.0 * pc; +AMREX_GPU_MANAGED Real Sigma_star = 208 * Msun/pc/pc; +AMREX_GPU_MANAGED Real rho_dm = 2.4e-2 * Msun/pc/pc/pc; +AMREX_GPU_MANAGED Real R0 = 4.e3 * pc; +AMREX_GPU_MANAGED Real ks_sigma_sfr = 1.3857971188361884e-54; +AMREX_GPU_MANAGED Real hscale= 150. * pc; +AMREX_GPU_MANAGED Real sigma1 = 7.* kmps; +AMREX_GPU_MANAGED Real sigma2 = 10. * 7.* kmps; +AMREX_GPU_MANAGED Real rho01 =1.668*Const_mH; +AMREX_GPU_MANAGED Real rho02 = 1.e-5 * 1.668*Const_mH;; \ No newline at end of file diff --git a/src/problems/MetalAdvectionProblem/test_sne.cpp b/src/problems/MetalAdvectionProblem/test_sne.cpp new file mode 100644 index 000000000..967edbc4e --- /dev/null +++ b/src/problems/MetalAdvectionProblem/test_sne.cpp @@ -0,0 +1,550 @@ + +//============================================================================== +// TwoMomentRad - a radiation transport library for patch-based AMR codes +// Copyright 2020 Benjamin Wibking. +// Released under the MIT license. See LICENSE file included in the GitHub repo. +//============================================================================== +/// \file test_hydro3d_blast.cpp +/// \brief Defines a test problem for a 3D explosion. +/// + +#include +#include +#include +#include + +#include "AMReX.H" +#include "AMReX_BC_TYPES.H" +#include "AMReX_BLassert.H" +#include "AMReX_Config.H" +#include "AMReX_FabArrayUtility.H" +#include "AMReX_GpuDevice.H" +#include "AMReX_MultiFab.H" +#include "AMReX_ParallelContext.H" +#include "AMReX_ParallelDescriptor.H" +#include "AMReX_ParmParse.H" +#include "AMReX_Print.H" +#include "AMReX_Random.H" +#include "AMReX_RandomEngine.H" +#include "AMReX_SPACE.H" +#include "AMReX_TableData.H" + +#include "AMReX_TableData.H" +#include "QuokkaSimulation.hpp" +#include "hydro/hydro_system.hpp" +#include "math/quadrature.hpp" +#include "radiation/radiation_system.hpp" +#include "test_sne.hpp" + +// global variables needed for Dirichlet boundary condition and initial conditions +#if 0 // workaround AMDGPU compiler bug +namespace +{ +#endif // +AMREX_GPU_MANAGED amrex::GpuArray logphi_data{ + 0.83638381, 4.50705067, 5.10271383, 5.45268878, 5.70140736, 5.89447928, 6.05229308, 6.1857468, 6.30135334, 6.40331919, 6.49451684, 6.57699822, 6.65227803, + 6.72150811, 6.78558362, 6.84521466, 6.90097381, 6.95332936, 7.00266958, 7.04931937, 7.09355406, 7.1356083, 7.17568429, 7.21395706, 7.25057937, 7.28568514, + 7.31939257, 7.35180662, 7.38302073, 7.41311872, 7.44217572, 7.47025967, 7.49743193, 7.52374829, 7.54925951, 7.57401195, 7.5980481, 7.62140688, 7.64412421, + 7.66623309, 7.68776413, 7.70874559, 7.72920373, 7.74916295, 7.76864596, 7.78767394, 7.80626666, 7.82444266, 7.84221924, 7.85961268, 7.87663823, 7.89331028, + 7.90964238, 7.92564731, 7.94133718, 7.95672339, 7.9718168, 7.98662762, 8.00116561, 8.01543998, 8.02945951, 8.04323254, 8.056767, 8.07007045}; +AMREX_GPU_MANAGED amrex::GpuArray logg_data{ + -11.24587667, -10.68390187, -10.45405356, -10.31158848, -10.21082455, -10.13474898, -10.0750678, -10.02714739, -9.98806155, -9.95589857, -9.92931504, + -9.90722171, -9.88888792, -9.87375864, -9.86122421, -9.85088355, -9.84247276, -9.83565399, -9.83007266, -9.8256095, -9.82214665, -9.81939435, + -9.81722964, -9.81559821, -9.81433642, -9.8133616, -9.81265555, -9.81214602, -9.81177066, -9.81149643, -9.81129813, -9.81115635, -9.8110556, + -9.81098525, -9.81093618, -9.81090247, -9.81087945, -9.81086374, -9.81085319, -9.81084604, -9.8108412, -9.81083792, -9.81083566, -9.81083407, + -9.81083293, -9.81083206, -9.81083139, -9.81083083, -9.81083035, -9.81082991, -9.81082951, -9.81082913, -9.81082876, -9.8108284, -9.81082804, + -9.81082769, -9.81082735, -9.810827, -9.81082666, -9.81082632, -9.81082599, -9.81082565, -9.81082532, -9.81082499}; +AMREX_GPU_MANAGED amrex::GpuArray z_data{ + 1.50900000e+20, 5.55695238e+20, 9.60490476e+20, 1.36528571e+21, 1.77008095e+21, 2.17487619e+21, 2.57967143e+21, 2.98446667e+21, + 3.38926190e+21, 3.79405714e+21, 4.19885238e+21, 4.60364762e+21, 5.00844286e+21, 5.41323810e+21, 5.81803333e+21, 6.22282857e+21, + 6.62762381e+21, 7.03241905e+21, 7.43721429e+21, 7.84200952e+21, 8.24680476e+21, 8.65160000e+21, 9.05639524e+21, 9.46119048e+21, + 9.86598571e+21, 1.02707810e+22, 1.06755762e+22, 1.10803714e+22, 1.14851667e+22, 1.18899619e+22, 1.22947571e+22, 1.26995524e+22, + 1.31043476e+22, 1.35091429e+22, 1.39139381e+22, 1.43187333e+22, 1.47235286e+22, 1.51283238e+22, 1.55331190e+22, 1.59379143e+22, + 1.63427095e+22, 1.67475048e+22, 1.71523000e+22, 1.75570952e+22, 1.79618905e+22, 1.83666857e+22, 1.87714810e+22, 1.91762762e+22, + 1.95810714e+22, 1.99858667e+22, 2.03906619e+22, 2.07954571e+22, 2.12002524e+22, 2.16050476e+22, 2.20098429e+22, 2.24146381e+22, + 2.28194333e+22, 2.32242286e+22, 2.36290238e+22, 2.40338190e+22, 2.44386143e+22, 2.48434095e+22, 2.52482048e+22, 2.56530000e+22}; + +AMREX_GPU_MANAGED Real z_star = 245.0 * pc; +AMREX_GPU_MANAGED Real Sigma_star = 1.71 * Msun / pc / pc; +AMREX_GPU_MANAGED Real rho_dm = 1.4e-3 * Msun / pc / pc / pc; +AMREX_GPU_MANAGED Real R0 = 16.e3 * pc; +AMREX_GPU_MANAGED Real ks_sigma_sfr = 5.499927024044233e-57; +AMREX_GPU_MANAGED Real hscale = 300. * pc; +AMREX_GPU_MANAGED Real sigma1 = 11. * kmps; +AMREX_GPU_MANAGED Real sigma2 = 10. * 11. * kmps; +AMREX_GPU_MANAGED Real rho01 = 0.0268988 * Const_mH; +AMREX_GPU_MANAGED Real rho02 = 1.e-5 * 0.0268988 * Const_mH; +; +#if 0 // workaround AMDGPU compiler bug +}; // namespace +#endif + +struct NewProblem { +}; + +template <> struct HydroSystem_Traits { + static constexpr double gamma = 5. / 3.; + static constexpr bool reconstruct_eint = true; // Set to true - temperature +}; + +template <> struct quokka::EOS_Traits { + static constexpr double gamma = 5. / 3.; + static constexpr double mean_molecular_weight = C::m_u; + static constexpr double boltzmann_constant = C::k_B; +}; + +template <> struct Physics_Traits { + static constexpr bool is_hydro_enabled = true; + static constexpr bool is_radiation_enabled = false; + static constexpr bool is_chemistry_enabled = false; + static constexpr bool is_mhd_enabled = false; + static constexpr int numMassScalars = 0; // number of mass scalars + static constexpr int numPassiveScalars = 1; // number of passive scalars + static constexpr int nGroups = 1; // number of radiation groups +}; + +template <> struct SimulationData { + + std::unique_ptr> blast_x; + std::unique_ptr> blast_y; + std::unique_ptr> blast_z; + + int nblast = 0; + int SN_counter_cumulative = 0; + Real SN_rate_per_vol = NAN; // rate per unit time per unit volume + Real E_blast = 1.0e51; // ergs + Real M_ejecta = 5.0 * Msun; // 5.0 * Msun; // g + Real refine_threshold = 1.0; // gradient refinement threshold +}; + +template <> void QuokkaSimulation::setInitialConditionsOnGrid(quokka::grid const &grid_elem) +{ + + amrex::GpuArray dx = grid_elem.dx_; + amrex::GpuArray prob_lo = grid_elem.prob_lo_; + const amrex::Box &indexRange = grid_elem.indexRange_; + const amrex::Array4 &state_cc = grid_elem.array_; + + double vol = AMREX_D_TERM(dx[0], *dx[1], *dx[2]); + + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + amrex::Real const z = prob_lo[2] + (k + amrex::Real(0.5)) * dx[2]; + + // Calculate DM Potential + double prefac; + prefac = 2. * M_PI * Const_G * rho_dm * std::pow(R0, 2); + double Phidm = (prefac * std::log(1. + std::pow(z / R0, 2))); + + // Calculate Stellar Disk Potential + double prefac2; + prefac2 = 2. * M_PI * Const_G * Sigma_star * z_star; + double Phist = prefac2 * (std::pow(1. + z * z / z_star / z_star, 0.5) - 1.); + + // Calculate Gas Disk Potential + + double Phigas; + // Interpolate to find the accurate g-value from array-- because linterp doesn't work on Setonix + // TODO - AV to find out why linterp doesn't work + size_t ii = 0; + double x_interp = std::abs(z); + amrex::GpuArray xx = z_data; + amrex::GpuArray yy = logphi_data; + while (ii < xx.size() - 1 && x_interp > xx[ii + 1]) { + ii++; + } + + // Perform linear interpolation + const Real x1 = xx[ii]; + const Real x2 = xx[ii + 1]; + const Real y1 = yy[ii]; + const Real y2 = yy[ii + 1]; + amrex::Real phi_interp = (y1 + (y2 - y1) * (x_interp - x1) / (x2 - x1)); + Phigas = std::pow(10., phi_interp); + + double Phitot = Phist + Phidm + Phigas; + + double rho, rho_disk, rho_halo; + rho_disk = rho01 * std::exp(-Phitot / std::pow(sigma1, 2.0)); + rho_halo = rho02 * std::exp(-Phitot / std::pow(sigma2, 2.0)); // in g/cc + rho = (rho_disk + rho_halo); + + double P = rho_disk * std::pow(sigma1, 2.0) + rho_halo * std::pow(sigma2, 2.0); + + AMREX_ASSERT(!std::isnan(rho)); + + const auto gamma = HydroSystem::gamma_; + + state_cc(i, j, k, HydroSystem::density_index) = rho; + state_cc(i, j, k, HydroSystem::x1Momentum_index) = 0.0; + state_cc(i, j, k, HydroSystem::x2Momentum_index) = 0.0; + state_cc(i, j, k, HydroSystem::x3Momentum_index) = 0.0; + state_cc(i, j, k, HydroSystem::internalEnergy_index) = P / (gamma - 1.); + state_cc(i, j, k, HydroSystem::energy_index) = P / (gamma - 1.); + state_cc(i, j, k, Physics_Indices::pscalarFirstIndex) = 1.e-5 / vol; // Injected tracer + }); +} + +void AddSupernova(amrex::MultiFab &mf, amrex::GpuArray prob_lo, amrex::GpuArray prob_hi, + amrex::GpuArray dx, SimulationData const &userData, int level) +{ + // TODO for AV - ave (and restore) the RNG state in the metadata.yaml file + // inject energy into cells with stochastic sampling + BL_PROFILE("QuokkaSimulation::Addsupernova()") + + const Real cell_vol = AMREX_D_TERM(dx[0], *dx[1], *dx[2]); // cm^3 + const Real rho_eint_blast = userData.E_blast / cell_vol; // ergs cm^-3 + const Real rho_blast = userData.M_ejecta / cell_vol; // g cm^-3 + const Real scalar_blast = 1.e3 / cell_vol; // g cm^-3 + const int cum_sn = userData.SN_counter_cumulative; + + for (amrex::MFIter iter(mf); iter.isValid(); ++iter) { + const amrex::Box &box = iter.validbox(); + auto const &state = mf.array(iter); + auto const &px = userData.blast_x->table(); + auto const &py = userData.blast_y->table(); + auto const &pz = userData.blast_z->table(); + const int np = userData.nblast; + + amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const Real xc = prob_lo[0] + static_cast(i) * dx[0] + 0.5 * dx[0]; + const Real yc = prob_lo[1] + static_cast(j) * dx[1] + 0.5 * dx[1]; + const Real zc = prob_lo[2] + static_cast(k) * dx[2] + 0.5 * dx[2]; + + for (int n = 0; n < np; ++n) { + Real x0 = std::abs(xc - px(n)); + Real y0 = std::abs(yc - py(n)); + Real z0 = std::abs(zc - pz(n)); + + if (x0 < 0.5 * dx[0] && y0 < 0.5 * dx[1] && z0 < 0.5 * dx[2]) { + state(i, j, k, HydroSystem::density_index) += rho_blast; + state(i, j, k, HydroSystem::energy_index) += rho_eint_blast; + state(i, j, k, HydroSystem::internalEnergy_index) += rho_eint_blast; + state(i, j, k, Physics_Indices::pscalarFirstIndex) += scalar_blast; + + printf("The total number of SN gone off=%d\n", cum_sn); + Real Rpds = 14. * std::pow(state(i, j, k, HydroSystem::density_index) / Const_mH, -3. / 7.); + printf("Rpds = %.2e pc\n", Rpds); + } + } + }); + } +} + +template <> void QuokkaSimulation::computeBeforeTimestep() +{ + // compute how many (and where) SNe will go off on the this coarse timestep + // sample from Poisson distribution + + const Real dt_coarse = dt_[0]; + const Real domain_area = geom[0].ProbLength(0) * geom[0].ProbLength(1); + const Real mean = 0.0; + const Real stddev = hscale / geom[0].ProbLength(2) / 2.; + + const Real expectation_value = ks_sigma_sfr * domain_area * dt_coarse; + + const int count = static_cast(amrex::RandomPoisson(expectation_value)); + + // resize particle arrays + amrex::Array const lo{0}; + amrex::Array const hi{count}; + userData_.blast_x = std::make_unique>(lo, hi, amrex::The_Pinned_Arena()); + userData_.blast_y = std::make_unique>(lo, hi, amrex::The_Pinned_Arena()); + userData_.blast_z = std::make_unique>(lo, hi, amrex::The_Pinned_Arena()); + userData_.nblast = count; + userData_.SN_counter_cumulative += count; + + // for each, sample location at random + auto const &px = userData_.blast_x->table(); + auto const &py = userData_.blast_y->table(); + auto const &pz = userData_.blast_z->table(); + for (int i = 0; i < count; ++i) { + px(i) = geom[0].ProbLength(0) * amrex::Random(); + py(i) = geom[0].ProbLength(1) * amrex::Random(); + pz(i) = geom[0].ProbLength(2) * amrex::RandomNormal(mean, stddev); + } +} + +template <> void QuokkaSimulation::computeAfterLevelAdvance(int lev, amrex::Real time, amrex::Real dt_lev, int ncycle) +{ + amrex::GpuArray prob_lo = geom[lev].ProbLoArray(); + amrex::GpuArray prob_hi = geom[lev].ProbHiArray(); + amrex::GpuArray const &dx = geom[lev].CellSizeArray(); + + AddSupernova(state_new_cc_[lev], prob_lo, prob_hi, dx, userData_, lev); +} + +template <> +AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto +HydroSystem::GetGradFixedPotential(amrex::GpuArray posvec) -> amrex::GpuArray +{ + + amrex::GpuArray grad_potential; +#if (AMREX_SPACEDIM >= 3) + double z = posvec[2]; + + // Interpolate to find the accurate g-value from array-- because linterp doesn't work on Setonix + size_t i = 0; + double x_interp = std::abs(z); + amrex::GpuArray xx = z_data; + amrex::GpuArray yy = logg_data; + while (i < xx.size() - 1 && x_interp > xx[i + 1]) { + i++; + } + + // Perform linear interpolation + const Real x1 = xx[i]; + const Real x2 = xx[i + 1]; + const Real y1 = yy[i]; + const Real y2 = yy[i + 1]; + + amrex::Real ginterp = (y1 + (y2 - y1) * (x_interp - x1) / (x2 - x1)); + + grad_potential[2] = 2. * M_PI * Const_G * rho_dm * std::pow(R0, 2) * (2. * z / std::pow(R0, 2)) / (1. + std::pow(z, 2) / std::pow(R0, 2)); + grad_potential[2] += 2. * M_PI * Const_G * Sigma_star * (z / z_star) * (std::pow(1. + z * z / (z_star * z_star), -0.5)); + grad_potential[2] += (z / std::abs(z)) * std::pow(10., ginterp); + +#endif + + return grad_potential; +} + +// Add Strang Split Source Term for External Fixed Potential Here +template <> void QuokkaSimulation::addStrangSplitSources(amrex::MultiFab &mf, int lev, amrex::Real time, amrex::Real dt_lev) +{ + amrex::GpuArray prob_lo = geom[lev].ProbLoArray(); + amrex::GpuArray const &dx = geom[lev].CellSizeArray(); + const Real dt = dt_lev; + + for (amrex::MFIter iter(mf); iter.isValid(); ++iter) { + const amrex::Box &indexRange = iter.validbox(); + auto const &state = mf.array(iter); + + amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + amrex::GpuArray posvec, GradPhi; + double x1mom_new, x2mom_new, x3mom_new; + + const Real rho = state(i, j, k, HydroSystem::density_index); + const Real x1mom = state(i, j, k, HydroSystem::x1Momentum_index); + const Real x2mom = state(i, j, k, HydroSystem::x2Momentum_index); + const Real x3mom = state(i, j, k, HydroSystem::x3Momentum_index); + const Real Egas = state(i, j, k, HydroSystem::energy_index); + + Real Eint = RadSystem::ComputeEintFromEgas(rho, x1mom, x2mom, x3mom, Egas); + + posvec[0] = prob_lo[0] + (i + 0.5) * dx[0]; + +#if (AMREX_SPACEDIM >= 2) + posvec[1] = prob_lo[1] + (j + 0.5) * dx[1]; +#endif + +#if (AMREX_SPACEDIM >= 3) + posvec[2] = prob_lo[2] + (k + 0.5) * dx[2]; +#endif + + GradPhi = HydroSystem::GetGradFixedPotential(posvec); + + x1mom_new = state(i, j, k, x1mom) + dt * (-rho * GradPhi[0]); + x2mom_new = state(i, j, k, x2mom) + dt * (-rho * GradPhi[1]); + x3mom_new = state(i, j, k, x3mom) + dt * (-rho * GradPhi[2]); + + state(i, j, k, x1mom) = x1mom_new; + state(i, j, k, x2mom) = x2mom_new; + state(i, j, k, x3mom) = x3mom_new; + + state(i, j, k, HydroSystem::energy_index) = + RadSystem::ComputeEgasFromEint(rho, x1mom_new, x2mom_new, x3mom_new, Eint); + }); + } +} + +// Code for producing in-situ Projection plots +template <> auto QuokkaSimulation::ComputeProjections(const int dir) const -> std::unordered_map> +{ + // compute density projection + std::unordered_map> proj; + + proj["mass_outflow"] = computePlaneProjection( + [=] AMREX_GPU_DEVICE(int i, int j, int k, amrex::Array4 const &state) noexcept { + // int nmscalars = Physics_Traits::numMassScalars; + Real const rho = state(i, j, k, HydroSystem::density_index); + Real const vx3 = state(i, j, k, HydroSystem::x3Momentum_index) / rho; + return (rho * vx3); + }, + dir); + + proj["hot_mass_outflow"] = computePlaneProjection( + [=] AMREX_GPU_DEVICE(int i, int j, int k, amrex::Array4 const &state) noexcept { + double flux; + Real const rho = state(i, j, k, HydroSystem::density_index); + Real const vx3 = state(i, j, k, HydroSystem::x3Momentum_index) / rho; + Real const Eint = state(i, j, k, HydroSystem::internalEnergy_index); + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(state, i, j, k); + Real const primTemp = quokka::EOS::ComputeTgasFromEint(rho, Eint, massScalars); + if (primTemp > 1.e6) { + flux = rho * vx3; + } else { + flux = 0.0; + } + return flux; + }, + dir); + + proj["warm_mass_outflow"] = computePlaneProjection( + [=] AMREX_GPU_DEVICE(int i, int j, int k, amrex::Array4 const &state) noexcept { + double flux; + Real const rho = state(i, j, k, HydroSystem::density_index); + Real const vx3 = state(i, j, k, HydroSystem::x3Momentum_index) / rho; + Real const Eint = state(i, j, k, HydroSystem::internalEnergy_index); + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(state, i, j, k); + Real const primTemp = quokka::EOS::ComputeTgasFromEint(rho, Eint, massScalars); + if (primTemp < 2.e4) { + flux = rho * vx3; + } else { + flux = 0.0; + } + return flux; + }, + dir); + + proj["scalar_outflow"] = computePlaneProjection( + [=] AMREX_GPU_DEVICE(int i, int j, int k, amrex::Array4 const &state) noexcept { + Real const rho = state(i, j, k, HydroSystem::density_index); + Real const rhoZ = state(i, j, k, Physics_Indices::pscalarFirstIndex); + Real const vz = state(i, j, k, HydroSystem::x3Momentum_index) / rho; + return (rhoZ * vz); + }, + dir); + + proj["warm_scalar_outflow"] = computePlaneProjection( + [=] AMREX_GPU_DEVICE(int i, int j, int k, amrex::Array4 const &state) noexcept { + double flux; + Real const rho = state(i, j, k, HydroSystem::density_index); + Real const rhoZ = state(i, j, k, Physics_Indices::pscalarFirstIndex); + Real const vx3 = state(i, j, k, HydroSystem::x3Momentum_index) / rho; + Real const Eint = state(i, j, k, HydroSystem::internalEnergy_index); + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(state, i, j, k); + Real const primTemp = quokka::EOS::ComputeTgasFromEint(rho, Eint, massScalars); + if (primTemp < 2.e4) { + flux = rhoZ * vx3; + } else { + flux = 0.0; + } + return flux; + }, + dir); + + proj["hot_scalar_outflow"] = computePlaneProjection( + [=] AMREX_GPU_DEVICE(int i, int j, int k, amrex::Array4 const &state) noexcept { + double flux; + Real const rho = state(i, j, k, HydroSystem::density_index); + Real const rhoZ = state(i, j, k, Physics_Indices::pscalarFirstIndex); + Real const vx3 = state(i, j, k, HydroSystem::x3Momentum_index) / rho; + Real const Eint = state(i, j, k, HydroSystem::internalEnergy_index); + amrex::GpuArray massScalars = RadSystem::ComputeMassScalars(state, i, j, k); + Real const primTemp = quokka::EOS::ComputeTgasFromEint(rho, Eint, massScalars); + if (primTemp > 1.e6) { + flux = rhoZ * vx3; + } else { + flux = 0.0; + } + return flux; + }, + dir); + + proj["rho"] = computePlaneProjection( + [=] AMREX_GPU_DEVICE(int i, int j, int k, amrex::Array4 const &state) noexcept { + Real const rho = state(i, j, k, HydroSystem::density_index); + return (rho); + }, + dir); + + proj["scalar"] = computePlaneProjection( + [=] AMREX_GPU_DEVICE(int i, int j, int k, amrex::Array4 const &state) noexcept { + Real const rhoZ = state(i, j, k, Physics_Indices::pscalarFirstIndex); + return (rhoZ); + }, + dir); + return proj; +} + +// Implement User-defined diode BC +template <> +AMREX_GPU_DEVICE AMREX_FORCE_INLINE void AMRSimulation::setCustomBoundaryConditions(const amrex::IntVect &iv, amrex::Array4 const &consVar, + int /*dcomp*/, int /*numcomp*/, amrex::GeometryData const &geom, + const Real /*time*/, const amrex::BCRec * /*bcr*/, + int /*bcomp*/, int /*orig_comp*/) +{ + auto [i, j, k] = iv.dim3(); + amrex::Box const &box = geom.Domain(); + const auto &domain_lo = box.loVect3d(); + const auto &domain_hi = box.hiVect3d(); + const int klo = domain_lo[2]; + const int khi = domain_hi[2]; + int kedge = 0; + int normal = 0; + + if (k < klo) { + kedge = klo; + normal = -1; + } else if (k > khi) { + kedge = khi; + normal = 1.0; + } + + const double rho_edge = consVar(i, j, kedge, HydroSystem::density_index); + const double x1Mom_edge = consVar(i, j, kedge, HydroSystem::x1Momentum_index); + const double x2Mom_edge = consVar(i, j, kedge, HydroSystem::x2Momentum_index); + double x3Mom_edge = consVar(i, j, kedge, HydroSystem::x3Momentum_index); + const double etot_edge = consVar(i, j, kedge, HydroSystem::energy_index); + const double eint_edge = consVar(i, j, kedge, HydroSystem::internalEnergy_index); + + if ((x3Mom_edge * normal) < 0) { // gas is inflowing + x3Mom_edge = -1. * consVar(i, j, kedge, HydroSystem::x3Momentum_index); + } + + consVar(i, j, k, HydroSystem::density_index) = rho_edge; + consVar(i, j, k, HydroSystem::x1Momentum_index) = x1Mom_edge; + consVar(i, j, k, HydroSystem::x2Momentum_index) = x2Mom_edge; + consVar(i, j, k, HydroSystem::x3Momentum_index) = x3Mom_edge; + consVar(i, j, k, HydroSystem::energy_index) = etot_edge; + consVar(i, j, k, HydroSystem::internalEnergy_index) = eint_edge; +} + +auto problem_main() -> int +{ + + const int ncomp_cc = Physics_Indices::nvarTotal_cc; + amrex::Vector BCs_cc(ncomp_cc); + + for (int n = 0; n < ncomp_cc; ++n) { + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + // diode boundary conditions + if (i == 2) { + BCs_cc[n].setLo(i, amrex::BCType::ext_dir); + BCs_cc[n].setHi(i, amrex::BCType::ext_dir); + } else { + BCs_cc[n].setLo(i, amrex::BCType::int_dir); // periodic + BCs_cc[n].setHi(i, amrex::BCType::int_dir); // periodic + } + } + } + + // set random state + const int seed = 42; + amrex::InitRandom(seed, 1); // all ranks should produce the same values + + // Problem initialization + QuokkaSimulation sim(BCs_cc); + + sim.reconstructionOrder_ = 3; // 2=PLM, 3=PPM + sim.cflNumber_ = 0.3; // *must* be less than 1/3 in 3D! + + sim.setInitialConditions(); + + // evolve + sim.evolve(); + + // Cleanup and exit + amrex::Print() << "Finished." << std::endl; + return 0; +} diff --git a/src/problems/MetalAdvectionProblem/test_sne.hpp b/src/problems/MetalAdvectionProblem/test_sne.hpp new file mode 100644 index 000000000..8671594b9 --- /dev/null +++ b/src/problems/MetalAdvectionProblem/test_sne.hpp @@ -0,0 +1,34 @@ +#ifndef TEST_SNE_ // NOLINT +#define TEST_SNE_ +//============================================================================== +// TwoMomentRad - a radiation transport library for patch-based AMR codes +// Copyright 2020 Benjamin Wibking. +// Released under the MIT license. See LICENSE file included in the GitHub repo. +//============================================================================== +/// \file test_hydro3d_blast.hpp +/// \brief Defines a test problem for a 3D explosion. +/// + +// external headers +#include +#include + +// TODO: replace with Microphysics constant + +constexpr double Msun = 2.e33; +constexpr double Const_G = 6.67e-8; +constexpr double yr_to_s = 3.154e7; +constexpr double Myr = 1.e6 * yr_to_s; +constexpr double pc = 3.018e18; +constexpr double kpc = 1.e3 * pc; +constexpr double Mu = 0.6; +constexpr double kmps = 1.e5; +constexpr double Const_mH = 1.67e-24; +constexpr double kb = 1.3807e-16; +constexpr double sqrtpi = 1.772453; + +// internal headers +#include "hydro/hydro_system.hpp" +#include "math/interpolate.hpp" + +#endif // TEST_HYDRO3D_BLAST_HPP_ diff --git a/tests/metal_problem_regression.in b/tests/metal_problem_regression.in new file mode 100644 index 000000000..6542f47e5 --- /dev/null +++ b/tests/metal_problem_regression.in @@ -0,0 +1,38 @@ +# ***************************************************************** +# Problem size and geometry +# ***************************************************************** +geometry.prob_lo = 0.0 0.0 -2.4144e+22 +geometry.prob_hi = 3.018e21 3.018e21 2.4144e+22 +geometry.is_periodic = 1 1 0 + +# ***************************************************************** +# VERBOSITY +# ***************************************************************** +amr.v = 1 # verbosity in Amr + +# ***************************************************************** +# Resolution and refinement +# ***************************************************************** +amr.n_cell = 32 32 256 +amr.max_level = 0 # number of levels = max_level + 1 +amr.max_grid_size = 128 # at least 128 for GPUs +amr.blocking_factor = 32 # grid size must be divisible by this +amr.n_error_buf = 3 # minimum 3 cell buffer around tagged cells +amr.grid_eff = 0.7 # default + +do_reflux = 1 +do_subcycle = 1 +stop_time = 9.462e+15 +max_timesteps = 1000 +plotfile_interval = 500 +checkpoint_interval = 500 + +density_floor = 1.67e-30 #in g cm^-3 +temperature_floor = 10.0 #K + +cooling.enabled = 1 +cooling.cooling_table_type = "grackle" +cooling.hdf5_data_file = "CloudyData_UVB=HM2012.h5" + +projection_interval = 500 +projection.dirs = y