diff --git a/Docs/source/AMR.rst b/Docs/source/AMR.rst index 7658fb2a57..19f64193d7 100644 --- a/Docs/source/AMR.rst +++ b/Docs/source/AMR.rst @@ -148,6 +148,8 @@ fine cell data above it.) The synchronization consists of two parts: +.. index:: castro.do_reflux + - Step 1: Hyperbolic reflux In the hyperbolic reflux step, we update the conserved variables with @@ -159,7 +161,7 @@ The synchronization consists of two parts: where :math:`V` is the volume of the cell and the correction from :math:`\delta\Fb` is supported only on coarse cells adjacent to fine grids. - Note: this can be enabled/disabled via castro.do_reflux. Generally, + Note: this can be enabled/disabled via ``castro.do_reflux``. Generally, it should be enabled (1). Also note that for axisymmetric or 1D spherical coordinates, the diff --git a/Docs/source/Hydrodynamics.rst b/Docs/source/Hydrodynamics.rst index 1f9dc0234f..8d76bd5fdf 100644 --- a/Docs/source/Hydrodynamics.rst +++ b/Docs/source/Hydrodynamics.rst @@ -1020,9 +1020,13 @@ parameters apply: - 1: the Colella & Glaz solver - - 2: the HLLC solver. Note: this should only be used with Cartesian - geometries because it relies on the pressure term being part of the flux - in the momentum equation. + - 2: the HLLC solver. + + .. note:: + + HLLC should only be used with Cartesian + geometries because it relies on the pressure term being part of the flux + in the momentum equation. The default is to use the solver based on an unpublished Colella, Glaz, & Ferguson manuscript (it also appears in :cite:`pember:1996`), diff --git a/Exec/radiation_tests/Rad2Tshock/ci-benchmarks/Rad2TShock-1d.out b/Exec/radiation_tests/Rad2Tshock/ci-benchmarks/Rad2TShock-1d.out index 239b56319f..4507245d57 100644 --- a/Exec/radiation_tests/Rad2Tshock/ci-benchmarks/Rad2TShock-1d.out +++ b/Exec/radiation_tests/Rad2Tshock/ci-benchmarks/Rad2TShock-1d.out @@ -1,40 +1,40 @@ plotfile = plt_00025 - time = 0.00026084061812361401 + time = 0.00026083835470476599 variables minimum value maximum value - density 5.4596904436e-13 1.2717694091e-12 - xmom 1.2328047524e-07 1.4591226172e-07 + density 5.4596904436e-13 1.2720098871e-12 + xmom 1.2327065163e-07 1.4587103381e-07 ymom 0 0 zmom 0 0 - rho_E 0.021940622285 0.039194728059 - rho_e 0.0068091599514 0.032333527949 - Temp 100.00003115 217.29660609 - rho_X 5.4596904436e-13 1.2717694091e-12 - rad 7.5662675907e-07 1.4083418723e-05 - pressure 0.0045394399678 0.0215556853 - kineng 0.0064473434191 0.015131462334 - soundspeed 117717.62829 173527.33922 + rho_E 0.021940622286 0.039199658391 + rho_e 0.0068091599527 0.032333529333 + Temp 100.00003116 217.3564598 + rho_X 5.4596904436e-13 1.2720098871e-12 + rad 7.5662679486e-07 1.4083425328e-05 + pressure 0.0045394399687 0.021555686223 + kineng 0.0064467244856 0.015131462333 + soundspeed 117717.6283 173551.23637 Gamma_1 1.6666666667 1.6666666667 - MachNumber 0.6068973935 1.9999997207 - uplusc 269518.17809 362247.33487 - uminusc -66804.446805 117717.59541 - entropy 2458725989.8 2507184972.2 + MachNumber 0.60689740202 1.9999997205 + uplusc 269483.97988 362284.15407 + uminusc -66810.704217 117717.5954 + entropy 2458725989.8 2507243428 magvort 0 0 - divu -20283.836079 245.94831837 - eint_E 12471696009 27100563708 - eint_e 12471696009 27100563708 - logden -12.26283198 -11.895591626 - StateErr_0 5.4596904436e-13 1.2717694091e-12 - StateErr_1 100.00003115 217.29660609 + divu -20286.061564 245.50876324 + eint_E 12471696011 27108028479 + eint_e 12471696011 27108028479 + logden -12.26283198 -11.895509513 + StateErr_0 5.4596904436e-13 1.2720098871e-12 + StateErr_1 100.00003116 217.3564598 StateErr_2 1 1 X(X) 1 1 abar 1 1 - x_velocity 102299.57948 235435.22371 + x_velocity 102258.95554 235435.2237 y_velocity 0 0 z_velocity 0 0 - magvel 102299.57948 235435.22371 - radvel -235435.22371 114042.93363 + magvel 102258.95554 235435.2237 + radvel -235435.2237 114038.34336 circvel 0 0.005524271728 - magmom 1.2328047524e-07 1.4591226172e-07 + magmom 1.2327065163e-07 1.4587103381e-07 angular_momentum_x 0 0 angular_momentum_y 0 0 angular_momentum_z 0 0 diff --git a/Exec/science/Detonation/ci-benchmarks/sdc_det_plt00040_extrema.out b/Exec/science/Detonation/ci-benchmarks/sdc_det_plt00040_extrema.out index 29337a5500..9fd16bf092 100644 --- a/Exec/science/Detonation/ci-benchmarks/sdc_det_plt00040_extrema.out +++ b/Exec/science/Detonation/ci-benchmarks/sdc_det_plt00040_extrema.out @@ -1,79 +1,79 @@ plotfile = det_x_plt00040 time = 5.1558159137010903e-06 variables minimum value maximum value - density 185259703.62 216632855.61 - xmom -73614997193 2.950511888e+16 + density 185260168 216581718.87 + xmom -65665876144 2.9498376337e+16 ymom 0 0 zmom 0 0 - rho_E 1.3062473822e+26 2.789305721e+26 - rho_e 1.3062473822e+26 2.7750368401e+26 - Temp 50000002.242 7845855083.7 - rho_H1 2.1206512613e-22 0.020000133523 - rho_He3 0.0017224949189 0.021023375462 - rho_He4 94357639.072 200001777.65 - rho_C12 0.020000000214 21232035.46 - rho_N14 1.9999997735e-22 1309.9957212 - rho_O16 0.019999999998 19209.484633 - rho_Ne20 0.019999999998 15403.34656 - rho_Mg24 0.019999999998 23316.245425 - rho_Si28 0.019999999998 2015224.2113 - rho_S32 0.019999999998 1656287.9261 - rho_Ar36 0.019999999998 821891.34664 - rho_Ca40 0.019999999998 725035.59799 - rho_Ti44 0.019999999998 34177.948134 - rho_Cr48 0.019999999998 78136.778748 - rho_Fe52 0.019999999998 278009.71311 - rho_Fe54 0.019999999998 95011242.573 - rho_Ni56 0.019999999998 2237420.8112 - rho_n 2.1206512613e-22 234580.45592 - rho_p 0.019999995436 3609780.2473 - rho_enuc -4.6854492557e+29 3.5779455282e+32 - pressure 5.5236728662e+25 1.1610551634e+26 - kineng 0 2.0451461945e+24 - soundspeed 612864631.28 895226108.54 - Gamma_1 1.3599756153 1.3820354139 - MachNumber 0 0.16107073838 - uplusc 612864631.28 999603676.63 - uminusc -895226447.4 -612860134.18 - entropy 98214769.614 336273476.62 + rho_E 1.3062473822e+26 2.7889587452e+26 + rho_e 1.3062473822e+26 2.7747032576e+26 + Temp 50000002.242 7845854636.1 + rho_H1 2.1207474382e-22 0.020000122856 + rho_He3 0.0017224302211 0.021022965925 + rho_He4 94383054.619 200001670.89 + rho_C12 0.020000000214 21356780.225 + rho_N14 1.9999989146e-22 0.020000168049 + rho_O16 0.019999999998 19209.480543 + rho_Ne20 0.019999999998 13162.840422 + rho_Mg24 0.019999999998 23319.285605 + rho_Si28 0.019999999998 2015271.217 + rho_S32 0.019999999998 1656229.9122 + rho_Ar36 0.019999999998 821818.64951 + rho_Ca40 0.019999999998 725327.27962 + rho_Ti44 0.019999999998 34171.893679 + rho_Cr48 0.019999999998 78189.409032 + rho_Fe52 0.019999999998 278027.80691 + rho_Fe54 0.019999999998 95004483.767 + rho_Ni56 0.019999999998 2238843.5743 + rho_n 2.1207474382e-22 234580.14578 + rho_p 0.019999995436 3609938.2338 + rho_enuc -4.6757680073e+29 3.5780501818e+32 + pressure 5.5236728662e+25 1.1610545851e+26 + kineng 0 2.0436963139e+24 + soundspeed 612864631.28 895226072.17 + Gamma_1 1.3599756231 1.3820776319 + MachNumber 0 0.16098334941 + uplusc 612864631.28 999588584.12 + uminusc -895226385.45 -612860404.21 + entropy 98214769.614 336273430.26 magvort 0 0 - divu -97913.435105 33932.952247 - eint_E 6.5312369112e+17 1.3804344965e+18 - eint_e 6.5312369112e+17 1.3804344965e+18 - logden 8.2677809649 8.3357243245 - StateErr_0 185259703.62 216632855.61 - StateErr_1 50000002.242 7845855083.7 + divu -97898.205698 33936.615428 + eint_E 6.5312369112e+17 1.3804343907e+18 + eint_e 6.5312369112e+17 1.3804343907e+18 + logden 8.2677820536 8.3356217961 + StateErr_0 185260168 216581718.87 + StateErr_1 50000002.242 7845854636.1 StateErr_2 1e-30 9.9999779314e-11 X(H1) 1e-30 9.9999779314e-11 - X(He3) 8.9854753938e-12 9.9999601244e-11 - X(He4) 0.48402991788 0.9999999982 - X(C12) 1.0000000107e-10 0.10012035382 - X(N14) 1e-30 6.1773274327e-06 - X(O16) 9.999999999e-11 9.6047421179e-05 - X(Ne20) 9.999999999e-11 7.2634981719e-05 - X(Mg24) 9.999999999e-11 0.00011116058012 - X(Si28) 9.999999999e-11 0.010055718145 - X(S32) 9.999999999e-11 0.0083470010273 - X(Ar36) 9.999999999e-11 0.004165653613 - X(Ca40) 9.999999999e-11 0.0037097977227 - X(Ti44) 9.999999999e-11 0.00017298351185 - X(Cr48) 9.999999999e-11 0.00040082116268 - X(Fe52) 9.999999999e-11 0.0013384275423 - X(Fe54) 9.999999999e-11 0.4649087176 - X(Ni56) 9.999999999e-11 0.010328169312 - X(n) 1e-30 0.0011729013039 - X(p) 9.9999977059e-11 0.017597003408 - abar 4.000000001 6.7314197717 - Ye 0.49998669277 0.50001557826 - x_velocity -368.07483814 138629924.04 + X(He3) 8.9861891466e-12 9.9999601244e-11 + X(He4) 0.48404046077 0.9999999982 + X(C12) 1.0000000107e-10 0.10070402463 + X(N14) 1e-30 1.0000000461e-10 + X(O16) 9.999999999e-11 9.6047399871e-05 + X(Ne20) 9.999999999e-11 6.2066987254e-05 + X(Mg24) 9.999999999e-11 0.00011117102392 + X(Si28) 9.999999999e-11 0.010061112286 + X(S32) 9.999999999e-11 0.0083507383612 + X(Ar36) 9.999999999e-11 0.0041671868568 + X(Ca40) 9.999999999e-11 0.0037113184679 + X(Ti44) 9.999999999e-11 0.00017306364204 + X(Cr48) 9.999999999e-11 0.00040099187008 + X(Fe52) 9.999999999e-11 0.0013389369622 + X(Fe54) 9.999999999e-11 0.46488111991 + X(Ni56) 9.999999999e-11 0.010337177053 + X(n) 1e-30 0.0011729003504 + X(p) 9.9999977036e-11 0.017600501144 + abar 4.000000001 6.7311470542 + Ye 0.49998666383 0.5000155877 + x_velocity -328.32930659 138563308.75 y_velocity 0 0 z_velocity 0 0 - t_sound_t_enuc 3.441191753e-13 0.97463512872 - enuc -2.4912109618e+21 1.6516172111e+24 - magvel 0 138629924.04 - radvel -368.07483814 138629924.04 + t_sound_t_enuc 3.4411947755e-13 0.97502907407 + enuc -2.4860875728e+21 1.6520554923e+24 + magvel 0 138563308.75 + radvel -328.32930659 138563308.75 circvel 0 2 - magmom 0 2.950511888e+16 + magmom 0 2.9498376337e+16 angular_momentum_x 0 0 angular_momentum_y 0 0 angular_momentum_z 0 0 diff --git a/Exec/unit_tests/diffusion_test/Prob.cpp b/Exec/unit_tests/diffusion_test/Prob.cpp index 5700979d86..f588c6f94b 100644 --- a/Exec/unit_tests/diffusion_test/Prob.cpp +++ b/Exec/unit_tests/diffusion_test/Prob.cpp @@ -21,8 +21,8 @@ void Castro::problem_post_simulation(Vector >& amr_lev Castro& castro = dynamic_cast(*amr_level[n]); Real time = castro.get_state_data(State_Type).curTime(); - const int* domain_lo = castro.geom.Domain().loVect(); - const int* domain_hi = castro.geom.Domain().hiVect(); + auto domain_lo = castro.geom.Domain().loVect3d(); + auto domain_hi = castro.geom.Domain().hiVect3d(); // the state data MultiFab& S = castro.get_new_data(State_Type); @@ -33,12 +33,18 @@ void Castro::problem_post_simulation(Vector >& amr_lev #ifdef TRUE_SDC // if we are fourth-order, we need to convert to averages if (sdc_order == 4) { - for (MFIter mfi(*analytic); mfi.isValid(); ++mfi) { + FArrayBox tmp; - const Box& gbx = mfi.growntilebox(1); - ca_make_fourth_in_place(AMREX_INT_ANYD(gbx.loVect()), AMREX_INT_ANYD(gbx.hiVect()), - BL_TO_FORTRAN_FAB((*analytic)[mfi]), - AMREX_INT_ANYD(domain_lo), AMREX_INT_ANYD(domain_hi)); + for (MFIter mfi(*analytic); mfi.isValid(); ++mfi) { + + const Box& bx = mfi.tilebox(); + tmp.resize(bx, 1); + auto tmp_arr = tmp.array(); + + castro.make_fourth_in_place(bx, + analytic->array(mfi), + tmp_arr, + domain_lo, domain_hi); } } diff --git a/Exec/unit_tests/diffusion_test/README.md b/Exec/unit_tests/diffusion_test/README.md index 5b57e4951b..c5e7e084d4 100644 --- a/Exec/unit_tests/diffusion_test/README.md +++ b/Exec/unit_tests/diffusion_test/README.md @@ -69,10 +69,10 @@ analytic solution, giving: A convergence test of the 4th-order SDC algorithm can be run as: ``` -./Castro1d.gnu.ex inputs.1d castro.time_integration_method=2 castro.sdc_order=4 amr.n_cell=64 -./Castro1d.gnu.ex inputs.1d castro.time_integration_method=2 castro.sdc_order=4 amr.n_cell=128 -./Castro1d.gnu.ex inputs.1d castro.time_integration_method=2 castro.sdc_order=4 amr.n_cell=256 -./Castro1d.gnu.ex inputs.1d castro.time_integration_method=2 castro.sdc_order=4 amr.n_cell=512 +./Castro1d.gnu.TRUESDC.ex inputs.1d castro.time_integration_method=2 castro.sdc_order=4 amr.n_cell=64 +./Castro1d.gnu.TRUESDC.ex inputs.1d castro.time_integration_method=2 castro.sdc_order=4 amr.n_cell=128 +./Castro1d.gnu.TRUESDC.ex inputs.1d castro.time_integration_method=2 castro.sdc_order=4 amr.n_cell=256 +./Castro1d.gnu.TRUESDC.ex inputs.1d castro.time_integration_method=2 castro.sdc_order=4 amr.n_cell=512 ``` Note: this is Cartesian, not spherical (since we don't have spherical diff --git a/Source/driver/Castro_util.H b/Source/driver/Castro_util.H index a2cb5b1b6b..91de1e146b 100644 --- a/Source/driver/Castro_util.H +++ b/Source/driver/Castro_util.H @@ -69,18 +69,14 @@ bool mom_flux_has_p (const int mom_dir, const int flux_dir, const int coord) if (mom_dir == 0 && flux_dir == 0) { - if (AMREX_SPACEDIM == 1 || (AMREX_SPACEDIM == 2 && coord == CoordSys::RZ)) { + if (coord != CoordSys::cartesian) { return false; - } - else { + } else { return true; } - } - else { - + } else { return (mom_dir == flux_dir); - } #else // AMREX_SPACEDIM == 3; Cartesian diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp index fe1a510c02..c159e83e43 100644 --- a/Source/hydro/Castro_ctu_hydro.cpp +++ b/Source/hydro/Castro_ctu_hydro.cpp @@ -40,7 +40,7 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) GeometryData geomdata = geom.data(); #endif -#if AMREX_SPACEDIM == 2 +#if AMREX_SPACEDIM <= 2 int coord = geom.Coord(); #endif @@ -1314,12 +1314,7 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) // get the scaled radial pressure -- we need to treat this specially #if AMREX_SPACEDIM <= 2 - -#if AMREX_SPACEDIM == 1 - if (!Geom().IsCartesian()) { -#elif AMREX_SPACEDIM == 2 if (!mom_flux_has_p(0, 0, coord)) { -#endif amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -1363,12 +1358,7 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) #endif #if AMREX_SPACEDIM <= 2 - -#if AMREX_SPACEDIM == 1 - if (idir == 0 && !Geom().IsCartesian()) { -#elif AMREX_SPACEDIM == 2 if (idir == 0 && !mom_flux_has_p(0, 0, coord)) { -#endif Array4 pradial_fab = pradial.array(); Array4 P_radial_fab = P_radial.array(mfi); diff --git a/Source/hydro/Castro_mol.cpp b/Source/hydro/Castro_mol.cpp index 907be84528..685a4e800f 100644 --- a/Source/hydro/Castro_mol.cpp +++ b/Source/hydro/Castro_mol.cpp @@ -234,9 +234,6 @@ Castro::mol_consup(const Box& bx, #if AMREX_SPACEDIM <= 2 const auto dx = geom.CellSizeArray(); -#endif - -#if AMREX_SPACEDIM == 2 auto coord = geom.Coord(); #endif @@ -257,22 +254,14 @@ Castro::mol_consup(const Box& bx, flux2(i,j,k,n) * area2(i,j,k) - flux2(i,j,k+1,n) * area2(i,j,k+1) ) / vol(i,j,k); #endif -#if AMREX_SPACEDIM == 1 - if (do_hydro == 1) { - if (n == UMX) { - update(i,j,k,UMX) -= (q0(i+1,j,k,GDPRES) - q0(i,j,k,GDPRES)) / dx[0]; - } - } -#endif +#if AMREX_SPACEDIM <= 2 + if (n == UMX && do_hydro == 1) { + // Add gradp term to momentum equation -- only for axisymmetric + // coords (and only for the radial flux). -#if AMREX_SPACEDIM == 2 - if (do_hydro == 1) { - if (n == UMX) { - // add the pressure source term for axisymmetry - if (coord > 0) { - update(i,j,k,n) -= (q0(i+1,j,k,GDPRES) - q0(i,j,k,GDPRES)) / dx[0]; + if (!mom_flux_has_p(0, 0, coord)) { + update(i,j,k,UMX) -= (q0(i+1,j,k,GDPRES) - q0(i,j,k,GDPRES)) / dx[0]; } - } } #endif diff --git a/Source/hydro/Castro_mol_hydro.cpp b/Source/hydro/Castro_mol_hydro.cpp index 2a7e4e0370..6422ae8c64 100644 --- a/Source/hydro/Castro_mol_hydro.cpp +++ b/Source/hydro/Castro_mol_hydro.cpp @@ -693,17 +693,8 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) Array4 const qex_fab = qe[idir].array(); const int prescomp = GDPRES; -#if AMREX_SPACEDIM == 1 - if (!Geom().IsCartesian()) { - amrex::ParallelFor(nbx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - pradial_fab(i,j,k) = qex_fab(i,j,k,prescomp) * dt; - }); - } -#endif -#if AMREX_SPACEDIM == 2 +#if AMREX_SPACEDIM <= 2 if (!mom_flux_has_p(0, 0, coord)) { amrex::ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept diff --git a/Source/hydro/fourth_order.cpp b/Source/hydro/fourth_order.cpp index 44c721d4f1..6f134967bc 100644 --- a/Source/hydro/fourth_order.cpp +++ b/Source/hydro/fourth_order.cpp @@ -883,40 +883,40 @@ Castro::fourth_add_diffusive_flux(const Box& bx, if (idir == 0) { if (is_avg) { - // we are working with the cell-center state - dTdx = (-q_arr(i+1,j,k,temp_comp) + 27*q_arr(i,j,k,temp_comp) - - 27*q_arr(i-1,j,k,temp_comp) + q_arr(i-2,j,k,temp_comp)) / (24.0_rt * dx[0]); - - } else { // we are working with the cell-average state dTdx = (-q_arr(i+1,j,k,temp_comp) + 15*q_arr(i,j,k,temp_comp) - 15*q_arr(i-1,j,k,temp_comp) + q_arr(i-2,j,k,temp_comp)) / (12.0_rt * dx[0]); + + } else { + // we are working with the cell-center state + dTdx = (-q_arr(i+1,j,k,temp_comp) + 27*q_arr(i,j,k,temp_comp) - + 27*q_arr(i-1,j,k,temp_comp) + q_arr(i-2,j,k,temp_comp)) / (24.0_rt * dx[0]); } } else if (idir == 1) { if (is_avg) { - // we are working with the cell-center state - dTdx = (-q_arr(i,j+1,k,temp_comp) + 27*q_arr(i,j,k,temp_comp) - - 27*q_arr(i,j-1,k,temp_comp) + q_arr(i,j-2,k,temp_comp)) / (24.0_rt * dx[1]); - - } else { // we are working with the cell-average state dTdx = (-q_arr(i,j+1,k,temp_comp) + 15*q_arr(i,j,k,temp_comp) - 15*q_arr(i,j-1,k,temp_comp) + q_arr(i,j-2,k,temp_comp)) / (12.0_rt * dx[1]); + + } else { + // we are working with the cell-center state + dTdx = (-q_arr(i,j+1,k,temp_comp) + 27*q_arr(i,j,k,temp_comp) - + 27*q_arr(i,j-1,k,temp_comp) + q_arr(i,j-2,k,temp_comp)) / (24.0_rt * dx[1]); } } else { if (is_avg) { - // we are working with the cell-center state - dTdx = (-q_arr(i,j,k+1,temp_comp) + 27*q_arr(i,j,k,temp_comp) - - 27*q_arr(i,j,k-1,temp_comp) + q_arr(i,j,k-2,temp_comp)) / (24.0_rt * dx[2]); - - } else { // we are working with the cell-average state dTdx = (-q_arr(i,j,k+1,temp_comp) + 15*q_arr(i,j,k,temp_comp) - 15*q_arr(i,j,k-1,temp_comp) + q_arr(i,j,k-2,temp_comp)) / (12.0_rt * dx[2]); + + } else { + // we are working with the cell-center state + dTdx = (-q_arr(i,j,k+1,temp_comp) + 27*q_arr(i,j,k,temp_comp) - + 27*q_arr(i,j,k-1,temp_comp) + q_arr(i,j,k-2,temp_comp)) / (24.0_rt * dx[2]); } }