From 5835ef88049b3625adb68790852be78b050bc6cc Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Wed, 11 Dec 2024 09:21:53 +1100 Subject: [PATCH] Add radiation-driven shell to regression tests (#810) ### Description Adds this problem to the nightly regression test suite. There was also a bug in the initial conditions in the radiation flux in converting from spherical to Cartesian flux components (now fixed). At late times, the simulation still has a slightly discrepancy w.r.t. the semi-analytic solution (similar to that shown in the original code paper). This is most likely due to initial shell thickness being quite significant. (There is not an easy solution to this, since there is no steady-state solution that satisfies the M1 closure for sufficiently thin shells.) ![shell_velocity](https://github.com/user-attachments/assets/0473f97c-3506-4d40-94c5-a9b84e2900f0) https://github.com/user-attachments/assets/8851d7cd-9d92-4de6-ac7e-ab5fef035af5 ### Related issues Closes https://github.com/quokka-astro/quokka/issues/743. ### Checklist _Before this pull request can be reviewed, all of these tasks should be completed. Denote completed tasks with an `x` inside the square brackets `[ ]` in the Markdown source below:_ - [x] I have added a description (see above). - [x] I have added a link to any related issues (if applicable; see above). - [x] I have read the [Contributing Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md). - [ ] I have added tests for any new physics that this PR adds to the code. - [x] *(For quokka-astro org members)* I have manually triggered the GPU tests with the magic comment `/azp run`. --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- extern/dust_shell/initial_conditions.txt | 82 ++++---- extern/dust_shell/matplotlibrc | 2 +- extern/dust_shell/solution.py | 7 +- regression/quokka-tests.ini | 16 ++ src/problems/RadhydroShell/CMakeLists.txt | 2 + .../RadhydroShell/test_radhydro_shell.cpp | 193 +++++------------- tests/radhydro_shell.in | 17 +- tests/radhydro_shell_regression.in | 33 +++ 8 files changed, 163 insertions(+), 189 deletions(-) create mode 100644 tests/radhydro_shell_regression.in diff --git a/extern/dust_shell/initial_conditions.txt b/extern/dust_shell/initial_conditions.txt index f3ace09f0..b4f4b7e37 100644 --- a/extern/dust_shell/initial_conditions.txt +++ b/extern/dust_shell/initial_conditions.txt @@ -1,32 +1,32 @@ # r_over_r0 reduced_flux Erad_cgs Frad_cgs -1.000000000000000048e-04 3.755101031908026650e-05 5.849260643223253650e-07 6.584810833497746207e-01 -1.441718602973596219e-03 5.413927440706735881e-04 5.849081618804471932e-07 9.493378914228637555e+00 -2.704299640581599162e-02 1.013562762624721503e-02 5.846128128060447862e-07 1.776395559784417060e+02 -5.264427420865838703e-02 1.961668425674432600e-02 5.840360439316570708e-07 3.434677422808198344e+02 -7.824555201150078243e-02 2.888029728914901992e-02 5.830889736094660357e-07 5.048439908936409779e+02 -1.038468298143431778e-01 3.783067358743209818e-02 5.817195410834095839e-07 6.597485279428594822e+02 -1.779956739478935068e-01 6.114121918945376183e-02 5.759656028151289446e-07 1.055726310902443856e+03 -2.521445180814438358e-01 7.899856349142272971e-02 5.696285693863173242e-07 1.349061225503343167e+03 -3.262933622149941648e-01 9.023184801610512062e-02 5.638115733303117799e-07 1.525156961640026111e+03 -4.062237502412902046e-01 9.486076267381476912e-02 5.588646660393681869e-07 1.589329583582680470e+03 -4.861541382675862444e-01 9.293464799427836320e-02 5.553646422860039258e-07 1.547307347592608494e+03 -5.660845262938822842e-01 8.651743901096156653e-02 5.525405597518683318e-07 1.433139683424766645e+03 -6.460149143201783239e-01 7.823584123424670578e-02 5.460810899423697447e-07 1.280806719571427266e+03 -6.807459332032566035e-01 7.481212301298670320e-02 5.394891460892043165e-07 1.209972204515023122e+03 -7.154769520863348831e-01 7.186691853201436497e-02 5.286291820213997449e-07 1.138940038957400475e+03 -7.502079709694131626e-01 6.968198530201698271e-02 5.117644526806205409e-07 1.069082781484670022e+03 -7.849389898524914422e-01 6.855134780844877873e-02 4.872960447401361441e-07 1.001450729546252319e+03 -8.196700087355697217e-01 6.879633527401403381e-02 4.542155802142597972e-07 9.368024854275187181e+02 -8.521778939234854189e-01 7.061262999568297671e-02 4.154382001870662751e-07 8.794466892024265690e+02 -8.846857791114011160e-01 7.438479746469706111e-02 3.701317120995621678e-07 8.253937649803749537e+02 -9.171936642993168132e-01 8.066352339760292367e-02 3.203808828403107721e-07 7.747551733282971327e+02 -9.497015494872325103e-01 9.020744304431384253e-02 2.690248161789397429e-07 7.275375597170433366e+02 -1.003030219610614449e+00 1.158923399574568469e-01 1.891775328524417536e-07 6.572717886984859206e+02 -1.024463928995228423e+00 1.311152943576374130e-01 1.606398288646572669e-07 6.314330214299012596e+02 -1.045897638379842398e+00 1.499304699665272855e-01 1.350219548090641344e-07 6.068970081536025418e+02 -1.067331347764456373e+00 1.729284538875688892e-01 1.125730104992235057e-07 5.836082760506781142e+02 -1.088765057149070348e+00 2.007054354675900376e-01 9.332050127159230878e-08 5.615092306215482267e+02 -1.118146872845352524e+00 2.474951844047469651e-01 7.184122515985754712e-08 5.330417009732285578e+02 +1.000000000000000048e-04 3.755101032144143522e-05 5.849260642855459238e-07 6.584810833497746207e-01 +1.441718602968267149e-03 5.413927440846509056e-04 5.849081618524834593e-07 9.493378914019864112e+00 +2.704299640580912212e-02 1.013562762624422263e-02 5.846128128060968788e-07 1.776395559784050988e+02 +5.264427420864997709e-02 1.961668425674107166e-02 5.840360439316672352e-07 3.434677422807688458e+02 +7.824555201149083206e-02 2.888029728914537353e-02 5.830889736094605300e-07 5.048439908935724816e+02 +1.038468298143316870e-01 3.783067358742812220e-02 5.817195410834083134e-07 6.597485279427886553e+02 +1.779956739478792960e-01 6.114121918944982748e-02 5.759656028151295798e-07 1.055726310902377008e+03 +2.521445180814269049e-01 7.899856349141945455e-02 5.696285693863187006e-07 1.349061225503290416e+03 +3.262933622149745139e-01 9.023184801610305283e-02 5.638115733303131563e-07 1.525156961639994961e+03 +4.062237502412764378e-01 9.486076267381457483e-02 5.588646660393688221e-07 1.589329583582679106e+03 +4.861541382675783618e-01 9.293464799427880729e-02 5.553646422860042435e-07 1.547307347592616679e+03 +5.660845262938802858e-01 8.651743901096177469e-02 5.525405597518682259e-07 1.433139683424769828e+03 +6.460149143201822097e-01 7.823584123424630332e-02 5.460810899423690035e-07 1.280806719571418853e+03 +6.807459332032604893e-01 7.481212301298635625e-02 5.394891460892032577e-07 1.209972204515015164e+03 +7.154769520863387688e-01 7.186691853201411517e-02 5.286291820213980508e-07 1.138940038957392744e+03 +7.502079709694170484e-01 6.968198530201683005e-02 5.117644526806179998e-07 1.069082781484662291e+03 +7.849389898524953280e-01 6.855134780844875098e-02 4.872960447401328619e-07 1.001450729546245043e+03 +8.196700087355736075e-01 6.879633527401417259e-02 4.542155802142554032e-07 9.368024854275116695e+02 +8.521778939234886385e-01 7.061262999568328202e-02 4.154382001870618282e-07 8.794466892024209983e+02 +8.846857791114036695e-01 7.438479746469750520e-02 3.701317120995581444e-07 8.253937649803708609e+02 +9.171936642993187006e-01 8.066352339760342327e-02 3.203808828403075957e-07 7.747551733282942905e+02 +9.497015494872337316e-01 9.020744304431432825e-02 2.690248161789377312e-07 7.275375597170417450e+02 +1.003030219610615559e+00 1.158923399574575547e-01 1.891775328524401654e-07 6.572717886984843290e+02 +1.024463928995229312e+00 1.311152943576380792e-01 1.606398288646561817e-07 6.314330214299002364e+02 +1.045897638379843064e+00 1.499304699665278962e-01 1.350219548090634198e-07 6.068970081536018597e+02 +1.067331347764456817e+00 1.729284538875693888e-01 1.125730104992230822e-07 5.836082760506776594e+02 +1.088765057149070570e+00 2.007054354675903429e-01 9.332050127159211025e-08 5.615092306215478857e+02 +1.118146872845352524e+00 2.474951844047470206e-01 7.184122515985753388e-08 5.330417009732285578e+02 1.147446050287565411e+00 3.040530980337316502e-01 5.557934245197780048e-08 5.066214110713187893e+02 1.176745227729778298e+00 3.688774428023082197e-01 4.358752202844885829e-08 4.820199144752725715e+02 1.206044405171991185e+00 4.385500681459409411e-01 3.491927786461010067e-08 4.590977238732691035e+02 @@ -58,14 +58,20 @@ 1.546576247924941638e+00 8.636240495720141075e-01 1.079436353580702643e-08 2.794746821838933784e+02 1.596404594890389994e+00 8.791200480310583654e-01 9.952498041173527408e-09 2.623016290403775201e+02 1.646232941855838350e+00 8.916142219547158465e-01 9.227997027531266418e-09 2.466636400124249064e+02 -1.734375178853998056e+00 9.091912467741695147e-01 8.153162862422567086e-09 2.222296828330790390e+02 -1.822517415852157763e+00 9.224110306819669480e-01 7.277792911337121628e-09 2.012541684514484359e+02 +1.734375178853998056e+00 9.091912467741694037e-01 8.153162862422568740e-09 2.222296828330790390e+02 +1.822517415852157763e+00 9.224110306819671701e-01 7.277792911337119974e-09 2.012541684514484359e+02 1.910659652850317469e+00 9.323283560477764187e-01 6.551368291360504820e-09 1.831140256663184118e+02 -1.998801889848477176e+00 9.400305847321852593e-01 5.937260295798599041e-09 1.673203545558839949e+02 -2.086944126846637104e+00 9.463952883490297507e-01 5.409702260886093429e-09 1.534852463072307103e+02 -2.272433441672774723e+00 9.563968807764520186e-01 4.514889045582422516e-09 1.294511568363405729e+02 -2.457922756498912342e+00 9.638255901339627396e-01 3.829417015942469982e-09 1.106501019861635768e+02 -2.643412071325049961e+00 9.695659269448310136e-01 3.291247133925901363e-09 9.566620401136817975e+01 -2.828901386151187580e+00 9.740199695318354500e-01 2.860645873254432982e-09 8.353195821790863818e+01 -3.195695155395610776e+00 9.804830097508788755e-01 2.226880004038710651e-09 6.545722516738128149e+01 -3.500000000000000000e+00 9.840687773348448797e-01 1.849720182475050407e-09 5.456977847985999830e+01 +1.998801889848477176e+00 9.400305847321853703e-01 5.937260295798599041e-09 1.673203545558839949e+02 +2.086944126846637104e+00 9.463952883490298618e-01 5.409702260886092602e-09 1.534852463072307103e+02 +2.272433441672772059e+00 9.563968807762568414e-01 4.514889045583355576e-09 1.294511568363408855e+02 +2.457922756498907013e+00 9.638255901347262400e-01 3.829417015939454081e-09 1.106501019861640600e+02 +2.643412071325041968e+00 9.695659269473928532e-01 3.291247133917225065e-09 9.566620401136876239e+01 +2.828901386151176922e+00 9.740199695356472898e-01 2.860645873243259426e-09 8.353195821790926345e+01 +3.195695155818024880e+00 9.804830097539137812e-01 2.226880003443110004e-09 6.545722515007672371e+01 +3.562488925484872837e+00 9.847259406867946430e-01 1.784206642841335117e-09 5.267217275611908178e+01 +3.929282695151720795e+00 9.876556558038823708e-01 1.462296600552641392e-09 4.329739107983610324e+01 +4.570182323792682411e+00 9.908450535180021168e-01 1.077444177193692332e-09 3.200525022874398928e+01 +5.211081952433644027e+00 9.929357294822037661e-01 8.269716695772742945e-10 2.461684964939163933e+01 +5.851981581074605643e+00 9.944685441777345591e-01 6.547426396019746781e-10 1.952011529386182431e+01 +6.492881209715567259e+00 9.956025175470230026e-01 5.312594591625017127e-10 1.585672027136650897e+01 +7.000000000000000000e+00 9.962532121817468944e-01 4.567744063764533480e-10 1.364244461996499957e+01 diff --git a/extern/dust_shell/matplotlibrc b/extern/dust_shell/matplotlibrc index dd7cd9c2c..1652be0d4 100644 --- a/extern/dust_shell/matplotlibrc +++ b/extern/dust_shell/matplotlibrc @@ -1,7 +1,7 @@ # Set the backend, otherwise the figure won't show up. Note that this will # depend on your system setup; to see which backend is the default, # run "matplotlib.get_backend()" in the Python interpreter. -backend : GTK3Agg +backend : Agg # Increase the default DPI, and change the file type from png to pdf savefig.dpi : 300 diff --git a/extern/dust_shell/solution.py b/extern/dust_shell/solution.py index 102e916d3..67edadb80 100644 --- a/extern/dust_shell/solution.py +++ b/extern/dust_shell/solution.py @@ -112,7 +112,7 @@ def func_M1(t, y): # solve for critical point def g(r): - return 3.0*r_0*dlnF_dr(r*r_0) + (4.0/r) + 2.0*sqrt(3.0)*r_0*rho(r*r_0)*kappa0 + return float(3.0*r_0*dlnF_dr(r*r_0) + (4.0/r) + 2.0*sqrt(3.0)*r_0*rho(r*r_0)*kappa0) #plt.figure() #r_arr = np.linspace(1.0e-3, 2.0, 100) @@ -139,7 +139,7 @@ def solve_branch(f0, r0, r1): # solve right branch r_start = r_crit + eps - r_end = 3.5 # in units of r_0 + r_end = 7.0 # in units of r_0 f_start = f_crit + eps * df_dr_crit print(f"r_crit = {r_crit}; r_start = {r_start}") print(f"f_crit = {f_crit}; f_start = {f_start}") @@ -166,7 +166,8 @@ def solve_branch(f0, r0, r1): plt.xlabel('radius r (dimensionless)') plt.ylabel('reduced flux f (dimensionless)') - plt.xlim(0., 3.5) + #plt.xlim(0., 3.5) + plt.xlim(0., 7.0) plt.ylim(0., 1.) plt.savefig('solution.pdf') diff --git a/regression/quokka-tests.ini b/regression/quokka-tests.ini index cdee09f6e..972336c46 100644 --- a/regression/quokka-tests.ini +++ b/regression/quokka-tests.ini @@ -123,3 +123,19 @@ doVis = 1 visVar = temperature testSrcTree = . ignore_return_code = 0 + +[RadhydroShell-GPU] +buildDir = . +target = test_radhydro3d_shell +inputFile = tests/radhydro_shell_regression.in +link1File = extern/dust_shell/initial_conditions.txt +dim = 3 +cmakeSetupOpts = -DAMReX_GPU_BACKEND=CUDA +restartTest = 0 +useMPI = 1 +numprocs = 1 +compileTest = 0 +doVis = 1 +visVar = gasDensity +testSrcTree = . +ignore_return_code = 0 diff --git a/src/problems/RadhydroShell/CMakeLists.txt b/src/problems/RadhydroShell/CMakeLists.txt index dcccaa3db..1dc0b7278 100644 --- a/src/problems/RadhydroShell/CMakeLists.txt +++ b/src/problems/RadhydroShell/CMakeLists.txt @@ -4,4 +4,6 @@ if (AMReX_SPACEDIM EQUAL 3) if(AMReX_GPU_BACKEND MATCHES "CUDA") setup_target_for_cuda_compilation(test_radhydro3d_shell) endif() + + add_test(NAME RadhydroShell COMMAND test_radhydro3d_shell radhydro_shell.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests) endif() diff --git a/src/problems/RadhydroShell/test_radhydro_shell.cpp b/src/problems/RadhydroShell/test_radhydro_shell.cpp index 5e60a3592..55281f759 100644 --- a/src/problems/RadhydroShell/test_radhydro_shell.cpp +++ b/src/problems/RadhydroShell/test_radhydro_shell.cpp @@ -7,34 +7,24 @@ /// \brief Defines a test problem for a 3D radiation pressure-driven shell. /// -#include - -#include "AMReX_Arena.H" #include "AMReX_Array.H" #include "AMReX_Array4.H" #include "AMReX_BC_TYPES.H" #include "AMReX_BLassert.H" -#include "AMReX_Config.H" -#include "AMReX_Extension.H" -#include "AMReX_FArrayBox.H" -#include "AMReX_FabArrayUtility.H" -#include "AMReX_Loop.H" -#include "AMReX_ParallelDescriptor.H" #include "AMReX_ParmParse.H" -#include "AMReX_Print.H" #include "AMReX_REAL.H" #include "AMReX_Vector.H" #include "QuokkaSimulation.hpp" +#include "SimulationData.hpp" #include "hydro/hydro_system.hpp" #include "math/interpolate.hpp" #include "radiation/radiation_system.hpp" -#include "test_radhydro_shell.hpp" struct ShellProblem { }; // if false, use octant symmetry -constexpr bool simulate_full_box = true; +constexpr bool simulate_full_box = false; constexpr double a_rad = C::a_rad; constexpr double c = C::c_light; @@ -86,9 +76,7 @@ constexpr amrex::Real H_shell = 0.3 * r_0; // cm constexpr amrex::Real kappa0 = 20.0; // specific opacity [cm^2 g^-1] constexpr amrex::Real rho_0 = M_shell / ((4. / 3.) * M_PI * r_0 * r_0 * r_0); // g cm^-3 - -constexpr amrex::Real P_0 = gamma_gas * rho_0 * (a0 * a0); // erg cm^-3 -constexpr double c_v = k_B / ((2.2 * m_H) * (gamma_gas - 1.0)); +constexpr amrex::Real c_v = k_B / ((2.2 * m_H) * (gamma_gas - 1.0)); template <> void RadSystem::SetRadEnergySource(array_t &radEnergy, const amrex::Box &indexRange, amrex::GpuArray const &dx, @@ -133,19 +121,24 @@ template <> AMREX_GPU_HOST_DEVICE auto RadSystem::ComputeFluxMeanO return ComputePlanckOpacity(0.0, 0.0); } -// declare global variables -// initial conditions read from file -amrex::Gpu::HostVector r_arr; -amrex::Gpu::HostVector Erad_arr; -amrex::Gpu::HostVector Frad_arr; +template <> struct SimulationData { + // initial conditions read from file + amrex::Gpu::HostVector r_arr; + amrex::Gpu::HostVector Erad_arr; + amrex::Gpu::HostVector Frad_arr; -amrex::Gpu::DeviceVector r_arr_g; -amrex::Gpu::DeviceVector Erad_arr_g; -amrex::Gpu::DeviceVector Frad_arr_g; + amrex::Gpu::DeviceVector r_arr_g; + amrex::Gpu::DeviceVector Erad_arr_g; + amrex::Gpu::DeviceVector Frad_arr_g; +}; template <> void QuokkaSimulation::preCalculateInitialConditions() { - std::string filename = "./initial_conditions.txt"; + // get filename from input file + amrex::ParmParse pp("shell_problem"); + std::string filename{}; + pp.query("filename", filename); + std::ifstream fstream(filename, std::ios::in); AMREX_ALWAYS_ASSERT(fstream.is_open()); std::string header; @@ -161,19 +154,19 @@ template <> void QuokkaSimulation::preCalculateInitialConditions() auto Erad = values.at(2); // cgs auto Frad = values.at(3); // cgs - r_arr.push_back(r); - Erad_arr.push_back(Erad); - Frad_arr.push_back(Frad); + userData_.r_arr.push_back(r); + userData_.Erad_arr.push_back(Erad); + userData_.Frad_arr.push_back(Frad); } // copy to device - r_arr_g.resize(r_arr.size()); - Erad_arr_g.resize(Erad_arr.size()); - Frad_arr_g.resize(Frad_arr.size()); + userData_.r_arr_g.resize(userData_.r_arr.size()); + userData_.Erad_arr_g.resize(userData_.Erad_arr.size()); + userData_.Frad_arr_g.resize(userData_.Frad_arr.size()); - amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, r_arr.begin(), r_arr.end(), r_arr_g.begin()); - amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, Erad_arr.begin(), Erad_arr.end(), Erad_arr_g.begin()); - amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, Frad_arr.begin(), Frad_arr.end(), Frad_arr_g.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, userData_.r_arr.begin(), userData_.r_arr.end(), userData_.r_arr_g.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, userData_.Erad_arr.begin(), userData_.Erad_arr.end(), userData_.Erad_arr_g.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, userData_.Frad_arr.begin(), userData_.Frad_arr.end(), userData_.Frad_arr_g.begin()); amrex::Gpu::streamSynchronizeAll(); } @@ -199,10 +192,10 @@ template <> void QuokkaSimulation::setInitialConditionsOnGrid(quok z0 = 0.; } - auto const &r_ptr = r_arr_g.dataPtr(); - auto const &Erad_ptr = Erad_arr_g.dataPtr(); - auto const &Frad_ptr = Frad_arr_g.dataPtr(); - int r_size = static_cast(r_arr_g.size()); + auto const &r_ptr = userData_.r_arr_g.dataPtr(); + auto const &Erad_ptr = userData_.Erad_arr_g.dataPtr(); + auto const &Frad_ptr = userData_.Frad_arr_g.dataPtr(); + int r_size = static_cast(userData_.r_arr_g.size()); // loop over the grid and set the initial condition amrex::ParallelFor(indexRange, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { @@ -210,6 +203,9 @@ template <> void QuokkaSimulation::setInitialConditionsOnGrid(quok amrex::Real const y = prob_lo[1] + (j + amrex::Real(0.5)) * dx[1]; amrex::Real const z = prob_lo[2] + (k + amrex::Real(0.5)) * dx[2]; amrex::Real const r = std::sqrt(std::pow(x - x0, 2) + std::pow(y - y0, 2) + std::pow(z - z0, 2)); + amrex::Real const rhat_x = (x - x0) / r; + amrex::Real const rhat_y = (y - y0) / r; + amrex::Real const rhat_z = (z - z0) / r; double sigma_sh = H_shell / (2.0 * std::sqrt(2.0 * std::log(2.0))); double rho_norm = M_shell / (4.0 * M_PI * r * r * std::sqrt(2.0 * M_PI * sigma_sh * sigma_sh)); @@ -217,7 +213,7 @@ template <> void QuokkaSimulation::setInitialConditionsOnGrid(quok double rho = std::max(rho_shell, 1.0e-8 * rho_0); // interpolate Frad from table - const double Frad = interpolate_value(r, r_ptr, Frad_ptr, r_size); + const double Frad_r = interpolate_value(r, r_ptr, Frad_ptr, r_size); // interpolate Erad from table const double Erad = interpolate_value(r, r_ptr, Erad_ptr, r_size); @@ -236,86 +232,18 @@ template <> void QuokkaSimulation::setInitialConditionsOnGrid(quok state_cc(i, j, k, HydroSystem::x3Momentum_index) = 0; state_cc(i, j, k, HydroSystem::energy_index) = Eint; - const double Frad_xyz = Frad / std::sqrt(3.0); + const double Frad_x = Frad_r * rhat_x; + const double Frad_y = Frad_r * rhat_y; + const double Frad_z = Frad_r * rhat_z; + state_cc(i, j, k, RadSystem::gasInternalEnergy_index) = Eint; state_cc(i, j, k, RadSystem::radEnergy_index) = Erad; - state_cc(i, j, k, RadSystem::x1RadFlux_index) = Frad_xyz; - state_cc(i, j, k, RadSystem::x2RadFlux_index) = Frad_xyz; - state_cc(i, j, k, RadSystem::x3RadFlux_index) = Frad_xyz; + state_cc(i, j, k, RadSystem::x1RadFlux_index) = Frad_x; + state_cc(i, j, k, RadSystem::x2RadFlux_index) = Frad_y; + state_cc(i, j, k, RadSystem::x3RadFlux_index) = Frad_z; }); } -AMREX_GPU_DEVICE AMREX_FORCE_INLINE auto vec_dot_r(amrex::GpuArray vec, int i, int j, int k) -> amrex::Real -{ - // compute dot product of vec into rhat - amrex::Real xhat = (i + amrex::Real(0.5)); - amrex::Real yhat = (j + amrex::Real(0.5)); - amrex::Real zhat = (k + amrex::Real(0.5)); - amrex::Real const norminv = 1.0 / std::sqrt(xhat * xhat + yhat * yhat + zhat * zhat); - - xhat *= norminv; - yhat *= norminv; - zhat *= norminv; - - amrex::Real const dotproduct = vec[0] * xhat + vec[1] * yhat + vec[2] * zhat; - return dotproduct; -} - -#if 0 -template <> void QuokkaSimulation::computeAfterTimestep() { - // compute radial momentum for gas, radiation on level 0 - // (assuming octant symmetry) - - amrex::GpuArray const &dx0 = - geom[0].CellSizeArray(); - amrex::Real const vol = AMREX_D_TERM(dx0[0], *dx0[1], *dx0[2]); - auto const &state = state_new_cc_[0]; - - double radialMom = - vol * - amrex::ReduceSum( - state, 0, - [=] AMREX_GPU_DEVICE(amrex::Box const &bx, - amrex::Array4 const &arr) { - amrex::Real result = 0.; - amrex::Loop(bx, [&](int i, int j, int k) { - amrex::GpuArray vec{ - arr(i, j, k, RadSystem::x1GasMomentum_index), - arr(i, j, k, RadSystem::x2GasMomentum_index), - arr(i, j, k, RadSystem::x3GasMomentum_index)}; - result += vec_dot_r(vec, i, j, k); - }); - return result; - }); - - amrex::ParallelAllReduce::Sum(radialMom, - amrex::ParallelContext::CommunicatorSub()); - - double radialRadMom = - (vol / c) * - amrex::ReduceSum( - state, 0, - [=] AMREX_GPU_DEVICE(amrex::Box const &bx, - amrex::Array4 const &arr) { - amrex::Real result = 0.; - amrex::Loop(bx, [&](int i, int j, int k) { - amrex::GpuArray vec{ - arr(i, j, k, RadSystem::x1RadFlux_index), - arr(i, j, k, RadSystem::x2RadFlux_index), - arr(i, j, k, RadSystem::x3RadFlux_index)}; - result += vec_dot_r(vec, i, j, k); - }); - return result; - }); - - amrex::ParallelAllReduce::Sum(radialRadMom, - amrex::ParallelContext::CommunicatorSub()); - - amrex::Print() << "radial gas momentum = " << radialMom << std::endl; - amrex::Print() << "radial radiation momentum = " << radialRadMom << std::endl; -} -#endif - template <> void QuokkaSimulation::ErrorEst(int lev, amrex::TagBoxArray &tags, amrex::Real /*time*/, int /*ngrow*/) { // tag cells for refinement @@ -387,9 +315,12 @@ auto problem_main() -> int for (int n = 0; n < ncomp_cc; ++n) { for (int i = 0; i < AMREX_SPACEDIM; ++i) { if constexpr (simulate_full_box) { - // periodic boundaries - BCs_cc[n].setLo(i, amrex::BCType::int_dir); - BCs_cc[n].setHi(i, amrex::BCType::int_dir); + // periodic boundaries (not recommended for this problem) + // BCs_cc[n].setLo(i, amrex::BCType::int_dir); + // BCs_cc[n].setHi(i, amrex::BCType::int_dir); + // outflow boundaries + BCs_cc[n].setLo(i, amrex::BCType::foextrap); + BCs_cc[n].setHi(i, amrex::BCType::foextrap); } else { // reflecting boundaries, outflow boundaries if (isNormalComp(n, i)) { @@ -405,37 +336,13 @@ auto problem_main() -> int // Problem initialization QuokkaSimulation sim(BCs_cc); - - sim.cflNumber_ = 0.3; - sim.densityFloor_ = 1.0e-8 * rho_0; - sim.pressureFloor_ = 1.0e-8 * P_0; - // reconstructionOrder: 1 == donor cell, 2 == PLM, 3 == PPM (not recommended - // for this problem) - sim.reconstructionOrder_ = 2; - sim.radiationReconstructionOrder_ = 2; - sim.integratorOrder_ = 2; // RK2 - constexpr amrex::Real t0_hydro = r_0 / a0; // seconds - sim.stopTime_ = 0.125 * t0_hydro; // 0.124 * t0_hydro; - - // for production - // sim.checkpointInterval_ = 1000; - // sim.plotfileInterval_ = 100; - // sim.maxTimesteps_ = 5000; - - // for scaling tests - sim.checkpointInterval_ = -1; - sim.plotfileInterval_ = -1; - sim.maxTimesteps_ = 50; - - // initialize + // sim.densityFloor_ = 1.0e-8 * rho_0; + sim.stopTime_ = 0.125 * t0_hydro; sim.setInitialConditions(); - sim.computeAfterTimestep(); - // evolve + // run sim.evolve(); - // Cleanup and exit - amrex::Print() << "Finished." << std::endl; return 0; } diff --git a/tests/radhydro_shell.in b/tests/radhydro_shell.in index 1405a7ce2..3457c8d79 100644 --- a/tests/radhydro_shell.in +++ b/tests/radhydro_shell.in @@ -1,19 +1,19 @@ # ***************************************************************** # Problem size and geometry [length in cm] # ***************************************************************** -geometry.prob_lo = 0.0 0.0 0.0 +geometry.prob_lo = 0. 0. 0. geometry.prob_hi = 3.086e19 3.086e19 3.086e19 -geometry.is_periodic = 1 1 1 +geometry.is_periodic = 0 0 0 # ***************************************************************** # VERBOSITY # ***************************************************************** -amr.v = 0 # verbosity in Amr +amr.v = 1 # verbosity in Amr # ***************************************************************** # Resolution and refinement # ***************************************************************** -amr.n_cell = 128 128 128 +amr.n_cell = 64 64 64 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 @@ -22,3 +22,12 @@ amr.grid_eff = 0.7 # default do_reflux = 1 do_subcycle = 1 + +cfl = 0.2 + +checkpoint_interval = -1 +plotfile_interval = 5000 +max_timesteps = 5000 + +shell_problem.filename = "../extern/dust_shell/initial_conditions.txt" +density_floor = 1.0e-28 # g cm^{-3} diff --git a/tests/radhydro_shell_regression.in b/tests/radhydro_shell_regression.in new file mode 100644 index 000000000..0fa9d45c2 --- /dev/null +++ b/tests/radhydro_shell_regression.in @@ -0,0 +1,33 @@ +# ***************************************************************** +# Problem size and geometry [length in cm] +# ***************************************************************** +geometry.prob_lo = 0. 0. 0. +geometry.prob_hi = 3.086e19 3.086e19 3.086e19 +geometry.is_periodic = 0 0 0 + +# ***************************************************************** +# VERBOSITY +# ***************************************************************** +amr.v = 1 # verbosity in Amr + +# ***************************************************************** +# Resolution and refinement +# ***************************************************************** +amr.n_cell = 128 128 128 +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 + +cfl = 0.2 + +checkpoint_interval = -1 +plotfile_interval = 5000 +max_timesteps = 5000 + +shell_problem.filename = "./initial_conditions.txt" +density_floor = 1.0e-28 # g cm^{-3}