From 2c65efc908c230dff47388ad790ac4afb304deb7 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 5 Oct 2023 13:26:36 -0400 Subject: [PATCH] include pressure in x-momentum flux in 1-d Cartesian (#2468) previously, we were explicitly adding the pressure gradient in the conservative update to the momentum for all 1-d geometries. But we were only doing the flux correction on P_radial for non-Cartesian geometries. Now we treat all Cartesian geometries the same, regardless of dimensionality. This will likely change 1-d Cartesian answers. --- Docs/source/AMR.rst | 4 +- Docs/source/Hydrodynamics.rst | 10 +- .../ci-benchmarks/Rad2TShock-1d.out | 50 +++---- .../sdc_det_plt00040_extrema.out | 130 +++++++++--------- Source/driver/Castro_util.H | 10 +- Source/hydro/Castro_ctu_hydro.cpp | 12 +- Source/hydro/Castro_mol.cpp | 23 +--- Source/hydro/Castro_mol_hydro.cpp | 11 +- 8 files changed, 111 insertions(+), 139 deletions(-) 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/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 875cd724b2..79bff263e3 100644 --- a/Source/hydro/Castro_mol_hydro.cpp +++ b/Source/hydro/Castro_mol_hydro.cpp @@ -687,17 +687,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