diff --git a/.clang-tidy b/.clang-tidy index 36809183ac..391149e1da 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -7,7 +7,6 @@ Checks: > clang-diagnostic-*, cppcoreguidelines-*, -cppcoreguidelines-avoid-c-arrays, - -cppcoreguidelines-avoid-goto, -cppcoreguidelines-avoid-magic-numbers, -cppcoreguidelines-avoid-non-const-global-variables, -cppcoreguidelines-init-variables, @@ -17,11 +16,17 @@ Checks: > -cppcoreguidelines-non-private-member-variables-in-classes, -cppcoreguidelines-owning-memory, -cppcoreguidelines-pro-*, + misc-*, + -misc-const-correctness, + -misc-include-cleaner, + -misc-non-private-member-variables-in-classes, modernize-*, -modernize-avoid-c-arrays, -modernize-use-trailing-return-type, -modernize-use-using, performance-*, + -performance-avoid-endl, + portability-*, readability-*, -readability-avoid-const-params-in-decls, -readability-braces-around-statements, diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 0000000000..60a11755c8 --- /dev/null +++ b/.github/dependabot.yml @@ -0,0 +1,9 @@ +# Dependabot configuration +# ref: https://docs.github.com/en/code-security/supply-chain-security/keeping-your-dependencies-updated-automatically/configuration-options-for-dependency-updates +version: 2 +updates: + - package-ecosystem: "github-actions" + directory: "/" + schedule: + interval: "weekly" + target-branch: "development" diff --git a/.github/workflows/c-linter.yml b/.github/workflows/c-linter.yml index 435a7fa60a..d575a2f236 100644 --- a/.github/workflows/c-linter.yml +++ b/.github/workflows/c-linter.yml @@ -30,7 +30,7 @@ jobs: tar xfz v2.28.0.tar.gz cd hypre-2.28.0/src ./configure --with-cxxstandard=17 --without-MPI - make -j 2 + make -j 4 make install cd ../../ diff --git a/.github/workflows/check-ifdefs.yml b/.github/workflows/check-ifdefs.yml index 3f4ab8daae..9afa5f0167 100644 --- a/.github/workflows/check-ifdefs.yml +++ b/.github/workflows/check-ifdefs.yml @@ -14,7 +14,7 @@ jobs: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 with: fetch-depth: 0 diff --git a/.github/workflows/codespell.yml b/.github/workflows/codespell.yml index 6964398e22..d1454aceb2 100644 --- a/.github/workflows/codespell.yml +++ b/.github/workflows/codespell.yml @@ -19,7 +19,7 @@ jobs: fetch-depth: 0 - name: Setup Python - uses: actions/setup-python@v4 + uses: actions/setup-python@v5 with: python-version: '3.10' cache: "pip" diff --git a/.github/workflows/compiler-warnings.yml b/.github/workflows/compiler-warnings.yml index 0f6bd12fb1..e8092be3d3 100644 --- a/.github/workflows/compiler-warnings.yml +++ b/.github/workflows/compiler-warnings.yml @@ -26,14 +26,14 @@ jobs: - name: Compile Detonation run: | cd Exec/science/Detonation - make USE_MPI=FALSE USE_OMP=FALSE DEBUG=TRUE WARN_ALL=TRUE WARN_ERROR=TRUE -j 2 + make USE_MPI=FALSE USE_OMP=FALSE DEBUG=TRUE WARN_ALL=TRUE WARN_ERROR=TRUE -j 4 - name: Compile subchandra run: | cd Exec/science/subchandra - make USE_MPI=FALSE USE_OMP=FALSE DEBUG=TRUE WARN_ALL=TRUE WARN_ERROR=TRUE -j 2 + make USE_MPI=FALSE USE_OMP=FALSE DEBUG=TRUE WARN_ALL=TRUE WARN_ERROR=TRUE -j 4 - name: Compile wdmerger run: | cd Exec/science/wdmerger - make USE_MPI=FALSE USE_OMP=FALSE DEBUG=TRUE WARN_ALL=TRUE WARN_ERROR=TRUE -j 2 + make USE_MPI=FALSE USE_OMP=FALSE DEBUG=TRUE WARN_ALL=TRUE WARN_ERROR=TRUE -j 4 diff --git a/.github/workflows/detonation-sdc-compare.yml b/.github/workflows/detonation-sdc-compare.yml index be850b35cf..b8fe7e8568 100644 --- a/.github/workflows/detonation-sdc-compare.yml +++ b/.github/workflows/detonation-sdc-compare.yml @@ -26,7 +26,7 @@ jobs: - name: Compile Detonation run: | cd Exec/science/Detonation - make DIM=1 USE_MPI=FALSE -j 2 USE_SIMPLIFIED_SDC=TRUE + make DIM=1 USE_MPI=FALSE -j 4 USE_SIMPLIFIED_SDC=TRUE - name: Run Detonation run: | @@ -36,7 +36,7 @@ jobs: - name: Build the fextrema tool run: | cd external/amrex/Tools/Plotfile - make programs=fextrema -j 2 + make programs=fextrema -j 4 - name: Check the extrema run: | diff --git a/.github/workflows/good_defines.txt b/.github/workflows/good_defines.txt index 42cbeab4ad..c0f6c73287 100644 --- a/.github/workflows/good_defines.txt +++ b/.github/workflows/good_defines.txt @@ -3,6 +3,7 @@ AMREX_PARTICLES AMREX_SPACEDIM AMREX_USE_CUDA AMREX_USE_GPU +AMREX_USE_HIP AMREX_USE_OMP AUX_THERMO BL_FORT_USE_LOWERCASE diff --git a/.github/workflows/gpu_action.yml b/.github/workflows/gpu_action.yml index 6861161cbb..95d65496b0 100644 --- a/.github/workflows/gpu_action.yml +++ b/.github/workflows/gpu_action.yml @@ -39,7 +39,7 @@ jobs: tar xfz v2.28.0.tar.gz cd hypre-2.28.0/src CUDA_HOME=/usr/local/cuda HYPRE_CUDA_SM=60 ./configure --with-cxxstandard=17 --with-cuda --enable-unified-memory --without-MPI - make -j 2 + make -j 4 make install cd ../../ diff --git a/.github/workflows/hip.yml b/.github/workflows/hip.yml index f74b1e7c7e..1bd1aa4b3e 100644 --- a/.github/workflows/hip.yml +++ b/.github/workflows/hip.yml @@ -29,4 +29,4 @@ jobs: - name: compile flame_wave run: | cd Exec/science/flame_wave - make COMP=gnu USE_HIP=TRUE USE_MPI=FALSE USE_OMP=FALSE USE_CUDA=FALSE -j 2 + make COMP=gnu USE_HIP=TRUE USE_MPI=FALSE USE_OMP=FALSE USE_CUDA=FALSE -j 4 diff --git a/.github/workflows/mhd-compare.yml b/.github/workflows/mhd-compare.yml index 403ff79f5e..e5655ce236 100644 --- a/.github/workflows/mhd-compare.yml +++ b/.github/workflows/mhd-compare.yml @@ -26,7 +26,7 @@ jobs: - name: Compile OrszagTang run: | cd Exec/mhd_tests/OrszagTang - make USE_MPI=TRUE -j 2 + make USE_MPI=TRUE -j 4 - name: Run OrszagTang-3d run: | @@ -36,7 +36,7 @@ jobs: - name: Build the fextrema tool run: | cd external/amrex/Tools/Plotfile - make programs=fextrema -j 2 + make programs=fextrema -j 4 - name: Check the extrema run: | diff --git a/.github/workflows/rad-compare.yml b/.github/workflows/rad-compare.yml index 775024000d..8645eb9ed0 100644 --- a/.github/workflows/rad-compare.yml +++ b/.github/workflows/rad-compare.yml @@ -29,7 +29,7 @@ jobs: tar xfz v2.26.0.tar.gz cd hypre-2.26.0/src ./configure --with-cxxstandard=17 - make -j 2 + make -j 4 make install cd ../../ @@ -47,7 +47,7 @@ jobs: - name: Build the fextrema tool run: | cd external/amrex/Tools/Plotfile - make programs=fextrema -j 2 + make programs=fextrema -j 4 - name: Check the extrema run: | diff --git a/.github/workflows/reacting-convergence-true-sdc.yml b/.github/workflows/reacting-convergence-true-sdc.yml index bd810420f8..fbc3d48cb2 100644 --- a/.github/workflows/reacting-convergence-true-sdc.yml +++ b/.github/workflows/reacting-convergence-true-sdc.yml @@ -26,7 +26,7 @@ jobs: - name: Compile reacting_convergence run: | cd Exec/reacting_tests/reacting_convergence - make DIM=2 USE_MPI=FALSE -j 2 USE_TRUE_SDC=TRUE + make DIM=2 USE_MPI=FALSE -j 4 USE_TRUE_SDC=TRUE - name: Run reacting_convergence run: | @@ -36,7 +36,7 @@ jobs: - name: Build the fextrema tool run: | cd external/amrex/Tools/Plotfile - make programs=fextrema -j 2 + make programs=fextrema -j 4 - name: Check the extrema run: | diff --git a/.github/workflows/uniform_cube.yml b/.github/workflows/uniform_cube.yml index edff5a8494..952d133213 100644 --- a/.github/workflows/uniform_cube.yml +++ b/.github/workflows/uniform_cube.yml @@ -26,7 +26,7 @@ jobs: - name: Compile uniform_cube run: | cd Exec/gravity_tests/uniform_cube - make USE_MPI=FALSE -j 2 + make USE_MPI=FALSE -j 4 - name: Run uniform_cube run: | diff --git a/.github/workflows/uniform_sphere.yml b/.github/workflows/uniform_sphere.yml index e9ad3b9e4e..d8d066c34a 100644 --- a/.github/workflows/uniform_sphere.yml +++ b/.github/workflows/uniform_sphere.yml @@ -26,7 +26,7 @@ jobs: - name: Compile uniform_sphere run: | cd Exec/gravity_tests/uniform_sphere - make USE_MPI=FALSE -j 2 + make USE_MPI=FALSE -j 4 - name: Run uniform_sphere run: | diff --git a/CHANGES.md b/CHANGES.md index 0bb87dc598..12be083506 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -1,3 +1,15 @@ +# 24.02 + + * Lot's of code fixes from coverity and clang-tidy (#2736, #2734, + #2735, #2731, #2732, #2733) + + * Fix the boundary condition logic at a wall for Detonation (#2722) + + * Reimplement the shock detection algorithm to account for sources + and do a better job in multidimensions (#2711, #2710, #2709, #2704) + + * Start the process of moving the runtime parameters to structs (#2688) + # 24.01 * An option for unlimited PPM reconstruction was added (#2670) diff --git a/Diagnostics/DustCollapse/main.cpp b/Diagnostics/DustCollapse/main.cpp index 7e2a0763de..6fad04b96b 100644 --- a/Diagnostics/DustCollapse/main.cpp +++ b/Diagnostics/DustCollapse/main.cpp @@ -102,18 +102,18 @@ int main(int argc, char* argv[]) #if (AMREX_SPACEDIM == 1) int nbins = domain.length(0); #elif (AMREX_SPACEDIM == 2) - auto x_maxdist = fmax(fabs(problo[0]), fabs(probhi[0])); - auto y_maxdist = fmax(fabs(problo[1] - yctr), fabs(probhi[1] - yctr)); + auto x_maxdist = std::max(std::abs(problo[0]), std::abs(probhi[0])); + auto y_maxdist = std::max(std::abs(problo[1] - yctr), std::abs(probhi[1] - yctr)); - auto max_dist = sqrt(x_maxdist*x_maxdist + y_maxdist*y_maxdist); + auto max_dist = std::sqrt(x_maxdist*x_maxdist + y_maxdist*y_maxdist); int nbins = int(max_dist / dx_fine); #else - auto x_maxdist = fmax(fabs(problo[0] - xctr), fabs(probhi[0] - xctr)); - auto y_maxdist = fmax(fabs(problo[1] - yctr), fabs(probhi[1] - yctr)); - auto z_maxdist = fmax(fabs(problo[2] - zctr), fabs(probhi[2] - zctr)); + auto x_maxdist = std::max(std::abs(problo[0] - xctr), std::abs(probhi[0] - xctr)); + auto y_maxdist = std::max(std::abs(problo[1] - yctr), std::abs(probhi[1] - yctr)); + auto z_maxdist = std::max(std::abs(problo[2] - zctr), std::abs(probhi[2] - zctr)); - auto max_dist = sqrt(x_maxdist*x_maxdist + y_maxdist*y_maxdist + + auto max_dist = std::sqrt(x_maxdist*x_maxdist + y_maxdist*y_maxdist + z_maxdist*z_maxdist); int nbins = int(max_dist / dx_fine); @@ -147,7 +147,7 @@ int main(int argc, char* argv[]) // over levels, we will compare to the finest level index space to // determine if we've already output here int mask_size = domain.length().max(); - Vector imask(pow(mask_size, AMREX_SPACEDIM), 1); + Vector imask(std::pow(mask_size, AMREX_SPACEDIM), 1); // counter int cnt = 0; @@ -172,17 +172,17 @@ int main(int argc, char* argv[]) const Box& bx = mfi.tilebox(); #if (AMREX_SPACEDIM == 1) - fdustcollapse1d(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), + fdustcollapse1d(AMREX_ARLIM_3D(bx.loVect()), AMREX_ARLIM_3D(bx.hiVect()), BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), nbins, dens.dataPtr(), imask.dataPtr(), mask_size, r1, dens_comp, &cnt); #elif (AMREX_SPACEDIM == 2) - fdustcollapse2d(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), + fdustcollapse2d(AMREX_ARLIM_3D(bx.loVect()), AMREX_ARLIM_3D(bx.hiVect()), BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), nbins, dens.dataPtr(), volcount.dataPtr(), imask.dataPtr(), mask_size, r1, - ZFILL(level_dx), dx_fine, yctr, dens_comp); + AMREX_ZFILL(level_dx), dx_fine, yctr, dens_comp); #else fdustcollapse3d(bx.loVect(), bx.hiVect(), BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), @@ -210,21 +210,21 @@ int main(int argc, char* argv[]) // analytic expression for the radius as a function of time t = 0.00 Real max_dens = 1.e9; - if (fabs(data.Time()) <= 1.e-8) + if (std::abs(data.Time()) <= 1.e-8) max_dens = 1.e9; - else if (fabs(data.Time() - 0.01) <= 1.e-8) + else if (std::abs(data.Time() - 0.01) <= 1.e-8) max_dens = 1.043345e9; - else if (fabs(data.Time() - 0.02) <= 1.e-8) + else if (std::abs(data.Time() - 0.02) <= 1.e-8) max_dens = 1.192524e9; - else if (fabs(data.Time() - 0.03) <= 1.e-8) + else if (std::abs(data.Time() - 0.03) <= 1.e-8) max_dens = 1.527201e9; - else if (fabs(data.Time() - 0.04) <= 1.e-8) + else if (std::abs(data.Time() - 0.04) <= 1.e-8) max_dens = 2.312884e9; - else if (fabs(data.Time() - 0.05) <= 1.e-8) + else if (std::abs(data.Time() - 0.05) <= 1.e-8) max_dens = 4.779133e9; - else if (fabs(data.Time() - 0.06) <= 1.e-8) + else if (std::abs(data.Time() - 0.06) <= 1.e-8) max_dens = 24.472425e9; - else if (fabs(data.Time() - 0.065) <= 1.e-8) + else if (std::abs(data.Time() - 0.065) <= 1.e-8) max_dens = 423.447291e9; else { Print() << "Dont know the maximum density at this time: " << data.Time() < r, slicefile << std::setw(w) << r[i]; for (auto it=vars.begin(); it!=vars.end(); ++it) { - if(fabs((*it)[i]) < SMALL) (*it)[i] = 0.0; + if(std::abs((*it)[i]) < SMALL) (*it)[i] = 0.0; slicefile << std::setw(w) << (*it)[i]; } diff --git a/Diagnostics/Radiation/gaussian_pulse.cpp b/Diagnostics/Radiation/gaussian_pulse.cpp index b04020a567..a934fa7db4 100644 --- a/Diagnostics/Radiation/gaussian_pulse.cpp +++ b/Diagnostics/Radiation/gaussian_pulse.cpp @@ -64,9 +64,9 @@ int main(int argc, char* argv[]) // compute the size of the radially-binned array -- we'll do it to // the furtherest corner of the domain - double x_maxdist = max(fabs(probhi[0] - xctr), fabs(problo[0] - xctr)); - double y_maxdist = max(fabs(probhi[1] - yctr), fabs(problo[1] - yctr)); - double maxdist = sqrt(x_maxdist*x_maxdist + y_maxdist*y_maxdist); + double x_maxdist = std::max(std::abs(probhi[0] - xctr), std::abs(problo[0] - xctr)); + double y_maxdist = std::max(std::abs(probhi[1] - yctr), std::abs(problo[1] - yctr)); + double maxdist = std::sqrt(x_maxdist*x_maxdist + y_maxdist*y_maxdist); double dx_fine = *(std::min_element(dx.begin(), dx.end())); @@ -102,7 +102,7 @@ int main(int argc, char* argv[]) // over levels, we will compare to the finest level index space to // determine if we've already output here int mask_size = domain.length().max(); - Vector imask(pow(mask_size, AMREX_SPACEDIM), 1); + Vector imask(std::pow(mask_size, AMREX_SPACEDIM), 1); // loop over the data, starting at the finest grid, and if we haven't // already stored data in that grid location (according to imask), @@ -123,11 +123,11 @@ int main(int argc, char* argv[]) for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) { const Box& bx = mfi.tilebox(); - fgaussian_pulse(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), + fgaussian_pulse(AMREX_ARLIM_3D(bx.loVect()), AMREX_ARLIM_3D(bx.hiVect()), BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), nbins, rad_bin.dataPtr(), ncount.dataPtr(), imask.dataPtr(), mask_size, r1, - rad_comp, ZFILL(dx), dx_fine, xctr, yctr); + rad_comp, AMREX_ZFILL(dx), dx_fine, xctr, yctr); } diff --git a/Diagnostics/Radiation/lgt_frnt1d.cpp b/Diagnostics/Radiation/lgt_frnt1d.cpp index 55f61a5ffe..cf9ebfa5b5 100644 --- a/Diagnostics/Radiation/lgt_frnt1d.cpp +++ b/Diagnostics/Radiation/lgt_frnt1d.cpp @@ -53,7 +53,7 @@ int main(int argc, char* argv[]) // compute the size of the radially-binned array -- we'll do it to // the furtherest corner of the domain - double maxdist = fabs(probhi[0] - problo[0]); + double maxdist = std::abs(probhi[0] - problo[0]); double dx_fine = *(std::min_element(dx.begin(), dx.end())); int nbins = int(maxdist / dx_fine); @@ -93,7 +93,7 @@ int main(int argc, char* argv[]) // over levels, we will compare to the finest level index space to // determine if we've already output here int mask_size = domain.length().max(); - Vector imask(pow(mask_size, AMREX_SPACEDIM), 1); + Vector imask(std::pow(mask_size, AMREX_SPACEDIM), 1); // loop over the data, starting at the finest grid, and if we haven't // already stored data in that grid location (according to imask), @@ -114,13 +114,13 @@ int main(int argc, char* argv[]) for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) { const Box& bx = mfi.tilebox(); - flgt_frnt1d(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), + flgt_frnt1d(AMREX_ARLIM_3D(bx.loVect()), AMREX_ARLIM_3D(bx.hiVect()), BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), nbins, dens_bin.dataPtr(), vel_bin.dataPtr(), pres_bin.dataPtr(), rad_bin.dataPtr(), imask.dataPtr(), mask_size, r1, dens_comp, xmom_comp, pres_comp, rad_comp, - ZFILL(dx), dx_fine); + AMREX_ZFILL(dx), dx_fine); } // adjust r1 for the next lowest level diff --git a/Diagnostics/Radiation/rad_shock.cpp b/Diagnostics/Radiation/rad_shock.cpp index 43176c4e27..1707bb7fcc 100644 --- a/Diagnostics/Radiation/rad_shock.cpp +++ b/Diagnostics/Radiation/rad_shock.cpp @@ -106,7 +106,7 @@ int main(int argc, char* argv[]) // over levels, we will compare to the finest level index space to // determine if we've already output here int mask_size = domain.length().max(); - Vector imask(pow(mask_size, AMREX_SPACEDIM), 1); + Vector imask(std::pow(mask_size, AMREX_SPACEDIM), 1); // counter auto cnt = 0; @@ -135,7 +135,7 @@ int main(int argc, char* argv[]) for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) { const Box& bx = mfi.tilebox(); - fradshock(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), + fradshock(AMREX_ARLIM_3D(bx.loVect()), AMREX_ARLIM_3D(bx.hiVect()), problo.dataPtr(), probhi.dataPtr(), BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), nbins, vars_bin.dataPtr(), @@ -214,7 +214,7 @@ int main(int argc, char* argv[]) slicefile << std::setw(w) << vars_bin[isv[i]+((*it)+1)*nbins]; auto it = comps.end()-1; - slicefile << std::setw(w) << pow(vars_bin[isv[i]+((*it)+1)*nbins] / arad, 0.25); + slicefile << std::setw(w) << std::pow(vars_bin[isv[i]+((*it)+1)*nbins] / arad, 0.25); slicefile << std::endl; } diff --git a/Diagnostics/Radiation/rad_source.cpp b/Diagnostics/Radiation/rad_source.cpp index fae6885035..20260dc530 100644 --- a/Diagnostics/Radiation/rad_source.cpp +++ b/Diagnostics/Radiation/rad_source.cpp @@ -81,7 +81,7 @@ int main(int argc, char* argv[]) Real rhoe, rad; - fradsource(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), + fradsource(AMREX_ARLIM_3D(bx.loVect()), AMREX_ARLIM_3D(bx.hiVect()), BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), &rhoe, &rad, rhoe_comp, rad_comp); diff --git a/Diagnostics/Radiation/rad_sphere.cpp b/Diagnostics/Radiation/rad_sphere.cpp index 8eda551489..a60f03486c 100644 --- a/Diagnostics/Radiation/rad_sphere.cpp +++ b/Diagnostics/Radiation/rad_sphere.cpp @@ -142,7 +142,7 @@ int main(int argc, char* argv[]) // over levels, we will compare to the finest level index space to // determine if we've already output here int mask_size = domain.length()[0]; - Vector imask(pow(mask_size, AMREX_SPACEDIM), 1); + Vector imask(std::pow(mask_size, AMREX_SPACEDIM), 1); // counter int cnt = 0; @@ -163,12 +163,12 @@ int main(int argc, char* argv[]) for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) { const Box& bx = mfi.tilebox(); - fradsphere(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), - ZFILL(problo), ZFILL(probhi), + fradsphere(AMREX_ARLIM_3D(bx.loVect()), AMREX_ARLIM_3D(bx.hiVect()), + AMREX_ZFILL(problo), AMREX_ZFILL(probhi), BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), nbins, vars_bin.dataPtr(), imask.dataPtr(), mask_size, r1, - ZFILL(dx), &cnt); + AMREX_ZFILL(dx), &cnt); } // adjust r1 for the next lowest level diff --git a/Diagnostics/Radiation/rhd_shocktube.cpp b/Diagnostics/Radiation/rhd_shocktube.cpp index 9ed79f00ca..c81529fd1b 100644 --- a/Diagnostics/Radiation/rhd_shocktube.cpp +++ b/Diagnostics/Radiation/rhd_shocktube.cpp @@ -150,7 +150,7 @@ int main(int argc, char* argv[]) for (MFIter mfi(lev_data_mf, true); mfi.isValid(); ++mfi) { const Box& bx = mfi.tilebox(); - frhdshocktube(ARLIM_3D(bx.loVect()), ARLIM_3D(bx.hiVect()), + frhdshocktube(AMREX_ARLIM_3D(bx.loVect()), AMREX_ARLIM_3D(bx.hiVect()), BL_TO_FORTRAN_FAB(lev_data_mf[mfi]), nbins, dens_bin.dataPtr(), vel_bin.dataPtr(), pres_bin.dataPtr(), rad_bin.dataPtr(), diff --git a/Diagnostics/Sedov/main.cpp b/Diagnostics/Sedov/main.cpp index f6153d6bc0..4a32f24b6b 100644 --- a/Diagnostics/Sedov/main.cpp +++ b/Diagnostics/Sedov/main.cpp @@ -148,7 +148,7 @@ int main(int argc, char* argv[]) #else double x_maxdist = amrex::max(std::abs(probhi[0] - xctr), - std::fabs(problo[0] - xctr)); + std::abs(problo[0] - xctr)); double y_maxdist = amrex::max(std::abs(probhi[1] - yctr), std::abs(problo[1] - yctr)); double z_maxdist = amrex::max(std::abs(probhi[2] - zctr), @@ -342,10 +342,10 @@ int main(int argc, char* argv[]) // write the data in columns const auto SMALL = 1.e-20; for (auto i = 0; i < nbins; i++) { - if (fabs(dens_bin[i]) < SMALL) dens_bin[i] = 0.0; - if (fabs( vel_bin[i]) < SMALL) vel_bin[i] = 0.0; - if (fabs(pres_bin[i]) < SMALL) pres_bin[i] = 0.0; - if (fabs( e_bin[i]) < SMALL) e_bin[i] = 0.0; + if (std::abs(dens_bin[i]) < SMALL) dens_bin[i] = 0.0; + if (std::abs( vel_bin[i]) < SMALL) vel_bin[i] = 0.0; + if (std::abs(pres_bin[i]) < SMALL) pres_bin[i] = 0.0; + if (std::abs( e_bin[i]) < SMALL) e_bin[i] = 0.0; slicefile << std::setw(w) << r[i] << std::setw(w) << dens_bin[i] << std::setw(w) << vel_bin[i] << std::setw(w) << pres_bin[i] << std::setw(w) << e_bin[i] << std::endl; } diff --git a/Docs/source/FlowChart.rst b/Docs/source/FlowChart.rst index 5c306258f8..f24c4cfa7b 100644 --- a/Docs/source/FlowChart.rst +++ b/Docs/source/FlowChart.rst @@ -304,52 +304,72 @@ In the code, the objective is to evolve the state from the old time, #. *Initialize* - A. In ``initialize_do_advance()``, create ``Sborder``, initialized from ``S_old`` + In ``initialize_do_advance()``: - B. Call ``clean_state()`` to make sure the thermodynamics are in sync, in particular, - compute the temperature. + A. Create ``Sborder``, initialized from ``S_old`` + B. Call ``clean_state()`` to make sure the thermodynamics are in + sync, in particular, compute the temperature. -#. *React* :math:`\Delta t/2` [``strang_react_first_half()`` ] + C. [``SHOCK_VAR``] zero out the shock flag. - Update the solution due to the effect of reactions over half a time - step. The integration method and system of equations used here is - determined by a host of runtime parameters that are part of the - Microphysics package. But the basic idea is to evolve the energy - release from the reactions, the species mass fractions, and - temperature through :math:`\Delta t/2`. + D. Create the source corrector (if ``castro.source_term_predictor`` = 1) - Using the notation above, we begin with the time-level :math:`n` state, - :math:`\Ub^n`, and produce a state that has evolved only due to reactions, - :math:`\Ub^\star`. +#. *Do the pre-advance operations.* - .. math:: + This is handled by ``pre_advance_operators()`` and the main thing + that it does is the first half of the Strang burn. + + The steps are: + + A. *React* :math:`\Delta t/2` [``do_old_reactions()`` ] + + Update the solution due to the effect of reactions over half a + time step. The integration method and system of equations used + here is determined by a host of runtime parameters that are part + of the Microphysics package. But the basic idea is to evolve the + energy release from the reactions, the species mass fractions, + and temperature through :math:`\Delta t/2`. + + Using the notation above, we begin with the time-level :math:`n` state, + :math:`\Ub^n`, and produce a state that has evolved only due to reactions, + :math:`\Ub^\star`. - \begin{aligned} + .. math:: + + \begin{aligned} (\rho e)^\star &= (\rho e)^n + \frac{\dt}{2} \rho H_\mathrm{nuc} \\ (\rho E)^\star &= (\rho E)^n + \frac{\dt}{2} \rho H_\mathrm{nuc} \\ (\rho X_k)^\star &= (\rho X_k)^n + \frac{\dt}{2}(\rho\omegadot_k). \end{aligned} - Here, :math:`H_\mathrm{nuc}` is the energy release (erg/g/s) over the - burn, and :math:`\omegadot_k` is the creation rate for species :math:`k`. + Here, :math:`H_\mathrm{nuc}` is the energy release (erg/g/s) over the + burn, and :math:`\omegadot_k` is the creation rate for species :math:`k`. + + After exiting the burner, we call the EOS with :math:`\rho^\star`, + :math:`e^\star`, and :math:`X_k^\star` to get the new temperature, :math:`T^\star`. + + .. note:: - After exiting the burner, we call the EOS with :math:`\rho^\star`, - :math:`e^\star`, and :math:`X_k^\star` to get the new temperature, :math:`T^\star`. + The density, :math:`\rho`, does not change via reactions in the + Strang-split formulation. - Note that the density, :math:`\rho`, does not change via reactions in the - Strang-split formulation. + The reaction data needs to be valid in the ghost cells, so the reactions + are applied to the entire patch, including ghost cells. - The reaction data needs to be valid in the ghost cells, so the reactions - are applied to the entire patch, including ghost cells. + After reactions, ``clean_state`` is called. - After reactions, ``clean_state`` is called. + B. *Construct the gravitational potential at time $n$.* + + This is done by calling ``construct_old_gravity()`` + + C. *Initialize ``S_new`` with the current state* (``Sborder``). At the end of this step, ``Sborder`` sees the effects of the reactions. #. *Construct time-level n sources and apply* - [``construct_old_gravity()``, ``do_old_sources()`` ] + [``do_old_sources()`` ] The time level :math:`n` sources are computed, and added to the StateData ``Source_Type``. @@ -426,80 +446,81 @@ In the code, the objective is to evolve the state from the old time, with Strang-splitting, since the hydro and sources takes place completely inside of the surrounding burn operations. - The old-time source terms are stored in ``old_source``. + The old-time source terms are stored in ``old_source`` (and a ghost + cell fill is performed). The sources are then applied to the state after the burn, :math:`\Ub^\star` with a full :math:`\Delta t` weighting (this will be corrected later). This produces the intermediate state, :math:`\Ub^{n+1,(a)}` (stored in ``S_new``). +#. *Do pre-hydro operations* [``pre_hydro_operators()``] + + For Strang+CTU, nothing is done here. + #. *Construct the hydro / MHD update* [``construct_ctu_hydro_source()``, ``construct_ctu_mhd_source()``] The goal is to advance our system considering only the advective terms (which in Cartesian coordinates can be written as the divergence of a flux). - We do the hydro update in two parts—first we construct the - advective update and store it in the hydro_source - MultiFab, then we do the conservative update in a separate step. This - separation allows us to use the advective update separately in more - complex time-integration schemes. + A. In the Strang-split formulation, we start the reconstruction + using the state after burning, :math:`\Ub^\star` (``Sborder``). + For the CTU method, we predict to the half-time (:math:`n+1/2`) + to get a second-order accurate method. Note: ``Sborder`` does + not know of any sources except for reactions. - In the Strang-split formulation, we start the reconstruction using - the state after burning, :math:`\Ub^\star` (``Sborder``). For the - CTU method, we predict to the half-time (:math:`n+1/2`) to get a - second-order accurate method. Note: ``Sborder`` does not know of - any sources except for reactions. + The method done here differs depending on whether we are doing hydro or MHD. - The method done here differs depending on whether we are doing hydro or MHD. + * hydrodynamics - A. hydrodynamics + The advection step is complicated, and more detail is given in + Section :ref:`Sec:Advection Step`. Here is the summarized version: - The advection step is complicated, and more detail is given in - Section :ref:`Sec:Advection Step`. Here is the summarized version: + i. Compute primitive variables. - i. Compute primitive variables. + ii. Convert the source terms to those acting on primitive variables - ii. Convert the source terms to those acting on primitive variables + iii. Predict primitive variables to time-centered edges. - iii. Predict primitive variables to time-centered edges. + iv. Solve the Riemann problem. - iv. Solve the Riemann problem. + v. Compute fluxes and advective term. - v. Compute fluxes and advective term. + * MHD - B. MHD + The MHD update is described in :ref:`ch:mhd`. - The MHD update is described in :ref:`ch:mhd`. + To start the hydrodynamics/MHD source construction, we need to know + the hydrodynamics source terms at time-level :math:`n`, since this + enters into the prediction to the interface states. This is + essentially the same vector that was computed in the previous step, + with a few modifications. The most important is that if we set + ``castro.source_term_predictor``, then we extrapolate the source + terms from :math:`n` to :math:`n+1/2`, using the change from the + previous step. - To start the hydrodynamics/MHD source construction, we need to know - the hydrodynamics source terms at time-level :math:`n`, since this - enters into the prediction to the interface states. This is - essentially the same vector that was computed in the previous step, - with a few modifications. The most important is that if we set - ``castro.source_term_predictor``, then we extrapolate the source - terms from :math:`n` to :math:`n+1/2`, using the change from the - previous step. + Note: we neglect the reaction source terms, since those are already + accounted for in the state directly, due to the Strang-splitting + nature of this method. - Note: we neglect the reaction source terms, since those are already - accounted for in the state directly, due to the Strang-splitting - nature of this method. + The update computed here is then immediately applied to + ``S_new``. - The update computed here is then immediately applied to - ``S_new``. + B. *Clean State and check for NaNs* [``clean_state()``] -#. *Clean State* [``clean_state()``] + This is done on ``S_new``. - This is done on ``S_new``. + C. *Update the center of mass for monopole gravity* - After these checks, we check the state for NaNs. + This quantities are computed using ``S_new``. -#. *Update radial data and center of mass for monopole gravity* +#. *Do post-hydro operations* [``post_hydro_operators()``] - These quantities are computed using ``S_new``. + This constructs the new gravitational potential. #. *Correct the source terms with the n+1 - contribution* [``construct_new_gravity()``, ``do_new_sources`` ] + contribution* [``do_new_sources`` ] If we are doing self-gravity, then we first compute the updated gravitational potential using the updated density from ``S_new``. @@ -521,24 +542,19 @@ In the code, the objective is to evolve the state from the old time, In the process of updating the sources, we update the temperature to make it consistent with the new state. -#. *React* :math:`\Delta t/2` [``strang_react_second_half()``] +#. *Do post advance operations* [``post_advance_operators()``] - We do the final :math:`\dt/2` reacting on the state, beginning with :math:`\Ub^{n+1,(c)}` to - give us the final state on this level, :math:`\Ub^{n+1}`. + This simply does the final :math:`\dt/2` reacting on the state, + beginning with :math:`\Ub^{n+1,(c)}` to give us the final state on + this level, :math:`\Ub^{n+1}`. This is largely the same as ``strang_react_first_half()``, but it does not currently fill the reactions in the ghost cells. #. *Finalize* [``finalize_do_advance()``] - Finalize does the following: - - A. for the momentum sources, we compute :math:`d\Sb/dt`, to use in the - source term prediction/extrapolation for the hydrodynamic - interface states during the next step. - - B. If we are doing the hybrid momentum algorithm, then we sync up - the hybrid and linear momenta + This checks to ensure that we didn't violate the CFL criteria + during the advance. A summary of which state is the input and which is updated for each of these processes is presented below: @@ -559,7 +575,7 @@ these processes is presented below: +--------------------+-----------+---------------------+---------------------+ | 5. clean | | | input / updated | +--------------------+-----------+---------------------+---------------------+ - | 6. radial / center | | | input | + | 6. center of mass | | | input | +--------------------+-----------+---------------------+---------------------+ | 7. correct sources | | | input / updated | +--------------------+-----------+---------------------+---------------------+ diff --git a/Docs/source/faq.rst b/Docs/source/faq.rst index 4e3d03daf3..c723b745db 100644 --- a/Docs/source/faq.rst +++ b/Docs/source/faq.rst @@ -106,10 +106,10 @@ Debugging :: - print_state(mf, IntVect(D_DECL(10, 20, 30))); + print_state(mf, IntVect(AMREX_D_DECL(10, 20, 30))); Here, the IntVect has the dimension that we were compiled with - (and this is handled through the preprocessor ``D_DECL``). In + (and this is handled through the preprocessor ``AMREX_D_DECL``). In this case, we are inspecting zone (10, 20, 30), in the global index space. Note that since a multifab exists only on a single level, the integer indices here refer to the global index space on that level. diff --git a/Docs/source/mpi_plus_x.rst b/Docs/source/mpi_plus_x.rst index 5906fa68f6..fe583bc8fb 100644 --- a/Docs/source/mpi_plus_x.rst +++ b/Docs/source/mpi_plus_x.rst @@ -95,6 +95,16 @@ To enable this, compile with:: USE_OMP = FALSE USE_CUDA = TRUE +.. note:: + + For recent GPUs, like the NVIDIA RTX 4090, you may need to change + the default CUDA architecture. This can be done by adding: + + .. code:: + + CUDA_ARCH=89 + + to the ``make`` line or ``GNUmakefile``. AMD GPUs -------- diff --git a/Exec/gravity_tests/DustCollapse/problem_bc_fill.H b/Exec/gravity_tests/DustCollapse/problem_bc_fill.H index 2a8df96b84..b4a937f816 100644 --- a/Exec/gravity_tests/DustCollapse/problem_bc_fill.H +++ b/Exec/gravity_tests/DustCollapse/problem_bc_fill.H @@ -5,8 +5,8 @@ AMREX_GPU_HOST_DEVICE AMREX_INLINE void problem_bc_fill(int i, int j, int k, Array4 const& state, Real time, - Array1D bcs, - GeometryData geomdata) + const Array1D& bcs, + const GeometryData& geomdata) { #if AMREX_SPACEDIM == 3 const Real* dx = geomdata.CellSize(); @@ -30,7 +30,7 @@ void problem_bc_fill(int i, int j, int k, Real r = std::sqrt(x * x + y * y + z * z); // XLO - if (bcs(0).lo(0) == FOEXTRAP && i < domlo[0]) { + if (bcs(0).lo(0) == amrex::BCType::foextrap && i < domlo[0]) { int ic = domlo[0]; Real xc = problo[0] + (static_cast(ic) + 0.5_rt) * dx[0] - problem::center[0]; @@ -50,7 +50,7 @@ void problem_bc_fill(int i, int j, int k, } // XHI - if (bcs(0).hi(0) == FOEXTRAP && i > domhi[0]) { + if (bcs(0).hi(0) == amrex::BCType::foextrap && i > domhi[0]) { int ic = domhi[0]; Real xc = problo[0] + (static_cast(ic) + 0.5_rt) * dx[0] - problem::center[0]; @@ -71,7 +71,7 @@ void problem_bc_fill(int i, int j, int k, #if AMREX_SPACEDIM >= 2 // YLO - if (bcs(0).lo(1) == FOEXTRAP && j < domlo[1]) { + if (bcs(0).lo(1) == amrex::BCType::foextrap && j < domlo[1]) { int jc = domlo[1]; Real yc = problo[1] + (static_cast(jc) + 0.5_rt) * dx[1] - problem::center[1]; @@ -91,7 +91,7 @@ void problem_bc_fill(int i, int j, int k, } // YHI - if (bcs(0).hi(1) == FOEXTRAP && j > domhi[1]) { + if (bcs(0).hi(1) == amrex::BCType::foextrap && j > domhi[1]) { int jc = domhi[1]; Real yc = problo[1] + (static_cast(jc) + 0.5_rt) * dx[1] - problem::center[1]; @@ -112,7 +112,7 @@ void problem_bc_fill(int i, int j, int k, #if AMREX_SPACEDIM == 3 // ZLO - if (bcs(0).lo(2) == FOEXTRAP && k < domlo[2]) { + if (bcs(0).lo(2) == amrex::BCType::foextrap && k < domlo[2]) { int kc = domlo[2]; Real zc = problo[2] + (static_cast(kc) + 0.5_rt) * dx[2] - problem::center[2]; @@ -132,7 +132,7 @@ void problem_bc_fill(int i, int j, int k, } // ZHI - if (bcs(0).hi(2) == FOEXTRAP && k > domhi[2]) { + if (bcs(0).hi(2) == amrex::BCType::foextrap && k > domhi[2]) { int kc = domhi[2]; Real zc = problo[2] + (static_cast(kc) + 0.5_rt) * dx[2] - problem::center[2]; diff --git a/Exec/gravity_tests/hydrostatic_adjust/problem_bc_fill.H b/Exec/gravity_tests/hydrostatic_adjust/problem_bc_fill.H index 96ab44b365..9a28bef1d3 100644 --- a/Exec/gravity_tests/hydrostatic_adjust/problem_bc_fill.H +++ b/Exec/gravity_tests/hydrostatic_adjust/problem_bc_fill.H @@ -5,8 +5,8 @@ AMREX_GPU_HOST_DEVICE AMREX_INLINE void problem_bc_fill(int i, int j, int k, Array4 const& state, Real time, - Array1D bcs, - GeometryData geomdata) + const Array1D& bcs, + const GeometryData& geomdata) { // Set the velocity to zero at the top physical boundary. @@ -14,7 +14,7 @@ void problem_bc_fill(int i, int j, int k, if (AMREX_SPACEDIM == 1) { - if (bcs(UMX).hi(0) == FOEXTRAP && i > domhi[0]) { + if (bcs(UMX).hi(0) == amrex::BCType::foextrap && i > domhi[0]) { Real vel = amrex::max(state(i,j,k,UMX) / state(i,j,k,URHO), 0.0); state(i,j,k,URHO) = problem::hse_rho_top; @@ -36,7 +36,7 @@ void problem_bc_fill(int i, int j, int k, } else if (AMREX_SPACEDIM == 2) { - if (bcs(UMY).hi(1) == FOEXTRAP && j > domhi[1]) { + if (bcs(UMY).hi(1) == amrex::BCType::foextrap && j > domhi[1]) { Real vel = amrex::max(state(i,j,k,UMY) / state(i,j,k,URHO), 0.0); state(i,j,k,URHO) = problem::hse_rho_top; @@ -58,7 +58,7 @@ void problem_bc_fill(int i, int j, int k, } else { - if (bcs(UMZ).hi(2) == FOEXTRAP && k > domhi[2]) { + if (bcs(UMZ).hi(2) == amrex::BCType::foextrap && k > domhi[2]) { Real vel = amrex::max(state(i,j,k,UMZ) / state(i,j,k,URHO), 0.0); state(i,j,k,URHO) = problem::hse_rho_top; diff --git a/Exec/hydro_tests/Noh/problem_bc_fill.H b/Exec/hydro_tests/Noh/problem_bc_fill.H index 3c79eb0741..3a6470744f 100644 --- a/Exec/hydro_tests/Noh/problem_bc_fill.H +++ b/Exec/hydro_tests/Noh/problem_bc_fill.H @@ -8,8 +8,8 @@ AMREX_GPU_HOST_DEVICE AMREX_INLINE void problem_bc_fill(int i, int j, int k, Array4 const& state, Real time, - Array1D bcs, - GeometryData geomdata) + const Array1D& bcs, + const GeometryData& geomdata) { const Real pres_init = 1.0e-6_rt; const Real rho_init = 1.0e0_rt; diff --git a/Exec/hydro_tests/double_mach_reflection/problem_bc_fill.H b/Exec/hydro_tests/double_mach_reflection/problem_bc_fill.H index c1edfe36c7..aaba204f55 100644 --- a/Exec/hydro_tests/double_mach_reflection/problem_bc_fill.H +++ b/Exec/hydro_tests/double_mach_reflection/problem_bc_fill.H @@ -5,8 +5,8 @@ AMREX_GPU_HOST_DEVICE AMREX_INLINE void problem_bc_fill(int i, int j, int k, Array4 const& state, Real time, - Array1D bcs, - GeometryData geomdata) + const Array1D& bcs, + const GeometryData& geomdata) { const int* domlo = geomdata.Domain().loVect(); const int* domhi = geomdata.Domain().hiVect(); @@ -18,7 +18,7 @@ void problem_bc_fill(int i, int j, int k, // x boundaries //------------------------------------------------------------------------- - if ((bcs(URHO).lo(0) == EXT_DIR || bcs(URHO).lo(0) == FOEXTRAP) && i < domlo[0]) { + if ((bcs(URHO).lo(0) == amrex::BCType::ext_dir || bcs(URHO).lo(0) == amrex::BCType::foextrap) && i < domlo[0]) { state(i,j,k,URHO) = problem::rho_l; state(i,j,k,UMX) = problem::rho_l * problem::u_l; state(i,j,k,UMY) = problem::rho_l * problem::v_l; @@ -43,7 +43,7 @@ void problem_bc_fill(int i, int j, int k, Real y = problo[1] + dx[1] * (static_cast(j) + 0.5e0_rt); // YLO - if ((bcs(URHO).lo(1) == EXT_DIR || bcs(URHO).lo(1) == FOEXTRAP) && j < domlo[1]) { + if ((bcs(URHO).lo(1) == amrex::BCType::ext_dir || bcs(URHO).lo(1) == amrex::BCType::foextrap) && j < domlo[1]) { if (x < 1.e0_rt / 6.e0_rt) { // ICs state(i,j,k,URHO) = problem::rho_l; @@ -71,7 +71,7 @@ void problem_bc_fill(int i, int j, int k, } // YHI - if ((bcs(URHO).hi(2) == EXT_DIR || bcs(URHO).hi(2) == FOEXTRAP) && j > domhi[1]) { + if ((bcs(URHO).hi(2) == amrex::BCType::ext_dir || bcs(URHO).hi(2) == amrex::BCType::foextrap) && j > domhi[1]) { state(i,j,k,URHO ) = 0.e0_rt; state(i,j,k,UMX ) = 0.e0_rt; state(i,j,k,UMY ) = 0.e0_rt; diff --git a/Exec/radiation_tests/Rad2Tshock/problem_bc_fill.H b/Exec/radiation_tests/Rad2Tshock/problem_bc_fill.H index bb2f45ce6d..d23413fef3 100644 --- a/Exec/radiation_tests/Rad2Tshock/problem_bc_fill.H +++ b/Exec/radiation_tests/Rad2Tshock/problem_bc_fill.H @@ -5,8 +5,8 @@ AMREX_GPU_HOST_DEVICE AMREX_INLINE void problem_bc_fill(int i, int j, int k, Array4 const& state, Real time, - Array1D bcs, - GeometryData geomdata) + const Array1D& bcs, + const GeometryData& geomdata) { // The strategy here is to set Dirichlet condition for inflow and // outflow boundaries, and let the Riemann solver sort out the @@ -18,7 +18,7 @@ void problem_bc_fill(int i, int j, int k, if (problem::idir == 1) { // XLO - if ((bcs(URHO).lo(0) == EXT_DIR || bcs(URHO).lo(0) == FOEXTRAP) && i < domlo[0]) { + if ((bcs(URHO).lo(0) == amrex::BCType::ext_dir || bcs(URHO).lo(0) == amrex::BCType::foextrap) && i < domlo[0]) { state(i,j,k,URHO) = problem::rho0; state(i,j,k,UMX) = problem::rho0 * problem::v0; state(i,j,k,UMY) = 0.e0_rt; @@ -33,7 +33,7 @@ void problem_bc_fill(int i, int j, int k, } // XHI - if ((bcs(URHO).hi(0) == EXT_DIR || bcs(URHO).hi(0) == FOEXTRAP) && i > domhi[0]) { + if ((bcs(URHO).hi(0) == amrex::BCType::ext_dir || bcs(URHO).hi(0) == amrex::BCType::foextrap) && i > domhi[0]) { state(i,j,k,URHO) = problem::rho1; state(i,j,k,UMX) = problem::rho1 * problem::v1; state(i,j,k,UMY) = 0.e0_rt; @@ -51,7 +51,7 @@ void problem_bc_fill(int i, int j, int k, else if (problem::idir == 2) { // YLO - if ((bcs(URHO).lo(1) == EXT_DIR || bcs(URHO).lo(1) == FOEXTRAP) && j < domlo[1]) { + if ((bcs(URHO).lo(1) == amrex::BCType::ext_dir || bcs(URHO).lo(1) == amrex::BCType::foextrap) && j < domlo[1]) { state(i,j,k,URHO) = problem::rho0; state(i,j,k,UMX) = 0.e0_rt; state(i,j,k,UMY) = problem::rho0 * problem::v0; @@ -66,7 +66,7 @@ void problem_bc_fill(int i, int j, int k, } // YHI - if ((bcs(URHO).hi(1) == EXT_DIR || bcs(URHO).hi(1) == FOEXTRAP) && j > domhi[1]) { + if ((bcs(URHO).hi(1) == amrex::BCType::ext_dir || bcs(URHO).hi(1) == amrex::BCType::foextrap) && j > domhi[1]) { state(i,j,k,URHO) = problem::rho1; state(i,j,k,UMX) = 0.e0_rt; state(i,j,k,UMY) = problem::rho1 * problem::v1; @@ -84,7 +84,7 @@ void problem_bc_fill(int i, int j, int k, else if (problem::idir == 3) { // ZLO - if ((bcs(URHO).lo(2) == EXT_DIR || bcs(URHO).lo(2) == FOEXTRAP) && k < domlo[2]) { + if ((bcs(URHO).lo(2) == amrex::BCType::ext_dir || bcs(URHO).lo(2) == amrex::BCType::foextrap) && k < domlo[2]) { state(i,j,k,URHO) = problem::rho0; state(i,j,k,UMX) = 0.e0_rt; state(i,j,k,UMY) = 0.e0_rt; @@ -99,7 +99,7 @@ void problem_bc_fill(int i, int j, int k, } // ZHI - if ((bcs(URHO).hi(2) == EXT_DIR || bcs(URHO).hi(2) == FOEXTRAP) && k > domhi[2]) { + if ((bcs(URHO).hi(2) == amrex::BCType::ext_dir || bcs(URHO).hi(2) == amrex::BCType::foextrap) && k > domhi[2]) { state(i,j,k,URHO) = problem::rho1; state(i,j,k,UMX) = 0.e0_rt; state(i,j,k,UMY) = 0.e0_rt; diff --git a/Exec/science/Detonation/problem_bc_fill.H b/Exec/science/Detonation/problem_bc_fill.H index 02e4089e05..d16f07bf0d 100644 --- a/Exec/science/Detonation/problem_bc_fill.H +++ b/Exec/science/Detonation/problem_bc_fill.H @@ -45,8 +45,8 @@ AMREX_GPU_HOST_DEVICE AMREX_INLINE void problem_bc_fill(int i, int j, int k, Array4 const& state, Real time, - Array1D bcs, - GeometryData geomdata) + const Array1D& bcs, + const GeometryData& geomdata) { // Override the generic routine at the physical boundaries by // setting the material to the ambient state. Note that we @@ -63,39 +63,39 @@ void problem_bc_fill(int i, int j, int k, // Note: in these checks we're only looking at the boundary // conditions on the first state component (density). - if (i < domlo[0] && (bcs(0).lo(0) != REFLECT_ODD && - bcs(0).lo(0) != INT_DIR && - bcs(0).lo(0) != REFLECT_EVEN)) { + if (i < domlo[0] && (bcs(0).lo(0) != amrex::BCType::reflect_odd && + bcs(0).lo(0) != amrex::BCType::int_dir && + bcs(0).lo(0) != amrex::BCType::reflect_even)) { do_ambient_fill = true; } - else if (i > domhi[0] && (bcs(0).hi(0) != REFLECT_ODD && - bcs(0).hi(0) != INT_DIR && - bcs(0).hi(0) != REFLECT_ODD)) { + else if (i > domhi[0] && (bcs(0).hi(0) != amrex::BCType::reflect_odd && + bcs(0).hi(0) != amrex::BCType::int_dir && + bcs(0).hi(0) != amrex::BCType::reflect_even)) { do_ambient_fill = true; } if (AMREX_SPACEDIM >= 2) { - if (j < domlo[1] && (bcs(0).lo(1) != REFLECT_ODD && - bcs(0).lo(1) != INT_DIR && - bcs(0).lo(1) != REFLECT_ODD)) { + if (j < domlo[1] && (bcs(0).lo(1) != amrex::BCType::reflect_odd && + bcs(0).lo(1) != amrex::BCType::int_dir && + bcs(0).lo(1) != amrex::BCType::reflect_even)) { do_ambient_fill = true; } - else if (j > domhi[1] && (bcs(0).hi(1) != REFLECT_ODD && - bcs(0).hi(1) != INT_DIR && - bcs(0).hi(1) != REFLECT_ODD)) { + else if (j > domhi[1] && (bcs(0).hi(1) != amrex::BCType::reflect_odd && + bcs(0).hi(1) != amrex::BCType::int_dir && + bcs(0).hi(1) != amrex::BCType::reflect_even)) { do_ambient_fill = true; } } - + if (AMREX_SPACEDIM == 3) { - if (k < domlo[2] && (bcs(0).lo(2) != REFLECT_ODD && - bcs(0).lo(2) != INT_DIR && - bcs(0).lo(2) != REFLECT_ODD)) { + if (k < domlo[2] && (bcs(0).lo(2) != amrex::BCType::reflect_odd && + bcs(0).lo(2) != amrex::BCType::int_dir && + bcs(0).lo(2) != amrex::BCType::reflect_even)) { do_ambient_fill = true; } - else if (k > domhi[2] && (bcs(0).hi(2) != REFLECT_ODD && - bcs(0).hi(2) != INT_DIR && - bcs(0).hi(2) != REFLECT_ODD)) { + else if (k > domhi[2] && (bcs(0).hi(2) != amrex::BCType::reflect_odd && + bcs(0).hi(2) != amrex::BCType::int_dir && + bcs(0).hi(2) != amrex::BCType::reflect_even)) { do_ambient_fill = true; } } diff --git a/Exec/science/massive_star/inputs_2d.nse b/Exec/science/massive_star/inputs_2d.nse index ed703795a0..38d15ec74a 100644 --- a/Exec/science/massive_star/inputs_2d.nse +++ b/Exec/science/massive_star/inputs_2d.nse @@ -58,11 +58,23 @@ amr.regrid_int = 10000 # how often to regrid amr.n_error_buf = 2 2 2 2 # number of buffer cells in error est amr.grid_eff = 0.7 # what constitutes an efficient grid -amr.check_file = chk # root name of checkpoint file -amr.check_int = 25 # number of timesteps between checkpoints -amr.plot_file = plt # root name of plot file +amr.check_file = massive_star_chk # root name of checkpoint file +amr.check_int = 50 # number of timesteps between checkpoints + +amr.plot_file = massive_star_plt # root name of plot file +amr.plot_per = 5.0 +amr.derive_plot_vars = ALL +castro.store_burn_weights = 1 + +amr.small_plot_file = massive_star_smallplt # root name of plot file +amr.small_plot_per = 0.5 +amr.small_plot_vars = density Temp +amr.derive_small_plot_vars = abar Ye enuc MachNumber magvel magvort in_nse + +fab.format = NATIVE_32 + castro.plot_per_is_exact = 0 -amr.plot_per = 0.25 + amr.max_grid_size = 512 # maximum grid size allowed -- used to control parallelism amr.blocking_factor = 32 # block factor in grid generation @@ -70,9 +82,6 @@ amr.blocking_factor = 32 # block factor in grid generation amr.v = 1 # control verbosity in Amr.cpp castro.v = 1 # control verbosity in Castro.cpp -amr.derive_plot_vars = ALL - -castro.store_burn_weights = 1 castro.small_dens = 1.0 castro.small_temp = 1.e6 @@ -94,7 +103,7 @@ problem.interpolate_pres = 1 castro.drive_initial_convection = 1 castro.drive_initial_convection_reinit_period = 2 -castro.drive_initial_convection_tmax = 250 +castro.drive_initial_convection_tmax = 100 # refinement diff --git a/Exec/science/nova/problem_initialize_state_data.H b/Exec/science/nova/problem_initialize_state_data.H index 37a7bb3894..4dea3d43eb 100644 --- a/Exec/science/nova/problem_initialize_state_data.H +++ b/Exec/science/nova/problem_initialize_state_data.H @@ -121,10 +121,10 @@ void problem_initialize_state_data (int i, int j, int k, if (problem::num_vortices % 2 == 0){ temp_pert = delta * (1.0_rt + std::sin(static_cast(problem::num_vortices)*M_PI*x/problem::L) )* - delta*exp(-pow(ydist / problem::width, 2)); + delta*std::exp(-pow(ydist / problem::width, 2)); } else { temp_pert = delta * (1.0_rt + std::cos(static_cast(problem::num_vortices)*M_PI*x/problem::L) )* - delta*exp(-pow(ydist / problem::width, 2)); + delta*std::exp(-pow(ydist / problem::width, 2)); } state(i,j,k,UTEMP) = state(i,j,k, UTEMP)*(1.0_rt + temp_pert); diff --git a/Exec/science/planet/HotJupiter.cpp b/Exec/science/planet/HotJupiter.cpp index 108db17923..059b255df5 100644 --- a/Exec/science/planet/HotJupiter.cpp +++ b/Exec/science/planet/HotJupiter.cpp @@ -19,10 +19,10 @@ long double Temp(long double P1,long double Pc,long double P_char,long double T_ long double Temp; if(Pc>P1){ - Temp = Td * pow(1.0e0 + (Lad / (Lin - Lad))*pow((P1/Pc),exp1),exp2); + Temp = Td * std::pow(1.0e0 + (Lad / (Lin - Lad)) * std::pow((P1/Pc),exp1),exp2); } else{ - Temp = T_char * pow(P1/P_char,Lad); + Temp = T_char * std::pow(P1/P_char,Lad); } return Temp; } @@ -95,9 +95,9 @@ int main(){ //Defining the properties of the planet atmosphere // the density limit at the top of the atmosphere P_top=den_top*R*Td ; // the pressure limit at the density limit - Tc = Td * pow(Lin/(Lin-Lad),exp2); //the temperature at the radiative-convective boundary (RCB) - Pc = P_char * pow(Tc/T_char,1.0/Lad); //the pressure at the RCB - Pd = pow((Lin-Lad)/(2.0*Lin-Lad),1.0/exp1)*Pc; //the characteristic pressure defining the radiative zone + Tc = Td * std::pow(Lin/(Lin-Lad),exp2); //the temperature at the radiative-convective boundary (RCB) + Pc = P_char * std::pow(Tc/T_char,1.0/Lad); //the pressure at the RCB + Pd = std::pow((Lin-Lad)/(2.0*Lin-Lad),1.0/exp1)*Pc; //the characteristic pressure defining the radiative zone @@ -112,7 +112,7 @@ int main(){ } buffer_height=z_top*2.25/4.0; // the height at which the buffer and the top of the atmosphere meet (not continuous at RCB).-> needs sponge - optical_depth_buffer=6.35e-3*T_buffer*pow(den_buffer,2.0)*(z_top-buffer_height);// the optical depth in the buffer region [dimensionless] + optical_depth_buffer=6.35e-3*T_buffer*std::pow(den_buffer,2.0)*(z_top-buffer_height);// the optical depth in the buffer region [dimensionless] @@ -157,14 +157,14 @@ int main(){ den2 = density(P2, T2, R); T3 = T1; den3 = den1; - dz2 = abs ((-P1+P2)* (1.0/grav)/((1.0/2.0)*(den1+den2))); - dz3 = abs ((-P1+P3)* (1.0/grav)/((1.0/2.0)*(den1+den3))); + dz2 = std::abs ((-P1+P2)* (1.0/grav)/((1.0/2.0)*(den1+den2))); + dz3 = std::abs ((-P1+P3)* (1.0/grav)/((1.0/2.0)*(den1+den3))); a: if(difference(dz2, delta_z)*difference(dz3, delta_z)<0.0){ P_sol = (P2+P3)/2.0; T_sol = Temp(P_sol, Pc, P_char, T_char, Td, Lin, Lad, exp1, exp2); den_sol = density(P_sol,T_sol, R); - dz_sol = abs ((-P1+P_sol)* (1.0/grav)/((1.0/2.0)*(den1+den_sol))); + dz_sol = std::abs ((-P1+P_sol)* (1.0/grav)/((1.0/2.0)*(den1+den_sol))); count2++; } else{ @@ -172,13 +172,13 @@ int main(){ P3 = P3 - P3/dP_factor/1.0; T3 = Temp(P3, Pc, P_char, T_char, Td, Lin, Lad, exp1, exp2); den3 = density(P3, T3, R); - dz3 = abs ((-P1+P3)* (1.0/grav)/((1.0/2.0)*(den1+den3))); + dz3 = std::abs ((-P1+P3)* (1.0/grav)/((1.0/2.0)*(den1+den3))); } else if(difference(dz2, delta_z)<0.0){ P2 = P2 + P2/dP_factor/1.0; T2 = Temp(P2, Pc, P_char, T_char, Td, Lin, Lad, exp1, exp2); den2 = density(P2,T2,R); - dz2 = abs ((-P1+P2)* (1.0/grav)/((1.0/2.0)*(den1+den2))); + dz2 = std::abs ((-P1+P2)* (1.0/grav)/((1.0/2.0)*(den1+den2))); } goto a; } @@ -192,11 +192,11 @@ int main(){ if(count2>100000){ goto b; } - if(abs(difference(dz_sol, delta_z)) >= error_expected){ + if(std::abs(difference(dz_sol, delta_z)) >= error_expected){ goto a; } - b: cout< >& amr_lev auto deltax = it->m_pos[0] - p.m_pos[0]; auto deltay = it->m_pos[1] - p.m_pos[1]; - auto delta = sqrt(deltax*deltax + deltay*deltay); + auto delta = std::sqrt(deltax*deltax + deltay*deltay); // quick nan check here if (delta == delta) diff --git a/README.md b/README.md index e52b0137be..79fe71126d 100644 --- a/README.md +++ b/README.md @@ -3,6 +3,7 @@ [![AMReX](https://amrex-codes.github.io/badges/powered%20by-AMReX-red.svg)](https://amrex-codes.github.io) [![yt-project](https://img.shields.io/static/v1?label="works%20with"&message="yt"&color="blueviolet")](https://yt-project.org) [![github pages](https://github.com/AMReX-Astro/Castro/workflows/github%20pages/badge.svg)](https://github.com/AMReX-Astro/Castro/actions?query=workflow%3A%22github+pages%22) +[![coverity](https://scan.coverity.com/projects/29689/badge.svg)](https://scan.coverity.com/projects/amrex-astro-castro) ![Castro](https://github.com/AMReX-Astro/Castro/blob/development/Util/logo/castro_logo_hot_200.png) diff --git a/Source/diffusion/Diffusion.cpp b/Source/diffusion/Diffusion.cpp index 3fc02063d8..2f478fdc39 100644 --- a/Source/diffusion/Diffusion.cpp +++ b/Source/diffusion/Diffusion.cpp @@ -105,12 +105,12 @@ Diffusion::make_mg_bc () mlmg_lobc[idim] = MLLinOp::BCType::Periodic; mlmg_hibc[idim] = MLLinOp::BCType::Periodic; } else { - if (phys_bc->lo(idim) == Inflow) { + if (phys_bc->lo(idim) == amrex::PhysBCType::inflow) { mlmg_lobc[idim] = MLLinOp::BCType::Dirichlet; } else { mlmg_lobc[idim] = MLLinOp::BCType::Neumann; } - if (phys_bc->hi(idim) == Symmetry) { + if (phys_bc->hi(idim) == amrex::PhysBCType::symmetry) { mlmg_hibc[idim] = MLLinOp::BCType::Dirichlet; } else { mlmg_hibc[idim] = MLLinOp::BCType::Neumann; diff --git a/Source/driver/Castro.H b/Source/driver/Castro.H index 9cc45996c3..11b028f5f9 100644 --- a/Source/driver/Castro.H +++ b/Source/driver/Castro.H @@ -1487,8 +1487,10 @@ protected: /// sdc /// int sdc_iteration; - int current_sdc_node; +#ifdef TRUE_SDC + int current_sdc_node; +#endif /* problem-specific includes */ diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp index 7d114996e8..a0217208cf 100644 --- a/Source/driver/Castro.cpp +++ b/Source/driver/Castro.cpp @@ -245,14 +245,14 @@ Castro::read_params () { if (dgeom.isPeriodic(dir)) { - if (lo_bc[dir] != Interior) + if (lo_bc[dir] != amrex::PhysBCType::interior) { std::cerr << "Castro::read_params:periodic in direction " << dir << " but low BC is not Interior\n"; amrex::Error(); } - if (hi_bc[dir] != Interior) + if (hi_bc[dir] != amrex::PhysBCType::interior) { std::cerr << "Castro::read_params:periodic in direction " << dir @@ -269,14 +269,14 @@ Castro::read_params () // for (int dir=0; dir HISTORY_SIZE) { amrex::Error("error in riemanncg: cg_maxiter > HISTORY_SIZE"); @@ -978,13 +984,13 @@ Castro::initData () // make sure dx = dy = dz -- that's all we guarantee to support #if (AMREX_SPACEDIM == 2) const Real SMALL = 1.e-13; - if (fabs(dx[0] - dx[1]) > SMALL*dx[0]) + if (std::abs(dx[0] - dx[1]) > SMALL*dx[0]) { amrex::Abort("We don't support dx != dy"); } #elif (AMREX_SPACEDIM == 3) const Real SMALL = 1.e-13; - if ( (fabs(dx[0] - dx[1]) > SMALL*dx[0]) || (fabs(dx[0] - dx[2]) > SMALL*dx[0]) ) + if ( (std::abs(dx[0] - dx[1]) > SMALL*dx[0]) || (std::abs(dx[0] - dx[2]) > SMALL*dx[0]) ) { amrex::Abort("We don't support dx != dy != dz"); } @@ -3545,12 +3551,12 @@ Castro::apply_tagging_restrictions(TagBoxArray& tags, [[maybe_unused]] Real time int boundary_buf = n_error_buf[dim] + blocking_factor[dim] / ref_ratio[dim]; - if ((physbc_lo[dim] != Symmetry && physbc_lo[dim] != Interior) && + if ((physbc_lo[dim] != amrex::PhysBCType::symmetry && physbc_lo[dim] != amrex::PhysBCType::interior) && (idx[dim] <= domlo[dim] + boundary_buf)) { outer_boundary_test[dim] = true; } - if ((physbc_hi[dim] != Symmetry && physbc_lo[dim] != Interior) && + if ((physbc_hi[dim] != amrex::PhysBCType::symmetry && physbc_lo[dim] != amrex::PhysBCType::interior) && (idx[dim] >= domhi[dim] - boundary_buf)) { outer_boundary_test[dim] = true; } @@ -4209,12 +4215,12 @@ Castro::get_numpts () #elif (AMREX_SPACEDIM == 2) long ny = bx.size()[1]; Real ndiagsq = Real(nx*nx + ny*ny); - numpts_1d = int(sqrt(ndiagsq))+2*NUM_GROW; + numpts_1d = int(std::sqrt(ndiagsq))+2*NUM_GROW; #elif (AMREX_SPACEDIM == 3) long ny = bx.size()[1]; long nz = bx.size()[2]; Real ndiagsq = Real(nx*nx + ny*ny + nz*nz); - numpts_1d = int(sqrt(ndiagsq))+2*NUM_GROW; + numpts_1d = int(std::sqrt(ndiagsq))+2*NUM_GROW; #endif if (verbose && ParallelDescriptor::IOProcessor()) { diff --git a/Source/driver/Castro_setup.cpp b/Source/driver/Castro_setup.cpp index 9506fdde2c..2e3a3c1bfc 100644 --- a/Source/driver/Castro_setup.cpp +++ b/Source/driver/Castro_setup.cpp @@ -31,23 +31,23 @@ typedef StateDescriptor::BndryFunc BndryFunc; // static int scalar_bc[] = { - INT_DIR, EXT_DIR, FOEXTRAP, REFLECT_EVEN, REFLECT_EVEN, REFLECT_EVEN + amrex::BCType::int_dir, amrex::BCType::ext_dir, amrex::BCType::foextrap, amrex::BCType::reflect_even, amrex::BCType::reflect_even, amrex::BCType::reflect_even }; static int norm_vel_bc[] = { - INT_DIR, EXT_DIR, FOEXTRAP, REFLECT_ODD, REFLECT_ODD, REFLECT_ODD + amrex::BCType::int_dir, amrex::BCType::ext_dir, amrex::BCType::foextrap, amrex::BCType::reflect_odd, amrex::BCType::reflect_odd, amrex::BCType::reflect_odd }; static int tang_vel_bc[] = { - INT_DIR, EXT_DIR, FOEXTRAP, REFLECT_EVEN, REFLECT_EVEN, REFLECT_EVEN + amrex::BCType::int_dir, amrex::BCType::ext_dir, amrex::BCType::foextrap, amrex::BCType::reflect_even, amrex::BCType::reflect_even, amrex::BCType::reflect_even }; #ifdef MHD static int mag_field_bc[] = { - INT_DIR, EXT_DIR, FOEXTRAP, REFLECT_EVEN, FOEXTRAP, HOEXTRAP + amrex::BCType::int_dir, amrex::BCType::ext_dir, amrex::BCType::foextrap, amrex::BCType::reflect_even, amrex::BCType::foextrap, amrex::BCType::hoextrap }; #endif @@ -146,11 +146,11 @@ void replace_inflow_bc (BCRec& bc) { for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { - if (bc.lo(dir) == EXT_DIR) { - bc.setLo(dir, FOEXTRAP); + if (bc.lo(dir) == amrex::BCType::ext_dir) { + bc.setLo(dir, amrex::BCType::foextrap); } - if (bc.hi(dir) == EXT_DIR) { - bc.setHi(dir, FOEXTRAP); + if (bc.hi(dir) == amrex::BCType::ext_dir) { + bc.setHi(dir, amrex::BCType::foextrap); } } } diff --git a/Source/driver/MGutils.cpp b/Source/driver/MGutils.cpp index 62e5bbd85c..2a7455a6fc 100644 --- a/Source/driver/MGutils.cpp +++ b/Source/driver/MGutils.cpp @@ -20,7 +20,7 @@ apply_metric(const Box& bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - IntVect idx(D_DECL(i, j, k)); + IntVect idx(AMREX_D_DECL(i, j, k)); // at centers if (rbx.contains(idx)) { diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters index f4b4c41f13..5fbba137f4 100644 --- a/Source/driver/_cpp_parameters +++ b/Source/driver/_cpp_parameters @@ -422,6 +422,9 @@ disable_shock_burning int 0 # shock detection threshold for grad{P} / P shock_detection_threshold Real 0.6666666666666666666666_rt +# do we subtract off the hydrostatic pressure when evaluating a shock? +shock_detection_include_sources int 1 + # initial guess for the temperature when inverting the EoS (e.g. when # calling eos_input_re) T_guess Real 1.e8 @@ -621,12 +624,12 @@ max_tagging_radius real 10.0e0 (v, verbose) int 0 # do we dump the old state into the checkpoint files too? -dump_old bool false +dump_old int 0 # do we assume the domain is plane parallel when computing some of the derived # quantities (e.g. radial velocity). Note: this will always assume that the # last spatial dimension is vertical -domain_is_plane_parallel bool false +domain_is_plane_parallel int 0 # display information about updates to the state (how much mass, momentum, energy added) print_update_diagnostics int (0, 1) diff --git a/Source/driver/sum_integrated_quantities.cpp b/Source/driver/sum_integrated_quantities.cpp index 6dc1be4004..2398480765 100644 --- a/Source/driver/sum_integrated_quantities.cpp +++ b/Source/driver/sum_integrated_quantities.cpp @@ -131,7 +131,7 @@ Castro::sum_integrated_quantities () auto R = R_new[mfi].array(); #endif - auto level_mask = mask_available ? mask_mf[mfi].array() : Array4{}; + const auto & level_mask = mask_available ? mask_mf[mfi].array() : Array4{}; reduce_op.eval(bx, reduce_data, [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple diff --git a/Source/driver/sum_utils.cpp b/Source/driver/sum_utils.cpp index f0fa10e35e..438c744f0f 100644 --- a/Source/driver/sum_utils.cpp +++ b/Source/driver/sum_utils.cpp @@ -517,7 +517,7 @@ Castro::gwstrain (Real time, // Standard Kronecker delta. - Real delta[3][3] = {0.0}; + Real delta[3][3] = {{0.0}}; for (int i = 0; i < 3; ++i) { delta[i][i] = 1.0; @@ -539,7 +539,7 @@ Castro::gwstrain (Real time, // Projection operator onto the unit vector n. - Real proj[3][3][3][3] = {0.0}; + Real proj[3][3][3][3] = {{0.0}}; for (int l = 0; l < 3; ++l) { for (int k = 0; k < 3; ++k) { @@ -554,7 +554,7 @@ Castro::gwstrain (Real time, // Now we can calculate the strain tensor. - Real h[3][3] = {0.0}; + Real h[3][3] = {{0.0}}; for (int l = 0; l < 3; ++l) { for (int k = 0; k < 3; ++k) { diff --git a/Source/driver/timestep.cpp b/Source/driver/timestep.cpp index 42ac0325b2..a7b5936e4e 100644 --- a/Source/driver/timestep.cpp +++ b/Source/driver/timestep.cpp @@ -46,7 +46,7 @@ Castro::estdt_cfl (int is_new) Array4 const& u = ua[box_no]; - IntVect idx(D_DECL(i,j,k)); + IntVect idx(AMREX_D_DECL(i,j,k)); Real rhoInv = 1.0_rt / u(i,j,k,URHO); @@ -150,7 +150,7 @@ Castro::estdt_mhd (int is_new) Array4 const& by_arr = bya[box_no]; Array4 const& bz_arr = bza[box_no]; - IntVect idx(D_DECL(i,j,k)); + IntVect idx(AMREX_D_DECL(i,j,k)); Real rhoInv = 1.0_rt / u_arr(i,j,k,URHO); Real bcx = 0.5_rt * (bx_arr(i+1,j,k) + bx_arr(i,j,k)); @@ -253,7 +253,7 @@ Castro::estdt_temp_diffusion (int is_new) Array4 const& ustate = ua[box_no]; - IntVect idx(D_DECL(i,j,k)); + IntVect idx(AMREX_D_DECL(i,j,k)); if (ustate(i,j,k,URHO) > ldiffuse_cutoff_density) { @@ -315,7 +315,7 @@ Castro::estdt_burning (int is_new) { if (castro::dtnuc_e > 1.e199_rt && castro::dtnuc_X > 1.e199_rt) { - IntVect idx(D_DECL(0,0,0)); + IntVect idx(AMREX_D_DECL(0,0,0)); return {ValLocPair{1.e200_rt, idx}}; } @@ -331,7 +331,7 @@ Castro::estdt_burning (int is_new) Array4 const& S = ua[box_no]; - IntVect idx(D_DECL(i,j,k)); + IntVect idx(AMREX_D_DECL(i,j,k)); // Set a floor on the minimum size of a derivative. This floor // is small enough such that it will result in no timestep limiting. @@ -368,7 +368,7 @@ Castro::estdt_burning (int is_new) #if AMREX_SPACEDIM == 1 burn_state.dx = dx[0]; #else - burn_state.dx = amrex::min(D_DECL(dx[0], dx[1], dx[2])); + burn_state.dx = amrex::min(AMREX_D_DECL(dx[0], dx[1], dx[2])); #endif burn_state.rho = S(i,j,k,URHO); diff --git a/Source/gravity/Castro_gravity.cpp b/Source/gravity/Castro_gravity.cpp index 46eccda99f..7b300310ec 100644 --- a/Source/gravity/Castro_gravity.cpp +++ b/Source/gravity/Castro_gravity.cpp @@ -386,7 +386,7 @@ void Castro::construct_old_gravity_source(MultiFab& source, MultiFab& state_in, } #endif - Real SrE; + Real SrE{}; if (castro::grav_source_type == 1 || castro::grav_source_type == 2) { @@ -584,7 +584,7 @@ void Castro::construct_new_gravity_source(MultiFab& source, MultiFab& state_old, // Correct energy - Real SrEcorr; + Real SrEcorr{}; if (castro::grav_source_type == 1) { diff --git a/Source/gravity/Gravity.H b/Source/gravity/Gravity.H index c2c4ca8ebf..af7ff177fa 100644 --- a/Source/gravity/Gravity.H +++ b/Source/gravity/Gravity.H @@ -449,7 +449,6 @@ public: amrex::Vector area; int Density; - int finest_level; int finest_level_allocated; amrex::BCRec* phys_bc; diff --git a/Source/gravity/Gravity.cpp b/Source/gravity/Gravity.cpp index 6665b0ef38..8ef59f0e46 100644 --- a/Source/gravity/Gravity.cpp +++ b/Source/gravity/Gravity.cpp @@ -97,6 +97,7 @@ Gravity::Gravity(Amr* Parent, int _finest_level, BCRec* _phys_bc, int _density) init_multipole_grav(); } max_rhs = 0.0; + numpts_at_level = -1; } Gravity::~Gravity() = default; @@ -1569,7 +1570,7 @@ Gravity::compute_radial_mass(const Box& bx, r = std::sqrt(xxsq + yysq + zzsq); index = static_cast(r * drinv); - Real vol_frac; + Real vol_frac{}; if (coord_type == 0) { @@ -1655,7 +1656,7 @@ Gravity::init_multipole_grav() for (int b = 0; b < AMREX_SPACEDIM; ++b) { - if ((lo_bc[b] == Symmetry) && (parent->Geom(0).Coord() == 0)) { + if ((lo_bc[b] == amrex::PhysBCType::symmetry) && (parent->Geom(0).Coord() == 0)) { if (std::abs(problem::center[b] - problo[b]) < edgeTolerance) { multipole::volumeFactor *= 2.0_rt; multipole::doReflectionLo(b) = true; @@ -1666,7 +1667,7 @@ Gravity::init_multipole_grav() } } - if ((hi_bc[b] == Symmetry) && (parent->Geom(0).Coord() == 0)) { + if ((hi_bc[b] == amrex::PhysBCType::symmetry) && (parent->Geom(0).Coord() == 0)) { if (std::abs(problem::center[b] - probhi[b]) < edgeTolerance) { multipole::volumeFactor *= 2.0_rt; multipole::doReflectionHi(b) = true; @@ -1800,9 +1801,9 @@ Gravity::fill_multipole_BCs(int crse_level, int fine_level, const Vectorlo(idim) == Symmetry) { + if (phys_bc->lo(idim) == amrex::PhysBCType::symmetry) { mlmg_lobc[idim] = MLLinOp::BCType::Neumann; } else { mlmg_lobc[idim] = MLLinOp::BCType::Dirichlet; } - if (phys_bc->hi(idim) == Symmetry) { + if (phys_bc->hi(idim) == amrex::PhysBCType::symmetry) { mlmg_hibc[idim] = MLLinOp::BCType::Neumann; } else { mlmg_hibc[idim] = MLLinOp::BCType::Dirichlet; @@ -3524,10 +3525,10 @@ Gravity::sanity_check (int level) { if (!geom.isPeriodic(dir)) { - if (phys_bc->lo(dir) != Symmetry) { + if (phys_bc->lo(dir) != amrex::PhysBCType::symmetry) { shrunk_domain.growLo(dir,-1); } - if (phys_bc->hi(dir) != Symmetry) { + if (phys_bc->hi(dir) != amrex::PhysBCType::symmetry) { shrunk_domain.growHi(dir,-1); } } diff --git a/Source/hydro/Castro_ctu_hydro.cpp b/Source/hydro/Castro_ctu_hydro.cpp index 96d64c4f55..c1a29ae0c2 100644 --- a/Source/hydro/Castro_ctu_hydro.cpp +++ b/Source/hydro/Castro_ctu_hydro.cpp @@ -283,7 +283,7 @@ Castro::construct_ctu_hydro_source(Real time, Real dt) }); if (hybrid_riemann == 1 || compute_shock) { - shock(obx, q_arr, src_q_arr, shk_arr); + shock(obx, q_arr, old_src_arr, shk_arr); } else { amrex::ParallelFor(obx, diff --git a/Source/hydro/Castro_hydro.H b/Source/hydro/Castro_hydro.H index 647abc10ac..94093c6392 100644 --- a/Source/hydro/Castro_hydro.H +++ b/Source/hydro/Castro_hydro.H @@ -822,7 +822,7 @@ void do_enforce_minimum_density(const amrex::Box& bx, amrex::Array4 const& state_arr, - const int verbose); + const int verbose_warnings); #ifdef HYBRID_MOMENTUM diff --git a/Source/hydro/Castro_mol.cpp b/Source/hydro/Castro_mol.cpp index 1ea775f2cc..23a004c719 100644 --- a/Source/hydro/Castro_mol.cpp +++ b/Source/hydro/Castro_mol.cpp @@ -29,8 +29,8 @@ Castro::mol_plm_reconstruct(const Box& bx, const int* lo_bc = phys_bc.lo(); const int* hi_bc = phys_bc.hi(); - bool lo_symm = lo_bc[idir] == Symmetry; - bool hi_symm = hi_bc[idir] == Symmetry; + bool lo_symm = lo_bc[idir] == amrex::PhysBCType::symmetry; + bool hi_symm = hi_bc[idir] == amrex::PhysBCType::symmetry; const auto domlo = geom.Domain().loVect3d(); const auto domhi = geom.Domain().hiVect3d(); diff --git a/Source/hydro/Castro_mol_hydro.cpp b/Source/hydro/Castro_mol_hydro.cpp index aadc6a41d8..b73e6b7047 100644 --- a/Source/hydro/Castro_mol_hydro.cpp +++ b/Source/hydro/Castro_mol_hydro.cpp @@ -143,7 +143,7 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) src_q.resize(qbx, NQSRC); Array4 const src_q_arr = src_q.array(); - if (hybrid_riemann == 1 || compute_shock || (sdc_order == 2)) { + if (sdc_order == 2) { amrex::ParallelFor(qbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -162,7 +162,7 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) if (hybrid_riemann == 1 || compute_shock) { - shock(obx, q_arr, src_q_arr, shk_arr); + shock(obx, q_arr, source_in_arr, shk_arr); } else { amrex::ParallelFor(obx, @@ -314,8 +314,8 @@ Castro::construct_mol_hydro_source(Real time, Real dt, MultiFab& A_update) GpuArray lo_periodic; GpuArray hi_periodic; for (int idir = 0; idir < AMREX_SPACEDIM; idir++) { - lo_periodic[idir] = lo_bc[idir] == Interior; - hi_periodic[idir] = hi_bc[idir] == Interior; + lo_periodic[idir] = lo_bc[idir] == amrex::PhysBCType::interior; + hi_periodic[idir] = hi_bc[idir] == amrex::PhysBCType::interior; } amrex::ParallelFor(nbx, NQ, diff --git a/Source/hydro/advection_util.cpp b/Source/hydro/advection_util.cpp index c75bf51268..ec7505a438 100644 --- a/Source/hydro/advection_util.cpp +++ b/Source/hydro/advection_util.cpp @@ -61,7 +61,7 @@ Castro::ctoprim(const Box& bx, void Castro::shock(const Box& bx, Array4 const& q_arr, - Array4 const& q_src_arr, + Array4 const& U_src_arr, Array4 const& shk) { // This is a basic multi-dimensional shock detection algorithm. @@ -136,14 +136,23 @@ Castro::shock(const Box& bx, // We subtract off the hydrostatic force, since the pressure that // balances that is not available to make a shock. // We'll use a centered diff for the pressure gradient. - Real dP_x = 0.5_rt * (q_arr(i+1,j,k,QPRES) - q_arr(i-1,j,k,QPRES)) - dx[0] * q_src_arr(i,j,k,QU); + Real dP_x = 0.5_rt * (q_arr(i+1,j,k,QPRES) - q_arr(i-1,j,k,QPRES)); + if (shock_detection_include_sources == 1) { + dP_x += -0.25_rt * dx[0] * (U_src_arr(i+1,j,k,UMX) + 2.0_rt * U_src_arr(i,j,k,UMX) + U_src_arr(i-1,j,k,UMX)); + } Real dP_y = 0.0_rt; Real dP_z = 0.0_rt; #if AMREX_SPACEDIM >= 2 - dP_y = 0.5_rt * (q_arr(i,j+1,k,QPRES) - q_arr(i,j-1,k,QPRES)) - dx[1] * q_src_arr(i,j,k,QV); + dP_y = 0.5_rt * (q_arr(i,j+1,k,QPRES) - q_arr(i,j-1,k,QPRES)); + if (shock_detection_include_sources == 1) { + dP_y += -0.25_rt * dx[1] * (U_src_arr(i,j+1,k,UMY) + 2.0_rt * U_src_arr(i,j,k,UMY) + U_src_arr(i,j-1,k,UMY)); + } #endif #if AMREX_SPACEDIM == 3 - dP_z = 0.5_rt * (q_arr(i,j,k+1,QPRES) - q_arr(i,j,k-1,QPRES)) - dx[2] * q_src_arr(i,j,k,QW); + dP_z = 0.5_rt * (q_arr(i,j,k+1,QPRES) - q_arr(i,j,k-1,QPRES)); + if (shock_detection_include_sources == 1) { + dP_z += -0.25_rt * dx[2] * (U_src_arr(i,j,k+1,UMZ) + 2.0_rt * U_src_arr(i,j,k,UMZ) + U_src_arr(i,j,k-1,UMZ)); + } #endif //Real gradPdx_over_P = std::sqrt(dP_x * dP_x + dP_y * dP_y + dP_z * dP_z) / q_arr(i,j,k,QPRES); @@ -151,10 +160,13 @@ Castro::shock(const Box& bx, q_arr(i,j,k,QV) * q_arr(i,j,k,QV) + q_arr(i,j,k,QW) * q_arr(i,j,k,QW)); - Real gradPdx_over_P = std::abs(dP_x * q_arr(i,j,k,QU) + - dP_y * q_arr(i,j,k,QV) + - dP_z * q_arr(i,j,k,QW)) / vel; - gradPdx_over_P /= (q_arr(i,j,k,QPRES) / std::max(dx[0], std::max(dx[1], dx[2]))); + Real gradPdx_over_P{0.0_rt}; + if (vel != 0.0) { + gradPdx_over_P = std::abs(dP_x * q_arr(i,j,k,QU) + + dP_y * q_arr(i,j,k,QV) + + dP_z * q_arr(i,j,k,QW)) / vel; + } + gradPdx_over_P /= q_arr(i,j,k,QPRES); if (gradPdx_over_P > castro::shock_detection_threshold && div_u < 0.0_rt) { shk(i,j,k) = 1.0_rt; @@ -576,13 +588,13 @@ Castro::limit_hydro_fluxes_on_small_dens(const Box& bx, void // NOLINTNEXTLINE(readability-convert-member-functions-to-static) Castro::do_enforce_minimum_density(const Box& bx, Array4 const& state_arr, - const int verbose) { + const int verbose_warnings) { #ifdef HYBRID_MOMENTUM GeometryData geomdata = geom.data(); #endif - amrex::ignore_unused(verbose); + amrex::ignore_unused(verbose_warnings); amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -591,8 +603,8 @@ Castro::do_enforce_minimum_density(const Box& bx, if (state_arr(i,j,k,URHO) < small_dens) { #ifndef AMREX_USE_GPU - if (verbose > 1 || - (verbose > 0 && state_arr(i,j,k,URHO) > castro::retry_small_density_cutoff)) { + if (verbose_warnings > 1 || + (verbose_warnings > 0 && state_arr(i,j,k,URHO) > castro::retry_small_density_cutoff)) { std::cout << " " << std::endl; if (state_arr(i,j,k,URHO) < 0.0_rt) { std::cout << ">>> RESETTING NEG. DENSITY AT " << i << ", " << j << ", " << k << std::endl; @@ -682,8 +694,8 @@ Castro::enforce_reflect_states(const Box& bx, const int idir, const auto domlo = geom.Domain().loVect3d(); const auto domhi = geom.Domain().hiVect3d(); - bool lo_bc_test = lo_bc[idir] == Symmetry; - bool hi_bc_test = hi_bc[idir] == Symmetry; + bool lo_bc_test = lo_bc[idir] == amrex::PhysBCType::symmetry; + bool hi_bc_test = hi_bc[idir] == amrex::PhysBCType::symmetry; // normal velocity const int QUN = QU + idir; diff --git a/Source/hydro/fourth_center_average.cpp b/Source/hydro/fourth_center_average.cpp index 5c4566b615..070fcd8f3e 100644 --- a/Source/hydro/fourth_center_average.cpp +++ b/Source/hydro/fourth_center_average.cpp @@ -35,8 +35,8 @@ Castro::make_cell_center(const Box& bx, GpuArray lo_periodic; GpuArray hi_periodic; for (int idir = 0; idir < AMREX_SPACEDIM; idir++) { - lo_periodic[idir] = lo_bc[idir] == Interior; - hi_periodic[idir] = hi_bc[idir] == Interior; + lo_periodic[idir] = lo_bc[idir] == amrex::PhysBCType::interior; + hi_periodic[idir] = hi_bc[idir] == amrex::PhysBCType::interior; } amrex::ParallelFor(bx, U.nComp(), @@ -68,8 +68,8 @@ Castro::make_cell_center_in_place(const Box& bx, GpuArray lo_periodic; GpuArray hi_periodic; for (int idir = 0; idir < AMREX_SPACEDIM; idir++) { - lo_periodic[idir] = lo_bc[idir] == Interior; - hi_periodic[idir] = hi_bc[idir] == Interior; + lo_periodic[idir] = lo_bc[idir] == amrex::PhysBCType::interior; + hi_periodic[idir] = hi_bc[idir] == amrex::PhysBCType::interior; } for (int n = 0; n < U.nComp(); n++) { @@ -105,8 +105,8 @@ Castro::compute_lap_term(const Box& bx, GpuArray lo_periodic; GpuArray hi_periodic; for (int idir = 0; idir < AMREX_SPACEDIM; idir++) { - lo_periodic[idir] = lo_bc[idir] == Interior; - hi_periodic[idir] = hi_bc[idir] == Interior; + lo_periodic[idir] = lo_bc[idir] == amrex::PhysBCType::interior; + hi_periodic[idir] = hi_bc[idir] == amrex::PhysBCType::interior; } amrex::ParallelFor(bx, @@ -136,8 +136,8 @@ Castro::make_fourth_average(const Box& bx, GpuArray lo_periodic; GpuArray hi_periodic; for (int idir = 0; idir < AMREX_SPACEDIM; idir++) { - lo_periodic[idir] = lo_bc[idir] == Interior; - hi_periodic[idir] = hi_bc[idir] == Interior; + lo_periodic[idir] = lo_bc[idir] == amrex::PhysBCType::interior; + hi_periodic[idir] = hi_bc[idir] == amrex::PhysBCType::interior; } amrex::ParallelFor(bx, q.nComp(), @@ -186,8 +186,8 @@ Castro::make_fourth_in_place_n(const Box& bx, GpuArray lo_periodic; GpuArray hi_periodic; for (int idir = 0; idir < AMREX_SPACEDIM; idir++) { - lo_periodic[idir] = lo_bc[idir] == Interior; - hi_periodic[idir] = hi_bc[idir] == Interior; + lo_periodic[idir] = lo_bc[idir] == amrex::PhysBCType::interior; + hi_periodic[idir] = hi_bc[idir] == amrex::PhysBCType::interior; } amrex::ParallelFor(bx, diff --git a/Source/hydro/fourth_order.cpp b/Source/hydro/fourth_order.cpp index e8f563fe21..89183a5f27 100644 --- a/Source/hydro/fourth_order.cpp +++ b/Source/hydro/fourth_order.cpp @@ -37,25 +37,25 @@ Castro::fourth_interfaces(const Box& bx, // interpolate to the edges -- this is a_{i-1/2} // note for non-periodic physical boundaries, we use a special stencil - if (i == domlo[0]+1 && lo_bc[0] == Symmetry) { + if (i == domlo[0]+1 && lo_bc[0] == amrex::PhysBCType::symmetry) { // use a stencil for the interface that is one zone // from the left physical boundary, MC Eq. 22 a_int(i,j,k) = (1.0_rt/12.0_rt)*(3.0_rt*a(i-1,j,k,ncomp) + 13.0_rt*a(i,j,k,ncomp) - 5.0_rt*a(i+1,j,k,ncomp) + a(i+2,j,k,ncomp)); - } else if (i == domlo[0] && lo_bc[0] == Symmetry) { + } else if (i == domlo[0] && lo_bc[0] == amrex::PhysBCType::symmetry) { // use a stencil for when the interface is on the // left physical boundary MC Eq. 21 a_int(i,j,k) = (1.0_rt/12.0_rt)*(25.0_rt*a(i,j,k,ncomp) - 23.0_rt*a(i+1,j,k,ncomp) + 13.0_rt*a(i+2,j,k,ncomp) - 3.0_rt*a(i+3,j,k,ncomp)); - } else if (i == domhi[0] && hi_bc[0] == Symmetry) { + } else if (i == domhi[0] && hi_bc[0] == amrex::PhysBCType::symmetry) { // use a stencil for the interface that is one zone // from the right physical boundary, MC Eq. 22 a_int(i,j,k) = (1.0_rt/12.0_rt)*(3.0_rt*a(i,j,k,ncomp) + 13.0_rt*a(i-1,j,k,ncomp) - 5.0_rt*a(i-2,j,k,ncomp) + a(i-3,j,k,ncomp)); - } else if (i == domhi[0]+1 && hi_bc[0] == Symmetry) { + } else if (i == domhi[0]+1 && hi_bc[0] == amrex::PhysBCType::symmetry) { // use a stencil for when the interface is on the // right physical boundary MC Eq. 21 a_int(i,j,k) = (1.0_rt/12.0_rt)*(25.0_rt*a(i-1,j,k,ncomp) - 23.0_rt*a(i-2,j,k,ncomp) + @@ -78,25 +78,25 @@ Castro::fourth_interfaces(const Box& bx, // interpolate to the edges - if (j == domlo[1]+1 && lo_bc[1] == Symmetry) { + if (j == domlo[1]+1 && lo_bc[1] == amrex::PhysBCType::symmetry) { // use a stencil for the interface that is one zone // from the left physical boundary, MC Eq. 22 a_int(i,j,k) = (1.0_rt/12.0_rt)*(3.0_rt*a(i,j-1,k,ncomp) + 13.0_rt*a(i,j,k,ncomp) - 5.0_rt*a(i,j+1,k,ncomp) + a(i,j+2,k,ncomp)); - } else if (j == domlo[1] && lo_bc[1] == Symmetry) { + } else if (j == domlo[1] && lo_bc[1] == amrex::PhysBCType::symmetry) { // use a stencil for when the interface is on the // left physical boundary MC Eq. 21 a_int(i,j,k) = (1.0_rt/12.0_rt)*(25.0_rt*a(i,j,k,ncomp) - 23.0_rt*a(i,j+1,k,ncomp) + 13.0_rt*a(i,j+2,k,ncomp) - 3.0_rt*a(i,j+3,k,ncomp)); - } else if (j == domhi[1] && hi_bc[1] == Symmetry) { + } else if (j == domhi[1] && hi_bc[1] == amrex::PhysBCType::symmetry) { // use a stencil for the interface that is one zone // from the right physical boundary, MC Eq. 22 a_int(i,j,k) = (1.0_rt/12.0_rt)*(3.0_rt*a(i,j,k,ncomp) + 13.0_rt*a(i,j-1,k,ncomp) - 5.0_rt*a(i,j-2,k,ncomp) + a(i,j-3,k,ncomp)); - } else if (j == domhi[1]+1 && hi_bc[1] == Symmetry) { + } else if (j == domhi[1]+1 && hi_bc[1] == amrex::PhysBCType::symmetry) { // use a stencil for when the interface is on the // right physical boundary MC Eq. 21 a_int(i,j,k) = (1.0_rt/12.0_rt)*(25.0_rt*a(i,j-1,k,ncomp) - 23.0_rt*a(i,j-2,k,ncomp) + @@ -120,25 +120,25 @@ Castro::fourth_interfaces(const Box& bx, // interpolate to the edges - if (k == domlo[2]+1 && lo_bc[2] == Symmetry) { + if (k == domlo[2]+1 && lo_bc[2] == amrex::PhysBCType::symmetry) { // use a stencil for the interface that is one zone // from the left physical boundary, MC Eq. 22 a_int(i,j,k) = (1.0_rt/12.0_rt)*(3.0_rt*a(i,j,k-1,ncomp) + 13.0_rt*a(i,j,k,ncomp) - 5.0_rt*a(i,j,k+1,ncomp) + a(i,j,k+2,ncomp)); - } else if (k == domlo[2] && lo_bc[2] == Symmetry) { + } else if (k == domlo[2] && lo_bc[2] == amrex::PhysBCType::symmetry) { // use a stencil for when the interface is on the // left physical boundary MC Eq. 21 a_int(i,j,k) = (1.0_rt/12.0_rt)*(25.0_rt*a(i,j,k,ncomp) - 23.0_rt*a(i,j,k+1,ncomp) + 13.0_rt*a(i,j,k+2,ncomp) - 3.0_rt*a(i,j,k+3,ncomp)); - } else if (k == domhi[2] && hi_bc[2] == Symmetry) { + } else if (k == domhi[2] && hi_bc[2] == amrex::PhysBCType::symmetry) { // use a stencil for the interface that is one zone // from the right physical boundary, MC Eq. 22 a_int(i,j,k) = (1.0_rt/12.0_rt)*(3.0_rt*a(i,j,k,ncomp) + 13.0_rt*a(i,j,k-1,ncomp) - 5.0_rt*a(i,j,k-2,ncomp) + a(i,j,k-3,ncomp)); - } else if (k == domhi[2]+1 && hi_bc[2] == Symmetry) { + } else if (k == domhi[2]+1 && hi_bc[2] == amrex::PhysBCType::symmetry) { // use a stencil for when the interface is on the // right physical boundary MC Eq. 21 a_int(i,j,k) = (1.0_rt/12.0_rt)*(25.0_rt*a(i,j,k-1,ncomp) - 23.0_rt*a(i,j,k-2,ncomp) + @@ -300,17 +300,17 @@ Castro::states(const Box& bx, // reset the left state at domlo[0] if needed -- it is outside the domain - if (lo_bc[0] == Outflow) { + if (lo_bc[0] == amrex::PhysBCType::outflow) { //al(domlo[0],j,k,ncomp) = ar(domlo[0],j,k,ncomp); - } else if (lo_bc[0] == Symmetry) { + } else if (lo_bc[0] == amrex::PhysBCType::symmetry) { if (ncomp == QU) { al(domlo[0],j,k,QU) = -ar(domlo[0],j,k,QU); } else { al(domlo[0],j,k,ncomp) = ar(domlo[0],j,k,ncomp); } - } else if (lo_bc[0] == Interior) { + } else if (lo_bc[0] == amrex::PhysBCType::interior) { // we don't need to do anything here } else { @@ -324,17 +324,17 @@ Castro::states(const Box& bx, // reset the right state at domhi[0]+1 if needed -- it is outside the domain - if (hi_bc[0] == Outflow) { + if (hi_bc[0] == amrex::PhysBCType::outflow) { //ar(domhi[0]+1,j,k,ncomp) = al(domhi[0]+1,j,k,ncomp); - } else if (hi_bc[0] == Symmetry) { + } else if (hi_bc[0] == amrex::PhysBCType::symmetry) { if (ncomp == QU) { ar(domhi[0]+1,j,k,QU) = -al(domhi[0]+1,j,k,QU); } else { ar(domhi[0]+1,j,k,ncomp) = al(domhi[0]+1,j,k,ncomp); } - } else if (hi_bc[0] == Interior) { + } else if (hi_bc[0] == amrex::PhysBCType::interior) { // we don't need to do anything here } else { @@ -465,17 +465,17 @@ Castro::states(const Box& bx, // reset the left state at domlo[1] if needed -- it is outside the domain - if (lo_bc[1] == Outflow) { + if (lo_bc[1] == amrex::PhysBCType::outflow) { //al(i,domlo[1],k,ncomp) = ar(i,domlo[1],k,ncomp); - } else if (lo_bc[1] == Symmetry) { + } else if (lo_bc[1] == amrex::PhysBCType::symmetry) { if (ncomp == QV) { al(i,domlo[1],k,QV) = -ar(i,domlo[1],k,QV); } else { al(i,domlo[1],k,ncomp) = ar(i,domlo[1],k,ncomp); } - } else if (lo_bc[1] == Interior) { + } else if (lo_bc[1] == amrex::PhysBCType::interior) { // we don't need to do anything here } else { @@ -489,17 +489,17 @@ Castro::states(const Box& bx, // reset the right state at domhi[1]+1 if needed -- it is outside the domain - if (hi_bc[1] == Outflow) { + if (hi_bc[1] == amrex::PhysBCType::outflow) { //ar(i,domhi[1]+1,k,ncomp) = al(i,domhi[1]+1,k,ncomp); - } else if (hi_bc[1] == Symmetry) { + } else if (hi_bc[1] == amrex::PhysBCType::symmetry) { if (ncomp == QV) { ar(i,domhi[1]+1,k,QV) = -al(i,domhi[1]+1,k,QV); } else { ar(i,domhi[1]+1,k,ncomp) = al(i,domhi[1]+1,k,ncomp); } - } else if (hi_bc[1] == Interior) { + } else if (hi_bc[1] == amrex::PhysBCType::interior) { // we don't need to do anything here } else { @@ -628,17 +628,17 @@ Castro::states(const Box& bx, if (k == domlo[2]) { // reset the left state at domlo[2] if needed -- it is outside the domain - if (lo_bc[2] == Outflow) { + if (lo_bc[2] == amrex::PhysBCType::outflow) { //al(i,j,domlo[2],ncomp) = ar(i,j,domlo[2],ncomp); - } else if (lo_bc[2] == Symmetry) { + } else if (lo_bc[2] == amrex::PhysBCType::symmetry) { if (ncomp == QW) { al(i,j,domlo[2],QW) = -ar(i,j,domlo[2],QW); } else { al(i,j,domlo[2],ncomp) = ar(i,j,domlo[2],ncomp); } - } else if (lo_bc[2] == Interior) { + } else if (lo_bc[2] == amrex::PhysBCType::interior) { // we don't need to do anything here } else { @@ -651,17 +651,17 @@ Castro::states(const Box& bx, if (k == domhi[2]+1) { // reset the right state at domhi[2]+1 if needed -- it is outside the domain - if (hi_bc[2] == Outflow) { + if (hi_bc[2] == amrex::PhysBCType::outflow) { //ar(i,j,domhi[2]+1,ncomp) = al(i,j,domhi[2]+1,ncomp); - } else if (hi_bc[2] == Symmetry) { + } else if (hi_bc[2] == amrex::PhysBCType::symmetry) { if (ncomp == QW) { ar(i,j,domhi[2]+1,QW) = -al(i,j,domhi[2]+1,QW); } else { ar(i,j,domhi[2]+1,ncomp) = al(i,j,domhi[2]+1,ncomp); } - } else if (lo_bc[2] == Interior) { + } else if (lo_bc[2] == amrex::PhysBCType::interior) { // we don't need to do anything here } else { diff --git a/Source/hydro/riemann.cpp b/Source/hydro/riemann.cpp index 8a7a1a7c57..96d2c2e6f0 100644 --- a/Source/hydro/riemann.cpp +++ b/Source/hydro/riemann.cpp @@ -51,12 +51,12 @@ Castro::cmpflx_plus_godunov(const Box& bx, const int* hi_bc = phys_bc.hi(); // do we want to force the flux to zero at the boundary? - const bool special_bnd_lo = (lo_bc[idir] == Symmetry || - lo_bc[idir] == SlipWall || - lo_bc[idir] == NoSlipWall); - const bool special_bnd_hi = (hi_bc[idir] == Symmetry || - hi_bc[idir] == SlipWall || - hi_bc[idir] == NoSlipWall); + const bool special_bnd_lo = (lo_bc[idir] == amrex::PhysBCType::symmetry || + lo_bc[idir] == amrex::PhysBCType::slipwall || + lo_bc[idir] == amrex::PhysBCType::noslipwall); + const bool special_bnd_hi = (hi_bc[idir] == amrex::PhysBCType::symmetry || + hi_bc[idir] == amrex::PhysBCType::slipwall || + hi_bc[idir] == amrex::PhysBCType::noslipwall); auto coord = geom.Coord(); diff --git a/Source/hydro/trace_plm.cpp b/Source/hydro/trace_plm.cpp index 9a1e563e6f..8b075b47eb 100644 --- a/Source/hydro/trace_plm.cpp +++ b/Source/hydro/trace_plm.cpp @@ -37,8 +37,8 @@ Castro::trace_plm(const Box& bx, const int idir, const int* lo_bc = phys_bc.lo(); const int* hi_bc = phys_bc.hi(); - bool lo_symm = lo_bc[idir] == Symmetry; - bool hi_symm = hi_bc[idir] == Symmetry; + bool lo_symm = lo_bc[idir] == amrex::PhysBCType::symmetry; + bool hi_symm = hi_bc[idir] == amrex::PhysBCType::symmetry; const auto domlo = geom.Domain().loVect3d(); const auto domhi = geom.Domain().hiVect3d(); diff --git a/Source/mhd/mhd_plm.cpp b/Source/mhd/mhd_plm.cpp index 99c6dc64e7..2747812513 100644 --- a/Source/mhd/mhd_plm.cpp +++ b/Source/mhd/mhd_plm.cpp @@ -26,8 +26,8 @@ Castro::plm(const Box& bx, const int* lo_bc = phys_bc.lo(); const int* hi_bc = phys_bc.hi(); - bool lo_symm = lo_bc[idir] == Symmetry; - bool hi_symm = hi_bc[idir] == Symmetry; + bool lo_symm = lo_bc[idir] == amrex::PhysBCType::symmetry; + bool hi_symm = hi_bc[idir] == amrex::PhysBCType::symmetry; const auto domlo = geom.Domain().loVect3d(); const auto domhi = geom.Domain().hiVect3d(); diff --git a/Source/particles/CastroParticles.cpp b/Source/particles/CastroParticles.cpp index 10e98b078c..5c3adb1152 100644 --- a/Source/particles/CastroParticles.cpp +++ b/Source/particles/CastroParticles.cpp @@ -137,7 +137,7 @@ Castro::ParticleDerive(const std::string& name, // auto derive_dat = ParticleDerive("particle_count",time,ngrow); - IntVect trr(D_DECL(1,1,1)); + IntVect trr(AMREX_D_DECL(1,1,1)); for (int lev = level+1; lev <= parent->finestLevel(); lev++) { diff --git a/Source/problems/Castro_bc_fill_nd.cpp b/Source/problems/Castro_bc_fill_nd.cpp index fda142da03..f06a7f4236 100644 --- a/Source/problems/Castro_bc_fill_nd.cpp +++ b/Source/problems/Castro_bc_fill_nd.cpp @@ -25,11 +25,11 @@ void ca_statefill(Box const& bx, FArrayBox& data, Vector bcr_noinflow{bcr}; for (auto & bc : bcr_noinflow) { for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { - if (bc.lo(dir) == EXT_DIR) { - bc.setLo(dir, FOEXTRAP); + if (bc.lo(dir) == amrex::BCType::ext_dir) { + bc.setLo(dir, amrex::BCType::foextrap); } - if (bc.hi(dir) == EXT_DIR) { - bc.setHi(dir, FOEXTRAP); + if (bc.hi(dir) == amrex::BCType::ext_dir) { + bc.setHi(dir, amrex::BCType::foextrap); } } } @@ -62,39 +62,39 @@ void ca_statefill(Box const& bx, FArrayBox& data, // corners. #if AMREX_SPACEDIM == 2 - if ((bcr[URHO].lo(0) == EXT_DIR && bcr[URHO].lo(1) == EXT_DIR) || - (bcr[URHO].lo(0) == EXT_DIR && bcr[URHO].hi(1) == EXT_DIR) || - (bcr[URHO].hi(0) == EXT_DIR && bcr[URHO].lo(1) == EXT_DIR) || - (bcr[URHO].hi(0) == EXT_DIR && bcr[URHO].hi(1) == EXT_DIR)) { + if ((bcr[URHO].lo(0) == amrex::BCType::ext_dir && bcr[URHO].lo(1) == amrex::BCType::ext_dir) || + (bcr[URHO].lo(0) == amrex::BCType::ext_dir && bcr[URHO].hi(1) == amrex::BCType::ext_dir) || + (bcr[URHO].hi(0) == amrex::BCType::ext_dir && bcr[URHO].lo(1) == amrex::BCType::ext_dir) || + (bcr[URHO].hi(0) == amrex::BCType::ext_dir && bcr[URHO].hi(1) == amrex::BCType::ext_dir)) { amrex::Error("Error: external boundaries meeting at a corner not supported"); } #endif #if AMREX_SPACEDIM == 3 - if ((bcr[URHO].lo(0) == EXT_DIR && // xl, yl, zl corner - (bcr[URHO].lo(1) == EXT_DIR || bcr[URHO].lo(2) == EXT_DIR)) || - (bcr[URHO].lo(1) == EXT_DIR && bcr[URHO].lo(2) == EXT_DIR) || - (bcr[URHO].lo(0) == EXT_DIR && // xl, yr, zl corner - (bcr[URHO].hi(1) == EXT_DIR || bcr[URHO].lo(2) == EXT_DIR)) || - (bcr[URHO].hi(1) == EXT_DIR && bcr[URHO].lo(2) == EXT_DIR) || - (bcr[URHO].lo(0) == EXT_DIR && // xl, yl, zr corner - (bcr[URHO].lo(1) == EXT_DIR || bcr[URHO].hi(2) == EXT_DIR)) || - (bcr[URHO].lo(1) == EXT_DIR && bcr[URHO].hi(2) == EXT_DIR) || - (bcr[URHO].lo(0) == EXT_DIR && // xl, yr, zr corner - (bcr[URHO].hi(1) == EXT_DIR || bcr[URHO].hi(2) == EXT_DIR)) || - (bcr[URHO].hi(1) == EXT_DIR && bcr[URHO].hi(2) == EXT_DIR) || - (bcr[URHO].hi(0) == EXT_DIR && // xr, yl, zl corner - (bcr[URHO].lo(1) == EXT_DIR || bcr[URHO].lo(2) == EXT_DIR)) || - (bcr[URHO].lo(1) == EXT_DIR && bcr[URHO].lo(2) == EXT_DIR) || - (bcr[URHO].hi(0) == EXT_DIR && // xr, yr, zl corner - (bcr[URHO].hi(1) == EXT_DIR || bcr[URHO].lo(2) == EXT_DIR)) || - (bcr[URHO].hi(1) == EXT_DIR && bcr[URHO].lo(2) == EXT_DIR) || - (bcr[URHO].hi(0) == EXT_DIR && // xr, yl, zr corner - (bcr[URHO].lo(1) == EXT_DIR || bcr[URHO].hi(2) == EXT_DIR)) || - (bcr[URHO].lo(1) == EXT_DIR && bcr[URHO].hi(2) == EXT_DIR) || - (bcr[URHO].hi(0) == EXT_DIR && // xr, yr, zr corner - (bcr[URHO].hi(1) == EXT_DIR || bcr[URHO].hi(2) == EXT_DIR)) || - (bcr[URHO].hi(1) == EXT_DIR && bcr[URHO].hi(2) == EXT_DIR)) { + if ((bcr[URHO].lo(0) == amrex::BCType::ext_dir && // xl, yl, zl corner + (bcr[URHO].lo(1) == amrex::BCType::ext_dir || bcr[URHO].lo(2) == amrex::BCType::ext_dir)) || + (bcr[URHO].lo(1) == amrex::BCType::ext_dir && bcr[URHO].lo(2) == amrex::BCType::ext_dir) || + (bcr[URHO].lo(0) == amrex::BCType::ext_dir && // xl, yr, zl corner + (bcr[URHO].hi(1) == amrex::BCType::ext_dir || bcr[URHO].lo(2) == amrex::BCType::ext_dir)) || + (bcr[URHO].hi(1) == amrex::BCType::ext_dir && bcr[URHO].lo(2) == amrex::BCType::ext_dir) || + (bcr[URHO].lo(0) == amrex::BCType::ext_dir && // xl, yl, zr corner + (bcr[URHO].lo(1) == amrex::BCType::ext_dir || bcr[URHO].hi(2) == amrex::BCType::ext_dir)) || + (bcr[URHO].lo(1) == amrex::BCType::ext_dir && bcr[URHO].hi(2) == amrex::BCType::ext_dir) || + (bcr[URHO].lo(0) == amrex::BCType::ext_dir && // xl, yr, zr corner + (bcr[URHO].hi(1) == amrex::BCType::ext_dir || bcr[URHO].hi(2) == amrex::BCType::ext_dir)) || + (bcr[URHO].hi(1) == amrex::BCType::ext_dir && bcr[URHO].hi(2) == amrex::BCType::ext_dir) || + (bcr[URHO].hi(0) == amrex::BCType::ext_dir && // xr, yl, zl corner + (bcr[URHO].lo(1) == amrex::BCType::ext_dir || bcr[URHO].lo(2) == amrex::BCType::ext_dir)) || + (bcr[URHO].lo(1) == amrex::BCType::ext_dir && bcr[URHO].lo(2) == amrex::BCType::ext_dir) || + (bcr[URHO].hi(0) == amrex::BCType::ext_dir && // xr, yr, zl corner + (bcr[URHO].hi(1) == amrex::BCType::ext_dir || bcr[URHO].lo(2) == amrex::BCType::ext_dir)) || + (bcr[URHO].hi(1) == amrex::BCType::ext_dir && bcr[URHO].lo(2) == amrex::BCType::ext_dir) || + (bcr[URHO].hi(0) == amrex::BCType::ext_dir && // xr, yl, zr corner + (bcr[URHO].lo(1) == amrex::BCType::ext_dir || bcr[URHO].hi(2) == amrex::BCType::ext_dir)) || + (bcr[URHO].lo(1) == amrex::BCType::ext_dir && bcr[URHO].hi(2) == amrex::BCType::ext_dir) || + (bcr[URHO].hi(0) == amrex::BCType::ext_dir && // xr, yr, zr corner + (bcr[URHO].hi(1) == amrex::BCType::ext_dir || bcr[URHO].hi(2) == amrex::BCType::ext_dir)) || + (bcr[URHO].hi(1) == amrex::BCType::ext_dir && bcr[URHO].hi(2) == amrex::BCType::ext_dir)) { amrex::Error("Error: external boundaries meeting at a corner not supported"); } #endif diff --git a/Source/problems/ambient_fill.cpp b/Source/problems/ambient_fill.cpp index 14f30ab302..c5cb7cc701 100644 --- a/Source/problems/ambient_fill.cpp +++ b/Source/problems/ambient_fill.cpp @@ -17,22 +17,22 @@ ambient_denfill(const Box& bx, Array4 const& state, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { bool ambient_x_lo = (castro::ambient_fill_dir == 0 || castro::ambient_fill_dir == -1) && - (bc.lo(0) == FOEXTRAP || bc.lo(0) == HOEXTRAP); + (bc.lo(0) == amrex::BCType::foextrap || bc.lo(0) == amrex::BCType::hoextrap); bool ambient_x_hi = (castro::ambient_fill_dir == 0 || castro::ambient_fill_dir == -1) && - (bc.hi(0) == FOEXTRAP || bc.hi(0) == HOEXTRAP); + (bc.hi(0) == amrex::BCType::foextrap || bc.hi(0) == amrex::BCType::hoextrap); #if AMREX_SPACEDIM >= 2 bool ambient_y_lo = (castro::ambient_fill_dir == 1 || castro::ambient_fill_dir == -1) && - (bc.lo(1) == FOEXTRAP || bc.lo(1) == HOEXTRAP); + (bc.lo(1) == amrex::BCType::foextrap || bc.lo(1) == amrex::BCType::hoextrap); bool ambient_y_hi = (castro::ambient_fill_dir == 1 || castro::ambient_fill_dir == -1) && - (bc.hi(1) == FOEXTRAP || bc.hi(1) == HOEXTRAP); + (bc.hi(1) == amrex::BCType::foextrap || bc.hi(1) == amrex::BCType::hoextrap); #endif #if AMREX_SPACEDIM == 3 bool ambient_z_lo = (castro::ambient_fill_dir == 2 || castro::ambient_fill_dir == -1) && - (bc.lo(2) == FOEXTRAP || bc.lo(2) == HOEXTRAP); + (bc.lo(2) == amrex::BCType::foextrap || bc.lo(2) == amrex::BCType::hoextrap); bool ambient_z_hi = (castro::ambient_fill_dir == 2 || castro::ambient_fill_dir == -1) && - (bc.hi(2) == FOEXTRAP || bc.hi(2) == HOEXTRAP); + (bc.hi(2) == amrex::BCType::foextrap || bc.hi(2) == amrex::BCType::hoextrap); #endif if (castro::fill_ambient_bc == 1) { @@ -76,22 +76,22 @@ ambient_fill(const Box& bx, Array4 const& state, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { bool ambient_x_lo = (castro::ambient_fill_dir == 0 || castro::ambient_fill_dir == -1) && - (bcs(URHO).lo(0) == FOEXTRAP || bcs(URHO).lo(0) == HOEXTRAP); + (bcs(URHO).lo(0) == amrex::BCType::foextrap || bcs(URHO).lo(0) == amrex::BCType::hoextrap); bool ambient_x_hi = (castro::ambient_fill_dir == 0 || castro::ambient_fill_dir == -1) && - (bcs(URHO).hi(0) == FOEXTRAP || bcs(URHO).hi(0) == HOEXTRAP); + (bcs(URHO).hi(0) == amrex::BCType::foextrap || bcs(URHO).hi(0) == amrex::BCType::hoextrap); #if AMREX_SPACEDIM >= 2 bool ambient_y_lo = (castro::ambient_fill_dir == 1 || castro::ambient_fill_dir == -1) && - (bcs(URHO).lo(1) == FOEXTRAP || bcs(URHO).lo(1) == HOEXTRAP); + (bcs(URHO).lo(1) == amrex::BCType::foextrap || bcs(URHO).lo(1) == amrex::BCType::hoextrap); bool ambient_y_hi = (castro::ambient_fill_dir == 1 || castro::ambient_fill_dir == -1) && - (bcs(URHO).hi(1) == FOEXTRAP || bcs(URHO).hi(1) == HOEXTRAP); + (bcs(URHO).hi(1) == amrex::BCType::foextrap || bcs(URHO).hi(1) == amrex::BCType::hoextrap); #endif #if AMREX_SPACEDIM == 3 bool ambient_z_lo = (castro::ambient_fill_dir == 2 || castro::ambient_fill_dir == -1) && - (bcs(URHO).lo(2) == FOEXTRAP || bcs(URHO).lo(2) == HOEXTRAP); + (bcs(URHO).lo(2) == amrex::BCType::foextrap || bcs(URHO).lo(2) == amrex::BCType::hoextrap); bool ambient_z_hi = (castro::ambient_fill_dir == 2 || castro::ambient_fill_dir == -1) && - (bcs(URHO).hi(2) == FOEXTRAP || bcs(URHO).hi(2) == HOEXTRAP); + (bcs(URHO).hi(2) == amrex::BCType::foextrap || bcs(URHO).hi(2) == amrex::BCType::hoextrap); #endif if (castro::fill_ambient_bc == 1) { diff --git a/Source/problems/hse_fill.cpp b/Source/problems/hse_fill.cpp index 6d210fdd34..b9d2c6ab77 100644 --- a/Source/problems/hse_fill.cpp +++ b/Source/problems/hse_fill.cpp @@ -36,7 +36,7 @@ hse_fill(const Box& bx, Array4 const& adv, // XLO - if (bcr[URHO].lo(0) == EXT_DIR && lo[0] < domlo[0]) { + if (bcr[URHO].lo(0) == amrex::BCType::ext_dir && lo[0] < domlo[0]) { if (xl_ext_bc_type == EXT_HSE) { @@ -44,8 +44,8 @@ hse_fill(const Box& bx, Array4 const& adv, // to domlo[0]-1, but we want that to be handled by a // single thread on the GPU - Box gbx(IntVect(D_DECL(domlo[0]-1, lo[1], lo[2])), - IntVect(D_DECL(domlo[0]-1, hi[1], hi[2]))); + Box gbx(IntVect(AMREX_D_DECL(domlo[0]-1, lo[1], lo[2])), + IntVect(AMREX_D_DECL(domlo[0]-1, hi[1], hi[2]))); amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -236,12 +236,12 @@ hse_fill(const Box& bx, Array4 const& adv, // XHI - if (bcr[URHO].hi(0) == EXT_DIR && hi[0] > domhi[0]) { + if (bcr[URHO].hi(0) == amrex::BCType::ext_dir && hi[0] > domhi[0]) { if (xr_ext_bc_type == EXT_HSE) { - Box gbx(IntVect(D_DECL(domhi[0]+1, lo[1], lo[2])), - IntVect(D_DECL(domhi[0]+1, hi[1], hi[2]))); + Box gbx(IntVect(AMREX_D_DECL(domhi[0]+1, lo[1], lo[2])), + IntVect(AMREX_D_DECL(domhi[0]+1, hi[1], hi[2]))); amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -432,12 +432,12 @@ hse_fill(const Box& bx, Array4 const& adv, // YLO - if (bcr[URHO].lo(1) == EXT_DIR && lo[1] < domlo[1]) { + if (bcr[URHO].lo(1) == amrex::BCType::ext_dir && lo[1] < domlo[1]) { if (yl_ext_bc_type == EXT_HSE) { - Box gbx(IntVect(D_DECL(lo[0], domlo[1]-1, lo[2])), - IntVect(D_DECL(hi[0], domlo[1]-1, hi[2]))); + Box gbx(IntVect(AMREX_D_DECL(lo[0], domlo[1]-1, lo[2])), + IntVect(AMREX_D_DECL(hi[0], domlo[1]-1, hi[2]))); amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -626,12 +626,12 @@ hse_fill(const Box& bx, Array4 const& adv, // YHI - if (bcr[URHO].hi(1) == EXT_DIR && hi[1] > domhi[1]) { + if (bcr[URHO].hi(1) == amrex::BCType::ext_dir && hi[1] > domhi[1]) { if (yr_ext_bc_type == EXT_HSE) { - Box gbx(IntVect(D_DECL(lo[0], domhi[1]+1, lo[2])), - IntVect(D_DECL(hi[0], domhi[1]+1, hi[2]))); + Box gbx(IntVect(AMREX_D_DECL(lo[0], domhi[1]+1, lo[2])), + IntVect(AMREX_D_DECL(hi[0], domhi[1]+1, hi[2]))); amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -821,12 +821,12 @@ hse_fill(const Box& bx, Array4 const& adv, // ZLO - if (bcr[URHO].lo(2) == EXT_DIR && lo[2] < domlo[2]) { + if (bcr[URHO].lo(2) == amrex::BCType::ext_dir && lo[2] < domlo[2]) { if (zl_ext_bc_type == EXT_HSE) { - Box gbx(IntVect(D_DECL(lo[0], lo[1], domlo[2]-1)), - IntVect(D_DECL(hi[0], hi[1], domlo[2]-1))); + Box gbx(IntVect(AMREX_D_DECL(lo[0], lo[1], domlo[2]-1)), + IntVect(AMREX_D_DECL(hi[0], hi[1], domlo[2]-1))); amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -1012,7 +1012,7 @@ hse_fill(const Box& bx, Array4 const& adv, // ZHI - if (bcr[URHO].hi(2) == EXT_DIR && hi[2] > domhi[2]) { + if (bcr[URHO].hi(2) == amrex::BCType::ext_dir && hi[2] > domhi[2]) { if (zr_ext_bc_type == EXT_HSE) { #ifndef AMREX_USE_GPU diff --git a/Source/problems/problem_bc_fill.H b/Source/problems/problem_bc_fill.H index 29c6815b46..1eeb1589a4 100644 --- a/Source/problems/problem_bc_fill.H +++ b/Source/problems/problem_bc_fill.H @@ -15,8 +15,8 @@ AMREX_GPU_HOST_DEVICE AMREX_INLINE void problem_bc_fill(int i, int j, int k, Array4 const& state, - Real time, Array1D bcs, - GeometryData geomdata) { + Real time, const Array1D& bcs, + const GeometryData& geomdata) { amrex::ignore_unused(i); amrex::ignore_unused(j); diff --git a/Source/radiation/HABEC.H b/Source/radiation/HABEC.H index 3039e759d3..2bd6e7aabb 100644 --- a/Source/radiation/HABEC.H +++ b/Source/radiation/HABEC.H @@ -89,7 +89,7 @@ namespace HABEC Real bfv, bfm, bfm2; - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { if (bho >= 1) { Real h2 = 0.5e0_rt * h; Real th2 = 3.e0_rt * h2; @@ -306,7 +306,7 @@ namespace HABEC else { bct = bctype; } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { if (bho >= 1) { h2 = 0.5e0_rt * h; th2 = 3.e0_rt * h2; @@ -319,12 +319,12 @@ namespace HABEC bfm = bfv; } } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfv = beta * r; bfm = 0.e0_rt; bfm2 = 0.e0_rt; } - else if (bct == LO_MARSHAK) { + else if (bct == AMREX_LO_MARSHAK) { bfv = 2.e0_rt * beta * r; if (bho >= 1) { bfm = 0.375e0_rt * c * bfv; @@ -334,7 +334,7 @@ namespace HABEC bfm = 0.25e0_rt * c * bfv; } } - else if (bct == LO_SANCHEZ_POMRANING) { + else if (bct == AMREX_LO_SANCHEZ_POMRANING) { bfv = 2.e0_rt * beta * r; if (bho >= 1) { bfm = 1.5e0_rt * spa(i,j,k) * c * bfv; @@ -457,7 +457,7 @@ namespace HABEC { if (mask.contains(i+icp,j+jcp,k+kcp)) { if (mask(i+icp,j+jcp,k+kcp) > 0) { - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { Real d_sign = 1.0_rt; if (iep != 0 || jep != 0 || kep != 0) { // right edge @@ -568,7 +568,7 @@ namespace HABEC else { bct = bctype; } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { Real d_sign = 1.0_rt; if (iep != 0 || jep != 0 || kep != 0) { // right edge @@ -576,7 +576,7 @@ namespace HABEC } dterm(i+iep,j+jep,k+kep) = d(i+iep,j+jep,k+kep) * d_sign * (er(i,j,k) - bcval(i+icp,j+jcp,k+kcp)) / (0.5_rt * h + bcl); } - else if (bct == LO_NEUMANN && bcval(i+icp,j+jcp,k+kcp) == 0.0_rt) { + else if (bct == AMREX_LO_NEUMANN && bcval(i+icp,j+jcp,k+kcp) == 0.0_rt) { dterm(i+iep,j+jep,k+kep) = 0.0_rt; } } diff --git a/Source/radiation/HypreABec.cpp b/Source/radiation/HypreABec.cpp index 638d4715bb..d8ab7d5a4a 100644 --- a/Source/radiation/HypreABec.cpp +++ b/Source/radiation/HypreABec.cpp @@ -507,11 +507,11 @@ void HypreABec::hbmat (const Box& bx, Real bfv, bfm; - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { bfv = fac * h / (0.5e0_rt * h + bcl); bfm = bfv - fac; } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfv = beta / h; bfm = -fac; } @@ -691,18 +691,18 @@ void HypreABec::hbmat3 (const Box& bx, bct = bctype; } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { bfv = fac * h / (0.5e0_rt * h + bcl); bfm = bfv * b(i,j,k); } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfm = 0.e0_rt; } - else if (bct == LO_MARSHAK) { + else if (bct == AMREX_LO_MARSHAK) { bfv = 2.e0_rt * beta * r / h; bfm = 0.25e0_rt * c * bfv; } - else if (bct == LO_SANCHEZ_POMRANING) { + else if (bct == AMREX_LO_SANCHEZ_POMRANING) { bfv = 2.e0_rt * beta * r / h; bfm = spa(i,j,k) * c * bfv; } @@ -729,18 +729,18 @@ void HypreABec::hbmat3 (const Box& bx, bct = bctype; } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { bfv = fac * h / (0.5e0_rt * h + bcl); bfm = bfv * b(i+1,j,k); } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfm = 0.e0_rt; } - else if (bct == LO_MARSHAK) { + else if (bct == AMREX_LO_MARSHAK) { bfv = 2.e0_rt * beta * r / h; bfm = 0.25e0_rt * c * bfv; } - else if (bct == LO_SANCHEZ_POMRANING) { + else if (bct == AMREX_LO_SANCHEZ_POMRANING) { bfv = 2.e0_rt * beta * r / h; bfm = spa(i,j,k) * c * bfv; } @@ -766,18 +766,18 @@ void HypreABec::hbmat3 (const Box& bx, bct = bctype; } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { bfv = fac * h / (0.5e0_rt * h + bcl); bfm = bfv * b(i,j,k); } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfm = 0.e0_rt; } - else if (bct == LO_MARSHAK) { + else if (bct == AMREX_LO_MARSHAK) { bfv = 2.e0_rt * beta * r / h; bfm = 0.25e0_rt * c * bfv; } - else if (bct == LO_SANCHEZ_POMRANING) { + else if (bct == AMREX_LO_SANCHEZ_POMRANING) { bfv = 2.e0_rt * beta * r / h; bfm = spa(i,j,k) * c * bfv; } @@ -804,18 +804,18 @@ void HypreABec::hbmat3 (const Box& bx, bct = bctype; } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { bfv = fac * h / (0.5e0_rt * h + bcl); bfm = bfv * b(i,j+1,k); } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfm = 0.e0_rt; } - else if (bct == LO_MARSHAK) { + else if (bct == AMREX_LO_MARSHAK) { bfv = 2.e0_rt * beta * r / h; bfm = 0.25e0_rt * c * bfv; } - else if (bct == LO_SANCHEZ_POMRANING) { + else if (bct == AMREX_LO_SANCHEZ_POMRANING) { bfv = 2.e0_rt * beta * r / h; bfm = spa(i,j,k) * c * bfv; } @@ -842,18 +842,18 @@ void HypreABec::hbmat3 (const Box& bx, bct = bctype; } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { bfv = fac * h / (0.5e0_rt * h + bcl); bfm = bfv * b(i,j,k); } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfm = 0.e0_rt; } - else if (bct == LO_MARSHAK) { + else if (bct == AMREX_LO_MARSHAK) { bfv = 2.e0_rt * beta * r / h; bfm = 0.25e0_rt * c * bfv; } - else if (bct == LO_SANCHEZ_POMRANING) { + else if (bct == AMREX_LO_SANCHEZ_POMRANING) { bfv = 2.e0_rt * beta * r / h; bfm = spa(i,j,k) * c * bfv; } @@ -880,18 +880,18 @@ void HypreABec::hbmat3 (const Box& bx, bct = bctype; } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { bfv = fac * h / (0.5e0_rt * h + bcl); bfm = bfv * b(i,j,k+1); } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfm = 0.e0_rt; } - else if (bct == LO_MARSHAK) { + else if (bct == AMREX_LO_MARSHAK) { bfv = 2.e0_rt * beta * r / h; bfm = 0.25e0_rt * c * bfv; } - else if (bct == LO_SANCHEZ_POMRANING) { + else if (bct == AMREX_LO_SANCHEZ_POMRANING) { bfv = 2.e0_rt * beta * r / h; bfm = spa(i,j,k) * c * bfv; } @@ -1275,7 +1275,7 @@ void HypreABec::hbvec3 (const Box& bx, bct = bctype; } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { if (bho >= 1) { h2 = 0.5e0_rt * h; th2 = 3.e0_rt * h2; @@ -1287,10 +1287,10 @@ void HypreABec::hbvec3 (const Box& bx, bfv = bfv * b(i,j,k); } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfv = beta * r / h; } - else if (bct == LO_MARSHAK || bct == LO_SANCHEZ_POMRANING) { + else if (bct == AMREX_LO_MARSHAK || bct == AMREX_LO_SANCHEZ_POMRANING) { bfv = 2.e0_rt * beta * r / h; } #ifndef AMREX_USE_GPU @@ -1315,7 +1315,7 @@ void HypreABec::hbvec3 (const Box& bx, bct = bctype; } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { if (bho >= 1) { h2 = 0.5e0_rt * h; th2 = 3.e0_rt * h2; @@ -1327,10 +1327,10 @@ void HypreABec::hbvec3 (const Box& bx, bfv = bfv * b(i+1,j,k); } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfv = beta * r / h; } - else if (bct == LO_MARSHAK || bct == LO_SANCHEZ_POMRANING) { + else if (bct == AMREX_LO_MARSHAK || bct == AMREX_LO_SANCHEZ_POMRANING) { bfv = 2.e0_rt * beta * r / h; } #ifndef AMREX_USE_GPU @@ -1355,7 +1355,7 @@ void HypreABec::hbvec3 (const Box& bx, bct = bctype; } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { if (bho >= 1) { h2 = 0.5e0_rt * h; th2 = 3.e0_rt * h2; @@ -1367,10 +1367,10 @@ void HypreABec::hbvec3 (const Box& bx, bfv = bfv * b(i,j,k); } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfv = beta * r / h; } - else if (bct == LO_MARSHAK || bct == LO_SANCHEZ_POMRANING) { + else if (bct == AMREX_LO_MARSHAK || bct == AMREX_LO_SANCHEZ_POMRANING) { bfv = 2.e0_rt * beta * r / h; } #ifndef AMREX_USE_GPU @@ -1395,7 +1395,7 @@ void HypreABec::hbvec3 (const Box& bx, bct = bctype; } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { if (bho >= 1) { h2 = 0.5e0_rt * h; th2 = 3.e0_rt * h2; @@ -1407,10 +1407,10 @@ void HypreABec::hbvec3 (const Box& bx, bfv = bfv * b(i,j+1,k); } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfv = beta * r / h; } - else if (bct == LO_MARSHAK || bct == LO_SANCHEZ_POMRANING) { + else if (bct == AMREX_LO_MARSHAK || bct == AMREX_LO_SANCHEZ_POMRANING) { bfv = 2.e0_rt * beta * r / h; } #ifndef AMREX_USE_GPU @@ -1435,7 +1435,7 @@ void HypreABec::hbvec3 (const Box& bx, bct = bctype; } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { if (bho >= 1) { h2 = 0.5e0_rt * h; th2 = 3.e0_rt * h2; @@ -1447,10 +1447,10 @@ void HypreABec::hbvec3 (const Box& bx, bfv = bfv * b(i,j,k); } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfv = beta * r / h; } - else if (bct == LO_MARSHAK || bct == LO_SANCHEZ_POMRANING) { + else if (bct == AMREX_LO_MARSHAK || bct == AMREX_LO_SANCHEZ_POMRANING) { bfv = 2.e0_rt * beta * r / h; } #ifndef AMREX_USE_GPU @@ -1475,7 +1475,7 @@ void HypreABec::hbvec3 (const Box& bx, bct = bctype; } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { if (bho >= 1) { h2 = 0.5e0_rt * h; th2 = 3.e0_rt * h2; @@ -1487,10 +1487,10 @@ void HypreABec::hbvec3 (const Box& bx, bfv = bfv * b(i,j,k+1); } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfv = beta * r / h; } - else if (bct == LO_MARSHAK || bct == LO_SANCHEZ_POMRANING) { + else if (bct == AMREX_LO_MARSHAK || bct == AMREX_LO_SANCHEZ_POMRANING) { bfv = 2.e0_rt * beta * r / h; } #ifndef AMREX_USE_GPU @@ -1591,7 +1591,7 @@ void HypreABec::hbvec (const Box& bx, } } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { if (bho >= 1) { h2 = 0.5e0_rt * h; th2 = 3.e0_rt * h2; @@ -1601,7 +1601,7 @@ void HypreABec::hbvec (const Box& bx, bfv = (beta / h) / (0.5e0_rt * h + bcl); } } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfv = beta / h; } else { @@ -1779,7 +1779,7 @@ void HypreABec::solve(MultiFab& dest, int icomp, MultiFab& rhs, BC_Mode inhom) Real bnorm; bnorm = hypre_StructInnerProd((hypre_StructVector *) b, (hypre_StructVector *) b); - bnorm = sqrt(bnorm); + bnorm = std::sqrt(bnorm); const BoxArray& grids = acoefs->boxArray(); Real volume = 0.0; @@ -1788,7 +1788,7 @@ void HypreABec::solve(MultiFab& dest, int icomp, MultiFab& rhs, BC_Mode inhom) } Real reltol_new = (bnorm > 0.0 - ? abstol / bnorm * sqrt(volume) + ? abstol / bnorm * std::sqrt(volume) : reltol); if (reltol_new > reltol) { @@ -1898,7 +1898,7 @@ Real HypreABec::getAbsoluteResidual() Real bnorm; bnorm = hypre_StructInnerProd((hypre_StructVector *) b, (hypre_StructVector *) b); - bnorm = sqrt(bnorm); + bnorm = std::sqrt(bnorm); Real res; if (solver_flag == 0) { @@ -1925,5 +1925,5 @@ Real HypreABec::getAbsoluteResidual() volume += grids[i].numPts(); } - return bnorm * res / sqrt(volume); + return bnorm * res / std::sqrt(volume); } diff --git a/Source/radiation/HypreExtMultiABec.cpp b/Source/radiation/HypreExtMultiABec.cpp index a65292751a..0930849b3e 100644 --- a/Source/radiation/HypreExtMultiABec.cpp +++ b/Source/radiation/HypreExtMultiABec.cpp @@ -423,7 +423,7 @@ void HypreExtMultiABec::loadMatrix() if (reg[ori] == domain[ori] && bd[level]->mixedBndry(ori)) { bct = (*(bd[level]->bndryTypes(ori)[i]))(v+vin); } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { if (bho == 1) { evalue(ori,i)(v+vin).push(level, v, (bcl * th2) / (h * (bcl + h2))); @@ -478,10 +478,10 @@ void HypreExtMultiABec::loadMatrix() entry(ori,i)(v).push(&(*ederiv[level])(ori,i)(v+vin), fac); } } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { // no more action required here } - else if (bct == LO_MARSHAK || bct == LO_SANCHEZ_POMRANING) { + else if (bct == AMREX_LO_MARSHAK || bct == AMREX_LO_SANCHEZ_POMRANING) { if (bho == 1) { evalue(ori,i)(v+vin).push(level, v, 1.5); evalue(ori,i)(v+vin).push(level, v-vin, -0.5); @@ -937,7 +937,7 @@ void HypreExtMultiABec::loadLevelVectorB(int level, if (reg[ori] == domain[ori] && bd[level]->mixedBndry(ori)) { bct = (*(bd[level]->bndryTypes(ori)[i]))(v+vin); } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { Real dfac, vfac; if (bho == 1) { dfac = 1.0 / ((bcl + h2) * (bcl + th2)); @@ -991,7 +991,7 @@ void HypreExtMultiABec::loadLevelVectorB(int level, HYPRE_SStructVectorAddToValues(b, part, getV1(v), 0, &tmp); } } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { // cmult should be c for photons, 1 for neutrinos Real cmult = 1.0; Real xi = 0.0; // xi should be passed in through RadBndry? @@ -1000,7 +1000,7 @@ void HypreExtMultiABec::loadLevelVectorB(int level, cmult * fs(v+vin,bdcomp)); HYPRE_SStructVectorAddToValues(b, part, getV1(v), 0, &tmp); } - else if (bct == LO_MARSHAK || bct == LO_SANCHEZ_POMRANING) { + else if (bct == AMREX_LO_MARSHAK || bct == AMREX_LO_SANCHEZ_POMRANING) { // cmult should be c for photons, 1 for neutrinos Real cmult = 1.0; Real xi = 0.0; // xi should be passed in through RadBndry? diff --git a/Source/radiation/HypreMultiABec.cpp b/Source/radiation/HypreMultiABec.cpp index 311c78c0a4..aa27b902b1 100644 --- a/Source/radiation/HypreMultiABec.cpp +++ b/Source/radiation/HypreMultiABec.cpp @@ -785,9 +785,9 @@ void HypreMultiABec::addLevel(int level, static void TransverseInterpolant(AuxVarBox& cintrp, const Mask& msk, const Box& reg, const Box& creg, - D_DECL(const IntVect& rat, const IntVect& vj1, const IntVect& vk1), - D_DECL(const IntVect& ve, const IntVect& vjr, const IntVect& vkr), - D_DECL(int idir, int jdir, int kdir), + AMREX_D_DECL(const IntVect& rat, const IntVect& vj1, const IntVect& vk1), + AMREX_D_DECL(const IntVect& ve, const IntVect& vjr, const IntVect& vkr), + AMREX_D_DECL(int idir, int jdir, int kdir), int clevel) { for (IntVect vc = creg.smallEnd(); vc <= creg.bigEnd(); creg.next(vc)) { @@ -1070,9 +1070,9 @@ void HypreMultiABec::buildMatrixStructure() const Mask &msk = bd[level]->bndryMasks(ori,i); TransverseInterpolant((*cintrp[level])(ori,i), msk, reg, creg, - D_DECL(rat, vj1, vk1), - D_DECL(ve, vjr, vkr), - D_DECL(idir, jdir, kdir), + AMREX_D_DECL(rat, vj1, vk1), + AMREX_D_DECL(ve, vjr, vkr), + AMREX_D_DECL(idir, jdir, kdir), level-1); NormalDerivative((*ederiv[level])(ori,i), @@ -1192,9 +1192,9 @@ void HypreMultiABec::buildMatrixStructure() const Mask &msk = c_cintrp[level]->mask(ori,i,j); // fine mask TransverseInterpolant((*c_cintrp[level])(ori,i,j), msk, reg, creg, - D_DECL(rat, vj1, vk1), - D_DECL(ve, vjr, vkr), - D_DECL(idir, jdir, kdir), + AMREX_D_DECL(rat, vj1, vk1), + AMREX_D_DECL(ve, vjr, vkr), + AMREX_D_DECL(idir, jdir, kdir), level-1); NormalDerivative((*c_ederiv[level])(ori,i,j), @@ -1502,7 +1502,7 @@ HypreMultiABec::hmmat (const Box& bx, Real bfm, bfv; Real bfm2, h2, th2; - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { if (bho >= 1) { @@ -1520,7 +1520,7 @@ HypreMultiABec::hmmat (const Box& bx, } } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfm = -fac; bfm2 = 0.e0_rt; @@ -1739,7 +1739,7 @@ HypreMultiABec::hmmat3 (const Box& bx, bct = bctype; } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { if (bho >= 1) { h2 = 0.5e0_rt * h; @@ -1753,13 +1753,13 @@ HypreMultiABec::hmmat3 (const Box& bx, } } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfm = 0.e0_rt; bfm2 = 0.e0_rt; } - else if (bct == LO_MARSHAK) { + else if (bct == AMREX_LO_MARSHAK) { bfv = 2.e0_rt * c * beta * r / h; @@ -1772,7 +1772,7 @@ HypreMultiABec::hmmat3 (const Box& bx, } } - else if (bct == LO_SANCHEZ_POMRANING) { + else if (bct == AMREX_LO_SANCHEZ_POMRANING) { bfv = 2.e0_rt * c * beta * r / h; @@ -1813,7 +1813,7 @@ HypreMultiABec::hmmat3 (const Box& bx, bct = bctype; } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { if (bho >= 1) { h2 = 0.5e0_rt * h; @@ -1827,13 +1827,13 @@ HypreMultiABec::hmmat3 (const Box& bx, } } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfm = 0.e0_rt; bfm2 = 0.e0_rt; } - else if (bct == LO_MARSHAK) { + else if (bct == AMREX_LO_MARSHAK) { bfv = 2.e0_rt * c * beta * r / h; @@ -1846,7 +1846,7 @@ HypreMultiABec::hmmat3 (const Box& bx, } } - else if (bct == LO_SANCHEZ_POMRANING) { + else if (bct == AMREX_LO_SANCHEZ_POMRANING) { bfv = 2.e0_rt * c * beta * r / h; @@ -1887,7 +1887,7 @@ HypreMultiABec::hmmat3 (const Box& bx, bct = bctype; } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { if (bho >= 1) { h2 = 0.5e0_rt * h; @@ -1901,13 +1901,13 @@ HypreMultiABec::hmmat3 (const Box& bx, } } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfm = 0.e0_rt; bfm2 = 0.e0_rt; } - else if (bct == LO_MARSHAK) { + else if (bct == AMREX_LO_MARSHAK) { bfv = 2.e0_rt * c * beta * r / h; @@ -1920,7 +1920,7 @@ HypreMultiABec::hmmat3 (const Box& bx, } } - else if (bct == LO_SANCHEZ_POMRANING) { + else if (bct == AMREX_LO_SANCHEZ_POMRANING) { bfv = 2.e0_rt * c * beta * r / h; @@ -1961,7 +1961,7 @@ HypreMultiABec::hmmat3 (const Box& bx, bct = bctype; } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { if (bho >= 1) { h2 = 0.5e0_rt * h; @@ -1975,13 +1975,13 @@ HypreMultiABec::hmmat3 (const Box& bx, } } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfm = 0.e0_rt; bfm2 = 0.e0_rt; } - else if (bct == LO_MARSHAK) { + else if (bct == AMREX_LO_MARSHAK) { bfv = 2.e0_rt * c * beta * r / h; @@ -1994,7 +1994,7 @@ HypreMultiABec::hmmat3 (const Box& bx, } } - else if (bct == LO_SANCHEZ_POMRANING) { + else if (bct == AMREX_LO_SANCHEZ_POMRANING) { bfv = 2.e0_rt * c * beta * r / h; @@ -2035,7 +2035,7 @@ HypreMultiABec::hmmat3 (const Box& bx, bct = bctype; } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { if (bho >= 1) { h2 = 0.5e0_rt * h; @@ -2049,13 +2049,13 @@ HypreMultiABec::hmmat3 (const Box& bx, } } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfm = 0.e0_rt; bfm2 = 0.e0_rt; } - else if (bct == LO_MARSHAK) { + else if (bct == AMREX_LO_MARSHAK) { bfv = 2.e0_rt * c * beta * r / h; @@ -2068,7 +2068,7 @@ HypreMultiABec::hmmat3 (const Box& bx, } } - else if (bct == LO_SANCHEZ_POMRANING) { + else if (bct == AMREX_LO_SANCHEZ_POMRANING) { bfv = 2.e0_rt * c * beta * r / h; @@ -2109,7 +2109,7 @@ HypreMultiABec::hmmat3 (const Box& bx, bct = bctype; } - if (bct == LO_DIRICHLET) { + if (bct == AMREX_LO_DIRICHLET) { if (bho >= 1) { h2 = 0.5e0_rt * h; @@ -2123,13 +2123,13 @@ HypreMultiABec::hmmat3 (const Box& bx, } } - else if (bct == LO_NEUMANN) { + else if (bct == AMREX_LO_NEUMANN) { bfm = 0.e0_rt; bfm2 = 0.e0_rt; } - else if (bct == LO_MARSHAK) { + else if (bct == AMREX_LO_MARSHAK) { bfv = 2.e0_rt * c * beta * r / h; @@ -2142,7 +2142,7 @@ HypreMultiABec::hmmat3 (const Box& bx, } } - else if (bct == LO_SANCHEZ_POMRANING) { + else if (bct == AMREX_LO_SANCHEZ_POMRANING) { bfv = 2.e0_rt * c * beta * r / h; @@ -2286,7 +2286,7 @@ void HypreMultiABec::loadMatrix() // level in the current linear system. Zero out the interior // stencil using Neumann BC: - const RadBoundCond bct_coarse = LO_NEUMANN; + const RadBoundCond bct_coarse = AMREX_LO_NEUMANN; hmmat(reg, matfab.array(), cdir, bct_coarse, bho, bcl, msk.array(), (*bcoefs[level])[idim][mfi].array(), beta, geom[level].CellSize()); @@ -3476,7 +3476,7 @@ void HypreMultiABec::solve() hypre_SStructInnerProd((hypre_SStructVector *) b, (hypre_SStructVector *) b, &bnorm); - bnorm = sqrt(bnorm); + bnorm = std::sqrt(bnorm); Real volume = 0.0; for (int level = crse_level; level <= fine_level; level++) { @@ -3486,7 +3486,7 @@ void HypreMultiABec::solve() } Real reltol_new = (bnorm > 0.0 - ? abstol / bnorm * sqrt(volume) + ? abstol / bnorm * std::sqrt(volume) : reltol); if (reltol_new > reltol) { @@ -3830,7 +3830,7 @@ Real HypreMultiABec::getAbsoluteResidual() hypre_SStructInnerProd((hypre_SStructVector *) b, (hypre_SStructVector *) b, &bnorm); - bnorm = sqrt(bnorm); + bnorm = std::sqrt(bnorm); Real res; if (solver_flag == 100) { @@ -3891,7 +3891,7 @@ Real HypreMultiABec::getAbsoluteResidual() } } - return bnorm * res / sqrt(volume); + return bnorm * res / std::sqrt(volume); } void HypreMultiABec::boundaryFlux(int level, diff --git a/Source/radiation/MGFLDRadSolver.cpp b/Source/radiation/MGFLDRadSolver.cpp index 250de9b1ae..3bb6187791 100644 --- a/Source/radiation/MGFLDRadSolver.cpp +++ b/Source/radiation/MGFLDRadSolver.cpp @@ -147,8 +147,8 @@ void Radiation::MGFLD_implicit_update(int level, int iteration, int ncycle) for (int idim=0; idim(value_nu[igroup], igroup); } @@ -240,7 +240,7 @@ void MGRadBndry::setBndryFluxConds(const BCRec& bc, const BC_Mode phys_bc_mode) BaseFab& tfab = *(bctypearray[face][i]); FORT_RADBNDRY2(BL_TO_FORTRAN(bnd_fab), - tfab.dataPtr(), ARLIM(domain.loVect()), ARLIM(domain.hiVect()), dx, xlo, time); + tfab.dataPtr(), AMREX_ARLIM(domain.loVect()), AMREX_ARLIM(domain.hiVect()), dx, xlo, time); #endif if (p_bcflag == 0) { // Homogeneous case. We called RADBNDRY2 only to set tfab right. diff --git a/Source/radiation/RadBndry.cpp b/Source/radiation/RadBndry.cpp index aa742febef..24079fd85f 100644 --- a/Source/radiation/RadBndry.cpp +++ b/Source/radiation/RadBndry.cpp @@ -46,7 +46,7 @@ RadBndry::RadBndry(const BoxArray& _grids, const DistributionMapping& _dmap, #if 0 FORT_RADBNDRY2(BL_TO_FORTRAN(bndry[face][bi]), bctypearray[face][igrid]->dataPtr(), - ARLIM(domain.loVect()), ARLIM(domain.hiVect()), dx, xlo, time); + AMREX_ARLIM(domain.loVect()), AMREX_ARLIM(domain.hiVect()), dx, xlo, time); #endif } } @@ -143,25 +143,25 @@ void RadBndry::setBndryConds(const BCRec& bc, if (domain[face] == grd[face] && !geom.isPeriodic(dir)) { /* // All physical bc values are located on face - if (p_bc == EXT_DIR) { - bctag[i] = LO_DIRICHLET; + if (p_bc == amrex::BCType::ext_dir) { + bctag[i] = AMREX_LO_DIRICHLET; bloc[i] = 0.; } - else if (p_bc == EXTRAP || p_bc == HOEXTRAP || p_bc == REFLECT_EVEN) { - bctag[i] = LO_NEUMANN; + else if (p_bc == EXTRAP || p_bc == amrex::BCType::hoextrap || p_bc == amrex::BCType::reflect_even) { + bctag[i] = AMREX_LO_NEUMANN; bloc[i] = 0.; } - else if (p_bc == REFLECT_ODD) { - bctag[i] = LO_REFLECT_ODD; + else if (p_bc == amrex::BCType::reflect_odd) { + bctag[i] = AMREX_LO_REFLECT_ODD; bloc[i] = 0.; } */ - if (p_bc == LO_DIRICHLET || p_bc == LO_NEUMANN || - p_bc == LO_REFLECT_ODD) { + if (p_bc == AMREX_LO_DIRICHLET || p_bc == AMREX_LO_NEUMANN || + p_bc == AMREX_LO_REFLECT_ODD) { bctag[i] = p_bc; bloc[i] = 0.; } - else if (p_bc == LO_MARSHAK || p_bc == LO_SANCHEZ_POMRANING) { + else if (p_bc == AMREX_LO_MARSHAK || p_bc == AMREX_LO_SANCHEZ_POMRANING) { bctag[i] = p_bc; //gives asymmetric, second-order version of Marshak b.c. // (worked for bbmg, works with nonsymmetric hypre solvers): @@ -177,7 +177,7 @@ void RadBndry::setBndryConds(const BCRec& bc, } else { // internal bndry - bctag[i] = LO_DIRICHLET; + bctag[i] = AMREX_LO_DIRICHLET; bloc[i] = 0.5*delta; } } @@ -211,8 +211,8 @@ void RadBndry::setBndryFluxConds(const BCRec& bc, const BC_Mode phys_bc_mode) if (domain[face] == grd[face] && !geom.isPeriodic(dir)) { if (bcflag[face] <= 1) { - if (p_bc == LO_MARSHAK || p_bc == LO_SANCHEZ_POMRANING || - p_bc == LO_DIRICHLET || p_bc == LO_NEUMANN) { + if (p_bc == AMREX_LO_MARSHAK || p_bc == AMREX_LO_SANCHEZ_POMRANING || + p_bc == AMREX_LO_DIRICHLET || p_bc == AMREX_LO_NEUMANN) { setValue(face, i, value); } } @@ -227,7 +227,7 @@ void RadBndry::setBndryFluxConds(const BCRec& bc, const BC_Mode phys_bc_mode) BaseFab& tfab = bctypearray[face][i]; FORT_RADBNDRY2(BL_TO_FORTRAN(bnd_fab), tfab.dataPtr(), - ARLIM(domain.loVect()), ARLIM(domain.hiVect()), dx, xlo, time); + AMREX_ARLIM(domain.loVect()), AMREX_ARLIM(domain.hiVect()), dx, xlo, time); if (p_bcflag == 0) { // Homogeneous case. We called RADBNDRY2 only to set tfab right. setValue(face, i, value); diff --git a/Source/radiation/RadMultiGroup.cpp b/Source/radiation/RadMultiGroup.cpp index 83a1811ffc..4a8e682d78 100644 --- a/Source/radiation/RadMultiGroup.cpp +++ b/Source/radiation/RadMultiGroup.cpp @@ -43,37 +43,37 @@ void Radiation::get_groups(int verbose) if (lowest == 0.0) { nugroup[0] = 0.5*dnugroup[0]; - dlognugroup[0] = 2.0 * (log(xnu[1]) - log(nugroup[0])); + dlognugroup[0] = 2.0 * (std::log(xnu[1]) - std::log(nugroup[0])); } else { - nugroup[0] = sqrt(xnu[0]*xnu[1]); - dlognugroup[0] = log(xnu[1]) - log(xnu[0]); + nugroup[0] = std::sqrt(xnu[0]*xnu[1]); + dlognugroup[0] = std::log(xnu[1]) - std::log(xnu[0]); } for (int i=1; i& lambda, int bool nexttoboundary=false; for (int idim=0; idimGeom(level).CellSize(); - const Real volume = D_TERM(dx[0], * dx[1], * dx[2]); + const Real volume = AMREX_D_TERM(dx[0], * dx[1], * dx[2]); if (flux_in) { for (int n = 0; n < AMREX_SPACEDIM; n++) { diff --git a/Source/radiation/SGRadSolver.cpp b/Source/radiation/SGRadSolver.cpp index eca2c5db37..0cf0b76904 100644 --- a/Source/radiation/SGRadSolver.cpp +++ b/Source/radiation/SGRadSolver.cpp @@ -129,8 +129,8 @@ void Radiation::single_group_update(int level, int iteration, int ncycle) for (int idim=0; idim tmplen) { + if (std::pow(mxlen,(float)AMREX_SPACEDIM-1) > tmplen) { delete [] derives; #if (AMREX_SPACEDIM == 1) derives = new Real[1]; diff --git a/Source/reactions/Castro_react.cpp b/Source/reactions/Castro_react.cpp index e050e2008f..6859b891a5 100644 --- a/Source/reactions/Castro_react.cpp +++ b/Source/reactions/Castro_react.cpp @@ -186,12 +186,14 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra MultiFab tmp_mask_mf; const MultiFab& mask_mf = mask_covered_zones ? getLevel(level+1).build_fine_mask() : tmp_mask_mf; - ReduceOps reduce_op; - ReduceData reduce_data(reduce_op); - using ReduceTuple = typename decltype(reduce_data)::Type; +#if defined(AMREX_USE_GPU) + Gpu::Buffer d_num_failed({0}); + auto* p_num_failed = d_num_failed.data(); +#endif + int num_failed = 0; #ifdef _OPENMP -#pragma omp parallel +#pragma omp parallel reduction(+:num_failed) #endif for (MFIter mfi(s, TilingIfNotGPU()); mfi.isValid(); ++mfi) { @@ -201,15 +203,18 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra auto U = s.array(mfi); auto reactions = r.array(mfi); auto weights = store_burn_weights ? burn_weights.array(mfi) : Array4{}; - auto mask = mask_covered_zones ? mask_mf.array(mfi) : Array4{}; + const auto mask = mask_covered_zones ? mask_mf.array(mfi) : Array4{}; const auto dx = geom.CellSizeArray(); #ifdef MODEL_PARSER const auto problo = geom.ProbLoArray(); #endif - reduce_op.eval(bx, reduce_data, - [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) -> ReduceTuple +#if defined(AMREX_USE_GPU) + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) +#else + LoopOnCpu(bx, [&] (int i, int j, int k) mutable +#endif { burn_t burn_state; @@ -223,14 +228,14 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra #if AMREX_SPACEDIM == 1 burn_state.dx = dx[0]; #else - burn_state.dx = amrex::min(D_DECL(dx[0], dx[1], dx[2])); + burn_state.dx = amrex::min(AMREX_D_DECL(dx[0], dx[1], dx[2])); #endif // Initialize some data for later. bool do_burn = true; burn_state.success = true; - Real burn_failed = 0.0_rt; + int burn_failed = 0; // Don't burn on zones inside shock regions, if the relevant option is set. @@ -329,7 +334,7 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra // If we were unsuccessful, update the failure count. if (!burn_state.success) { - burn_failed = 1.0_rt; + burn_failed = 1; } // Add burning rates to reactions MultiFab, but be @@ -399,19 +404,25 @@ Castro::react_state(MultiFab& s, MultiFab& r, Real time, Real dt, const int stra } - - return {burn_failed}; - +#if defined(AMREX_USE_GPU) + if (burn_failed) { + Gpu::Atomic::Add(p_num_failed, burn_failed); + } +#else + num_failed += burn_failed; +#endif }); +#if defined(AMREX_USE_HIP) + Gpu::streamSynchronize(); // otherwise HIP may faile to allocate the necessary resources. +#endif } - ReduceTuple hv = reduce_data.value(); - Real burn_failed = amrex::get<0>(hv); +#if defined(AMREX_USE_GPU) + num_failed = *(d_num_failed.copyToHost()); +#endif - if (burn_failed != 0.0) { - burn_success = 0; - } + burn_success = !num_failed; ParallelDescriptor::ReduceIntMin(burn_success); @@ -516,11 +527,13 @@ Castro::react_state(Real time, Real dt) int burn_success = 1; - ReduceOps reduce_op; - ReduceData reduce_data(reduce_op); - - using ReduceTuple = typename decltype(reduce_data)::Type; +#if defined(AMREX_USE_GPU) + Gpu::Buffer d_num_failed({0}); + auto* p_num_failed = d_num_failed.data(); +#endif + int num_failed = 0; + // why no omp here? for (MFIter mfi(S_new, TilingIfNotGPU()); mfi.isValid(); ++mfi) { const Box& bx = mfi.growntilebox(ng); @@ -535,22 +548,25 @@ Castro::react_state(Real time, Real dt) auto I = SDC_react.array(mfi); auto react_src = reactions.array(mfi); auto weights = store_burn_weights ? burn_weights.array(mfi) : Array4{}; - auto mask = mask_covered_zones ? mask_mf.array(mfi) : Array4{}; + const auto mask = mask_covered_zones ? mask_mf.array(mfi) : Array4{}; int lsdc_iteration = sdc_iteration; const auto dx = geom.CellSizeArray(); const auto problo = geom.ProbLoArray(); - reduce_op.eval(bx, reduce_data, - [=] AMREX_GPU_HOST_DEVICE (int i, int j, int k) -> ReduceTuple +#if defined(AMREX_USE_GPU) + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) +#else + LoopOnCpu(bx, [&] (int i, int j, int k) mutable +#endif { burn_t burn_state; #if AMREX_SPACEDIM == 1 burn_state.dx = dx[0]; #else - burn_state.dx = amrex::min(D_DECL(dx[0], dx[1], dx[2])); + burn_state.dx = amrex::min(AMREX_D_DECL(dx[0], dx[1], dx[2])); #endif #ifdef NSE_NET @@ -563,7 +579,7 @@ Castro::react_state(Real time, Real dt) bool do_burn = true; burn_state.success = true; - Real burn_failed = 0.0_rt; + int burn_failed = 0; // Don't burn on zones inside shock regions, if the // relevant option is set. @@ -687,7 +703,7 @@ Castro::react_state(Real time, Real dt) // If we were unsuccessful, update the failure count. if (!burn_state.success) { - burn_failed = 1.0_rt; + burn_failed = 1; } // update the state data. @@ -780,16 +796,25 @@ Castro::react_state(Real time, Real dt) } } - return {burn_failed}; +#if defined(AMREX_USE_GPU) + if (burn_failed) { + Gpu::Atomic::Add(p_num_failed, burn_failed); + } +#else + num_failed += burn_failed; +#endif }); + +#if defined(AMREX_USE_HIP) + Gpu::streamSynchronize(); // otherwise HIP may faile to allocate the necessary resources. +#endif } - ReduceTuple hv = reduce_data.value(); - Real burn_failed = amrex::get<0>(hv); +#if defined(AMREX_USE_GPU) + num_failed = *(d_num_failed.copyToHost()); +#endif - if (burn_failed != 0.0) { - burn_success = 0; - } + burn_success = !num_failed; ParallelDescriptor::ReduceIntMin(burn_success); diff --git a/Source/scf/scf_relax.cpp b/Source/scf/scf_relax.cpp index 1add2074e3..ec55704406 100644 --- a/Source/scf/scf_relax.cpp +++ b/Source/scf/scf_relax.cpp @@ -338,7 +338,7 @@ Castro::do_hscf_solve() amrex::Error(); } - Real omega = sqrt(omegasq); + Real omega = std::sqrt(omegasq); // Rotational period is 2 pi / omega. // Let's also be sure not to let the period diff --git a/Util/ConvertCheckpoint/Embiggen.cpp b/Util/ConvertCheckpoint/Embiggen.cpp index f8c2d43b97..32893be58c 100644 --- a/Util/ConvertCheckpoint/Embiggen.cpp +++ b/Util/ConvertCheckpoint/Embiggen.cpp @@ -726,7 +726,7 @@ static void ConvertData() { { int shiftx(jx*dlenx); int shifty(jy*dleny); - IntVect shift_box(D_DECL(shiftx,shifty,shiftz)); + IntVect shift_box(AMREX_D_DECL(shiftx,shifty,shiftz)); Box bx(falRef0.grids[n]); bx.shift(shift_box); @@ -762,7 +762,7 @@ static void ConvertData() { int shiftz(jz*dlenz - dlenz/2); #endif #endif - IntVect shift_box(D_DECL(shiftx,shifty,shiftz)); + IntVect shift_box(AMREX_D_DECL(shiftx,shifty,shiftz)); Box bx(falRef0.grids[n]); bx.shift(shift_box); @@ -792,7 +792,7 @@ static void ConvertData() { int shiftz(jz*dlenz); #endif #endif - IntVect shift_box(D_DECL(shiftx,shifty,shiftz)); + IntVect shift_box(AMREX_D_DECL(shiftx,shifty,shiftz)); Box bx(falRef0.grids[n]); bx.shift(shift_box); diff --git a/Util/exact_riemann/exact_riemann.H b/Util/exact_riemann/exact_riemann.H index 021120a322..d7e7c4831b 100644 --- a/Util/exact_riemann/exact_riemann.H +++ b/Util/exact_riemann/exact_riemann.H @@ -1,12 +1,16 @@ #ifndef EXACT_RIEMANN_H #define EXACT_RIEMANN_H +#include + #include #include #include #include +using namespace amrex::literals; + AMREX_INLINE void exact_riemann() { @@ -15,7 +19,7 @@ exact_riemann() { // exploring composition jumps here. We'll take a constant // composition. - Real xn[NumSpec] = {0.0}; + amrex::Real xn[NumSpec] = {0.0}; xn[0] = 1.0_rt; @@ -44,12 +48,12 @@ exact_riemann() { } if (problem::co_moving_frame) { - Real W_avg = 0.5_rt * (problem::u_l + problem::u_r); + amrex::Real W_avg = 0.5_rt * (problem::u_l + problem::u_r); problem::u_l -= W_avg; problem::u_r -= W_avg; } - Real ustar, pstar, W_l, W_r; + amrex::Real ustar, pstar, W_l, W_r; riemann_star_state(problem::rho_l, problem::u_l, problem::p_l, xn, problem::rho_r, problem::u_r, problem::p_r, xn, @@ -62,7 +66,7 @@ exact_riemann() { // This follows from the discussion around C&G Eq. 15 - Real dx = (problem::xmax - problem::xmin) / static_cast(problem::npts); + amrex::Real dx = (problem::xmax - problem::xmin) / static_cast(problem::npts); // loop over xi space @@ -72,9 +76,9 @@ exact_riemann() { // compute xi = x/t -- this is the similarity variable for the // solution - Real x = problem::xmin + (static_cast(i) - 0.5_rt) * dx; + amrex::Real x = problem::xmin + (static_cast(i) - 0.5_rt) * dx; - Real rho, u, p, xn_s[NumSpec]; + amrex::Real rho, u, p, xn_s[NumSpec]; riemann_sample(problem::rho_l, problem::u_l, problem::p_l, xn, problem::rho_r, problem::u_r, problem::p_r, xn, @@ -83,7 +87,7 @@ exact_riemann() { rho, u, p, xn_s); if (problem::co_moving_frame) { - Real W_avg = 0.5_rt * (problem::u_l + problem::u_r); + amrex::Real W_avg = 0.5_rt * (problem::u_l + problem::u_r); u += W_avg; x += problem::t * W_avg; } diff --git a/Util/exact_riemann/main.cpp b/Util/exact_riemann/main.cpp index 20d8447738..e91b337308 100644 --- a/Util/exact_riemann/main.cpp +++ b/Util/exact_riemann/main.cpp @@ -17,25 +17,25 @@ int main(int argc, char *argv[]) { amrex::Initialize(argc, argv); - std::cout << "starting the exact Riemann solver..." << std::endl; - std::cout << argv[1] << std::endl; - std::cout << strlen(argv[1]) << std::endl; + std::cout << "starting the exact Riemann solver..." << std::endl; + std::cout << argv[1] << std::endl; + std::cout << strlen(argv[1]) << std::endl; - // initialize the external runtime parameters in C++ + // initialize the external runtime parameters in C++ - init_prob_parameters(); + init_prob_parameters(); - init_extern_parameters(); + init_extern_parameters(); - // now initialize the C++ Microphysics + // now initialize the C++ Microphysics #ifdef REACTIONS - network_init(); + network_init(); #endif - Real small_temp = 1.e-200; - Real small_dens = 1.e-200; - eos_init(small_temp, small_dens); + amrex::Real small_temp = 1.e-200; + amrex::Real small_dens = 1.e-200; + eos_init(small_temp, small_dens); - exact_riemann(); + exact_riemann(); } diff --git a/Util/exact_riemann/riemann_sample.H b/Util/exact_riemann/riemann_sample.H index dcadceb8a9..c2580dd258 100644 --- a/Util/exact_riemann/riemann_sample.H +++ b/Util/exact_riemann/riemann_sample.H @@ -1,16 +1,20 @@ #ifndef RIEMANN_SAMPLE_MODULE_H #define RIEMANN_SAMPLE_MODULE_H +#include + #include +using namespace amrex::literals; + AMREX_INLINE void -riemann_sample(const Real rho_l, const Real u_l, const Real p_l, const Real* xn_l, - const Real rho_r, const Real u_r, const Real p_r, const Real* xn_r, - const Real ustar, const Real pstar, - const Real W_l, const Real W_r, - const Real x, const Real xjump, const Real time, - Real& rho, Real& u, Real& p, Real* xn) { +riemann_sample(const amrex::Real rho_l, const amrex::Real u_l, const amrex::Real p_l, const amrex::Real* xn_l, + const amrex::Real rho_r, const amrex::Real u_r, const amrex::Real p_r, const amrex::Real* xn_r, + const amrex::Real ustar, const amrex::Real pstar, + const amrex::Real W_l, const amrex::Real W_r, + const amrex::Real x, const amrex::Real xjump, const amrex::Real time, + amrex::Real& rho, amrex::Real& u, amrex::Real& p, amrex::Real* xn) { // get the initial sound speeds @@ -24,7 +28,7 @@ riemann_sample(const Real rho_l, const Real u_l, const Real p_l, const Real* xn_ eos(eos_input_rp, eos_state); - Real cs_l = std::sqrt(eos_state.gam1 * p_l / rho_l); + amrex::Real cs_l = std::sqrt(eos_state.gam1 * p_l / rho_l); eos_state.rho = rho_r; eos_state.p = p_r; @@ -35,7 +39,7 @@ riemann_sample(const Real rho_l, const Real u_l, const Real p_l, const Real* xn_ eos(eos_input_rp, eos_state); - Real cs_r = std::sqrt(eos_state.gam1 * p_r / rho_r); + amrex::Real cs_r = std::sqrt(eos_state.gam1 * p_r / rho_r); // find the solution as a function of xi = x/t @@ -45,14 +49,14 @@ riemann_sample(const Real rho_l, const Real u_l, const Real p_l, const Real* xn_ // compute xi = x/t -- this is the similarity variable for the // solution - Real xi = (x - xjump) / time; + amrex::Real xi = (x - xjump) / time; // check which side of the contact we need to worry about - Real chi = std::copysign(1.0_rt, xi - ustar); + amrex::Real chi = std::copysign(1.0_rt, xi - ustar); - Real rho_s, u_s, p_s, W_s, cs_s; - Real uhat_s, xihat, uhat_star; + amrex::Real rho_s, u_s, p_s, W_s, cs_s; + amrex::Real uhat_s, xihat, uhat_star; if (chi == -1.0_rt) { rho_s = rho_l; @@ -97,7 +101,7 @@ riemann_sample(const Real rho_l, const Real u_l, const Real p_l, const Real* xn_ // are we a shock or rarefaction? - Real rhostar, Z_temp; + amrex::Real rhostar, Z_temp; if (pstar > p_s) { // shock @@ -107,7 +111,7 @@ riemann_sample(const Real rho_l, const Real u_l, const Real p_l, const Real* xn_ } else { // rarefaction. Here we need to integrate the Riemann // invariant curves to our pstar to get rhostar and cs_star - Real Z_temp, W_temp; + amrex::Real Z_temp, W_temp; if (chi == -1.0_rt) { rarefaction(pstar, rho_s, u_s, p_s, xn, 1, Z_temp, W_temp, rhostar); } else { @@ -128,11 +132,11 @@ riemann_sample(const Real rho_l, const Real u_l, const Real p_l, const Real* xn_ eos(eos_input_rp, eos_state); - Real cs_star = std::sqrt(eos_state.gam1 * pstar / rhostar); + amrex::Real cs_star = std::sqrt(eos_state.gam1 * pstar / rhostar); // now deal with the cases where we are not spanning a rarefaction - Real lambdahat_s, lambdahat_star; + amrex::Real lambdahat_s, lambdahat_star; if (pstar <= p_s) { lambdahat_s = uhat_s + cs_s; diff --git a/Util/exact_riemann/riemann_star_state.H b/Util/exact_riemann/riemann_star_state.H index 9cc23a8158..b29c81dd50 100644 --- a/Util/exact_riemann/riemann_star_state.H +++ b/Util/exact_riemann/riemann_star_state.H @@ -1,17 +1,21 @@ #ifndef RIEMANN_STAR_STATE_H #define RIEMANN_STAR_STATE_H +#include + #include +using namespace amrex::literals; + AMREX_INLINE void -riemann_star_state(const Real rho_l, const Real u_l, const Real p_l, const Real* xn_l, - const Real rho_r, const Real u_r, const Real p_r, const Real* xn_r, - Real& ustar, Real& pstar, Real& W_l, Real& W_r) { +riemann_star_state(const amrex::Real rho_l, const amrex::Real u_l, const amrex::Real p_l, const amrex::Real* xn_l, + const amrex::Real rho_r, const amrex::Real u_r, const amrex::Real p_r, const amrex::Real* xn_r, + amrex::Real& ustar, amrex::Real& pstar, amrex::Real& W_l, amrex::Real& W_r) { - const Real tol = 1.e-10_rt; - const Real smallp = 1.e-8_rt; - const Real SMALL = 1.e-13_rt; + const amrex::Real tol = 1.e-10_rt; + const amrex::Real smallp = 1.e-8_rt; + const amrex::Real SMALL = 1.e-13_rt; // get the initial sound speeds @@ -26,10 +30,10 @@ riemann_star_state(const Real rho_l, const Real u_l, const Real p_l, const Real* eos(eos_input_rp, eos_state); - Real cs_l = std::sqrt(eos_state.gam1 * p_l / rho_l); + amrex::Real cs_l = std::sqrt(eos_state.gam1 * p_l / rho_l); - Real gammaE_l = p_l / (rho_l * eos_state.e) + 1.0_rt; - Real gammaC_l = eos_state.gam1; + amrex::Real gammaE_l = p_l / (rho_l * eos_state.e) + 1.0_rt; + amrex::Real gammaC_l = eos_state.gam1; eos_state.rho = rho_r; eos_state.p = p_r; @@ -40,13 +44,13 @@ riemann_star_state(const Real rho_l, const Real u_l, const Real p_l, const Real* eos(eos_input_rp, eos_state); - Real cs_r = std::sqrt(eos_state.gam1 * p_r / rho_r); + amrex::Real cs_r = std::sqrt(eos_state.gam1 * p_r / rho_r); - Real gammaE_r = p_r / (rho_r * eos_state.e) + 1.0_rt; - Real gammaC_r = eos_state.gam1; + amrex::Real gammaE_r = p_r / (rho_r * eos_state.e) + 1.0_rt; + amrex::Real gammaC_r = eos_state.gam1; - Real gammaE_bar = 0.5_rt * (gammaE_l + gammaE_r); - Real gammaC_bar = 0.5_rt * (gammaC_l + gammaC_r); + amrex::Real gammaE_bar = 0.5_rt * (gammaE_l + gammaE_r); + amrex::Real gammaC_bar = 0.5_rt * (gammaC_l + gammaC_r); // create an initial guess for pstar using a primitive variable @@ -78,7 +82,7 @@ riemann_star_state(const Real rho_l, const Real u_l, const Real p_l, const Real* // this procedure follows directly from Colella & Glaz 1985, section 1 - Real ustar_l, ustar_r; + amrex::Real ustar_l, ustar_r; bool converged = false; int iter = 1; @@ -87,7 +91,7 @@ riemann_star_state(const Real rho_l, const Real u_l, const Real p_l, const Real* // compute Z_l and Z_r -- the form of these depend on whether the // wave is a shock or a rarefaction - Real Z_l, Z_r; + amrex::Real Z_l, Z_r; // left wave @@ -96,7 +100,7 @@ riemann_star_state(const Real rho_l, const Real u_l, const Real p_l, const Real* shock(pstar, rho_l, u_l, p_l, xn_l, gammaE_bar, gammaC_bar, Z_l, W_l); } else { // left rarefaction - Real rhostar_dummy; + amrex::Real rhostar_dummy; rarefaction(pstar, rho_l, u_l, p_l, xn_l, 1, Z_l, W_l, rhostar_dummy); } @@ -107,21 +111,21 @@ riemann_star_state(const Real rho_l, const Real u_l, const Real p_l, const Real* shock(pstar, rho_r, u_r, p_r, xn_r, gammaE_bar, gammaC_bar, Z_r, W_r); } else { // right rarefaction - Real rhostar_dummy; + amrex::Real rhostar_dummy; rarefaction(pstar, rho_r, u_r, p_r, xn_r, 3, Z_r, W_r, rhostar_dummy); } ustar_l = u_l - (pstar - p_l) / W_l; ustar_r = u_r + (pstar - p_r) / W_r; - Real pstar_new = pstar - Z_l * Z_r * (ustar_r - ustar_l) / (Z_l + Z_r); + amrex::Real pstar_new = pstar - Z_l * Z_r * (ustar_r - ustar_l) / (Z_l + Z_r); std::cout << "done with iteration " << iter << std::endl; std::cout << "ustar_l/r, pstar: " << ustar_l << " " << ustar_r << " " << pstar_new << std::endl; // estimate the error in the current star solution - Real err1 = std::abs(ustar_r - ustar_l); - Real err2 = pstar_new - pstar; + amrex::Real err1 = std::abs(ustar_r - ustar_l); + amrex::Real err2 = pstar_new - pstar; if (err1 < tol * amrex::max(std::abs(ustar_l), std::abs(ustar_r)) && err2 < tol * pstar) { diff --git a/Util/exact_riemann/riemann_support.H b/Util/exact_riemann/riemann_support.H index 002e183815..80f3391492 100644 --- a/Util/exact_riemann/riemann_support.H +++ b/Util/exact_riemann/riemann_support.H @@ -1,22 +1,26 @@ #ifndef RIEMANN_SUPPORT_H #define RIEMANN_SUPPORT_H +#include + #include #include -const Real smallrho = 1.e-5_rt; +using namespace amrex::literals; + +const amrex::Real smallrho = 1.e-5_rt; const int max_iters = 100; -const Real tol = 1.e-6_rt; +const amrex::Real tol = 1.e-6_rt; AMREX_INLINE void -W_s_shock(const Real W_s, const Real pstar, const Real rho_s, const Real p_s, const Real e_s, const Real* xn, - Real& rhostar_s, eos_t& eos_state, Real& f, Real& fprime) { +W_s_shock(const amrex::Real W_s, const amrex::Real pstar, const amrex::Real rho_s, const amrex::Real p_s, const amrex::Real e_s, const amrex::Real* xn, + amrex::Real& rhostar_s, eos_t& eos_state, amrex::Real& f, amrex::Real& fprime) { // we need rhostar -- get it from the R-H conditions - Real taustar_s = (1.0_rt / rho_s) - (pstar - p_s) / (W_s * W_s); + amrex::Real taustar_s = (1.0_rt / rho_s) - (pstar - p_s) / (W_s * W_s); rhostar_s = 1.0_rt / taustar_s; // get the thermodynamics @@ -36,7 +40,7 @@ W_s_shock(const Real W_s, const Real pstar, const Real rho_s, const Real p_s, co // we need de/drho at constant p -- this is not returned by the EOS - Real dedrho_p = eos_state.dedr - eos_state.dedT * eos_state.dpdr / eos_state.dpdT; + amrex::Real dedrho_p = eos_state.dedr - eos_state.dedT * eos_state.dpdr / eos_state.dpdT; fprime = 2.0_rt * W_s * (eos_state.e - e_s) - 2.0_rt * dedrho_p * (pstar - p_s) * std::pow(rhostar_s, 2) / W_s; @@ -45,8 +49,8 @@ W_s_shock(const Real W_s, const Real pstar, const Real rho_s, const Real p_s, co AMREX_INLINE void -newton_shock(Real& W_s, const Real pstar, const Real rho_s, const Real p_s, const Real e_s, const Real* xn, - Real* rhostar_hist, Real* Ws_hist, +newton_shock(amrex::Real& W_s, const amrex::Real pstar, const amrex::Real rho_s, const amrex::Real p_s, const amrex::Real e_s, const amrex::Real* xn, + amrex::Real* rhostar_hist, amrex::Real* Ws_hist, eos_t& eos_state, bool& converged) { // Newton iterations -- we are zeroing the energy R-H jump condition @@ -61,12 +65,12 @@ newton_shock(Real& W_s, const Real pstar, const Real rho_s, const Real p_s, cons int iter = 1; while (! converged && iter < max_iters) { - Real rhostar_s, f, fprime; + amrex::Real rhostar_s, f, fprime; W_s_shock(W_s, pstar, rho_s, p_s, e_s, xn, rhostar_s, eos_state, f, fprime); - Real dW = -f / fprime; + amrex::Real dW = -f / fprime; if (std::abs(dW) < tol * W_s) { converged = true; @@ -85,14 +89,14 @@ newton_shock(Real& W_s, const Real pstar, const Real rho_s, const Real p_s, cons AMREX_INLINE void -shock(const Real pstar, const Real rho_s, const Real u_s, const Real p_s, const Real* xn, - const Real gammaE_bar, const Real gammaC_bar, - Real& Z_s, Real& W_s) { +shock(const amrex::Real pstar, const amrex::Real rho_s, const amrex::Real u_s, const amrex::Real p_s, const amrex::Real* xn, + const amrex::Real gammaE_bar, const amrex::Real gammaC_bar, + amrex::Real& Z_s, amrex::Real& W_s) { - const Real tol_p = 1.e-6_rt; + const amrex::Real tol_p = 1.e-6_rt; - Real rhostar_hist[max_iters], Ws_hist[max_iters]; + amrex::Real rhostar_hist[max_iters], Ws_hist[max_iters]; // compute the Z_s function for a shock following C&G Eq. 20 and // 23. Here the "_s" variables are the state (L or R) that we are @@ -113,15 +117,15 @@ shock(const Real pstar, const Real rho_s, const Real u_s, const Real p_s, const eos(eos_input_rp, eos_state); - Real e_s = eos_state.e; + amrex::Real e_s = eos_state.e; // to kick things off, we need a guess for W_s. We'll use the // approximation from Colella & Glaz (Eq. 34), which in turn // makes an approximation for gammaE_star, using equation 31. - Real gammaE_s = p_s / (rho_s * e_s) + 1.0_rt; + amrex::Real gammaE_s = p_s / (rho_s * e_s) + 1.0_rt; - Real gammaE_star = gammaE_s + + amrex::Real gammaE_star = gammaE_s + 2.0_rt * (1.0_rt - gammaE_bar / gammaC_bar) * (gammaE_bar - 1.0_rt) * (pstar - p_s) / (pstar + p_s); @@ -143,15 +147,15 @@ shock(const Real pstar, const Real rho_s, const Real u_s, const Real p_s, const // we need rhostar -- get it from the R-H conditions - Real taustar_s = (1.0_rt / rho_s) - (pstar - p_s) / (W_s * W_s); - Real rhostar_s; + amrex::Real taustar_s = (1.0_rt / rho_s) - (pstar - p_s) / (W_s * W_s); + amrex::Real rhostar_s; if (taustar_s < 0.0_rt) { rhostar_s = smallrho; W_s = std::sqrt((pstar - p_s) / (1.0_rt / rho_s - 1.0_rt / rhostar_s)); } - Real W_s_guess = W_s; + amrex::Real W_s_guess = W_s; // newton @@ -181,17 +185,17 @@ shock(const Real pstar, const Real rho_s, const Real u_s, const Real p_s, const // next we compute the derivative dW_s/dpstar -- the paper gives // dW**2/dpstar (Eq. 23), so we take 1/2W of that - Real C = std::sqrt(eos_state.gam1 * pstar * rhostar_s); + amrex::Real C = std::sqrt(eos_state.gam1 * pstar * rhostar_s); - Real p_e = eos_state.dpdT / eos_state.dedT; - Real p_rho = eos_state.dpdr - eos_state.dpdT * eos_state.dedr / eos_state.dedT; + amrex::Real p_e = eos_state.dpdT / eos_state.dedT; + amrex::Real p_rho = eos_state.dpdr - eos_state.dpdT * eos_state.dedr / eos_state.dedT; - Real p_tau = -std::pow(rhostar_s, 2) * p_rho; + amrex::Real p_tau = -std::pow(rhostar_s, 2) * p_rho; - Real dW2dpstar = (C*C - W_s*W_s) * W_s * W_s / + amrex::Real dW2dpstar = (C*C - W_s*W_s) * W_s * W_s / ((0.5_rt * (pstar + p_s) * p_e - p_tau) * (pstar - p_s)); - Real dWdpstar = 0.5_rt * dW2dpstar / W_s; + amrex::Real dWdpstar = 0.5_rt * dW2dpstar / W_s; // finally compute Z_s @@ -202,8 +206,8 @@ shock(const Real pstar, const Real rho_s, const Real u_s, const Real p_s, const AMREX_INLINE void -riemann_invariant_rhs(const Real p, const Real tau, const Real u, const Real* xn, const int iwave, - Real& dtaudp, Real& dudp) { +riemann_invariant_rhs(const amrex::Real p, const amrex::Real tau, const amrex::Real u, const amrex::Real* xn, const int iwave, + amrex::Real& dtaudp, amrex::Real& dudp) { // here, p is out independent variable, and tau, u are the // dependent variables. We return the derivatives of these @@ -221,7 +225,7 @@ riemann_invariant_rhs(const Real p, const Real tau, const Real u, const Real* xn eos(eos_input_rp, eos_state); - Real C = std::sqrt(eos_state.gam1 * p / tau); + amrex::Real C = std::sqrt(eos_state.gam1 * p / tau); dtaudp = -1.0_rt / (C * C); @@ -236,8 +240,8 @@ riemann_invariant_rhs(const Real p, const Real tau, const Real u, const Real* xn AMREX_INLINE void -riemann_invariant_rhs2(const Real u, const Real tau, const Real p, const Real* xn, const int iwave, - Real& dtaudu, Real& dpdu) { +riemann_invariant_rhs2(const amrex::Real u, const amrex::Real tau, const amrex::Real p, const amrex::Real* xn, const int iwave, + amrex::Real& dtaudu, amrex::Real& dpdu) { // here, u is out independent variable, and tau, p are the // dependent variables. We return the derivatives of these @@ -255,7 +259,7 @@ riemann_invariant_rhs2(const Real u, const Real tau, const Real p, const Real* x eos(eos_input_rp, eos_state); - Real C = std::sqrt(eos_state.gam1 * p / tau); + amrex::Real C = std::sqrt(eos_state.gam1 * p / tau); if (iwave == 3) { dpdu = C; @@ -270,8 +274,8 @@ riemann_invariant_rhs2(const Real u, const Real tau, const Real p, const Real* x AMREX_INLINE void -rarefaction(const Real pstar, const Real rho_s, const Real u_s, const Real p_s, - const Real* xn, const int iwave, Real& Z_s, Real& W_s, Real& rhostar) { +rarefaction(const amrex::Real pstar, const amrex::Real rho_s, const amrex::Real u_s, const amrex::Real p_s, + const amrex::Real* xn, const int iwave, amrex::Real& Z_s, amrex::Real& W_s, amrex::Real& rhostar) { const int npts = 1000; @@ -279,27 +283,27 @@ rarefaction(const Real pstar, const Real rho_s, const Real u_s, const Real p_s, // region by integrating the Riemann invariant from p_s to pstar. // This means solving a system of ODEs. We use 4th-order R-K. - Real tau = 1.0_rt / rho_s; - Real u = u_s; - Real p = p_s; + amrex::Real tau = 1.0_rt / rho_s; + amrex::Real u = u_s; + amrex::Real p = p_s; - Real dp = (pstar - p_s) / static_cast(npts); - Real dp2 = 0.5_rt * dp; + amrex::Real dp = (pstar - p_s) / static_cast(npts); + amrex::Real dp2 = 0.5_rt * dp; for (int i = 1; i <= npts; ++i) { // do 4th-order RT - Real dtaudp1, dudp1; + amrex::Real dtaudp1, dudp1; riemann_invariant_rhs(p, tau, u, xn, iwave, dtaudp1, dudp1); - Real dtaudp2, dudp2; + amrex::Real dtaudp2, dudp2; riemann_invariant_rhs(p+dp2, tau+dp2*dtaudp1, u+dp2*dudp1, xn, iwave, dtaudp2, dudp2); - Real dtaudp3, dudp3; + amrex::Real dtaudp3, dudp3; riemann_invariant_rhs(p+dp2, tau+dp2*dtaudp2, u+dp2*dudp2, xn, iwave, dtaudp3, dudp3); - Real dtaudp4, dudp4; + amrex::Real dtaudp4, dudp4; riemann_invariant_rhs(p+dp, tau+dp*dtaudp3, u+dp*dudp3, xn, iwave, dtaudp4, dudp4); p += dp; @@ -336,8 +340,8 @@ rarefaction(const Real pstar, const Real rho_s, const Real u_s, const Real p_s, AMREX_INLINE void -rarefaction_to_u(const Real rho_s, const Real u_s, const Real p_s, const Real* xn, const int iwave, const Real xi, - Real& rho, Real& p, Real& u) { +rarefaction_to_u(const amrex::Real rho_s, const amrex::Real u_s, const amrex::Real p_s, const amrex::Real* xn, const int iwave, const amrex::Real xi, + amrex::Real& rho, amrex::Real& p, amrex::Real& u) { const int npts = 1000; @@ -355,7 +359,7 @@ rarefaction_to_u(const Real rho_s, const Real u_s, const Real p_s, const Real* x // dp/du = C; dtau/du = -1/C for the 1-wave // dp/du = -C; dtau/du = 1/C for the 3-wave - Real tau = 1.0_rt / rho_s; + amrex::Real tau = 1.0_rt / rho_s; u = u_s; p = p_s; @@ -372,37 +376,37 @@ rarefaction_to_u(const Real rho_s, const Real u_s, const Real p_s, const Real* x eos(eos_input_rp, eos_state); - Real c = std::sqrt(eos_state.gam1 * p * tau); + amrex::Real c = std::sqrt(eos_state.gam1 * p * tau); - Real ustop; + amrex::Real ustop; if (iwave == 1) { ustop = xi + c; } else if (iwave == 3) { ustop = xi - c; } - Real du = (ustop - u_s) / static_cast(npts); + amrex::Real du = (ustop - u_s) / static_cast(npts); bool finished = false; std::cout << "integrating from u: " << u << " " << ustop << " " << xi << " " << c << std::endl; - Real du2 = 0.5_rt * du; + amrex::Real du2 = 0.5_rt * du; while (! finished) { // do 4th-order RT - Real dtaudu1, dpdu1; + amrex::Real dtaudu1, dpdu1; riemann_invariant_rhs2(u, tau, p, xn, iwave, dtaudu1, dpdu1); - Real dtaudu2, dpdu2; + amrex::Real dtaudu2, dpdu2; riemann_invariant_rhs2(u+du2, tau+du2*dtaudu1, p+du2*dpdu1, xn, iwave, dtaudu2, dpdu2); - Real dtaudu3, dpdu3; + amrex::Real dtaudu3, dpdu3; riemann_invariant_rhs2(u+du2, tau+du2*dtaudu2, p+du2*dpdu2, xn, iwave, dtaudu3, dpdu3); - Real dtaudu4, dpdu4; + amrex::Real dtaudu4, dpdu4; riemann_invariant_rhs2(u+du, tau+du*dtaudu3, p+du*dpdu3, xn, iwave, dtaudu4, dpdu4); u += du; diff --git a/external/Microphysics b/external/Microphysics index 9d0655b75c..649063a33f 160000 --- a/external/Microphysics +++ b/external/Microphysics @@ -1 +1 @@ -Subproject commit 9d0655b75c2d7c0690fe637f8491fa572227a1dc +Subproject commit 649063a33f1ef83dd9acc457d7a79ceebcf3b60c diff --git a/external/amrex b/external/amrex index a068330e6c..52393d3faf 160000 --- a/external/amrex +++ b/external/amrex @@ -1 +1 @@ -Subproject commit a068330e6c66b5d9a7c6ca0e1c874f318e73f4cc +Subproject commit 52393d3fafc835fa379f0420aa2de4ccdddf155a