diff --git a/.clang-tidy b/.clang-tidy index 84742de4e..384e2db97 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -2,34 +2,58 @@ Checks: > '-*, bugprone-*, -bugprone-branch-clone, + -bugprone-casting-through-void, -bugprone-easily-swappable-parameters, -bugprone-exception-escape, -bugprone-implicit-widening-of-multiplication-result, + -bugprone-multi-level-implicit-pointer-conversion, clang-analyzer-*, -clang-analyzer-optin.mpi.MPI-Checker, clang-diagnostic-*, + -clang-diagnostic-deprecated-declarations, cppcoreguidelines-*, -cppcoreguidelines-avoid-c-arrays, - -cppcoreguidelines-avoid-goto, + -cppcoreguidelines-avoid-do-while, -cppcoreguidelines-avoid-magic-numbers, -cppcoreguidelines-avoid-non-const-global-variables, -cppcoreguidelines-init-variables, -cppcoreguidelines-interfaces-global-init, + -cppcoreguidelines-macro-to-enum, -cppcoreguidelines-macro-usage, -cppcoreguidelines-no-malloc, -cppcoreguidelines-non-private-member-variables-in-classes, -cppcoreguidelines-owning-memory, -cppcoreguidelines-pro-*, + -cppcoreguidelines-use-enum-class, + google-build-explicit-make-pair, + google-build-namespaces, + misc-*, + -misc-const-correctness, + -misc-include-cleaner, + -misc-non-private-member-variables-in-classes, + -misc-no-recursion, + -misc-use-internal-linkage, modernize-*, -modernize-avoid-c-arrays, -modernize-macro-to-enum, -modernize-return-braced-init-list, + -modernize-use-constraints, + -modernize-use-designated-initializers, + -modernize-use-integer-sign-comparison, + -modernize-use-ranges, + -modernize-use-starts-ends-with, + -modernize-use-std-numbers, -modernize-use-trailing-return-type, - -modernize-use-using, performance-*, + -performance-enum-size, + portability-*, + -portability-avoid-pragma-once, + -portability-template-virtual-member-function, readability-*, + -readability-avoid-nested-conditional-operator, + -readability-avoid-unconditional-preprocessor-if, -readability-braces-around-statements, - -readability-container-data-pointer, + -readability-container-contains, -readability-else-after-return, -readability-function-cognitive-complexity, -readability-function-size, @@ -37,10 +61,29 @@ Checks: > -readability-implicit-bool-conversion, -readability-isolate-declaration, -readability-magic-numbers, + -readability-math-missing-parentheses, -readability-named-parameter, + -readability-redundant-casting, + -readability-redundant-member-init, + -readability-static-accessed-through-instance, -readability-simplify-boolean-expr, + -readability-use-concise-preprocessor-directives, + -readability-use-std-min-max, mpi-* ' # Files not ending with nolint.H will be filtered in. HeaderFilterRegex: '([^n].....|[^o]....|[^l]...|[^i]..|[^n].|[^t])\.H$' + +HeaderFileExtensions: ['', "H", 'h', 'hh', 'hpp', 'hxx'] + +# We will try to modernize this in the future +# -modernize-use-constraints +# -modernize-use-designated-initializers +# -modernize-use-std-numbers +# -modernize-use-ranges +# -modernize-use-integer-sign-comparison +# -modernize-use-starts-ends-with +# -readability-container-contains +# -readability-braces-around-statements +# -google-global-names-in-headers diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 000000000..60a11755c --- /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/cuda.yml b/.github/workflows/cuda.yml index 6b5c0b2c2..320566c01 100644 --- a/.github/workflows/cuda.yml +++ b/.github/workflows/cuda.yml @@ -18,11 +18,11 @@ jobs: env: {CXXFLAGS: "-fno-operator-names -Werror -Wall -Wextra -Wpedantic -Wnull-dereference -Wfloat-conversion -Wshadow -Woverloaded-virtual -Wextra-semi -Wunreachable-code -Wnon-virtual-dtor -Wlogical-op -Wmisleading-indentation -Wduplicated-cond -Wduplicated-branches"} steps: - name: Checkout AMReX-Hydro - uses: actions/checkout@v5 + uses: actions/checkout@v6 with: path: AMReX-Hydro - name: Checkout AMReX - uses: actions/checkout@v5 + uses: actions/checkout@v6 with: repository: AMReX-Codes/amrex path: amrex @@ -92,11 +92,11 @@ jobs: env: {CXXFLAGS: "-fno-operator-names -Werror -Wall -Wextra -Wpedantic -Wnull-dereference -Wfloat-conversion -Wshadow -Woverloaded-virtual -Wextra-semi -Wunreachable-code -Wnon-virtual-dtor -Wlogical-op -Wmisleading-indentation -Wduplicated-cond -Wduplicated-branches"} steps: - name: Checkout AMReX-Hydro - uses: actions/checkout@v5 + uses: actions/checkout@v6 with: path: AMReX-Hydro - name: Checkout AMReX - uses: actions/checkout@v5 + uses: actions/checkout@v6 with: repository: AMReX-Codes/amrex path: amrex diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 8c2b03f78..1848f35a1 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -10,7 +10,7 @@ jobs: runs-on: ubuntu-latest steps: - name: Checkout - uses: actions/checkout@v5 + uses: actions/checkout@v6 with: persist-credentials: false @@ -19,9 +19,8 @@ jobs: .github/workflows/dependencies/documentation.sh echo "Installing python packages for docs..." - pip install breathe sphinx sphinx_rtd_theme sphinxcontrib.bibtex - pip install docutils - pip install sphinx-toolbox + python3 -m pip install --upgrade pip + python3 -m pip install sphinx==5.0.0 sphinx_rtd_theme sphinxcontrib.bibtex docutils sphinx-toolbox - name: Build Doxygen Docs run: | @@ -38,25 +37,25 @@ jobs: - name: Deploy if: github.event_name == 'push' && github.repository == 'AMReX-Fluids/AMReX-Hydro' && github.ref == 'refs/heads/development' - uses: JamesIves/github-pages-deploy-action@3.7.1 + uses: JamesIves/github-pages-deploy-action@v4.8.0 with: - GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} - ACCESS_TOKEN: ${{ secrets.DEPLOY_DOCS_HYDRO }} - REPOSITORY_NAME: AMReX-Fluids/AMReX-Fluids.github.io - BRANCH: main # The branch the action should deploy to. - FOLDER: Docs/build/html # The folder the action should deploy. - TARGET_FOLDER: amrex-hydro/docs_html # The folder the action should deploy to. - CLEAN: false # Do not remove existing files from the deploy target. + github-token: ${{ secrets.GITHUB_TOKEN }} + access-token: ${{ secrets.DEPLOY_DOCS_HYDRO }} + repository-name: AMReX-Fluids/AMReX-Fluids.github.io + branch: main # The branch the action should deploy to. + folder: Docs/build/html # The folder the action should deploy. + target-folder: amrex-hydro/docs_html # The folder the action should deploy to. + clean: false # Do not remove existing files from the deploy target. - name: Deploy Doxygen if: github.event_name == 'push' && github.repository == 'AMReX-Fluids/AMReX-Hydro' && github.ref == 'refs/heads/development' - uses: JamesIves/github-pages-deploy-action@3.7.1 + uses: JamesIves/github-pages-deploy-action@v4.8.0 with: - GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} - ACCESS_TOKEN: ${{ secrets.DEPLOY_DOCS_HYDRO }} - REPOSITORY_NAME: AMReX-Fluids/AMReX-Fluids.github.io - BRANCH: main # The branch the action should deploy to. - FOLDER: Docs/Doxygen/html # The folder the action should deploy. - TARGET_FOLDER: amrex-hydro/Doxygen/html # The folder the action should deploy to. - CLEAN: false # Do not remove existing files from the deploy target. + github-token: ${{ secrets.GITHUB_TOKEN }} + access-token: ${{ secrets.DEPLOY_DOCS_HYDRO }} + repository-name: AMReX-Fluids/AMReX-Fluids.github.io + branch: main # The branch the action should deploy to. + folder: Docs/Doxygen/html # The folder the action should deploy. + target-folder: amrex-hydro/Doxygen/html # The folder the action should deploy to. + clean: false # Do not remove existing files from the deploy target. diff --git a/.github/workflows/gcc.yml b/.github/workflows/gcc.yml index 9a8b6d49f..372f29730 100644 --- a/.github/workflows/gcc.yml +++ b/.github/workflows/gcc.yml @@ -18,11 +18,11 @@ jobs: env: {CXXFLAGS: "-fno-operator-names -Werror -Wall -Wextra -Wpedantic -Wnull-dereference -Wfloat-conversion -Wshadow -Woverloaded-virtual -Wnon-virtual-dtor -Wlogical-op -Wmisleading-indentation -Wduplicated-cond -Wduplicated-branches -Wmissing-include-dirs"} steps: - name: Checkout AMReX-Hydro - uses: actions/checkout@v5 + uses: actions/checkout@v6 with: path: AMReX-Hydro - name: Checkout AMReX - uses: actions/checkout@v5 + uses: actions/checkout@v6 with: repository: AMReX-Codes/amrex path: amrex @@ -36,7 +36,7 @@ jobs: - name: Load Dependencies run: | AMReX-Hydro/.github/workflows/dependencies/dependencies.sh - amrex/.github/workflows/dependencies/dependencies_clang-tidy.sh 14 + amrex/.github/workflows/dependencies/dependencies_clang-tidy-apt-llvm.sh 21 amrex/.github/workflows/dependencies/dependencies_ccache.sh - name: Build & Install run: | @@ -76,7 +76,7 @@ jobs: ${{github.workspace}}/amrex/Tools/C_scripts/mmclt.py --input ${{github.workspace}}/ccache.log.txt --identifier Godunov make -j2 -f clang-tidy-ccache-misses.mak \ - CLANG_TIDY=clang-tidy-14 \ + CLANG_TIDY=clang-tidy-21 \ CLANG_TIDY_ARGS="--config-file=${{github.workspace}}/AMReX-Hydro/.clang-tidy --warnings-as-errors=*" ccache -s @@ -88,11 +88,11 @@ jobs: env: {CXXFLAGS: "-fno-operator-names -Werror -Wall -Wextra -Wpedantic -Wnull-dereference -Wfloat-conversion -Wshadow -Woverloaded-virtual -Wnon-virtual-dtor -Wlogical-op -Wmisleading-indentation -Wduplicated-cond -Wduplicated-branches -Wmissing-include-dirs"} steps: - name: Checkout AMReX-Hydro - uses: actions/checkout@v5 + uses: actions/checkout@v6 with: path: AMReX-Hydro - name: Checkout AMReX - uses: actions/checkout@v5 + uses: actions/checkout@v6 with: repository: AMReX-Codes/amrex path: amrex @@ -106,7 +106,7 @@ jobs: - name: Load Dependencies run: | AMReX-Hydro/.github/workflows/dependencies/dependencies.sh - amrex/.github/workflows/dependencies/dependencies_clang-tidy.sh 14 + amrex/.github/workflows/dependencies/dependencies_clang-tidy-apt-llvm.sh 21 amrex/.github/workflows/dependencies/dependencies_ccache.sh - name: Build & Install run: | @@ -148,7 +148,7 @@ jobs: ${{github.workspace}}/amrex/Tools/C_scripts/mmclt.py --input ${{github.workspace}}/ccache.log.txt --identifier Godunov make -j2 -f clang-tidy-ccache-misses.mak \ - CLANG_TIDY=clang-tidy-14 \ + CLANG_TIDY=clang-tidy-21 \ CLANG_TIDY_ARGS="--config-file=${{github.workspace}}/AMReX-Hydro/.clang-tidy --warnings-as-errors=*" ccache -s diff --git a/.github/workflows/hip.yml b/.github/workflows/hip.yml index 52c46323b..4a872ed57 100644 --- a/.github/workflows/hip.yml +++ b/.github/workflows/hip.yml @@ -22,11 +22,11 @@ jobs: env: {CXXFLAGS: "-Werror -Wall -Wextra -Wpedantic -Wnull-dereference -Wfloat-conversion -Wshadow -Woverloaded-virtual -Wextra-semi -Wunreachable-code -Wnon-virtual-dtor -Wno-deprecated-declarations -Wno-gnu-zero-variadic-macro-arguments"} steps: - name: Checkout AMReX-Hydro - uses: actions/checkout@v5 + uses: actions/checkout@v6 with: path: AMReX-Hydro - name: Checkout AMReX - uses: actions/checkout@v5 + uses: actions/checkout@v6 with: repository: AMReX-Codes/amrex path: amrex diff --git a/.github/workflows/sycl.yml b/.github/workflows/sycl.yml index 57456e3a1..375c037f4 100644 --- a/.github/workflows/sycl.yml +++ b/.github/workflows/sycl.yml @@ -18,11 +18,11 @@ jobs: env: {CXXFLAGS: "-fno-operator-names -Werror -Wall -Wextra -Wpedantic -Wnull-dereference -Wfloat-conversion -Wshadow -Woverloaded-virtual -Wextra-semi -Wunreachable-code -Wnon-virtual-dtor -Wno-sign-compare -Wno-missing-braces -Wno-unused-variable -Wno-shadow"} steps: - name: Checkout AMReX-Hydro - uses: actions/checkout@v5 + uses: actions/checkout@v6 with: path: AMReX-Hydro - name: Checkout AMReX - uses: actions/checkout@v5 + uses: actions/checkout@v6 with: repository: AMReX-Codes/amrex path: amrex diff --git a/BDS/hydro_bds_edge_state_2D.cpp b/BDS/hydro_bds_edge_state_2D.cpp index 814803aa2..ecae4fc30 100644 --- a/BDS/hydro_bds_edge_state_2D.cpp +++ b/BDS/hydro_bds_edge_state_2D.cpp @@ -10,7 +10,7 @@ using namespace amrex; -constexpr amrex::Real eps = 1.0e-8; +constexpr amrex::Real eps = Real(1e-8); /** * Uses the Bell-Dawson-Shubin (BDS) algorithm, a higher order Godunov @@ -128,19 +128,19 @@ BDS::ComputeSlopes ( Box const& bx, // set node values equal to the average of the ghost cell values since they store the physical condition on the boundary if ( i<=dlo.x && lo_x_physbc ) { - sint(i,j,k) = 0.5*(s(dlo.x-1,j,k,icomp) + s(dlo.x-1,j-1,k,icomp)); + sint(i,j,k) = Real(0.5)*(s(dlo.x-1,j,k,icomp) + s(dlo.x-1,j-1,k,icomp)); return; } if ( i>=dhi.x+1 && hi_x_physbc ) { - sint(i,j,k) = 0.5*(s(dhi.x+1,j,k,icomp) + s(dhi.x+1,j-1,k,icomp)); + sint(i,j,k) = Real(0.5)*(s(dhi.x+1,j,k,icomp) + s(dhi.x+1,j-1,k,icomp)); return; } if ( j<=dlo.y && lo_y_physbc ) { - sint(i,j,k) = 0.5*(s(i,dlo.y-1,k,icomp) + s(i-1,dlo.y-1,k,icomp)); + sint(i,j,k) = Real(0.5)*(s(i,dlo.y-1,k,icomp) + s(i-1,dlo.y-1,k,icomp)); return; } if ( j>=dhi.y+1 && hi_y_physbc ) { - sint(i,j,k) = 0.5*(s(i,dhi.y+1,k,icomp) + s(i-1,dhi.y+1,k,icomp)); + sint(i,j,k) = Real(0.5)*(s(i,dhi.y+1,k,icomp) + s(i-1,dhi.y+1,k,icomp)); return; } @@ -150,14 +150,14 @@ BDS::ComputeSlopes ( Box const& bx, (j==dlo.y+1 && lo_y_physbc) || (j==dhi.y && hi_y_physbc) ) { - sint(i,j,k) = 0.25* (s(i,j,k,icomp) + s(i-1,j,k,icomp) + s(i,j-1,k,icomp) + s(i-1,j-1,k,icomp)); + sint(i,j,k) = Real(0.25)* (s(i,j,k,icomp) + s(i-1,j,k,icomp) + s(i,j-1,k,icomp) + s(i-1,j-1,k,icomp)); return; } sint(i,j,k) = (s(i-2,j-2,k,icomp) + s(i-2,j+1,k,icomp) + s(i+1,j-2,k,icomp) + s(i+1,j+1,k,icomp) - - 7.0*(s(i-2,j-1,k,icomp) + s(i-2,j ,k,icomp) + s(i-1,j-2,k,icomp) + s(i ,j-2,k,icomp) + - s(i-1,j+1,k,icomp) + s(i ,j+1,k,icomp) + s(i+1,j-1,k,icomp) + s(i+1,j ,k,icomp)) - + 49.0*(s(i-1,j-1,k,icomp) + s(i ,j-1,k,icomp) + s(i-1,j ,k,icomp) + s(i ,j ,k,icomp)) ) / 144.0; + - Real(7)*(s(i-2,j-1,k,icomp) + s(i-2,j ,k,icomp) + s(i-1,j-2,k,icomp) + s(i ,j-2,k,icomp) + + s(i-1,j+1,k,icomp) + s(i ,j+1,k,icomp) + s(i+1,j-1,k,icomp) + s(i+1,j ,k,icomp)) + + Real(49)*(s(i-1,j-1,k,icomp) + s(i ,j-1,k,icomp) + s(i-1,j ,k,icomp) + s(i ,j ,k,icomp)) ) / Real(144); }); ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k){ @@ -173,25 +173,25 @@ BDS::ComputeSlopes ( Box const& bx, // compute initial estimates of slopes from unlimited corner points // sx - slopes(i,j,k,0) = 0.5*(sint(i+1,j+1,k) + sint(i+1,j,k) - sint(i,j+1,k) - sint(i,j,k)) / hx; + slopes(i,j,k,0) = Real(0.5)*(sint(i+1,j+1,k) + sint(i+1,j,k) - sint(i,j+1,k) - sint(i,j,k)) / hx; // sy - slopes(i,j,k,1) = 0.5*(sint(i+1,j+1,k) - sint(i+1,j,k) + sint(i,j+1,k) - sint(i,j,k)) / hy; + slopes(i,j,k,1) = Real(0.5)*(sint(i+1,j+1,k) - sint(i+1,j,k) + sint(i,j+1,k) - sint(i,j,k)) / hy; // sxy slopes(i,j,k,2) = (sint(i+1,j+1,k) - sint(i+1,j,k) - sint(i,j+1,k) + sint(i,j,k)) / (hx*hy); if (bds_limiter_type > 0) { // ++ / sint(i+1,j+1) - sc(4) = s(i,j,k,icomp) + 0.5*(hx*slopes(i,j,k,0) + hy*slopes(i,j,k,1)) + 0.25*hx*hy*slopes(i,j,k,2); + sc(4) = s(i,j,k,icomp) + Real(0.5)*(hx*slopes(i,j,k,0) + hy*slopes(i,j,k,1)) + Real(0.25)*hx*hy*slopes(i,j,k,2); // +- / sint(i+1,j ) - sc(3) = s(i,j,k,icomp) + 0.5*(hx*slopes(i,j,k,0) - hy*slopes(i,j,k,1)) - 0.25*hx*hy*slopes(i,j,k,2); + sc(3) = s(i,j,k,icomp) + Real(0.5)*(hx*slopes(i,j,k,0) - hy*slopes(i,j,k,1)) - Real(0.25)*hx*hy*slopes(i,j,k,2); // -+ / sint(i ,j+1) - sc(2) = s(i,j,k,icomp) - 0.5*(hx*slopes(i,j,k,0) - hy*slopes(i,j,k,1)) - 0.25*hx*hy*slopes(i,j,k,2); + sc(2) = s(i,j,k,icomp) - Real(0.5)*(hx*slopes(i,j,k,0) - hy*slopes(i,j,k,1)) - Real(0.25)*hx*hy*slopes(i,j,k,2); // -- / sint(i ,j ) - sc(1) = s(i,j,k,icomp) - 0.5*(hx*slopes(i,j,k,0) + hy*slopes(i,j,k,1)) + 0.25*hx*hy*slopes(i,j,k,2); + sc(1) = s(i,j,k,icomp) - Real(0.5)*(hx*slopes(i,j,k,0) + hy*slopes(i,j,k,1)) + Real(0.25)*hx*hy*slopes(i,j,k,2); // enforce max/min bounds smin(4) = amrex::min(s(i,j,k,icomp), s(i+1,j,k,icomp), s(i,j+1,k,icomp), s(i+1,j+1,k,icomp)); @@ -214,11 +214,11 @@ BDS::ComputeSlopes ( Box const& bx, for(int ll=1; ll<=3; ++ll){ // compute the amount by which the average of the nodal values differs from cell-center value - sumloc = 0.25*(sc(4) + sc(3) + sc(2) + sc(1)); - sumdif = (sumloc - s(i,j,k,icomp))*4.0; + sumloc = Real(0.25)*(sc(4) + sc(3) + sc(2) + sc(1)); + sumdif = (sumloc - s(i,j,k,icomp))*Real(4); // sgndif = +(-)1 if the node average is too large(small) - sgndif = std::copysign(1.0,sumdif); + sgndif = std::copysign(Real(1),sumdif); // compute how much each node is larger(smaller) than the cell-centered value for(int mm=1; mm<=4; ++mm){ @@ -241,7 +241,7 @@ BDS::ComputeSlopes ( Box const& bx, // how many node values are left to potentially adjust if (kdp<1) { - div = 1.0; + div = Real(1); } else { div = kdp; } @@ -295,9 +295,9 @@ BDS::ComputeSlopes ( Box const& bx, // final slopes // sx - slopes(i,j,k,0) = 0.5*( sc(4) + sc(3) -sc(1) - sc(2) )/hx; + slopes(i,j,k,0) = Real(0.5)*( sc(4) + sc(3) -sc(1) - sc(2) )/hx; // sy - slopes(i,j,k,1) = 0.5*( sc(4) + sc(2) -sc(1) - sc(3) )/hy; + slopes(i,j,k,1) = Real(0.5)*( sc(4) + sc(2) -sc(1) - sc(3) )/hy; // sxy slopes(i,j,k,2) = ( sc(1) + sc(4) -sc(2) - sc(3) )/(hx*hy); } @@ -376,8 +376,8 @@ BDS::ComputeConc (Box const& bx, Real hx = dx[0]; Real hy = dx[1]; - Real dt2 = dt/2.0; - Real dt3 = dt/3.0; + Real dt2 = dt/Real(2); + Real dt3 = dt/Real(3); Box const& domain = geom.Domain(); const auto dlo = amrex::lbound(domain); @@ -445,10 +445,10 @@ BDS::ComputeConc (Box const& bx, /////////////////////////////////////////////// if (umac(i,j,k) > 0) { - isign = 1.; + isign = Real(1); ioff = -1; } else { - isign = -1.; + isign = -Real(1); ioff = 0; } @@ -457,15 +457,15 @@ BDS::ComputeConc (Box const& bx, } // centroid of rectangular volume - del(1) = isign*0.5*hx - 0.5*umac(i,j,k)*dt; + del(1) = isign*Real(0.5)*hx - Real(0.5)*umac(i,j,k)*dt; del(2) = 0.; xedge_tmp = eval(s(i+ioff,j,k,icomp),slope_tmp,del); // source term if (iconserv[icomp]) { - xedge_tmp = xedge_tmp*(1. - dt2*ux(i+ioff,j,k)); + xedge_tmp = xedge_tmp*(Real(1)- dt2*ux(i+ioff,j,k)); } else { - xedge_tmp = xedge_tmp*(1. + dt2*vy(i+ioff,j,k)); + xedge_tmp = xedge_tmp*(Real(1)+ dt2*vy(i+ioff,j,k)); } if (force) { xedge_tmp += dt2*force(i+ioff,j,k,icomp); @@ -476,10 +476,10 @@ BDS::ComputeConc (Box const& bx, /////////////////////////////////////////////// if (vmac(i+ioff,j+1,k) > 0.) { - jsign = 1.; + jsign = Real(1); joff = 0; } else { - jsign = -1.; + jsign = -Real(1); joff = 1; } @@ -488,40 +488,40 @@ BDS::ComputeConc (Box const& bx, u = umac(i,j+joff,k); } - p1(1) = isign*0.5*hx; - p1(2) = jsign*0.5*hy; + p1(1) = isign*Real(0.5)*hx; + p1(2) = jsign*Real(0.5)*hy; - p2(1) = isign*0.5*hx - umac(i,j,k)*dt; - p2(2) = jsign*0.5*hy; + p2(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p2(2) = jsign*Real(0.5)*hy; - p3(1) = isign*0.5*hx - u*dt; - p3(2) = jsign*0.5*hy - vmac(i+ioff,j+1,k)*dt; + p3(1) = isign*Real(0.5)*hx - u*dt; + p3(2) = jsign*Real(0.5)*hy - vmac(i+ioff,j+1,k)*dt; for(int n=1; n<=3; ++n){ slope_tmp(n) = slopes(i+ioff,j+joff,k,n-1); } for (int ll=1; ll<=2; ++ll) { - del(ll) = (p2(ll)+p3(ll))/2.; + del(ll) = (p2(ll)+p3(ll))/Real(2); } val1 = eval(s(i+ioff,j+joff,k,icomp),slope_tmp,del); for (int ll=1; ll<=2; ++ll) { - del(ll) = (p1(ll)+p3(ll))/2.; + del(ll) = (p1(ll)+p3(ll))/Real(2); } val2 = eval(s(i+ioff,j+joff,k,icomp),slope_tmp,del); for (int ll=1; ll<=2; ++ll) { - del(ll) = (p1(ll)+p2(ll))/2.; + del(ll) = (p1(ll)+p2(ll))/Real(2); } val3 = eval(s(i+ioff,j+joff,k,icomp),slope_tmp,del); // average these centroid values to get the average value - gamma = (val1+val2+val3)/3.; + gamma = (val1+val2+val3)/Real(3); // source term if (iconserv[icomp]) { - gamma = gamma*(1. - dt3*divu(i+ioff,j+joff,k)); + gamma = gamma*(Real(1)- dt3*divu(i+ioff,j+joff,k)); } /////////////////////////////////////////////// @@ -530,17 +530,17 @@ BDS::ComputeConc (Box const& bx, gamma = gamma * vmac(i+ioff,j+1,k); - xedge_tmp = xedge_tmp - dt*gamma/(2.*hy); + xedge_tmp = xedge_tmp - dt*gamma/(Real(2)*hy); /////////////////////////////////////////////// // compute \Gamma^{y-} /////////////////////////////////////////////// if (vmac(i+ioff,j,k) > 0.) { - jsign = 1.; + jsign = Real(1); joff = -1; } else { - jsign = -1.; + jsign = -Real(1); joff = 0; } @@ -549,40 +549,40 @@ BDS::ComputeConc (Box const& bx, u = umac(i,j+joff,k); } - p1(1) = isign*0.5*hx; - p1(2) = jsign*0.5*hy; + p1(1) = isign*Real(0.5)*hx; + p1(2) = jsign*Real(0.5)*hy; - p2(1) = isign*0.5*hx - umac(i,j,k)*dt; - p2(2) = jsign*0.5*hy; + p2(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p2(2) = jsign*Real(0.5)*hy; - p3(1) = isign*0.5*hx - u*dt; - p3(2) = jsign*0.5*hy - vmac(i+ioff,j,k)*dt; + p3(1) = isign*Real(0.5)*hx - u*dt; + p3(2) = jsign*Real(0.5)*hy - vmac(i+ioff,j,k)*dt; for(int n=1; n<=3; ++n){ slope_tmp(n) = slopes(i+ioff,j+joff,k,n-1); } for (int ll=1; ll<=2; ++ll) { - del(ll) = (p2(ll)+p3(ll))/2.; + del(ll) = (p2(ll)+p3(ll))/Real(2); } val1 = eval(s(i+ioff,j+joff,k,icomp),slope_tmp,del); for (int ll=1; ll<=2; ++ll) { - del(ll) = (p1(ll)+p3(ll))/2.; + del(ll) = (p1(ll)+p3(ll))/Real(2); } val2 = eval(s(i+ioff,j+joff,k,icomp),slope_tmp,del); for (int ll=1; ll<=2; ++ll) { - del(ll) = (p1(ll)+p2(ll))/2.; + del(ll) = (p1(ll)+p2(ll))/Real(2); } val3 = eval(s(i+ioff,j+joff,k,icomp),slope_tmp,del); // average these centroid values to get the average value - gamma = (val1+val2+val3)/3.; + gamma = (val1+val2+val3)/Real(3); // source term if (iconserv[icomp]) { - gamma = gamma*(1. - dt3*divu(i+ioff,j+joff,k)); + gamma = gamma*(Real(1)- dt3*divu(i+ioff,j+joff,k)); } /////////////////////////////////////////////// @@ -590,7 +590,7 @@ BDS::ComputeConc (Box const& bx, /////////////////////////////////////////////// gamma = gamma * vmac(i+ioff,j,k); - sedgex(i,j,k,icomp) = xedge_tmp + dt*gamma/(2.*hy); + sedgex(i,j,k,icomp) = xedge_tmp + dt*gamma/(Real(2)*hy); }); // compute sedgey on y-faces @@ -640,10 +640,10 @@ BDS::ComputeConc (Box const& bx, // centroid of rectangular volume if (vmac(i,j,k) > 0.) { - jsign = 1.; + jsign = Real(1); joff = -1; } else { - jsign = -1.; + jsign = -Real(1); joff = 0; } @@ -652,14 +652,14 @@ BDS::ComputeConc (Box const& bx, } del(1) = 0.; - del(2) = jsign*0.5*hy - 0.5*vmac(i,j,k)*dt; + del(2) = jsign*Real(0.5)*hy - Real(0.5)*vmac(i,j,k)*dt; yedge_tmp = eval(s(i,j+joff,k,icomp),slope_tmp,del); // source term if (iconserv[icomp]) { - yedge_tmp = yedge_tmp*(1. - dt2*vy(i,j+joff,k)); + yedge_tmp = yedge_tmp*(Real(1)- dt2*vy(i,j+joff,k)); } else { - yedge_tmp = yedge_tmp*(1. + dt2*ux(i,j+joff,k)); + yedge_tmp = yedge_tmp*(Real(1)+ dt2*ux(i,j+joff,k)); } if (force) { yedge_tmp += dt2*force(i,j+joff,k,icomp); @@ -670,10 +670,10 @@ BDS::ComputeConc (Box const& bx, /////////////////////////////////////////////// if (umac(i+1,j+joff,k) > 0.) { - isign = 1.; + isign = Real(1); ioff = 0; } else { - isign = -1.; + isign = -Real(1); ioff = 1; } @@ -682,40 +682,40 @@ BDS::ComputeConc (Box const& bx, v = vmac(i+ioff,j,k); } - p1(1) = isign*0.5*hx; - p1(2) = jsign*0.5*hy; + p1(1) = isign*Real(0.5)*hx; + p1(2) = jsign*Real(0.5)*hy; - p2(1) = isign*0.5*hx; - p2(2) = jsign*0.5*hy - vmac(i,j,k)*dt; + p2(1) = isign*Real(0.5)*hx; + p2(2) = jsign*Real(0.5)*hy - vmac(i,j,k)*dt; - p3(1) = isign*0.5*hx - umac(i+1,j+joff,k)*dt; - p3(2) = jsign*0.5*hy - v*dt; + p3(1) = isign*Real(0.5)*hx - umac(i+1,j+joff,k)*dt; + p3(2) = jsign*Real(0.5)*hy - v*dt; for(int n=1; n<=3; ++n){ slope_tmp(n) = slopes(i+ioff,j+joff,k,n-1); } for (int ll=1; ll<=2; ++ll) { - del(ll) = (p2(ll)+p3(ll))/2.; + del(ll) = (p2(ll)+p3(ll))/Real(2); } val1 = eval(s(i+ioff,j+joff,k,icomp),slope_tmp,del); for (int ll=1; ll<=2; ++ll) { - del(ll) = (p1(ll)+p3(ll))/2.; + del(ll) = (p1(ll)+p3(ll))/Real(2); } val2 = eval(s(i+ioff,j+joff,k,icomp),slope_tmp,del); for (int ll=1; ll<=2; ++ll) { - del(ll) = (p1(ll)+p2(ll))/2.; + del(ll) = (p1(ll)+p2(ll))/Real(2); } val3 = eval(s(i+ioff,j+joff,k,icomp),slope_tmp,del); // average these centroid values to get the average value - gamma = (val1+val2+val3)/3.; + gamma = (val1+val2+val3)/Real(3); // source term if (iconserv[icomp]) { - gamma = gamma*(1. - dt3*divu(i+ioff,j+joff,k)); + gamma = gamma*(Real(1)- dt3*divu(i+ioff,j+joff,k)); } /////////////////////////////////////////////// @@ -723,17 +723,17 @@ BDS::ComputeConc (Box const& bx, /////////////////////////////////////////////// gamma = gamma * umac(i+1,j+joff,k); - yedge_tmp = yedge_tmp - dt*gamma/(2.*hx); + yedge_tmp = yedge_tmp - dt*gamma/(Real(2)*hx); /////////////////////////////////////////////// // compute \Gamma^{x-} /////////////////////////////////////////////// if (umac(i,j+joff,k) > 0.) { - isign = 1.; + isign = Real(1); ioff = -1; } else { - isign = -1.; + isign = -Real(1); ioff = 0; } @@ -742,40 +742,40 @@ BDS::ComputeConc (Box const& bx, v = vmac(i+ioff,j,k); } - p1(1) = isign*0.5*hx; - p1(2) = jsign*0.5*hy; + p1(1) = isign*Real(0.5)*hx; + p1(2) = jsign*Real(0.5)*hy; - p2(1) = isign*0.5*hx; - p2(2) = jsign*0.5*hy - vmac(i,j,k)*dt; + p2(1) = isign*Real(0.5)*hx; + p2(2) = jsign*Real(0.5)*hy - vmac(i,j,k)*dt; - p3(1) = isign*0.5*hx - umac(i,j+joff,k)*dt; - p3(2) = jsign*0.5*hy - v*dt; + p3(1) = isign*Real(0.5)*hx - umac(i,j+joff,k)*dt; + p3(2) = jsign*Real(0.5)*hy - v*dt; for(int n=1; n<=3; ++n){ slope_tmp(n) = slopes(i+ioff,j+joff,k,n-1); } for (int ll=1; ll<=2; ++ll) { - del(ll) = (p2(ll)+p3(ll))/2.; + del(ll) = (p2(ll)+p3(ll))/Real(2); } val1 = eval(s(i+ioff,j+joff,k,icomp),slope_tmp,del); for (int ll=1; ll<=2; ++ll) { - del(ll) = (p1(ll)+p3(ll))/2.; + del(ll) = (p1(ll)+p3(ll))/Real(2); } val2 = eval(s(i+ioff,j+joff,k,icomp),slope_tmp,del); for (int ll=1; ll<=2; ++ll) { - del(ll) = (p1(ll)+p2(ll))/2.; + del(ll) = (p1(ll)+p2(ll))/Real(2); } val3 = eval(s(i+ioff,j+joff,k,icomp),slope_tmp,del); // average these centroid values to get the average value - gamma = (val1+val2+val3)/3.; + gamma = (val1+val2+val3)/Real(3); // source term if (iconserv[icomp]) { - gamma = gamma*(1. - dt3*divu(i+ioff,j+joff,k)); + gamma = gamma*(Real(1)- dt3*divu(i+ioff,j+joff,k)); } /////////////////////////////////////////////// @@ -783,7 +783,7 @@ BDS::ComputeConc (Box const& bx, /////////////////////////////////////////////// gamma = gamma * umac(i,j+joff,k); - sedgey(i,j,k,icomp) = yedge_tmp + dt*gamma/(2.*hx); + sedgey(i,j,k,icomp) = yedge_tmp + dt*gamma/(Real(2)*hx); }); } diff --git a/BDS/hydro_bds_edge_state_3D.cpp b/BDS/hydro_bds_edge_state_3D.cpp index 15aafa190..5baedcc1a 100644 --- a/BDS/hydro_bds_edge_state_3D.cpp +++ b/BDS/hydro_bds_edge_state_3D.cpp @@ -214,39 +214,39 @@ BDS::ComputeSlopes ( Box const& bx, // compute initial estimates of slopes from unlimited corner points // sx - slopes(i,j,k,0) = 0.25*(( sint(i+1,j ,k ) + sint(i+1,j+1,k ) - +sint(i+1,j ,k+1) + sint(i+1,j+1,k+1) ) - -( sint(i ,j ,k ) + sint(i ,j+1,k ) - +sint(i ,j ,k+1) + sint(i ,j+1,k+1) )) / hx; + slopes(i,j,k,0) = Real(0.25)*(( sint(i+1,j ,k ) + sint(i+1,j+1,k ) + +sint(i+1,j ,k+1) + sint(i+1,j+1,k+1) ) + -( sint(i ,j ,k ) + sint(i ,j+1,k ) + +sint(i ,j ,k+1) + sint(i ,j+1,k+1) )) / hx; // sy - slopes(i,j,k,1) = 0.25*(( sint(i ,j+1,k ) + sint(i+1,j+1,k ) - +sint(i ,j+1,k+1) + sint(i+1,j+1,k+1) ) - -( sint(i ,j ,k ) + sint(i+1,j ,k ) - +sint(i ,j ,k+1) + sint(i+1,j ,k+1) )) / hy; + slopes(i,j,k,1) = Real(0.25)*(( sint(i ,j+1,k ) + sint(i+1,j+1,k ) + +sint(i ,j+1,k+1) + sint(i+1,j+1,k+1) ) + -( sint(i ,j ,k ) + sint(i+1,j ,k ) + +sint(i ,j ,k+1) + sint(i+1,j ,k+1) )) / hy; // sz - slopes(i,j,k,2) = 0.25*(( sint(i ,j ,k+1) + sint(i+1,j ,k+1) - +sint(i ,j+1,k+1) + sint(i+1,j+1,k+1) ) - -( sint(i ,j ,k ) + sint(i+1,j ,k ) - +sint(i ,j+1,k ) + sint(i+1,j+1,k ) )) / hz; + slopes(i,j,k,2) = Real(0.25)*(( sint(i ,j ,k+1) + sint(i+1,j ,k+1) + +sint(i ,j+1,k+1) + sint(i+1,j+1,k+1) ) + -( sint(i ,j ,k ) + sint(i+1,j ,k ) + +sint(i ,j+1,k ) + sint(i+1,j+1,k ) )) / hz; // sxy - slopes(i,j,k,3) = 0.5*( ( sint(i ,j ,k ) + sint(i ,j ,k+1) - +sint(i+1,j+1,k ) + sint(i+1,j+1,k+1) ) - -( sint(i+1,j ,k ) + sint(i+1,j ,k+1) - +sint(i ,j+1,k ) + sint(i ,j+1,k+1) )) / (hx*hy); + slopes(i,j,k,3) = Real(0.5)*( ( sint(i ,j ,k ) + sint(i ,j ,k+1) + +sint(i+1,j+1,k ) + sint(i+1,j+1,k+1) ) + -( sint(i+1,j ,k ) + sint(i+1,j ,k+1) + +sint(i ,j+1,k ) + sint(i ,j+1,k+1) )) / (hx*hy); // sxz - slopes(i,j,k,4) = 0.5*( ( sint(i ,j ,k ) + sint(i ,j+1,k ) - +sint(i+1,j ,k+1) + sint(i+1,j+1,k+1) ) - -( sint(i+1,j ,k ) + sint(i+1,j+1,k ) - +sint(i ,j ,k+1) + sint(i ,j+1,k+1) )) / (hx*hz); + slopes(i,j,k,4) = Real(0.5)*( ( sint(i ,j ,k ) + sint(i ,j+1,k ) + +sint(i+1,j ,k+1) + sint(i+1,j+1,k+1) ) + -( sint(i+1,j ,k ) + sint(i+1,j+1,k ) + +sint(i ,j ,k+1) + sint(i ,j+1,k+1) )) / (hx*hz); // syz - slopes(i,j,k,5) = 0.5*( ( sint(i ,j ,k ) + sint(i+1,j ,k ) - +sint(i ,j+1,k+1) + sint(i+1,j+1,k+1) ) - -( sint(i ,j ,k+1) + sint(i+1,j ,k+1) - +sint(i ,j+1,k ) + sint(i+1,j+1,k ) )) / (hy*hz); + slopes(i,j,k,5) = Real(0.5)*( ( sint(i ,j ,k ) + sint(i+1,j ,k ) + +sint(i ,j+1,k+1) + sint(i+1,j+1,k+1) ) + -( sint(i ,j ,k+1) + sint(i+1,j ,k+1) + +sint(i ,j+1,k ) + sint(i+1,j+1,k ) )) / (hy*hz); // sxyz slopes(i,j,k,6) = (-sint(i ,j ,k ) + sint(i+1,j ,k ) + sint(i ,j+1,k ) @@ -257,51 +257,51 @@ BDS::ComputeSlopes ( Box const& bx, // +++ / sint(i+1,j+1,k+1) sc(8) = s(i,j,k,icomp) - +0.5 *( hx*slopes(i,j,k,0)+ hy*slopes(i,j,k,1)+ hz*slopes(i,j,k,2)) - +0.25 *( hx*hy*slopes(i,j,k,3)+hx*hz*slopes(i,j,k,4)+hy*hz*slopes(i,j,k,5)) - +0.125*hx*hy*hz*slopes(i,j,k,6); + +Real(0.5 )*( hx*slopes(i,j,k,0)+ hy*slopes(i,j,k,1)+ hz*slopes(i,j,k,2)) + +Real(0.25 )*( hx*hy*slopes(i,j,k,3)+hx*hz*slopes(i,j,k,4)+hy*hz*slopes(i,j,k,5)) + +Real(0.125)*hx*hy*hz*slopes(i,j,k,6); // ++- / sint(i+1,j+1,k ) sc(7) = s(i,j,k,icomp) - +0.5 *( hx*slopes(i,j,k,0)+ hy*slopes(i,j,k,1)- hz*slopes(i,j,k,2)) - +0.25 *( hx*hy*slopes(i,j,k,3)-hx*hz*slopes(i,j,k,4)-hy*hz*slopes(i,j,k,5)) - -0.125*hx*hy*hz*slopes(i,j,k,6); + +Real(0.5 )*( hx*slopes(i,j,k,0)+ hy*slopes(i,j,k,1)- hz*slopes(i,j,k,2)) + +Real(0.25 )*( hx*hy*slopes(i,j,k,3)-hx*hz*slopes(i,j,k,4)-hy*hz*slopes(i,j,k,5)) + -Real(0.125)*hx*hy*hz*slopes(i,j,k,6); // +-+ / sint(i+1,j ,k+1) sc(6) = s(i,j,k,icomp) - +0.5 *( hx*slopes(i,j,k,0)- hy*slopes(i,j,k,1)+ hz*slopes(i,j,k,2)) - +0.25 *( -hx*hy*slopes(i,j,k,3)+hx*hz*slopes(i,j,k,4)-hy*hz*slopes(i,j,k,5)) - -0.125*hx*hy*hz*slopes(i,j,k,6); + +Real(0.5 )*( hx*slopes(i,j,k,0)- hy*slopes(i,j,k,1)+ hz*slopes(i,j,k,2)) + +Real(0.25 )*( -hx*hy*slopes(i,j,k,3)+hx*hz*slopes(i,j,k,4)-hy*hz*slopes(i,j,k,5)) + -Real(0.125)*hx*hy*hz*slopes(i,j,k,6); // +-- / sint(i+1,j ,k ) sc(5) = s(i,j,k,icomp) - +0.5 *( hx*slopes(i,j,k,0)- hy*slopes(i,j,k,1)- hz*slopes(i,j,k,2)) - +0.25 *( -hx*hy*slopes(i,j,k,3)-hx*hz*slopes(i,j,k,4)+hy*hz*slopes(i,j,k,5)) - +0.125*hx*hy*hz*slopes(i,j,k,6); + +Real(0.5 )*( hx*slopes(i,j,k,0)- hy*slopes(i,j,k,1)- hz*slopes(i,j,k,2)) + +Real(0.25 )*( -hx*hy*slopes(i,j,k,3)-hx*hz*slopes(i,j,k,4)+hy*hz*slopes(i,j,k,5)) + +Real(0.125)*hx*hy*hz*slopes(i,j,k,6); // -++ / sint(i ,j+1,k+1) sc(4) = s(i,j,k,icomp) - +0.5 *( -hx*slopes(i,j,k,0)+ hy*slopes(i,j,k,1)+ hz*slopes(i,j,k,2)) - +0.25 *( -hx*hy*slopes(i,j,k,3)-hx*hz*slopes(i,j,k,4)+hy*hz*slopes(i,j,k,5)) - -0.125*hx*hy*hz*slopes(i,j,k,6); + +Real(0.5 )*( -hx*slopes(i,j,k,0)+ hy*slopes(i,j,k,1)+ hz*slopes(i,j,k,2)) + +Real(0.25 )*( -hx*hy*slopes(i,j,k,3)-hx*hz*slopes(i,j,k,4)+hy*hz*slopes(i,j,k,5)) + -Real(0.125)*hx*hy*hz*slopes(i,j,k,6); // -+- / sint(i ,j+1,k ) sc(3) = s(i,j,k,icomp) - +0.5 *( -hx*slopes(i,j,k,0)+ hy*slopes(i,j,k,1)- hz*slopes(i,j,k,2)) - +0.25 *( -hx*hy*slopes(i,j,k,3)+hx*hz*slopes(i,j,k,4)-hy*hz*slopes(i,j,k,5)) - +0.125*hx*hy*hz*slopes(i,j,k,6); + +Real(0.5 )*( -hx*slopes(i,j,k,0)+ hy*slopes(i,j,k,1)- hz*slopes(i,j,k,2)) + +Real(0.25 )*( -hx*hy*slopes(i,j,k,3)+hx*hz*slopes(i,j,k,4)-hy*hz*slopes(i,j,k,5)) + +Real(0.125)*hx*hy*hz*slopes(i,j,k,6); // --+ / sint(i ,j ,k+1) sc(2) = s(i,j,k,icomp) - +0.5 *( -hx*slopes(i,j,k,0)- hy*slopes(i,j,k,1)+ hz*slopes(i,j,k,2)) - +0.25 *( hx*hy*slopes(i,j,k,3)-hx*hz*slopes(i,j,k,4)-hy*hz*slopes(i,j,k,5)) - +0.125*hx*hy*hz*slopes(i,j,k,6); + +Real(0.5 )*( -hx*slopes(i,j,k,0)- hy*slopes(i,j,k,1)+ hz*slopes(i,j,k,2)) + +Real(0.25 )*( hx*hy*slopes(i,j,k,3)-hx*hz*slopes(i,j,k,4)-hy*hz*slopes(i,j,k,5)) + +Real(0.125)*hx*hy*hz*slopes(i,j,k,6); // ---/ sint(i ,j ,k ) sc(1) = s(i,j,k,icomp) - +0.5 *( -hx*slopes(i,j,k,0)- hy*slopes(i,j,k,1)- hz*slopes(i,j,k,2)) - +0.25 *( hx*hy*slopes(i,j,k,3)+hx*hz*slopes(i,j,k,4)+hy*hz*slopes(i,j,k,5)) - -0.125*hx*hy*hz*slopes(i,j,k,6); + +Real(0.5 )*( -hx*slopes(i,j,k,0)- hy*slopes(i,j,k,1)- hz*slopes(i,j,k,2)) + +Real(0.25 )*( hx*hy*slopes(i,j,k,3)+hx*hz*slopes(i,j,k,4)+hy*hz*slopes(i,j,k,5)) + -Real(0.125)*hx*hy*hz*slopes(i,j,k,6); // enforce max/min bounds smin(8) = min(s(i ,j ,k ,icomp),s(i+1,j ,k ,icomp),s(i ,j+1,k ,icomp),s(i ,j ,k+1,icomp), @@ -434,40 +434,40 @@ BDS::ComputeSlopes ( Box const& bx, // final slopes // sx - slopes(i,j,k,0) = 0.25*( ( sc(5) + sc(7) - +sc(6) + sc(8)) - -( sc(1) + sc(3) - +sc(2) + sc(4)) ) / hx; + slopes(i,j,k,0) = Real(0.25)*( ( sc(5) + sc(7) + +sc(6) + sc(8)) + -( sc(1) + sc(3) + +sc(2) + sc(4)) ) / hx; // sy - slopes(i,j,k,1) = 0.25*( ( sc(3) + sc(7) - +sc(4) + sc(8)) - -( sc(1) + sc(5) - +sc(2) + sc(6)) ) / hy; + slopes(i,j,k,1) = Real(0.25)*( ( sc(3) + sc(7) + +sc(4) + sc(8)) + -( sc(1) + sc(5) + +sc(2) + sc(6)) ) / hy; // sz - slopes(i,j,k,2) = 0.25*( ( sc(2) + sc(6) - +sc(4) + sc(8)) - -( sc(1) + sc(5) - +sc(3) + sc(7)) ) / hz; + slopes(i,j,k,2) = Real(0.25)*( ( sc(2) + sc(6) + +sc(4) + sc(8)) + -( sc(1) + sc(5) + +sc(3) + sc(7)) ) / hz; // sxy - slopes(i,j,k,3) = 0.5*( ( sc(1) + sc(2) - +sc(7) + sc(8)) - -( sc(5) + sc(6) - +sc(3) + sc(4)) ) / (hx*hy); + slopes(i,j,k,3) = Real(0.5)*( ( sc(1) + sc(2) + +sc(7) + sc(8)) + -( sc(5) + sc(6) + +sc(3) + sc(4)) ) / (hx*hy); // sxz - slopes(i,j,k,4) = 0.5*( ( sc(1) + sc(3) - +sc(6) + sc(8)) - -( sc(5) + sc(7) - +sc(2) + sc(4)) ) / (hx*hz); + slopes(i,j,k,4) = Real(0.5)*( ( sc(1) + sc(3) + +sc(6) + sc(8)) + -( sc(5) + sc(7) + +sc(2) + sc(4)) ) / (hx*hz); // syz - slopes(i,j,k,5) = 0.5*( ( sc(1) + sc(5) - +sc(4) + sc(8)) - -( sc(2) + sc(6) - +sc(3) + sc(7)) ) / (hy*hz); + slopes(i,j,k,5) = Real(0.5)*( ( sc(1) + sc(5) + +sc(4) + sc(8)) + -( sc(2) + sc(6) + +sc(3) + sc(7)) ) / (hy*hz); // sxyz slopes(i,j,k,6) = (-sc(1) + sc(5) + sc(3) @@ -650,7 +650,7 @@ BDS::ComputeConc (Box const& bx, // centroid of rectangular volume - del(1) = isign*0.5*hx - 0.5*umac(i,j,k)*dt; + del(1) = isign*Real(0.5)*hx - Real(0.5)*umac(i,j,k)*dt; del(2) = 0.0; del(3) = 0.0; xedge_tmp = eval(s(i+ioff,j,k,icomp),slope_tmp,del); @@ -683,16 +683,16 @@ BDS::ComputeConc (Box const& bx, u = umac(i,j+joff,k); } - p1(1) = isign*0.5*hx; - p1(2) = jsign*0.5*hy; + p1(1) = isign*Real(0.5)*hx; + p1(2) = jsign*Real(0.5)*hy; p1(3) = 0.0; - p2(1) = isign*0.5*hx - umac(i,j,k)*dt; - p2(2) = jsign*0.5*hy; + p2(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p2(2) = jsign*Real(0.5)*hy; p2(3) = 0.0; - p3(1) = isign*0.5*hx - u*dt; - p3(2) = jsign*0.5*hy - vmac(i+ioff,j+1,k)*dt; + p3(1) = isign*Real(0.5)*hx - u*dt; + p3(2) = jsign*Real(0.5)*hy - vmac(i+ioff,j+1,k)*dt; p3(3) = 0.0; for(int n=1; n<=7; ++n){ @@ -746,52 +746,52 @@ BDS::ComputeConc (Box const& bx, vv = vmac(i+ioff,j+1,k+koff); } - p1(1) = isign*0.5*hx; - p1(2) = jsign*0.5*hy; - p1(3) = ksign*0.5*hz; + p1(1) = isign*Real(0.5)*hx; + p1(2) = jsign*Real(0.5)*hy; + p1(3) = ksign*Real(0.5)*hz; - p2(1) = isign*0.5*hx - umac(i,j,k)*dt; - p2(2) = jsign*0.5*hy; - p2(3) = ksign*0.5*hz; + p2(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p2(2) = jsign*Real(0.5)*hy; + p2(3) = ksign*Real(0.5)*hz; - p3(1) = isign*0.5*hx - umac(i,j,k)*dt; - p3(2) = jsign*0.5*hy - vmac(i+ioff,j+1,k)*dt; - p3(3) = ksign*0.5*hz; + p3(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p3(2) = jsign*Real(0.5)*hy - vmac(i+ioff,j+1,k)*dt; + p3(3) = ksign*Real(0.5)*hz; - p4(1) = isign*0.5*hx - uu*dt; - p4(2) = jsign*0.5*hy - vv*dt; - p4(3) = ksign*0.5*hz - wmac(i+ioff,j+joff,k+1)*dt; + p4(1) = isign*Real(0.5)*hx - uu*dt; + p4(2) = jsign*Real(0.5)*hy - vv*dt; + p4(3) = ksign*Real(0.5)*hz - wmac(i+ioff,j+joff,k+1)*dt; for(int n=1; n<=7; ++n){ slope_tmp(n) = slopes(i+ioff,j+joff,k+koff,n-1); } for(int ll=1; ll<=3; ++ll ){ - del(ll) = (p1(ll)+p2(ll)+p3(ll)+p4(ll))/4.0; + del(ll) = (p1(ll)+p2(ll)+p3(ll)+p4(ll))/Real(4); } val1 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p1(ll) + sixth*(p2(ll)+p3(ll)+p4(ll)); + del(ll) = Real(0.5)*p1(ll) + sixth*(p2(ll)+p3(ll)+p4(ll)); } val2 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p2(ll) + sixth*(p1(ll)+p3(ll)+p4(ll)); + del(ll) = Real(0.5)*p2(ll) + sixth*(p1(ll)+p3(ll)+p4(ll)); } val3 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p3(ll) + sixth*(p2(ll)+p1(ll)+p4(ll)); + del(ll) = Real(0.5)*p3(ll) + sixth*(p2(ll)+p1(ll)+p4(ll)); } val4 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p4(ll) + sixth*(p2(ll)+p3(ll)+p1(ll)); + del(ll) = Real(0.5)*p4(ll) + sixth*(p2(ll)+p3(ll)+p1(ll)); } val5 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); - gamma2 = -0.8*val1 + 0.45*(val2+val3+val4+val5); + gamma2 = -Real(0.8)*val1 + Real(0.45)*(val2+val3+val4+val5); // divu source term if (iconserv[icomp]) { @@ -800,76 +800,76 @@ BDS::ComputeConc (Box const& bx, gamma2 = gamma2 * wmac(i+ioff,j+joff,k+1); - gamma = gamma - dt*gamma2/(3.0*hz); + gamma = gamma - dt*gamma2/(Real(3)*hz); //////////////////////////////////////////////// // correct \Gamma^{y+} with \Gamma^{y+,z-} //////////////////////////////////////////////// - if (wmac(i+ioff,j+joff,k) > 0.0) { - ksign = 1.0; + if (wmac(i+ioff,j+joff,k) > Real(0)) { + ksign = Real(1); koff = -1; } else { - ksign = -1.0; + ksign = -Real(1); koff = 0; } - uu = 0.0; - if (umac(i,j,k)*umac(i,j+joff,k+koff) > 0.0) { + uu = Real(0); + if (umac(i,j,k)*umac(i,j+joff,k+koff) > Real(0)) { uu = umac(i,j+joff,k+koff); } - vv = 0.0; - if (vmac(i+ioff,j+1,k)*vmac(i+ioff,j+1,k+koff) > 0.0) { + vv = Real(0); + if (vmac(i+ioff,j+1,k)*vmac(i+ioff,j+1,k+koff) > Real(0)) { vv = vmac(i+ioff,j+1,k+koff); } - p1(1) = isign*0.5*hx; - p1(2) = jsign*0.5*hy; - p1(3) = ksign*0.5*hz; + p1(1) = isign*Real(0.5)*hx; + p1(2) = jsign*Real(0.5)*hy; + p1(3) = ksign*Real(0.5)*hz; - p2(1) = isign*0.5*hx - umac(i,j,k)*dt; - p2(2) = jsign*0.5*hy; - p2(3) = ksign*0.5*hz; + p2(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p2(2) = jsign*Real(0.5)*hy; + p2(3) = ksign*Real(0.5)*hz; - p3(1) = isign*0.5*hx - umac(i,j,k)*dt; - p3(2) = jsign*0.5*hy - vmac(i+ioff,j+1,k)*dt; - p3(3) = ksign*0.5*hz; + p3(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p3(2) = jsign*Real(0.5)*hy - vmac(i+ioff,j+1,k)*dt; + p3(3) = ksign*Real(0.5)*hz; - p4(1) = isign*0.5*hx - uu*dt; - p4(2) = jsign*0.5*hy - vv*dt; - p4(3) = ksign*0.5*hz - wmac(i+ioff,j+joff,k)*dt; + p4(1) = isign*Real(0.5)*hx - uu*dt; + p4(2) = jsign*Real(0.5)*hy - vv*dt; + p4(3) = ksign*Real(0.5)*hz - wmac(i+ioff,j+joff,k)*dt; for(int n=1; n<=7; ++n){ slope_tmp(n) = slopes(i+ioff,j+joff,k+koff,n-1); } for(int ll=1; ll<=3; ++ll ){ - del(ll) = (p1(ll)+p2(ll)+p3(ll)+p4(ll))/4.0; + del(ll) = (p1(ll)+p2(ll)+p3(ll)+p4(ll))/Real(4); } val1 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p1(ll) + sixth*(p2(ll)+p3(ll)+p4(ll)); + del(ll) = Real(0.5)*p1(ll) + sixth*(p2(ll)+p3(ll)+p4(ll)); } val2 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p2(ll) + sixth*(p1(ll)+p3(ll)+p4(ll)); + del(ll) = Real(0.5)*p2(ll) + sixth*(p1(ll)+p3(ll)+p4(ll)); } val3 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p3(ll) + sixth*(p2(ll)+p1(ll)+p4(ll)); + del(ll) = Real(0.5)*p3(ll) + sixth*(p2(ll)+p1(ll)+p4(ll)); } val4 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p4(ll) + sixth*(p2(ll)+p3(ll)+p1(ll)); + del(ll) = Real(0.5)*p4(ll) + sixth*(p2(ll)+p3(ll)+p1(ll)); } val5 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); - gamma2 = -0.8*val1 + 0.45*(val2+val3+val4+val5); + gamma2 = -Real(0.8)*val1 + Real(0.45)*(val2+val3+val4+val5); // divu source term if (iconserv[icomp]) { @@ -878,65 +878,65 @@ BDS::ComputeConc (Box const& bx, gamma2 = gamma2 * wmac(i+ioff,j+joff,k); - gamma = gamma + dt*gamma2/(3.0*hz); + gamma = gamma + dt*gamma2/(Real(3)*hz); //////////////////////////////////////////////// // correct sedgex with \Gamma^{y+} //////////////////////////////////////////////// gamma = gamma * vmac(i+ioff,j+1,k); - xedge_tmp = xedge_tmp - dt*gamma/(2.0*hy); + xedge_tmp = xedge_tmp - dt*gamma/(Real(2)*hy); //////////////////////////////////////////////// // compute \Gamma^{y-} without corner corrections //////////////////////////////////////////////// - if (vmac(i+ioff,j,k) > 0.0) { - jsign = 1.0; + if (vmac(i+ioff,j,k) > Real(0)) { + jsign = Real(1); joff = -1; } else { - jsign = -1.0; + jsign = -Real(1); joff = 0; } - u = 0.0; - if (umac(i,j,k)*umac(i,j+joff,k) > 0.0) { + u = Real(0); + if (umac(i,j,k)*umac(i,j+joff,k) > Real(0)) { u = umac(i,j+joff,k); } - p1(1) = isign*0.5*hx; - p1(2) = jsign*0.5*hy; - p1(3) = 0.0; + p1(1) = isign*Real(0.5)*hx; + p1(2) = jsign*Real(0.5)*hy; + p1(3) = Real(0); - p2(1) = isign*0.5*hx - umac(i,j,k)*dt; - p2(2) = jsign*0.5*hy; - p2(3) = 0.0; + p2(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p2(2) = jsign*Real(0.5)*hy; + p2(3) = Real(0); - p3(1) = isign*0.5*hx - u*dt; - p3(2) = jsign*0.5*hy - vmac(i+ioff,j,k)*dt; - p3(3) = 0.0; + p3(1) = isign*Real(0.5)*hx - u*dt; + p3(2) = jsign*Real(0.5)*hy - vmac(i+ioff,j,k)*dt; + p3(3) = Real(0); for(int n=1; n<=7; ++n){ slope_tmp(n) = slopes(i+ioff,j+joff,k,n-1); } for(int ll=1; ll<=3; ++ll ){ - del(ll) = (p2(ll)+p3(ll))/2.0; + del(ll) = (p2(ll)+p3(ll))*Real(0.5); } val1 = eval(s(i+ioff,j+joff,k,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = (p1(ll)+p3(ll))/2.0; + del(ll) = (p1(ll)+p3(ll))*Real(0.5); } val2 = eval(s(i+ioff,j+joff,k,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = (p1(ll)+p2(ll))/2.0; + del(ll) = (p1(ll)+p2(ll))*Real(0.5); } val3 = eval(s(i+ioff,j+joff,k,icomp),slope_tmp,del); // average these centroid values to get the average value - gamma = (val1+val2+val3)/3.0; + gamma = (val1+val2+val3)/Real(3); // source term if (iconserv[icomp]) { @@ -949,70 +949,70 @@ BDS::ComputeConc (Box const& bx, // correct \Gamma^{y-} with \Gamma^{y-,z+} //////////////////////////////////////////////// - if (wmac(i+ioff,j+joff,k+1) > 0.0) { - ksign = 1.0; + if (wmac(i+ioff,j+joff,k+1) > Real(0)) { + ksign = Real(1); koff = 0; } else { - ksign = -1.0; + ksign = -Real(1); koff = 1; } - uu = 0.0; - if (umac(i,j,k)*umac(i,j+joff,k+koff) > 0.0) { + uu = Real(0); + if (umac(i,j,k)*umac(i,j+joff,k+koff) > Real(0)) { uu = umac(i,j+joff,k+koff); } - vv = 0.0; - if (vmac(i+ioff,j,k)*vmac(i+ioff,j,k+koff) > 0.0) { + vv = Real(0); + if (vmac(i+ioff,j,k)*vmac(i+ioff,j,k+koff) > Real(0)) { vv = vmac(i+ioff,j,k+koff); } - p1(1) = isign*0.5*hx; - p1(2) = jsign*0.5*hy; - p1(3) = ksign*0.5*hz; + p1(1) = isign*Real(0.5)*hx; + p1(2) = jsign*Real(0.5)*hy; + p1(3) = ksign*Real(0.5)*hz; - p2(1) = isign*0.5*hx - umac(i,j,k)*dt; - p2(2) = jsign*0.5*hy; - p2(3) = ksign*0.5*hz; + p2(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p2(2) = jsign*Real(0.5)*hy; + p2(3) = ksign*Real(0.5)*hz; - p3(1) = isign*0.5*hx - umac(i,j,k)*dt; - p3(2) = jsign*0.5*hy - vmac(i+ioff,j,k)*dt; - p3(3) = ksign*0.5*hz; + p3(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p3(2) = jsign*Real(0.5)*hy - vmac(i+ioff,j,k)*dt; + p3(3) = ksign*Real(0.5)*hz; - p4(1) = isign*0.5*hx - uu*dt; - p4(2) = jsign*0.5*hy - vv*dt; - p4(3) = ksign*0.5*hz - wmac(i+ioff,j+joff,k+1)*dt; + p4(1) = isign*Real(0.5)*hx - uu*dt; + p4(2) = jsign*Real(0.5)*hy - vv*dt; + p4(3) = ksign*Real(0.5)*hz - wmac(i+ioff,j+joff,k+1)*dt; for(int n=1; n<=7; ++n){ slope_tmp(n) = slopes(i+ioff,j+joff,k+koff,n-1); } for(int ll=1; ll<=3; ++ll ){ - del(ll) = (p1(ll)+p2(ll)+p3(ll)+p4(ll))/4.0; + del(ll) = (p1(ll)+p2(ll)+p3(ll)+p4(ll))/Real(4); } val1 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p1(ll) + sixth*(p2(ll)+p3(ll)+p4(ll)); + del(ll) = Real(0.5)*p1(ll) + sixth*(p2(ll)+p3(ll)+p4(ll)); } val2 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p2(ll) + sixth*(p1(ll)+p3(ll)+p4(ll)); + del(ll) = Real(0.5)*p2(ll) + sixth*(p1(ll)+p3(ll)+p4(ll)); } val3 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p3(ll) + sixth*(p2(ll)+p1(ll)+p4(ll)); + del(ll) = Real(0.5)*p3(ll) + sixth*(p2(ll)+p1(ll)+p4(ll)); } val4 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p4(ll) + sixth*(p2(ll)+p3(ll)+p1(ll)); + del(ll) = Real(0.5)*p4(ll) + sixth*(p2(ll)+p3(ll)+p1(ll)); } val5 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); - gamma2 = -0.8*val1 + 0.45*(val2+val3+val4+val5); + gamma2 = -Real(0.8)*val1 + Real(0.45)*(val2+val3+val4+val5); // divu source term if (iconserv[icomp]) { @@ -1021,76 +1021,76 @@ BDS::ComputeConc (Box const& bx, gamma2 = gamma2 * wmac(i+ioff,j+joff,k+1); - gamma = gamma - dt*gamma2/(3.0*hz); + gamma = gamma - dt*gamma2/(Real(3)*hz); //////////////////////////////////////////////// // correct \Gamma^{y-} with \Gamma^{y-,z-} //////////////////////////////////////////////// - if (wmac(i+ioff,j+joff,k) > 0.0) { - ksign = 1.0; + if (wmac(i+ioff,j+joff,k) > Real(0)) { + ksign = Real(1); koff = -1; } else { - ksign = -1.0; + ksign = -Real(1); koff = 0; } - uu = 0.0; - if (umac(i,j,k)*umac(i,j+joff,k+koff) > 0.0) { + uu = Real(0); + if (umac(i,j,k)*umac(i,j+joff,k+koff) > Real(0)) { uu = umac(i,j+joff,k+koff); } - vv = 0.0; - if (vmac(i+ioff,j,k)*vmac(i+ioff,j,k+koff) > 0.0) { + vv = Real(0); + if (vmac(i+ioff,j,k)*vmac(i+ioff,j,k+koff) > Real(0)) { vv = vmac(i+ioff,j,k+koff); } - p1(1) = isign*0.5*hx; - p1(2) = jsign*0.5*hy; - p1(3) = ksign*0.5*hz; + p1(1) = isign*Real(0.5)*hx; + p1(2) = jsign*Real(0.5)*hy; + p1(3) = ksign*Real(0.5)*hz; - p2(1) = isign*0.5*hx - umac(i,j,k)*dt; - p2(2) = jsign*0.5*hy; - p2(3) = ksign*0.5*hz; + p2(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p2(2) = jsign*Real(0.5)*hy; + p2(3) = ksign*Real(0.5)*hz; - p3(1) = isign*0.5*hx - umac(i,j,k)*dt; - p3(2) = jsign*0.5*hy - vmac(i+ioff,j,k)*dt; - p3(3) = ksign*0.5*hz; + p3(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p3(2) = jsign*Real(0.5)*hy - vmac(i+ioff,j,k)*dt; + p3(3) = ksign*Real(0.5)*hz; - p4(1) = isign*0.5*hx - uu*dt; - p4(2) = jsign*0.5*hy - vv*dt; - p4(3) = ksign*0.5*hz - wmac(i+ioff,j+joff,k)*dt; + p4(1) = isign*Real(0.5)*hx - uu*dt; + p4(2) = jsign*Real(0.5)*hy - vv*dt; + p4(3) = ksign*Real(0.5)*hz - wmac(i+ioff,j+joff,k)*dt; for(int n=1; n<=7; ++n){ slope_tmp(n) = slopes(i+ioff,j+joff,k+koff,n-1); } for(int ll=1; ll<=3; ++ll ){ - del(ll) = (p1(ll)+p2(ll)+p3(ll)+p4(ll))/4.0; + del(ll) = (p1(ll)+p2(ll)+p3(ll)+p4(ll))*Real(0.25); } val1 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p1(ll) + sixth*(p2(ll)+p3(ll)+p4(ll)); + del(ll) = Real(0.5)*p1(ll) + sixth*(p2(ll)+p3(ll)+p4(ll)); } val2 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p2(ll) + sixth*(p1(ll)+p3(ll)+p4(ll)); + del(ll) = Real(0.5)*p2(ll) + sixth*(p1(ll)+p3(ll)+p4(ll)); } val3 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p3(ll) + sixth*(p2(ll)+p1(ll)+p4(ll)); + del(ll) = Real(0.5)*p3(ll) + sixth*(p2(ll)+p1(ll)+p4(ll)); } val4 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p4(ll) + sixth*(p2(ll)+p3(ll)+p1(ll)); + del(ll) = Real(0.5)*p4(ll) + sixth*(p2(ll)+p3(ll)+p1(ll)); } val5 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); - gamma2 = -0.8*val1 + 0.45*(val2+val3+val4+val5); + gamma2 = -Real(0.8)*val1 + Real(0.45)*(val2+val3+val4+val5); // divu source term if (iconserv[icomp]) { @@ -1099,141 +1099,141 @@ BDS::ComputeConc (Box const& bx, gamma2 = gamma2 * wmac(i+ioff,j+joff,k); - gamma = gamma + dt*gamma2/(3.0*hz); + gamma = gamma + dt*gamma2/(Real(3)*hz); //////////////////////////////////////////////// // correct sedgex with \Gamma^{y-} //////////////////////////////////////////////// gamma = gamma * vmac(i+ioff,j,k); - xedge_tmp = xedge_tmp + dt*gamma/(2.0*hy); + xedge_tmp = xedge_tmp + dt*gamma/(Real(2)*hy); //////////////////////////////////////////////// // compute \Gamma^{z+} without corner corrections //////////////////////////////////////////////// - if (wmac(i+ioff,j,k+1) > 0.0) { - ksign = 1.0; + if (wmac(i+ioff,j,k+1) > Real(0)) { + ksign = Real(1); koff = 0; } else { - ksign = -1.0; + ksign = -Real(1); koff = 1; } - u = 0.0; - if (umac(i,j,k)*umac(i,j,k+koff) > 0.0) { + u = Real(0); + if (umac(i,j,k)*umac(i,j,k+koff) > Real(0)) { u = umac(i,j,k+koff); } - p1(1) = isign*0.5*hx; - p1(2) = 0.0; - p1(3) = ksign*0.5*hz; + p1(1) = isign*Real(0.5)*hx; + p1(2) = Real(0); + p1(3) = ksign*Real(0.5)*hz; - p2(1) = isign*0.5*hx - umac(i,j,k)*dt; - p2(2) = 0.0; - p2(3) = ksign*0.5*hz; + p2(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p2(2) = Real(0); + p2(3) = ksign*Real(0.5)*hz; - p3(1) = isign*0.5*hx - u*dt; - p3(2) = 0.0; - p3(3) = ksign*0.5*hz - wmac(i+ioff,j,k+1)*dt; + p3(1) = isign*Real(0.5)*hx - u*dt; + p3(2) = Real(0); + p3(3) = ksign*Real(0.5)*hz - wmac(i+ioff,j,k+1)*dt; for(int n=1; n<=7; ++n){ slope_tmp(n) = slopes(i+ioff,j,k+koff,n-1); } for(int ll=1; ll<=3; ++ll ){ - del(ll) = (p2(ll)+p3(ll))/2.0; + del(ll) = (p2(ll)+p3(ll))*Real(0.5); } val1 = eval(s(i+ioff,j,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = (p1(ll)+p3(ll))/2.0; + del(ll) = (p1(ll)+p3(ll))*Real(0.5); } val2 = eval(s(i+ioff,j,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = (p1(ll)+p2(ll))/2.0; + del(ll) = (p1(ll)+p2(ll))*Real(0.5); } val3 = eval(s(i+ioff,j,k+koff,icomp),slope_tmp,del); // average these centroid values to get the average value - gamma = (val1+val2+val3)/3.0; + gamma = (val1+val2+val3)/Real(3); // source term if (iconserv[icomp]) { - gamma = gamma*(1. - dt3*(ux(i+ioff,j,k+koff)+wz(i+ioff,j,k+koff))); + gamma = gamma*(Real(1) - dt3*(ux(i+ioff,j,k+koff)+wz(i+ioff,j,k+koff))); } else { - gamma = gamma*(1. + dt3*vy(i+ioff,j,k+koff)); + gamma = gamma*(Real(1) + dt3*vy(i+ioff,j,k+koff)); } //////////////////////////////////////////////// // correct \Gamma^{z+} with \Gamma^{z+,y+} //////////////////////////////////////////////// - if (vmac(i+ioff,j+1,k+koff) > 0.0) { - jsign = 1.0; + if (vmac(i+ioff,j+1,k+koff) > Real(0)) { + jsign = Real(1); joff = 0; } else { - jsign = -1.0; + jsign = -Real(1); joff = 1; } - uu = 0.0; - if (umac(i,j,k)*umac(i,j+joff,k+koff) > 0.0) { + uu = Real(0); + if (umac(i,j,k)*umac(i,j+joff,k+koff) > Real(0)) { uu = umac(i,j+joff,k+koff); } - ww = 0.0; + ww = Real(0); if (wmac(i+ioff,j,k+1)*wmac(i+ioff,j+joff,k+1) > 0.0) { ww = wmac(i+ioff,j+joff,k+1); } - p1(1) = isign*0.5*hx; - p1(2) = jsign*0.5*hy; - p1(3) = ksign*0.5*hz; + p1(1) = isign*Real(0.5)*hx; + p1(2) = jsign*Real(0.5)*hy; + p1(3) = ksign*Real(0.5)*hz; - p2(1) = isign*0.5*hx - umac(i,j,k)*dt; - p2(2) = jsign*0.5*hy; - p2(3) = ksign*0.5*hz; + p2(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p2(2) = jsign*Real(0.5)*hy; + p2(3) = ksign*Real(0.5)*hz; - p3(1) = isign*0.5*hx - umac(i,j,k)*dt; - p3(2) = jsign*0.5*hy; - p3(3) = ksign*0.5*hz - wmac(i+ioff,j,k+1)*dt; + p3(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p3(2) = jsign*Real(0.5)*hy; + p3(3) = ksign*Real(0.5)*hz - wmac(i+ioff,j,k+1)*dt; - p4(1) = isign*0.5*hx - uu*dt; - p4(2) = jsign*0.5*hy - vmac(i+ioff,j+1,k+koff)*dt; - p4(3) = ksign*0.5*hz - ww*dt; + p4(1) = isign*Real(0.5)*hx - uu*dt; + p4(2) = jsign*Real(0.5)*hy - vmac(i+ioff,j+1,k+koff)*dt; + p4(3) = ksign*Real(0.5)*hz - ww*dt; for(int n=1; n<=7; ++n){ slope_tmp(n) = slopes(i+ioff,j+joff,k+koff,n-1); } for(int ll=1; ll<=3; ++ll ){ - del(ll) = (p1(ll)+p2(ll)+p3(ll)+p4(ll))/4.0; + del(ll) = (p1(ll)+p2(ll)+p3(ll)+p4(ll))*Real(0.25); } val1 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p1(ll) + sixth*(p2(ll)+p3(ll)+p4(ll)); + del(ll) = Real(0.5)*p1(ll) + sixth*(p2(ll)+p3(ll)+p4(ll)); } val2 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p2(ll) + sixth*(p1(ll)+p3(ll)+p4(ll)); + del(ll) = Real(0.5)*p2(ll) + sixth*(p1(ll)+p3(ll)+p4(ll)); } val3 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p3(ll) + sixth*(p2(ll)+p1(ll)+p4(ll)); + del(ll) = Real(0.5)*p3(ll) + sixth*(p2(ll)+p1(ll)+p4(ll)); } val4 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p4(ll) + sixth*(p2(ll)+p3(ll)+p1(ll)); + del(ll) = Real(0.5)*p4(ll) + sixth*(p2(ll)+p3(ll)+p1(ll)); } val5 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); - gamma2 = -0.8*val1 + 0.45*(val2+val3+val4+val5); + gamma2 = -Real(0.8)*val1 + Real(0.45)*(val2+val3+val4+val5); // divu source term if (iconserv[icomp]) { @@ -1242,76 +1242,76 @@ BDS::ComputeConc (Box const& bx, gamma2 = gamma2 * vmac(i+ioff,j+1,k+koff); - gamma = gamma - dt*gamma2/(3.0*hy); + gamma = gamma - dt*gamma2/(Real(3)*hy); //////////////////////////////////////////////// // correct \Gamma^{z+} with \Gamma^{z+,y-} //////////////////////////////////////////////// - if (vmac(i+ioff,j,k+koff) > 0.0) { - jsign = 1.0; + if (vmac(i+ioff,j,k+koff) > Real(0)) { + jsign = Real(1); joff = -1; } else { - jsign = -1.0; + jsign = -Real(1); joff = 0; } - uu = 0.0; - if (umac(i,j,k)*umac(i,j+joff,k+koff) > 0.0) { + uu = Real(0); + if (umac(i,j,k)*umac(i,j+joff,k+koff) > Real(0)) { uu = umac(i,j+joff,k+koff); } - ww = 0.0; - if (wmac(i+ioff,j,k+1)*wmac(i+ioff,j+joff,k+1) > 0.0) { + ww = Real(0); + if (wmac(i+ioff,j,k+1)*wmac(i+ioff,j+joff,k+1) > Real(0)) { ww = wmac(i+ioff,j+joff,k+1); } - p1(1) = isign*0.5*hx; - p1(2) = jsign*0.5*hy; - p1(3) = ksign*0.5*hz; + p1(1) = isign*Real(0.5)*hx; + p1(2) = jsign*Real(0.5)*hy; + p1(3) = ksign*Real(0.5)*hz; - p2(1) = isign*0.5*hx - umac(i,j,k)*dt; - p2(2) = jsign*0.5*hy; - p2(3) = ksign*0.5*hz; + p2(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p2(2) = jsign*Real(0.5)*hy; + p2(3) = ksign*Real(0.5)*hz; - p3(1) = isign*0.5*hx - umac(i,j,k)*dt; - p3(2) = jsign*0.5*hy; - p3(3) = ksign*0.5*hz - wmac(i+ioff,j,k+1)*dt; + p3(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p3(2) = jsign*Real(0.5)*hy; + p3(3) = ksign*Real(0.5)*hz - wmac(i+ioff,j,k+1)*dt; - p4(1) = isign*0.5*hx - uu*dt; - p4(2) = jsign*0.5*hy - vmac(i+ioff,j,k+koff)*dt; - p4(3) = ksign*0.5*hz - ww*dt; + p4(1) = isign*Real(0.5)*hx - uu*dt; + p4(2) = jsign*Real(0.5)*hy - vmac(i+ioff,j,k+koff)*dt; + p4(3) = ksign*Real(0.5)*hz - ww*dt; for(int n=1; n<=7; ++n){ slope_tmp(n) = slopes(i+ioff,j+joff,k+koff,n-1); } for(int ll=1; ll<=3; ++ll ){ - del(ll) = (p1(ll)+p2(ll)+p3(ll)+p4(ll))/4.0; + del(ll) = (p1(ll)+p2(ll)+p3(ll)+p4(ll))*Real(0.25); } val1 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p1(ll) + sixth*(p2(ll)+p3(ll)+p4(ll)); + del(ll) = Real(0.5)*p1(ll) + sixth*(p2(ll)+p3(ll)+p4(ll)); } val2 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p2(ll) + sixth*(p1(ll)+p3(ll)+p4(ll)); + del(ll) = Real(0.5)*p2(ll) + sixth*(p1(ll)+p3(ll)+p4(ll)); } val3 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p3(ll) + sixth*(p2(ll)+p1(ll)+p4(ll)); + del(ll) = Real(0.5)*p3(ll) + sixth*(p2(ll)+p1(ll)+p4(ll)); } val4 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p4(ll) + sixth*(p2(ll)+p3(ll)+p1(ll)); + del(ll) = Real(0.5)*p4(ll) + sixth*(p2(ll)+p3(ll)+p1(ll)); } val5 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); - gamma2 = -0.8*val1 + 0.45*(val2+val3+val4+val5); + gamma2 = -Real(0.8)*val1 + Real(0.45)*(val2+val3+val4+val5); // divu source term if (iconserv[icomp]) { @@ -1320,43 +1320,43 @@ BDS::ComputeConc (Box const& bx, gamma2 = gamma2 * vmac(i+ioff,j,k+koff); - gamma = gamma + dt*gamma2/(3.0*hy); + gamma = gamma + dt*gamma2/(Real(3)*hy); //////////////////////////////////////////////// // correct sedgex with \Gamma^{z+} //////////////////////////////////////////////// gamma = gamma * wmac(i+ioff,j,k+1); - xedge_tmp = xedge_tmp - dt*gamma/(2.0*hz); + xedge_tmp = xedge_tmp - dt*gamma/(Real(2)*hz); //////////////////////////////////////////////// // compute \Gamma^{z-} without corner corrections //////////////////////////////////////////////// - if (wmac(i+ioff,j,k) > 0.0) { - ksign = 1.0; + if (wmac(i+ioff,j,k) > Real(0)) { + ksign = Real(1); koff = -1; } else { - ksign = -1.0; + ksign = -Real(1); koff = 0; } - u = 0.0; - if (umac(i,j,k)*umac(i,j,k+koff) > 0.0) { + u = Real(0); + if (umac(i,j,k)*umac(i,j,k+koff) > Real(0)) { u = umac(i,j,k+koff); } - p1(1) = isign*0.5*hx; - p1(2) = 0.0; - p1(3) = ksign*0.5*hz; + p1(1) = isign*Real(0.5)*hx; + p1(2) = Real(0); + p1(3) = ksign*Real(0.5)*hz; - p2(1) = isign*0.5*hx - umac(i,j,k)*dt; - p2(2) = 0.0; - p2(3) = ksign*0.5*hz; + p2(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p2(2) = Real(0); + p2(3) = ksign*Real(0.5)*hz; - p3(1) = isign*0.5*hx - u*dt; - p3(2) = 0.0; - p3(3) = ksign*0.5*hz - wmac(i+ioff,j,k)*dt; + p3(1) = isign*Real(0.5)*hx - u*dt; + p3(2) = Real(0); + p3(3) = ksign*Real(0.5)*hz - wmac(i+ioff,j,k)*dt; for(int n=1; n<=7; ++n){ slope_tmp(n) = slopes(i+ioff,j,k+koff,n-1); @@ -1391,7 +1391,7 @@ BDS::ComputeConc (Box const& bx, // correct \Gamma^{z-} with \Gamma^{z-,y+} //////////////////////////////////////////////// - if (vmac(i+ioff,j+1,k+koff) > 0.0) { + if (vmac(i+ioff,j+1,k+koff) > Real(0)) { jsign = 1.0; joff = 0; } else { @@ -1399,31 +1399,31 @@ BDS::ComputeConc (Box const& bx, joff = 1; } - uu = 0.0; - if (umac(i,j,k)*umac(i,j+joff,k+koff) > 0.0) { + uu = Real(0); + if (umac(i,j,k)*umac(i,j+joff,k+koff) > Real(0)) { uu = umac(i,j+joff,k+koff); } - ww = 0.0; - if (wmac(i+ioff,j,k)*wmac(i+ioff,j+joff,k) > 0.0) { + ww = Real(0); + if (wmac(i+ioff,j,k)*wmac(i+ioff,j+joff,k) > Real(0)) { ww = wmac(i+ioff,j+joff,k); } - p1(1) = isign*0.5*hx; - p1(2) = jsign*0.5*hy; - p1(3) = ksign*0.5*hz; + p1(1) = isign*Real(0.5)*hx; + p1(2) = jsign*Real(0.5)*hy; + p1(3) = ksign*Real(0.5)*hz; - p2(1) = isign*0.5*hx - umac(i,j,k)*dt; - p2(2) = jsign*0.5*hy; - p2(3) = ksign*0.5*hz; + p2(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p2(2) = jsign*Real(0.5)*hy; + p2(3) = ksign*Real(0.5)*hz; - p3(1) = isign*0.5*hx - umac(i,j,k)*dt; - p3(2) = jsign*0.5*hy; - p3(3) = ksign*0.5*hz - wmac(i+ioff,j,k)*dt; + p3(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p3(2) = jsign*Real(0.5)*hy; + p3(3) = ksign*Real(0.5)*hz - wmac(i+ioff,j,k)*dt; - p4(1) = isign*0.5*hx - uu*dt; - p4(2) = jsign*0.5*hy - vmac(i+ioff,j+1,k+koff)*dt; - p4(3) = ksign*0.5*hz - ww*dt; + p4(1) = isign*Real(0.5)*hx - uu*dt; + p4(2) = jsign*Real(0.5)*hy - vmac(i+ioff,j+1,k+koff)*dt; + p4(3) = ksign*Real(0.5)*hz - ww*dt; for(int n=1; n<=7; ++n){ slope_tmp(n) = slopes(i+ioff,j+joff,k+koff,n-1); @@ -1435,26 +1435,26 @@ BDS::ComputeConc (Box const& bx, val1 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p1(ll) + sixth*(p2(ll)+p3(ll)+p4(ll)); + del(ll) = Real(0.5)*p1(ll) + sixth*(p2(ll)+p3(ll)+p4(ll)); } val2 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p2(ll) + sixth*(p1(ll)+p3(ll)+p4(ll)); + del(ll) = Real(0.5)*p2(ll) + sixth*(p1(ll)+p3(ll)+p4(ll)); } val3 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p3(ll) + sixth*(p2(ll)+p1(ll)+p4(ll)); + del(ll) = Real(0.5)*p3(ll) + sixth*(p2(ll)+p1(ll)+p4(ll)); } val4 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p4(ll) + sixth*(p2(ll)+p3(ll)+p1(ll)); + del(ll) = Real(0.5)*p4(ll) + sixth*(p2(ll)+p3(ll)+p1(ll)); } val5 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); - gamma2 = -0.8*val1 + 0.45*(val2+val3+val4+val5); + gamma2 = -Real(0.8)*val1 + Real(0.45)*(val2+val3+val4+val5); // divu source term if (iconserv[icomp]) { @@ -1469,7 +1469,7 @@ BDS::ComputeConc (Box const& bx, // correct \Gamma^{z-} with \Gamma^{z-,y-} //////////////////////////////////////////////// - if (vmac(i+ioff,j,k+koff) > 0.0) { + if (vmac(i+ioff,j,k+koff) > Real(0)) { jsign = 1.0; joff = -1; } else { @@ -1477,31 +1477,31 @@ BDS::ComputeConc (Box const& bx, joff = 0; } - uu = 0.0; - if (umac(i,j,k)*umac(i,j+joff,k+koff) > 0.0) { + uu = Real(0); + if (umac(i,j,k)*umac(i,j+joff,k+koff) > Real(0)) { uu = umac(i,j+joff,k+koff); } - ww = 0.0; - if (wmac(i+ioff,j,k)*wmac(i+ioff,j+joff,k) > 0.0) { + ww = Real(0); + if (wmac(i+ioff,j,k)*wmac(i+ioff,j+joff,k) > Real(0)) { ww = wmac(i+ioff,j+joff,k); } - p1(1) = isign*0.5*hx; - p1(2) = jsign*0.5*hy; - p1(3) = ksign*0.5*hz; + p1(1) = isign*Real(0.5)*hx; + p1(2) = jsign*Real(0.5)*hy; + p1(3) = ksign*Real(0.5)*hz; - p2(1) = isign*0.5*hx - umac(i,j,k)*dt; - p2(2) = jsign*0.5*hy; - p2(3) = ksign*0.5*hz; + p2(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p2(2) = jsign*Real(0.5)*hy; + p2(3) = ksign*Real(0.5)*hz; - p3(1) = isign*0.5*hx - umac(i,j,k)*dt; - p3(2) = jsign*0.5*hy; - p3(3) = ksign*0.5*hz - wmac(i+ioff,j,k)*dt; + p3(1) = isign*Real(0.5)*hx - umac(i,j,k)*dt; + p3(2) = jsign*Real(0.5)*hy; + p3(3) = ksign*Real(0.5)*hz - wmac(i+ioff,j,k)*dt; - p4(1) = isign*0.5*hx - uu*dt; - p4(2) = jsign*0.5*hy - vmac(i+ioff,j,k+koff)*dt; - p4(3) = ksign*0.5*hz - ww*dt; + p4(1) = isign*Real(0.5)*hx - uu*dt; + p4(2) = jsign*Real(0.5)*hy - vmac(i+ioff,j,k+koff)*dt; + p4(3) = ksign*Real(0.5)*hz - ww*dt; for(int n=1; n<=7; ++n){ slope_tmp(n) = slopes(i+ioff,j+joff,k+koff,n-1); @@ -1513,26 +1513,26 @@ BDS::ComputeConc (Box const& bx, val1 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p1(ll) + sixth*(p2(ll)+p3(ll)+p4(ll)); + del(ll) = Real(0.5)*p1(ll) + sixth*(p2(ll)+p3(ll)+p4(ll)); } val2 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p2(ll) + sixth*(p1(ll)+p3(ll)+p4(ll)); + del(ll) = Real(0.5)*p2(ll) + sixth*(p1(ll)+p3(ll)+p4(ll)); } val3 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p3(ll) + sixth*(p2(ll)+p1(ll)+p4(ll)); + del(ll) = Real(0.5)*p3(ll) + sixth*(p2(ll)+p1(ll)+p4(ll)); } val4 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p4(ll) + sixth*(p2(ll)+p3(ll)+p1(ll)); + del(ll) = Real(0.5)*p4(ll) + sixth*(p2(ll)+p3(ll)+p1(ll)); } val5 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); - gamma2 = -0.8*val1 + 0.45*(val2+val3+val4+val5); + gamma2 = -Real(0.8)*val1 + Real(0.45)*(val2+val3+val4+val5); // divu source term if (iconserv[icomp]) { @@ -1562,7 +1562,7 @@ BDS::ComputeConc (Box const& bx, sedgey(i,j,k,icomp) = s(i,j-1,k,icomp); if (is_velocity && icomp == YVEL && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) { // make sure velocity is not blowing inward - sedgey(i,j,k,icomp) = amrex::min(0._rt,sedgey(i,j,k,icomp)); + sedgey(i,j,k,icomp) = amrex::min(Real(0),sedgey(i,j,k,icomp)); } return; } @@ -1570,7 +1570,7 @@ BDS::ComputeConc (Box const& bx, sedgey(i,j,k,icomp) = s(i,j,k,icomp); if (is_velocity && icomp == YVEL && (bc.hi(1) == BCType::foextrap || bc.hi(1) == BCType::hoextrap) ) { // make sure velocity is not blowing inward - sedgey(i,j,k,icomp) = amrex::max(0._rt,sedgey(i,j,k,icomp)); + sedgey(i,j,k,icomp) = amrex::max(Real(0),sedgey(i,j,k,icomp)); } return; } @@ -1600,7 +1600,7 @@ BDS::ComputeConc (Box const& bx, //////////////////////////////////////////////// // centroid of rectangular volume - if (vmac(i,j,k) > 0.0) { + if (vmac(i,j,k) > Real(0)) { jsign = 1.0; joff = -1; } else { @@ -1608,9 +1608,9 @@ BDS::ComputeConc (Box const& bx, joff = 0; } - del(1) = 0.0; - del(2) = jsign*0.5*hy - 0.5*vmac(i,j,k)*dt; - del(3) = 0.0; + del(1) = Real(0); + del(2) = jsign*Real(0.5)*hy - Real(0.5)*vmac(i,j,k)*dt; + del(3) = Real(0); for(int n=1; n<=7; ++n){ slope_tmp(n) = slopes(i,j+joff,k,n-1); @@ -1632,7 +1632,7 @@ BDS::ComputeConc (Box const& bx, // compute \Gamma^{x+} without corner corrections //////////////////////////////////////////////// - if (umac(i+1,j+joff,k) > 0.0) { + if (umac(i+1,j+joff,k) > Real(0)) { isign = 1.0; ioff = 0; } else { @@ -1640,22 +1640,22 @@ BDS::ComputeConc (Box const& bx, ioff = 1; } - v = 0.0; - if (vmac(i,j,k)*vmac(i+ioff,j,k) > 0.0) { + v = Real(0); + if (vmac(i,j,k)*vmac(i+ioff,j,k) > Real(0)) { v = vmac(i+ioff,j,k); } - p1(1) = isign*0.5*hx; - p1(2) = jsign*0.5*hy; - p1(3) = 0.0; + p1(1) = isign*Real(0.5)*hx; + p1(2) = jsign*Real(0.5)*hy; + p1(3) = Real(0); - p2(1) = isign*0.5*hx; - p2(2) = jsign*0.5*hy - vmac(i,j,k)*dt; - p2(3) = 0.0; + p2(1) = isign*Real(0.5)*hx; + p2(2) = jsign*Real(0.5)*hy - vmac(i,j,k)*dt; + p2(3) = Real(0); - p3(1) = isign*0.5*hx - umac(i+1,j+joff,k)*dt; - p3(2) = jsign*0.5*hy - v*dt; - p3(3) = 0.0; + p3(1) = isign*Real(0.5)*hx - umac(i+1,j+joff,k)*dt; + p3(2) = jsign*Real(0.5)*hy - v*dt; + p3(3) = Real(0); for(int n=1; n<=7; ++n){ slope_tmp(n) = slopes(i+ioff,j+joff,k,n-1); @@ -1690,7 +1690,7 @@ BDS::ComputeConc (Box const& bx, // correct \Gamma^{x+} with \Gamma^{x+,z+} //////////////////////////////////////////////// - if (wmac(i+ioff,j+joff,k+1) > 0.0) { + if (wmac(i+ioff,j+joff,k+1) > Real(0)) { ksign = 1.0; koff = 0; } else { @@ -1698,31 +1698,31 @@ BDS::ComputeConc (Box const& bx, koff = 1; } - vv = 0.0; - if (vmac(i,j,k)*vmac(i+ioff,j,k+koff) > 0.0) { + vv = Real(0); + if (vmac(i,j,k)*vmac(i+ioff,j,k+koff) > Real(0)) { vv = vmac(i+ioff,j,k+koff); } - uu = 0.0; - if (umac(i+1,j+joff,k)*umac(i+1,j+joff,k+koff) > 0.0) { + uu = Real(0); + if (umac(i+1,j+joff,k)*umac(i+1,j+joff,k+koff) > Real(0)) { uu = umac(i+1,j+joff,k+koff); } - p1(1) = isign*0.5*hx; - p1(2) = jsign*0.5*hy; - p1(3) = ksign*0.5*hz; + p1(1) = isign*Real(0.5)*hx; + p1(2) = jsign*Real(0.5)*hy; + p1(3) = ksign*Real(0.5)*hz; - p2(1) = isign*0.5*hx; - p2(2) = jsign*0.5*hy - vmac(i,j,k)*dt; - p2(3) = ksign*0.5*hz; + p2(1) = isign*Real(0.5)*hx; + p2(2) = jsign*Real(0.5)*hy - vmac(i,j,k)*dt; + p2(3) = ksign*Real(0.5)*hz; - p3(1) = isign*0.5*hx - umac(i+1,j+joff,k)*dt; - p3(2) = jsign*0.5*hy - vmac(i,j,k)*dt; - p3(3) = ksign*0.5*hz; + p3(1) = isign*Real(0.5)*hx - umac(i+1,j+joff,k)*dt; + p3(2) = jsign*Real(0.5)*hy - vmac(i,j,k)*dt; + p3(3) = ksign*Real(0.5)*hz; - p4(1) = isign*0.5*hx - uu*dt; - p4(2) = jsign*0.5*hy - vv*dt; - p4(3) = ksign*0.5*hz - wmac(i+ioff,j+joff,k+1)*dt; + p4(1) = isign*Real(0.5)*hx - uu*dt; + p4(2) = jsign*Real(0.5)*hy - vv*dt; + p4(3) = ksign*Real(0.5)*hz - wmac(i+ioff,j+joff,k+1)*dt; for(int n=1; n<=7; ++n){ slope_tmp(n) = slopes(i+ioff,j+joff,k+koff,n-1); @@ -1734,26 +1734,26 @@ BDS::ComputeConc (Box const& bx, val1 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p1(ll) + sixth*(p2(ll)+p3(ll)+p4(ll)); + del(ll) = Real(0.5)*p1(ll) + sixth*(p2(ll)+p3(ll)+p4(ll)); } val2 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p2(ll) + sixth*(p1(ll)+p3(ll)+p4(ll)); + del(ll) = Real(0.5)*p2(ll) + sixth*(p1(ll)+p3(ll)+p4(ll)); } val3 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p3(ll) + sixth*(p2(ll)+p1(ll)+p4(ll)); + del(ll) = Real(0.5)*p3(ll) + sixth*(p2(ll)+p1(ll)+p4(ll)); } val4 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p4(ll) + sixth*(p2(ll)+p3(ll)+p1(ll)); + del(ll) = Real(0.5)*p4(ll) + sixth*(p2(ll)+p3(ll)+p1(ll)); } val5 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); - gamma2 = -0.8*val1 + 0.45*(val2+val3+val4+val5); + gamma2 = -Real(0.8)*val1 + Real(0.45)*(val2+val3+val4+val5); // divu source term if (iconserv[icomp]) { @@ -1768,7 +1768,7 @@ BDS::ComputeConc (Box const& bx, // correct \Gamma^{x+} with \Gamma^{x+,z-} //////////////////////////////////////////////// - if (wmac(i+ioff,j+joff,k) > 0.0) { + if (wmac(i+ioff,j+joff,k) > Real(0)) { ksign = 1.0; koff = -1; } else { @@ -1776,31 +1776,31 @@ BDS::ComputeConc (Box const& bx, koff = 0; } - vv = 0.0; - if (vmac(i,j,k)*vmac(i+ioff,j,k+koff) > 0.0) { + vv = Real(0); + if (vmac(i,j,k)*vmac(i+ioff,j,k+koff) > Real(0)) { vv = vmac(i+ioff,j,k+koff); } - uu = 0.0; + uu = Real(0); if (umac(i+1,j+joff,k)*umac(i+1,j+joff,k+koff) > 0.0) { uu = umac(i+1,j+joff,k+koff); } - p1(1) = isign*0.5*hx; - p1(2) = jsign*0.5*hy; - p1(3) = ksign*0.5*hz; + p1(1) = isign*Real(0.5)*hx; + p1(2) = jsign*Real(0.5)*hy; + p1(3) = ksign*Real(0.5)*hz; - p2(1) = isign*0.5*hx; - p2(2) = jsign*0.5*hy - vmac(i,j,k)*dt; - p2(3) = ksign*0.5*hz; + p2(1) = isign*Real(0.5)*hx; + p2(2) = jsign*Real(0.5)*hy - vmac(i,j,k)*dt; + p2(3) = ksign*Real(0.5)*hz; - p3(1) = isign*0.5*hx - umac(i+1,j+joff,k)*dt; - p3(2) = jsign*0.5*hy - vmac(i,j,k)*dt; - p3(3) = ksign*0.5*hz; + p3(1) = isign*Real(0.5)*hx - umac(i+1,j+joff,k)*dt; + p3(2) = jsign*Real(0.5)*hy - vmac(i,j,k)*dt; + p3(3) = ksign*Real(0.5)*hz; - p4(1) = isign*0.5*hx - uu*dt; - p4(2) = jsign*0.5*hy - vv*dt; - p4(3) = ksign*0.5*hz - wmac(i+ioff,j+joff,k)*dt; + p4(1) = isign*Real(0.5)*hx - uu*dt; + p4(2) = jsign*Real(0.5)*hy - vv*dt; + p4(3) = ksign*Real(0.5)*hz - wmac(i+ioff,j+joff,k)*dt; for(int n=1; n<=7; ++n){ slope_tmp(n) = slopes(i+ioff,j+joff,k+koff,n-1); @@ -1812,26 +1812,26 @@ BDS::ComputeConc (Box const& bx, val1 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p1(ll) + sixth*(p2(ll)+p3(ll)+p4(ll)); + del(ll) = Real(0.5)*p1(ll) + sixth*(p2(ll)+p3(ll)+p4(ll)); } val2 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p2(ll) + sixth*(p1(ll)+p3(ll)+p4(ll)); + del(ll) = Real(0.5)*p2(ll) + sixth*(p1(ll)+p3(ll)+p4(ll)); } val3 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p3(ll) + sixth*(p2(ll)+p1(ll)+p4(ll)); + del(ll) = Real(0.5)*p3(ll) + sixth*(p2(ll)+p1(ll)+p4(ll)); } val4 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); for(int ll=1; ll<=3; ++ll ){ - del(ll) = 0.5*p4(ll) + sixth*(p2(ll)+p3(ll)+p1(ll)); + del(ll) = Real(0.5)*p4(ll) + sixth*(p2(ll)+p3(ll)+p1(ll)); } val5 = eval(s(i+ioff,j+joff,k+koff,icomp),slope_tmp,del); - gamma2 = -0.8*val1 + 0.45*(val2+val3+val4+val5); + gamma2 = -Real(0.8)*val1 + Real(0.45)*(val2+val3+val4+val5); // divu source term if (iconserv[icomp]) { @@ -1853,30 +1853,30 @@ BDS::ComputeConc (Box const& bx, // compute \Gamma^{x-} without corner corrections //////////////////////////////////////////////// - if (umac(i,j+joff,k) > 0.0) { - isign = 1.0; + if (umac(i,j+joff,k) > Real(0)) { + isign = Real(1); ioff = -1; } else { - isign = -1.0; + isign = -Real(1); ioff = 0; } - v = 0.0; - if (vmac(i,j,k)*vmac(i+ioff,j,k) > 0.0) { + v = Real(0); + if (vmac(i,j,k)*vmac(i+ioff,j,k) > Real(0)) { v = vmac(i+ioff,j,k); } - p1(1) = isign*0.5*hx; - p1(2) = jsign*0.5*hy; - p1(3) = 0.0; + p1(1) = isign*Real(0.5)*hx; + p1(2) = jsign*Real(0.5)*hy; + p1(3) = Real(0); - p2(1) = isign*0.5*hx; - p2(2) = jsign*0.5*hy - vmac(i,j,k)*dt; - p2(3) = 0.0; + p2(1) = isign*Real(0.5)*hx; + p2(2) = jsign*Real(0.5)*hy - vmac(i,j,k)*dt; + p2(3) = Real(0); - p3(1) = isign*0.5*hx - umac(i,j+joff,k)*dt; - p3(2) = jsign*0.5*hy - v*dt; - p3(3) = 0.0; + p3(1) = isign*Real(0.5)*hx - umac(i,j+joff,k)*dt; + p3(2) = jsign*Real(0.5)*hy - v*dt; + p3(3) = Real(0); for(int n=1; n<=7; ++n){ slope_tmp(n) = slopes(i+ioff,j+joff,k,n-1); @@ -1911,7 +1911,7 @@ BDS::ComputeConc (Box const& bx, // correct \Gamma^{x-} with \Gamma^{x-,z+} //////////////////////////////////////////////// - if (wmac(i+ioff,j+joff,k+1) > 0.0) { + if (wmac(i+ioff,j+joff,k+1) > Real(0)) { ksign = 1.0; koff = 0; } else { @@ -1919,12 +1919,12 @@ BDS::ComputeConc (Box const& bx, koff = 1; } - vv = 0.0; - if (vmac(i,j,k)*vmac(i+ioff,j,k+koff) > 0.0) { + vv = Real(0); + if (vmac(i,j,k)*vmac(i+ioff,j,k+koff) > Real(0)) { vv = vmac(i+ioff,j,k+koff); } - uu = 0.0; + uu = Real(0); if (umac(i,j+joff,k)*umac(i,j+joff,k+koff) > 0.0) { uu = umac(i,j+joff,k+koff); } diff --git a/Docs/source/conf.py b/Docs/source/conf.py index af60791f6..0d9fc90c5 100644 --- a/Docs/source/conf.py +++ b/Docs/source/conf.py @@ -16,15 +16,6 @@ # add these directories to sys.path here. If the directory is relative to the # documentation root, use os.path.abspath to make it absolute, like shown here. -### Added to get Breathe working in Git workflow -import os -import sys -import breathe - -sys.path.insert(0, os.path.abspath('../../')) -sys.path.append(os.path.dirname(breathe.__file__)) -### end Breathe - # sys.path.insert(0, os.path.abspath('.')) import sphinx_rtd_theme @@ -49,8 +40,7 @@ def get_amrex_hydro_version(): 'sphinx.ext.viewcode', 'sphinx.ext.intersphinx', 'sphinx.ext.autosectionlabel', - 'sphinx_toolbox.collapse', - 'breathe'] + 'sphinx_toolbox.collapse'] # Allow for same subheading in multiple sections autosectionlabel_prefix_document = True @@ -212,11 +202,3 @@ def get_amrex_hydro_version(): author, 'AMReX-Hydro team', 'One line description of project.', 'Miscellaneous'), ] - - -# -- Breathe for Doxygen Conversion ------------------------------------------- -# see https://github.com/michaeljones/breathe - -breathe_projects = { "amrex-hydro": "../Doxygen/xml/"} - -breathe_default_project = "amrex-hydro" diff --git a/EBMOL/hydro_ebmol_edge_state_K.H b/EBMOL/hydro_ebmol_edge_state_K.H index 97bc8de78..45d0f4041 100644 --- a/EBMOL/hydro_ebmol_edge_state_K.H +++ b/EBMOL/hydro_ebmol_edge_state_K.H @@ -5,8 +5,6 @@ * */ -using namespace amrex::literals; - #ifndef HYDRO_EBMOL_EDGE_STATE_K_H_ #define HYDRO_EBMOL_EDGE_STATE_K_H_ @@ -21,9 +19,9 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real hydro_ebmol_xedge_state_extdir ( AMREX_D_DECL(int i, int j, int k), int n, amrex::Array4 const& q, amrex::Array4 const& umac, - AMREX_D_DECL(amrex::Array4 const& fcx, - amrex::Array4 const& fcy, - amrex::Array4 const& fcz), + AMREX_D_DECL(amrex::Array4 const& fcx, + amrex::Array4 const& fcy, + amrex::Array4 const& fcz), amrex::Array4 const& ccc, amrex::Array4 const& vfrac, amrex::Array4 const& flag, diff --git a/Godunov/hydro_godunov_edge_state_2D.cpp b/Godunov/hydro_godunov_edge_state_2D.cpp index 389148a4a..864f8f1ee 100644 --- a/Godunov/hydro_godunov_edge_state_2D.cpp +++ b/Godunov/hydro_godunov_edge_state_2D.cpp @@ -155,8 +155,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, if (use_forces_in_trans && fq) { - lo += 0.5*l_dt*fq(i-1,j,k,n); - hi += 0.5*l_dt*fq(i ,j,k,n); + lo += Real(0.5)*l_dt*fq(i-1,j,k,n); + hi += Real(0.5)*l_dt*fq(i ,j,k,n); } const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); @@ -174,8 +174,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, if (use_forces_in_trans && fq) { - lo += 0.5*l_dt*fq(i,j-1,k,n); - hi += 0.5*l_dt*fq(i,j ,k,n); + lo += Real(0.5)*l_dt*fq(i,j-1,k,n); + hi += Real(0.5)*l_dt*fq(i,j ,k,n); } const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); @@ -207,8 +207,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, HydroBC::SetEdgeBCsHi(1, i, j, k, n, qnph, l_yzlo, l_yzhi, vmac(i,j,k), bc.hi(1), dhi.y, is_velocity); Real st = (vad >= 0.) ? l_yzlo : l_yzhi; - Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0; - yzlo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_yzhi + l_yzlo); + Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : Real(1); + yzlo(i,j,k,n) = fu*st + (Real(1) - fu) * Real(0.5) * (l_yzhi + l_yzlo); }); // @@ -225,31 +225,31 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, // q + dx/2 q_x - dt/2 (u q_x + q u_x + (v q)_y) which is equivalent to // --> q + dx/2 q_x - dt/2 ( div (uvec q) ) Real quxl = (umac(i,j,k) - umac(i-1,j,k)) * q(i-1,j,k,n); - stl += ( - (0.5*dtdx) * quxl - - (0.5*dtdy) * (yzlo(i-1,j+1,k ,n)*vmac(i-1,j+1,k ) + stl += ( - (Real(0.5)*dtdx) * quxl + - (Real(0.5)*dtdy) * (yzlo(i-1,j+1,k ,n)*vmac(i-1,j+1,k ) -yzlo(i-1,j ,k ,n)*vmac(i-1,j ,k )) ); // Here we adjust for non-conservative by removing the q divu contribution to get // q + dx/2 q_x - dt/2 ( div (uvec q) - q divu ) which is equivalent to // --> q + dx/2 q_x - dt/2 ( uvec dot grad q) - stl += (!iconserv[n]) ? 0.5*l_dt* q(i-1,j,k,n)*divu(i-1,j,k) : 0.; + stl += (!iconserv[n]) ? Real(0.5)*l_dt* q(i-1,j,k,n)*divu(i-1,j,k) : Real(0); - stl += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i-1,j,k,n) : 0.; + stl += (!use_forces_in_trans && fq) ? Real(0.5)*l_dt*fq(i-1,j,k,n) : Real(0); // Here we add uq/r for RZ - stl += (is_rz) ? -0.25 * l_dt * q(i-1,j,k,n)*( umac(i,j,k) + umac(i-1,j,k) ) / ( dx*(amrex::Math::abs(Real(i)-0.5)) ) : 0.; + stl += (is_rz) ? -Real(0.25) * l_dt * q(i-1,j,k,n)*( umac(i,j,k) + umac(i-1,j,k) ) / ( dx*(amrex::Math::abs(Real(i)-Real(0.5))) ) : Real(0); // High side Real quxh = (umac(i+1,j,k) - umac(i,j,k)) * q(i,j,k,n); - sth += ( - (0.5*dtdx) * quxh - - (0.5*dtdy)*(yzlo(i,j+1,k,n)*vmac(i,j+1,k) + sth += ( - (Real(0.5)*dtdx) * quxh + - (Real(0.5)*dtdy)*(yzlo(i,j+1,k,n)*vmac(i,j+1,k) -yzlo(i,j ,k,n)*vmac(i,j ,k)) ); - sth += (!iconserv[n]) ? 0.5*l_dt* q(i ,j,k,n)*divu(i,j,k) : 0.; + sth += (!iconserv[n]) ? Real(0.5)*l_dt* q(i ,j,k,n)*divu(i,j,k) : Real(0); - sth += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i ,j,k,n) : 0.; + sth += (!use_forces_in_trans && fq) ? Real(0.5)*l_dt*fq(i ,j,k,n) : Real(0); - sth += (is_rz) ? -0.25 * l_dt * q(i,j,k,n)*( umac(i,j,k) + umac(i+1,j,k) ) / ( dx*(amrex::Math::abs(Real(i)+0.5)) ) : 0.; + sth += (is_rz) ? -Real(0.25) * l_dt * q(i,j,k,n)*( umac(i,j,k) + umac(i+1,j,k) ) / ( dx*(amrex::Math::abs(Real(i)+Real(0.5))) ) : Real(0); const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); @@ -270,7 +270,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, } Real temp = (umac(i,j,k) >= 0.) ? stl : sth; - temp = (amrex::Math::abs(umac(i,j,k)) < small_vel) ? 0.5*(stl + sth) : temp; + temp = (amrex::Math::abs(umac(i,j,k)) < small_vel) ? Real(0.5)*(stl + sth) : temp; xedge(i,j,k,n) = temp; }); @@ -294,8 +294,8 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, Real uad = umac(i,j,k); Real st = (uad >= 0.) ? l_xzlo : l_xzhi; - Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0; - xzlo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_xzhi + l_xzlo); + Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : Real(1); + xzlo(i,j,k,n) = fu*st + (Real(1) - fu) * Real(0.5) * (l_xzhi + l_xzlo); }); // @@ -313,31 +313,31 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, // q + dy/2 q_y - dt/2 (v q_y + q v_y + (u q)_x) which is equivalent to // --> q + dy/2 q_y - dt/2 ( div (uvec q) ) Real qvyl = (vmac(i,j,k) - vmac(i,j-1,k)) * q(i,j-1,k,n); - stl += ( - (0.5*dtdy)*qvyl - - (0.5*dtdx)*(xzlo(i+1,j-1,k ,n)*umac(i+1,j-1,k ) + stl += ( - (Real(0.5)*dtdy)*qvyl + - (Real(0.5)*dtdx)*(xzlo(i+1,j-1,k ,n)*umac(i+1,j-1,k ) -xzlo(i ,j-1,k ,n)*umac(i ,j-1,k )) ); // Here we adjust for non-conservative by removing the q divu contribution to get // q + dy/2 q_y - dt/2 ( div (uvec q) - q divu ) which is equivalent to // --> q + dy/2 q_y - dt/2 ( uvec dot grad q) - stl += (!iconserv[n]) ? 0.5*l_dt* q(i,j-1,k,n)*divu(i,j-1,k) : 0.; + stl += (!iconserv[n]) ? Real(0.5)*l_dt* q(i,j-1,k,n)*divu(i,j-1,k) : Real(0); - stl += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i,j-1,k,n) : 0.; + stl += (!use_forces_in_trans && fq) ? Real(0.5)*l_dt*fq(i,j-1,k,n) : Real(0); // Here we add uq/r for RZ - stl += (is_rz) ? -0.25 * l_dt * q(i,j-1,k,n)*( umac(i,j-1,k) + umac(i+1,j-1,k) ) / ( dx*(amrex::Math::abs(Real(i)+0.5)) ) : 0.; + stl += (is_rz) ? -Real(0.25) * l_dt * q(i,j-1,k,n)*( umac(i,j-1,k) + umac(i+1,j-1,k) ) / ( dx*(amrex::Math::abs(Real(i)+Real(0.5))) ) : Real(0); // High side Real qvyh = (vmac(i,j+1,k) - vmac(i,j,k)) * q(i,j,k,n); - sth += ( - (0.5*dtdy)*qvyh - - (0.5*dtdx)*(xzlo(i+1,j,k ,n)*umac(i+1,j,k ) + sth += ( - (Real(0.5)*dtdy)*qvyh + - (Real(0.5)*dtdx)*(xzlo(i+1,j,k ,n)*umac(i+1,j,k ) -xzlo(i ,j,k ,n)*umac(i ,j,k )) ); - sth += (!iconserv[n]) ? 0.5*l_dt* q(i,j,k,n)*divu(i,j,k) : 0.; + sth += (!iconserv[n]) ? Real(0.5)*l_dt* q(i,j,k,n)*divu(i,j,k) : Real(0); - sth += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i,j,k,n) : 0.; + sth += (!use_forces_in_trans && fq) ? Real(0.5)*l_dt*fq(i,j,k,n) : Real(0); - sth += (is_rz) ? -0.25 * l_dt * q(i,j,k,n)*( umac(i,j ,k) + umac(i+1,j ,k) ) / ( dx*(amrex::Math::abs(Real(i)+0.5)) ) : 0.; + sth += (is_rz) ? -Real(0.25) * l_dt * q(i,j,k,n)*( umac(i,j ,k) + umac(i+1,j ,k) ) / ( dx*(amrex::Math::abs(Real(i)+Real(0.5))) ) : Real(0); const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); @@ -347,18 +347,18 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, if (!allow_inflow_on_outflow) { if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) { - if ( vmac(i,j,k) >= 0. && n==YVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); + if ( vmac(i,j,k) >= Real(0) && n==YVEL && is_velocity ) sth = amrex::min(sth,Real(0)); stl = sth; } if ( (j==dhi.y+1) && (bc.hi(1) == BCType::foextrap || bc.hi(1) == BCType::hoextrap) ) { - if ( vmac(i,j,k) <= 0. && n==YVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); + if ( vmac(i,j,k) <= Real(0) && n==YVEL && is_velocity ) stl = amrex::max(stl,Real(0)); sth = stl; } } - Real temp = (vmac(i,j,k) >= 0.) ? stl : sth; - temp = (amrex::Math::abs(vmac(i,j,k)) < small_vel) ? 0.5*(stl + sth) : temp; + Real temp = (vmac(i,j,k) >= Real(0)) ? stl : sth; + temp = (amrex::Math::abs(vmac(i,j,k)) < small_vel) ? Real(0.5)*(stl + sth) : temp; yedge(i,j,k,n) = temp; }); diff --git a/Godunov/hydro_godunov_edge_state_3D.cpp b/Godunov/hydro_godunov_edge_state_3D.cpp index 2247c16e3..0a5b07854 100644 --- a/Godunov/hydro_godunov_edge_state_3D.cpp +++ b/Godunov/hydro_godunov_edge_state_3D.cpp @@ -170,15 +170,15 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, xebox, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { Real uad = umac(i,j,k); - Real fux = (amrex::Math::abs(uad) < small_vel)? 0. : 1.; - bool uval = uad >= 0.; + Real fux = (amrex::Math::abs(uad) < small_vel)? Real(0) : Real(1); + bool uval = uad >= Real(0); Real lo = Ipx(i-1,j,k,n); Real hi = Imx(i ,j,k,n); if (use_forces_in_trans && fq) { - lo += 0.5*l_dt*fq(i-1,j,k,n); - hi += 0.5*l_dt*fq(i ,j,k,n); + lo += Real(0.5)*l_dt*fq(i-1,j,k,n); + hi += Real(0.5)*l_dt*fq(i ,j,k,n); } const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); @@ -189,21 +189,21 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, xlo(i,j,k,n) = lo; xhi(i,j,k,n) = hi; Real st = (uval) ? lo : hi; - Imx(i,j,k,n) = fux*st + (1. - fux)*0.5*(hi + lo); + Imx(i,j,k,n) = fux*st + (Real(1) - fux)*Real(0.5)*(hi + lo); }, yebox, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { Real vad = vmac(i,j,k); - Real fuy = (amrex::Math::abs(vad) < small_vel)? 0. : 1.; - bool vval = vad >= 0.; + Real fuy = (amrex::Math::abs(vad) < small_vel)? Real(0) : Real(1); + bool vval = vad >= Real(0); Real lo = Ipy(i,j-1,k,n); Real hi = Imy(i,j ,k,n); if (use_forces_in_trans && fq) { - lo += 0.5*l_dt*fq(i,j-1,k,n); - hi += 0.5*l_dt*fq(i,j ,k,n); + lo += Real(0.5)*l_dt*fq(i,j-1,k,n); + hi += Real(0.5)*l_dt*fq(i,j ,k,n); } const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); @@ -214,20 +214,20 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, ylo(i,j,k,n) = lo; yhi(i,j,k,n) = hi; Real st = (vval) ? lo : hi; - Imy(i,j,k,n) = fuy*st + (1. - fuy)*0.5*(hi + lo); + Imy(i,j,k,n) = fuy*st + (Real(1) - fuy)*Real(0.5)*(hi + lo); }, zebox, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { Real wad = wmac(i,j,k); - Real fuz = (amrex::Math::abs(wad) < small_vel) ? 0. : 1.; - bool wval = wad >= 0.; + Real fuz = (amrex::Math::abs(wad) < small_vel) ? Real(0) : Real(1); + bool wval = wad >= Real(0); Real lo = Ipz(i,j,k-1,n); Real hi = Imz(i,j,k ,n); if (use_forces_in_trans && fq) { - lo += 0.5*l_dt*fq(i,j,k-1,n); - hi += 0.5*l_dt*fq(i,j,k ,n); + lo += Real(0.5)*l_dt*fq(i,j,k-1,n); + hi += Real(0.5)*l_dt*fq(i,j,k ,n); } const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); @@ -238,7 +238,7 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, zlo(i,j,k,n) = lo; zhi(i,j,k,n) = hi; Real st = (wval) ? lo : hi; - Imz(i,j,k,n) = fuz*st + (1. - fuz)*0.5*(hi + lo); + Imz(i,j,k,n) = fuz*st + (Real(1) - fuz)*Real(0.5)*(hi + lo); } ); @@ -265,9 +265,9 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, HydroBC::SetEdgeBCsHi(2, i, j, k, n, qnph, l_zylo, l_zyhi, wmac(i,j,k), bc.hi(2), dhi.z, is_velocity); Real wad = wmac(i,j,k); - Real st = (wad >= 0.) ? l_zylo : l_zyhi; - Real fu = (amrex::Math::abs(wad) < small_vel) ? 0.0 : 1.0; - zylo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_zyhi + l_zylo); + Real st = (wad >= Real(0)) ? l_zylo : l_zyhi; + Real fu = (amrex::Math::abs(wad) < small_vel) ? Real(0) : Real(1); + zylo(i,j,k,n) = fu*st + (Real(1) - fu) * Real(0.5) * (l_zyhi + l_zylo); }, Box(yzlo), ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept @@ -283,9 +283,9 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, HydroBC::SetEdgeBCsHi(1, i, j, k, n, qnph, l_yzlo, l_yzhi, vmac(i,j,k), bc.hi(1), dhi.y, is_velocity); Real vad = vmac(i,j,k); - Real st = (vad >= 0.) ? l_yzlo : l_yzhi; - Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0; - yzlo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_yzhi + l_yzlo); + Real st = (vad >= Real(0)) ? l_yzlo : l_yzhi; + Real fu = (amrex::Math::abs(vad) < small_vel) ? Real(0) : Real(1); + yzlo(i,j,k,n) = fu*st + (Real(1) - fu) * Real(0.5) * (l_yzhi + l_yzlo); }); @@ -302,30 +302,30 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, // q + dx/2 q_x - dt/2 (u q_x + q u_x + (v q)_y + (w q)_z) which is equivalent to // --> q + dx/2 q_x - dt/2 ( div (uvec q) ) Real quxl = (umac(i,j,k) - umac(i-1,j,k)) * q(i-1,j,k,n); - stl += ( - (0.5*dtdx) * quxl - - (0.5*dtdy)*(yzlo(i-1,j+1,k ,n)*vmac(i-1,j+1,k ) + stl += ( - (Real(0.5)*dtdx) * quxl + - (Real(0.5)*dtdy)*(yzlo(i-1,j+1,k ,n)*vmac(i-1,j+1,k ) - yzlo(i-1,j ,k ,n)*vmac(i-1,j ,k )) - - (0.5*dtdz)*(zylo(i-1,j ,k+1,n)*wmac(i-1,j ,k+1) + - (Real(0.5)*dtdz)*(zylo(i-1,j ,k+1,n)*wmac(i-1,j ,k+1) - zylo(i-1,j ,k ,n)*wmac(i-1,j ,k )) ); // Here we adjust for non-conservative by removing the q divu contribution to get // q + dx/2 q_x - dt/2 ( div (uvec q) - q divu ) which is equivalent to // --> q + dx/2 q_x - dt/2 ( uvec dot grad q) - stl += (!iconserv[n]) ? 0.5*l_dt* q(i-1,j,k,n)*divu(i-1,j,k) : 0.; + stl += (!iconserv[n]) ? Real(0.5)*l_dt* q(i-1,j,k,n)*divu(i-1,j,k) : Real(0); - stl += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i-1,j,k,n) : 0.; + stl += (!use_forces_in_trans && fq) ? Real(0.5)*l_dt*fq(i-1,j,k,n) : Real(0); // High side Real quxh = (umac(i+1,j,k) - umac(i,j,k)) * q(i,j,k,n); - sth += ( - (0.5*dtdx) * quxh - - (0.5*dtdy)*(yzlo(i,j+1,k ,n)*vmac(i,j+1,k ) + sth += ( - (Real(0.5)*dtdx) * quxh + - (Real(0.5)*dtdy)*(yzlo(i,j+1,k ,n)*vmac(i,j+1,k ) - yzlo(i,j ,k ,n)*vmac(i,j ,k )) - - (0.5*dtdz)*(zylo(i,j ,k+1,n)*wmac(i,j ,k+1) + - (Real(0.5)*dtdz)*(zylo(i,j ,k+1,n)*wmac(i,j ,k+1) - zylo(i,j ,k ,n)*wmac(i,j ,k )) ); - sth += (!iconserv[n]) ? 0.5*l_dt* q(i ,j,k,n)*divu(i,j,k) : 0.; + sth += (!iconserv[n]) ? Real(0.5)*l_dt* q(i ,j,k,n)*divu(i,j,k) : Real(0); - sth += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i ,j,k,n) : 0.; + sth += (!use_forces_in_trans && fq) ? Real(0.5)*l_dt*fq(i ,j,k,n) : Real(0); const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); @@ -335,18 +335,18 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, if (!allow_inflow_on_outflow) { if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) { - if ( umac(i,j,k) >= 0. && n==XVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); + if ( umac(i,j,k) >= Real(0) && n==XVEL && is_velocity ) sth = amrex::min(sth,Real(0)); stl = sth; } if ( (i==dhi.x+1) && (bc.hi(0) == BCType::foextrap || bc.hi(0) == BCType::hoextrap) ) { - if ( umac(i,j,k) <= 0. && n==XVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); + if ( umac(i,j,k) <= Real(0) && n==XVEL && is_velocity ) stl = amrex::max(stl,Real(0)); sth = stl; } } - Real temp = (umac(i,j,k) >= 0.) ? stl : sth; - temp = (amrex::Math::abs(umac(i,j,k)) < small_vel) ? 0.5*(stl + sth) : temp; + Real temp = (umac(i,j,k) >= Real(0)) ? stl : sth; + temp = (amrex::Math::abs(umac(i,j,k)) < small_vel) ? Real(0.5)*(stl + sth) : temp; xedge(i,j,k,n) = temp; }); @@ -371,9 +371,9 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, HydroBC::SetEdgeBCsHi(0, i, j, k, n, qnph, l_xzlo, l_xzhi, umac(i,j,k), bc.hi(0), dhi.x, is_velocity); Real uad = umac(i,j,k); - Real st = (uad >= 0.) ? l_xzlo : l_xzhi; - Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0; - xzlo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_xzhi + l_xzlo); + Real st = (uad >= Real(0)) ? l_xzlo : l_xzhi; + Real fu = (amrex::Math::abs(uad) < small_vel) ? Real(0) : Real(1); + xzlo(i,j,k,n) = fu*st + (Real(1) - fu) * Real(0.5) * (l_xzhi + l_xzlo); }, Box(zxlo), ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept @@ -389,9 +389,9 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, HydroBC::SetEdgeBCsHi(2, i, j, k, n, qnph, l_zxlo, l_zxhi, wmac(i,j,k), bc.hi(2), dhi.z, is_velocity); Real wad = wmac(i,j,k); - Real st = (wad >= 0.) ? l_zxlo : l_zxhi; - Real fu = (amrex::Math::abs(wad) < small_vel) ? 0.0 : 1.0; - zxlo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_zxhi + l_zxlo); + Real st = (wad >= Real(0)) ? l_zxlo : l_zxhi; + Real fu = (amrex::Math::abs(wad) < small_vel) ? Real(0) : Real(1); + zxlo(i,j,k,n) = fu*st + (Real(1) - fu) * Real(0.5) * (l_zxhi + l_zxlo); }); // @@ -406,30 +406,30 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, // q + dy/2 q_y - dt/2 (v q_y + q v_y + (u q)_x + (w q)_z) which is equivalent to // --> q + dy/2 q_y - dt/2 ( div (uvec q) ) Real qvyl = (vmac(i,j,k) - vmac(i,j-1,k)) * q(i,j-1,k,n); - stl += ( - (0.5*dtdy) * qvyl - - (0.5*dtdx)*(xzlo(i+1,j-1,k ,n)*umac(i+1,j-1,k ) + stl += ( - (Real(0.5)*dtdy) * qvyl + - (Real(0.5)*dtdx)*(xzlo(i+1,j-1,k ,n)*umac(i+1,j-1,k ) - xzlo(i ,j-1,k ,n)*umac(i ,j-1,k )) - - (0.5*dtdz)*(zxlo(i ,j-1,k+1,n)*wmac(i ,j-1,k+1) + - (Real(0.5)*dtdz)*(zxlo(i ,j-1,k+1,n)*wmac(i ,j-1,k+1) - zxlo(i ,j-1,k ,n)*wmac(i ,j-1,k )) ); // Here we adjust for non-conservative by removing the q divu contribution to get // q + dy/2 q_y - dt/2 ( div (uvec q) - q divu ) which is equivalent to // --> q + dy/2 q_y - dt/2 ( uvec dot grad q) - stl += (!iconserv[n]) ? 0.5*l_dt* q(i,j-1,k,n)*divu(i,j-1,k) : 0.; + stl += (!iconserv[n]) ? Real(0.5)*l_dt* q(i,j-1,k,n)*divu(i,j-1,k) : Real(0); - stl += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i,j-1,k,n) : 0.; + stl += (!use_forces_in_trans && fq) ? Real(0.5)*l_dt*fq(i,j-1,k,n) : Real(0); // High side Real qvyh = (vmac(i,j+1,k) - vmac(i,j,k)) * q(i,j,k,n); - sth += ( - (0.5*dtdy) * qvyh - - (0.5*dtdx)*(xzlo(i+1,j,k ,n)*umac(i+1,j,k ) + sth += ( - (Real(0.5)*dtdy) * qvyh + - (Real(0.5)*dtdx)*(xzlo(i+1,j,k ,n)*umac(i+1,j,k ) - xzlo(i ,j,k ,n)*umac(i ,j,k )) - - (0.5*dtdz)*(zxlo(i ,j,k+1,n)*wmac(i ,j,k+1) + - (Real(0.5)*dtdz)*(zxlo(i ,j,k+1,n)*wmac(i ,j,k+1) - zxlo(i ,j,k ,n)*wmac(i ,j,k )) ); - sth += (!iconserv[n]) ? 0.5*l_dt* q(i,j,k,n)*divu(i,j,k) : 0.; + sth += (!iconserv[n]) ? Real(0.5)*l_dt* q(i,j,k,n)*divu(i,j,k) : Real(0); - sth += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i,j,k,n) : 0.; + sth += (!use_forces_in_trans && fq) ? Real(0.5)*l_dt*fq(i,j,k,n) : Real(0); const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); @@ -439,18 +439,18 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, if (!allow_inflow_on_outflow) { if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) { - if ( vmac(i,j,k) >= 0. && n==YVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); + if ( vmac(i,j,k) >= Real(0) && n==YVEL && is_velocity ) sth = amrex::min(sth,Real(0)); stl = sth; } if ( (j==dhi.y+1) && (bc.hi(1) == BCType::foextrap || bc.hi(1) == BCType::hoextrap) ) { - if ( vmac(i,j,k) <= 0. && n==YVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); + if ( vmac(i,j,k) <= Real(0) && n==YVEL && is_velocity ) stl = amrex::max(stl,Real(0)); sth = stl; } } - Real temp = (vmac(i,j,k) >= 0.) ? stl : sth; - temp = (amrex::Math::abs(vmac(i,j,k)) < small_vel) ? 0.5*(stl + sth) : temp; + Real temp = (vmac(i,j,k) >= Real(0)) ? stl : sth; + temp = (amrex::Math::abs(vmac(i,j,k)) < small_vel) ? Real(0.5)*(stl + sth) : temp; yedge(i,j,k,n) = temp; }); @@ -475,9 +475,9 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, HydroBC::SetEdgeBCsHi(0, i, j, k, n, qnph, l_xylo, l_xyhi, umac(i,j,k), bc.hi(0), dhi.x, is_velocity); Real uad = umac(i,j,k); - Real st = (uad >= 0.) ? l_xylo : l_xyhi; - Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0; - xylo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_xyhi + l_xylo); + Real st = (uad >= Real(0)) ? l_xylo : l_xyhi; + Real fu = (amrex::Math::abs(uad) < small_vel) ? Real(0) : Real(1); + xylo(i,j,k,n) = fu*st + (Real(1) - fu) * Real(0.5) * (l_xyhi + l_xylo); }, Box(yxlo), ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept @@ -493,9 +493,9 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, HydroBC::SetEdgeBCsHi(1, i, j, k, n, qnph, l_yxlo, l_yxhi, vmac(i,j,k), bc.hi(1), dhi.y, is_velocity); Real vad = vmac(i,j,k); - Real st = (vad >= 0.) ? l_yxlo : l_yxhi; - Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0; - yxlo(i,j,k,n) = fu*st + (1.0 - fu) * 0.5 * (l_yxhi + l_yxlo); + Real st = (vad >= Real(0)) ? l_yxlo : l_yxhi; + Real fu = (amrex::Math::abs(vad) < small_vel) ? Real(0) : Real(1); + yxlo(i,j,k,n) = fu*st + (Real(1) - fu) * Real(0.5) * (l_yxhi + l_yxlo); }); // @@ -510,30 +510,30 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, // q + dz/2 q_z - dt/2 (w q_z + q w_z + (u q)_x + (v q)_y) which is equivalent to // --> q + dz/2 q_z - dt/2 ( div (uvec q) ) Real qwzl = (wmac(i,j,k) - wmac(i,j,k-1)) * q(i,j,k-1,n); - stl += ( - (0.5*dtdz) * qwzl - -(0.5*dtdx)*(xylo(i+1,j ,k-1,n)*umac(i+1,j ,k-1) + stl += ( -(Real(0.5)*dtdz) * qwzl + -(Real(0.5)*dtdx)*(xylo(i+1,j ,k-1,n)*umac(i+1,j ,k-1) -xylo(i ,j ,k-1,n)*umac(i ,j ,k-1)) - -(0.5*dtdy)*(yxlo(i ,j+1,k-1,n)*vmac(i ,j+1,k-1) + -(Real(0.5)*dtdy)*(yxlo(i ,j+1,k-1,n)*vmac(i ,j+1,k-1) -yxlo(i ,j ,k-1,n)*vmac(i ,j ,k-1)) ); // Here we adjust for non-conservative by removing the q divu contribution to get // q + dz/2 q_z - dt/2 ( div (uvec q) - q divu ) which is equivalent to // --> q + dz/2 q_z - dt/2 ( uvec dot grad q) - stl += (!iconserv[n]) ? 0.5*l_dt* q(i,j,k-1,n)*divu(i,j,k-1) : 0.; + stl += (!iconserv[n]) ? Real(0.5)*l_dt* q(i,j,k-1,n)*divu(i,j,k-1) : Real(0); - stl += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i,j,k-1,n) : 0.; + stl += (!use_forces_in_trans && fq) ? Real(0.5)*l_dt*fq(i,j,k-1,n) : Real(0); // High side Real qwzh = (wmac(i,j,k+1) - wmac(i,j,k)) * q(i,j,k,n); - sth += ( - (0.5*dtdz) * qwzh - - (0.5*dtdx)*(xylo(i+1,j ,k,n)*umac(i+1,j ,k) + sth += ( - (Real(0.5)*dtdz) * qwzh + - (Real(0.5)*dtdx)*(xylo(i+1,j ,k,n)*umac(i+1,j ,k) -xylo(i ,j ,k,n)*umac(i ,j ,k)) - - (0.5*dtdy)*(yxlo(i ,j+1,k,n)*vmac(i ,j+1,k) + - (Real(0.5)*dtdy)*(yxlo(i ,j+1,k,n)*vmac(i ,j+1,k) -yxlo(i ,j ,k,n)*vmac(i ,j ,k)) ); - sth += (!iconserv[n]) ? 0.5*l_dt* q(i,j,k,n)*divu(i,j,k) : 0.; + sth += (!iconserv[n]) ? Real(0.5)*l_dt* q(i,j,k,n)*divu(i,j,k) : Real(0); - sth += (!use_forces_in_trans && fq) ? 0.5*l_dt*fq(i,j,k,n) : 0.; + sth += (!use_forces_in_trans && fq) ? Real(0.5)*l_dt*fq(i,j,k,n) : Real(0); const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetEdgeBCsLo(2, i, j, k, n, qnph, stl, sth, wmac(i,j,k), bc.lo(2), dlo.z, is_velocity); @@ -542,18 +542,18 @@ Godunov::ComputeEdgeState (Box const& bx, int ncomp, if (!allow_inflow_on_outflow) { if ( (k==dlo.z) && (bc.lo(2) == BCType::foextrap || bc.lo(2) == BCType::hoextrap) ) { - if ( wmac(i,j,k) >= 0. && n==ZVEL && is_velocity ) sth = amrex::min(sth,0.0_rt); + if ( wmac(i,j,k) >= Real(0) && n==ZVEL && is_velocity ) sth = amrex::min(sth,Real(0)); stl = sth; } if ( (k==dhi.z+1) && (bc.hi(2) == BCType::foextrap || bc.hi(2) == BCType::hoextrap) ) { - if ( wmac(i,j,k) <= 0. && n==ZVEL && is_velocity ) stl = amrex::max(stl,0.0_rt); + if ( wmac(i,j,k) <= Real(0) && n==ZVEL && is_velocity ) stl = amrex::max(stl,Real(0)); sth = stl; } } - Real temp = (wmac(i,j,k) >= 0.) ? stl : sth; - temp = (amrex::Math::abs(wmac(i,j,k)) < small_vel) ? 0.5*(stl + sth) : temp; + Real temp = (wmac(i,j,k) >= Real(0)) ? stl : sth; + temp = (amrex::Math::abs(wmac(i,j,k)) < small_vel) ? Real(0.5)*(stl + sth) : temp; zedge(i,j,k,n) = temp; }); diff --git a/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp b/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp index 211593ad9..f99f0bdc4 100644 --- a/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp +++ b/Godunov/hydro_godunov_extrap_vel_to_faces_2D.cpp @@ -150,17 +150,17 @@ Godunov::ComputeAdvectiveVel ( Box const& xbx, if (l_use_forces_in_trans) { - lo += 0.5*l_dt*f(i-1,j,k,n); - hi += 0.5*l_dt*f(i ,j,k,n); + lo += Real(0.5)*l_dt*f(i-1,j,k,n); + hi += Real(0.5)*l_dt*f(i ,j,k,n); } const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetExtrapVelBCsLo(0, i, j, k, n, vel, lo, hi, bc.lo(0), dlo.x); HydroBC::SetExtrapVelBCsHi(0, i, j, k, n, vel, lo, hi, bc.hi(0), dhi.x); - Real st = ( (lo+hi) >= 0.) ? lo : hi; - bool ltm = ( (lo <= 0. && hi >= 0.) || (amrex::Math::abs(lo+hi) < small_vel) ); - u_ad(i,j,k) = ltm ? 0. : st; + Real st = ( (lo+hi) >= Real(0)) ? lo : hi; + bool ltm = ( (lo <= Real(0) && hi >= Real(0)) || (amrex::Math::abs(lo+hi) < small_vel) ); + u_ad(i,j,k) = ltm ? Real(0) : st; }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -172,17 +172,17 @@ Godunov::ComputeAdvectiveVel ( Box const& xbx, if (l_use_forces_in_trans) { - lo += 0.5*l_dt*f(i,j-1,k,n); - hi += 0.5*l_dt*f(i,j ,k,n); + lo += Real(0.5)*l_dt*f(i,j-1,k,n); + hi += Real(0.5)*l_dt*f(i,j ,k,n); } const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); HydroBC::SetExtrapVelBCsLo(1, i, j, k, n, vel, lo, hi, bc.lo(1), dlo.y); HydroBC::SetExtrapVelBCsHi(1, i, j, k, n, vel, lo, hi, bc.hi(1), dhi.y); - Real st = ( (lo+hi) >= 0.) ? lo : hi; - bool ltm = ( (lo <= 0. && hi >= 0.) || (amrex::Math::abs(lo+hi) < small_vel) ); - v_ad(i,j,k) = ltm ? 0. : st; + Real st = ( (lo+hi) >= Real(0)) ? lo : hi; + bool ltm = ( (lo <= Real(0) && hi >= Real(0)) || (amrex::Math::abs(lo+hi) < small_vel) ); + v_ad(i,j,k) = ltm ? Real(0) : st; } ); } @@ -240,8 +240,8 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp, if (l_use_forces_in_trans) { - lo += 0.5*l_dt*f(i-1,j,k,n); - hi += 0.5*l_dt*f(i ,j,k,n); + lo += Real(0.5)*l_dt*f(i-1,j,k,n); + hi += Real(0.5)*l_dt*f(i ,j,k,n); } const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); @@ -258,8 +258,8 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp, if (l_use_forces_in_trans) { - lo += 0.5*l_dt*f(i,j-1,k,n); - hi += 0.5*l_dt*f(i,j ,k,n); + lo += Real(0.5)*l_dt*f(i,j-1,k,n); + hi += Real(0.5)*l_dt*f(i,j ,k,n); } const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); @@ -293,22 +293,22 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp, Real vad = v_ad(i,j,k); Real st = (vad >= 0.) ? l_yzlo : l_yzhi; Real fu = (amrex::Math::abs(vad) < small_vel) ? 0.0 : 1.0; - yzlo(i,j,k) = fu*st + (1.0 - fu) * 0.5 * (l_yzhi + l_yzlo); + yzlo(i,j,k) = fu*st + (Real(1) - fu) * Real(0.5) * (l_yzhi + l_yzlo); }); amrex::ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { constexpr int n = 0; const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); - Real stl = xlo(i,j,k,n) - (0.25*l_dt/dy)*(v_ad(i-1,j+1,k )+v_ad(i-1,j,k))* - (yzlo(i-1,j+1,k )-yzlo(i-1,j,k)); - Real sth = xhi(i,j,k,n) - (0.25*l_dt/dy)*(v_ad(i ,j+1,k )+v_ad(i ,j,k))* - (yzlo(i ,j+1,k )-yzlo(i ,j,k)); + Real stl = xlo(i,j,k,n) - (Real(0.25)*l_dt/dy)*(v_ad(i-1,j+1,k )+v_ad(i-1,j,k))* + (yzlo(i-1,j+1,k )-yzlo(i-1,j,k)); + Real sth = xhi(i,j,k,n) - (Real(0.25)*l_dt/dy)*(v_ad(i ,j+1,k )+v_ad(i ,j,k))* + (yzlo(i ,j+1,k )-yzlo(i ,j,k)); if (!l_use_forces_in_trans) { - stl += 0.5 * l_dt * f(i-1,j,k,n); - sth += 0.5 * l_dt * f(i ,j,k,n); + stl += Real(0.5) * l_dt * f(i-1,j,k,n); + sth += Real(0.5) * l_dt * f(i ,j,k,n); } HydroBC::SetExtrapVelBCsLo(0, i, j, k, n, q, stl, sth, bc.lo(0), dlo.x); @@ -317,19 +317,19 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp, if (!allow_inflow_on_outflow) { if ( (i==dlo.x) && (bc.lo(0) == BCType::foextrap || bc.lo(0) == BCType::hoextrap) ) { - sth = amrex::min(sth,0.0_rt); + sth = amrex::min(sth,Real(0)); stl = sth; } if ( (i==dhi.x+1) && (bc.hi(0) == BCType::foextrap || bc.hi(0) == BCType::hoextrap) ) { - stl = amrex::max(stl,0.0_rt); + stl = amrex::max(stl,Real(0)); sth = stl; } } - Real st = ( (stl+sth) >= 0.) ? stl : sth; - bool ltm = ( (stl <= 0. && sth >= 0.) || (amrex::Math::abs(stl+sth) < small_vel) ); - qx(i,j,k) = ltm ? 0. : st; + Real st = ( (stl+sth) >= Real(0)) ? stl : sth; + bool ltm = ( (stl <= Real(0) && sth >= Real(0)) || (amrex::Math::abs(stl+sth) < small_vel) ); + qx(i,j,k) = ltm ? Real(0) : st; }); // @@ -352,8 +352,8 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp, Real uad = u_ad(i,j,k); Real st = (uad >= 0.) ? l_xzlo : l_xzhi; - Real fu = (amrex::Math::abs(uad) < small_vel) ? 0.0 : 1.0; - xzlo(i,j,k) = fu*st + (1.0 - fu) * 0.5 * (l_xzhi + l_xzlo); + Real fu = (amrex::Math::abs(uad) < small_vel) ? Real(0) : Real(1); + xzlo(i,j,k) = fu*st + (Real(1) - fu) * Real(0.5) * (l_xzhi + l_xzlo); }); amrex::ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept @@ -361,15 +361,15 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp, constexpr int n = 1; const auto bc = HydroBC::getBC(i, j, k, n, domain, pbc, bc_arr); - Real stl = ylo(i,j,k,n) - (0.25*l_dt/dx)*(u_ad(i+1,j-1,k )+u_ad(i,j-1,k))* - (xzlo(i+1,j-1,k )-xzlo(i,j-1,k)); - Real sth = yhi(i,j,k,n) - (0.25*l_dt/dx)*(u_ad(i+1,j ,k )+u_ad(i,j ,k))* - (xzlo(i+1,j ,k )-xzlo(i,j ,k)); + Real stl = ylo(i,j,k,n) - (Real(0.25)*l_dt/dx)*(u_ad(i+1,j-1,k )+u_ad(i,j-1,k))* + (xzlo(i+1,j-1,k )-xzlo(i,j-1,k)); + Real sth = yhi(i,j,k,n) - (Real(0.25)*l_dt/dx)*(u_ad(i+1,j ,k )+u_ad(i,j ,k))* + (xzlo(i+1,j ,k )-xzlo(i,j ,k)); if (!l_use_forces_in_trans) { - stl += 0.5 * l_dt * f(i,j-1,k,n); - sth += 0.5 * l_dt * f(i,j ,k,n); + stl += Real(0.5) * l_dt * f(i,j-1,k,n); + sth += Real(0.5) * l_dt * f(i,j ,k,n); } HydroBC::SetExtrapVelBCsLo(1, i, j, k, n, q, stl, sth, bc.lo(1), dlo.y); @@ -378,19 +378,19 @@ Godunov::ExtrapVelToFacesOnBox (Box const& bx, int ncomp, if (!allow_inflow_on_outflow) { if ( (j==dlo.y) && (bc.lo(1) == BCType::foextrap || bc.lo(1) == BCType::hoextrap) ) { - sth = amrex::min(sth,0.0_rt); + sth = amrex::min(sth,Real(0)); stl = sth; } if ( (j==dhi.y+1) && (bc.hi(1) == BCType::foextrap || bc.hi(1) == BCType::hoextrap) ) { - stl = amrex::max(stl,0.0_rt); + stl = amrex::max(stl,Real(0)); sth = stl; } } - Real st = ( (stl+sth) >= 0.) ? stl : sth; - bool ltm = ( (stl <= 0. && sth >= 0.) || (amrex::Math::abs(stl+sth) < small_vel) ); - qy(i,j,k) = ltm ? 0. : st; + Real st = ( (stl+sth) >= Real(0)) ? stl : sth; + bool ltm = ( (stl <= Real(0) && sth >= Real(0)) || (amrex::Math::abs(stl+sth) < small_vel) ); + qy(i,j,k) = ltm ? Real(0) : st; }); } /** @} */ diff --git a/Godunov/hydro_godunov_plm.cpp b/Godunov/hydro_godunov_plm.cpp index b33c67a02..cd8da1c29 100644 --- a/Godunov/hydro_godunov_plm.cpp +++ b/Godunov/hydro_godunov_plm.cpp @@ -68,9 +68,9 @@ PLM::PredictVelOnXFace ( Box const& xebox, int ncomp, int order = 4; - Real upls = q(i ,j,k,n) + 0.5 * (-1.0 - vcc(i ,j,k,0) * dtdx) * + Real upls = q(i ,j,k,n) + Real(0.5) * (-Real(1) - vcc(i ,j,k,0) * dtdx) * amrex_calc_xslope_extdir(i ,j,k,n,order,q,extdir_or_ho_ilo,extdir_or_ho_ihi,domain_ilo,domain_ihi); - Real umns = q(i-1,j,k,n) + 0.5 * ( 1.0 - vcc(i-1,j,k,0) * dtdx) * + Real umns = q(i-1,j,k,n) + Real(0.5) * ( Real(1) - vcc(i-1,j,k,0) * dtdx) * amrex_calc_xslope_extdir(i-1,j,k,n,order,q,extdir_or_ho_ilo,extdir_or_ho_ihi,domain_ilo,domain_ihi); Ipx(i-1,j,k,n) = umns; @@ -84,9 +84,9 @@ PLM::PredictVelOnXFace ( Box const& xebox, int ncomp, { int order = 4; - Real upls = q(i ,j,k,n) + 0.5 * (-1.0 - vcc(i ,j,k,0) * dtdx) * + Real upls = q(i ,j,k,n) + Real(0.5) * (-Real(1) - vcc(i ,j,k,0) * dtdx) * amrex_calc_xslope(i ,j,k,n,order,q); - Real umns = q(i-1,j,k,n) + 0.5 * ( 1.0 - vcc(i-1,j,k,0) * dtdx) * + Real umns = q(i-1,j,k,n) + Real(0.5) * ( Real(1) - vcc(i-1,j,k,0) * dtdx) * amrex_calc_xslope(i-1,j,k,n,order,q); Ipx(i-1,j,k,n) = umns; @@ -132,9 +132,9 @@ PLM::PredictVelOnYFace (Box const& yebox, int ncomp, int order = 4; - Real vpls = q(i,j ,k,n) + 0.5 * (-1.0 - vcc(i,j ,k,1) * dtdy) * + Real vpls = q(i,j ,k,n) + Real(0.5) * (-Real(1) - vcc(i,j ,k,1) * dtdy) * amrex_calc_yslope_extdir(i,j ,k,n,order,q,extdir_or_ho_jlo,extdir_or_ho_jhi,domain_jlo,domain_jhi); - Real vmns = q(i,j-1,k,n) + 0.5 * ( 1.0 - vcc(i,j-1,k,1) * dtdy) * + Real vmns = q(i,j-1,k,n) + Real(0.5) * ( Real(1) - vcc(i,j-1,k,1) * dtdy) * amrex_calc_yslope_extdir(i,j-1,k,n,order,q,extdir_or_ho_jlo,extdir_or_ho_jhi,domain_jlo,domain_jhi); Ipy(i,j-1,k,n) = vmns; @@ -148,9 +148,9 @@ PLM::PredictVelOnYFace (Box const& yebox, int ncomp, { int order = 4; - Real vpls = q(i,j ,k,n) + 0.5 * (-1.0 - vcc(i,j ,k,1) * dtdy) * + Real vpls = q(i,j ,k,n) + Real(0.5) * (-Real(1) - vcc(i,j ,k,1) * dtdy) * amrex_calc_yslope(i,j ,k,n,order,q); - Real vmns = q(i,j-1,k,n) + 0.5 * ( 1.0 - vcc(i,j-1,k,1) * dtdy) * + Real vmns = q(i,j-1,k,n) + Real(0.5) * ( Real(1) - vcc(i,j-1,k,1) * dtdy) * amrex_calc_yslope(i,j-1,k,n,order,q); Ipy(i,j-1,k,n) = vmns; @@ -197,9 +197,9 @@ PLM::PredictVelOnZFace ( Box const& zebox, int ncomp, int order = 4; - Real wpls = q(i,j,k ,n) + 0.5 * (-1.0 - vcc(i,j,k ,2) * dtdz) * + Real wpls = q(i,j,k ,n) + Real(0.5) * (-Real(1) - vcc(i,j,k ,2) * dtdz) * amrex_calc_zslope_extdir(i,j,k,n,order,q,extdir_or_ho_klo,extdir_or_ho_khi,domain_klo,domain_khi); - Real wmns = q(i,j,k-1,n) + 0.5 * ( 1.0 - vcc(i,j,k-1,2) * dtdz) * + Real wmns = q(i,j,k-1,n) + Real(0.5) * ( Real(1) - vcc(i,j,k-1,2) * dtdz) * amrex_calc_zslope_extdir(i,j,k-1,n,order,q,extdir_or_ho_klo,extdir_or_ho_khi,domain_klo,domain_khi); @@ -214,9 +214,9 @@ PLM::PredictVelOnZFace ( Box const& zebox, int ncomp, { int order = 4; - Real wpls = q(i,j,k ,n) + 0.5 * (-1.0 - vcc(i,j,k ,2) * dtdz) * + Real wpls = q(i,j,k ,n) + Real(0.5) * (-Real(1) - vcc(i,j,k ,2) * dtdz) * amrex_calc_zslope(i,j,k ,n,order,q); - Real wmns = q(i,j,k-1,n) + 0.5 * ( 1.0 - vcc(i,j,k-1,2) * dtdz) * + Real wmns = q(i,j,k-1,n) + Real(0.5) * ( Real(1) - vcc(i,j,k-1,2) * dtdz) * amrex_calc_zslope(i,j,k-1,n,order,q); Ipz(i,j,k-1,n) = wmns; diff --git a/Godunov/hydro_godunov_ppm.H b/Godunov/hydro_godunov_ppm.H index f955334c0..8f8eafe1c 100644 --- a/Godunov/hydro_godunov_ppm.H +++ b/Godunov/hydro_godunov_ppm.H @@ -101,7 +101,7 @@ struct nolimiter { amrex::Real d1 = amrex::Real(0.5) * (sp1 - sm1); amrex::Real d2 = amrex::Real(0.5) * (s0 - sm2); - amrex::Real sedge = amrex::Real(0.5)*(s0 + sm1) - amrex::Real(1./6.)*(d1 - d2); + amrex::Real sedge = amrex::Real(0.5)*(s0 + sm1) - (amrex::Real(1)/amrex::Real(6))*(d1 - d2); return sedge; } @@ -116,7 +116,7 @@ struct nolimiter { amrex::Real d1 = amrex::Real(0.5) * (sp2 - s0); amrex::Real d2 = amrex::Real(0.5) * (sp1 - sm1); - amrex::Real sedge = amrex::Real(0.5)*(sp1 + s0) - amrex::Real(1./6.)*(d1 - d2); + amrex::Real sedge = amrex::Real(0.5)*(sp1 + s0) - (amrex::Real(1)/amrex::Real(6))*(d1 - d2); return sedge; } @@ -159,7 +159,7 @@ struct vanleer { amrex::Real d1 = vanLeer(s0,sp1,sm1); amrex::Real d2 = vanLeer(sm1,s0,sm2); - amrex::Real sedge = amrex::Real(0.5)*(s0 + sm1) - amrex::Real(1./6.)*(d1 - d2); + amrex::Real sedge = amrex::Real(0.5)*(s0 + sm1) - (amrex::Real(1)/amrex::Real(6))*(d1 - d2); return amrex::min(amrex::max(sedge, amrex::min(s0, sm1)),amrex::max(s0,sm1)); } @@ -176,7 +176,7 @@ struct vanleer { amrex::Real d1 = vanLeer(sp1,sp2,s0); amrex::Real d2 = vanLeer(s0,sp1,sm1); - amrex::Real sedge = amrex::Real(0.5)*(sp1 + s0) - amrex::Real(1./6.)*(d1 - d2); + amrex::Real sedge = amrex::Real(0.5)*(sp1 + s0) - (amrex::Real(1)/amrex::Real(6))*(d1 - d2); return amrex::min(amrex::max(sedge, amrex::min(s0, sp1)),amrex::max(s0,sp1)); @@ -195,9 +195,9 @@ struct vanleer { sm_ = s0; sp_ = s0; } else if (amrex::Math::abs(sedge2-s0) >= amrex::Real(2.0)*amrex::Math::abs(sedge1-s0)) { - sp_ = amrex::Real(3.0)*s0 - amrex::Real(2.0)*sedge1; + sp_ = amrex::Real(3)*s0 - amrex::Real(2)*sedge1; } else if (amrex::Math::abs(sedge1-s0) >= amrex::Real(2.0)*amrex::Math::abs(sedge2-s0)) { - sm_ = amrex::Real(3.0)*s0 - amrex::Real(2.0)*sedge2; + sm_ = amrex::Real(3)*s0 - amrex::Real(2)*sedge2; } return amrex::makeTuple(sm_,sp_); @@ -209,10 +209,10 @@ struct vanleer { // amrex::Real sm_ = sedge1; // amrex::Real sp_ = sedge2; // if (amrex::Math::abs(sedge2-s0) >= amrex::Real(2.0)*amrex::Math::abs(sedge1-s0)) { -// sp_ = amrex::Real(3.0)*s0 - amrex::Real(2.0)*sedge1; +// sp_ = amrex::Real(3)*s0 - amrex::Real(2)*sedge1; // } // if (amrex::Math::abs(sedge1-s0) >= amrex::Real(2.0)*amrex::Math::abs(sedge2-s0)) { -// sm_ = amrex::Real(3.0)*s0 - amrex::Real(2.0)*sedge2; +// sm_ = amrex::Real(3)*s0 - amrex::Real(2)*sedge2; // } // return amrex::makeTuple(sm_,sp_); @@ -234,27 +234,27 @@ struct wenoz { constexpr auto eps = amrex::Real(1.0e-6); const amrex::Real beta1 = - amrex::Real(13.0) / amrex::Real(12.0) * (sm2 - amrex::Real(2.0) * sm1 + s0) * (sm2 - amrex::Real(2.0) * sm1 + s0) + - amrex::Real(0.25) * (sm2 - amrex::Real(4.0) * sm1 + amrex::Real(3.0) * s0) * (sm2 - amrex::Real(4.0) * sm1 + amrex::Real(3.0) * s0); + amrex::Real(13) / amrex::Real(12) * (sm2 - amrex::Real(2) * sm1 + s0) * (sm2 - amrex::Real(2) * sm1 + s0) + + amrex::Real(0.25) * (sm2 - amrex::Real(4) * sm1 + amrex::Real(3) * s0) * (sm2 - amrex::Real(4) * sm1 + amrex::Real(3) * s0); const amrex::Real beta2 = - amrex::Real(13.0) / amrex::Real(12.0) * (sm1 - amrex::Real(2.0) * s0 + sp1) * (sm1 - amrex::Real(2.0) * s0 + sp1) + + amrex::Real(13) / amrex::Real(12) * (sm1 - amrex::Real(2) * s0 + sp1) * (sm1 - amrex::Real(2) * s0 + sp1) + amrex::Real(0.25) * (sm1 - sp1) * (sm1 - sp1); const amrex::Real beta3 = - amrex::Real(13.0) / amrex::Real(12.0) * (s0 - amrex::Real(2.0) * sp1 + sp2) * (s0 - amrex::Real(2.0) * sp1 + sp2) + - amrex::Real(0.25) * (amrex::Real(3.0) * s0 - amrex::Real(4.0) * sp1 + sp2) * (amrex::Real(3.0) * s0 - amrex::Real(4.0) * sp1 + sp2); + amrex::Real(13) / amrex::Real(12) * (s0 - amrex::Real(2) * sp1 + sp2) * (s0 - amrex::Real(2) * sp1 + sp2) + + amrex::Real(0.25) * (amrex::Real(3) * s0 - amrex::Real(4) * sp1 + sp2) * (amrex::Real(3) * s0 - amrex::Real(4) * sp1 + sp2); const amrex::Real t5 = amrex::Math::abs(beta3 - beta1); - const amrex::Real omega1 = amrex::Real(0.1) * (amrex::Real(1.0) + t5 / (eps + beta1)); - const amrex::Real omega2 = amrex::Real(0.6) * (amrex::Real(1.0) + t5 / (eps + beta2)); - const amrex::Real omega3 = amrex::Real(0.3) * (amrex::Real(1.0) + t5 / (eps + beta3)); + const amrex::Real omega1 = amrex::Real(0.1) * (amrex::Real(1) + t5 / (eps + beta1)); + const amrex::Real omega2 = amrex::Real(0.6) * (amrex::Real(1) + t5 / (eps + beta2)); + const amrex::Real omega3 = amrex::Real(0.3) * (amrex::Real(1) + t5 / (eps + beta3)); const amrex::Real omega = omega1 + omega2 + omega3; - const amrex::Real v_1 = amrex::Real(2.0) * sm2 - amrex::Real(7.0) * sm1 + amrex::Real(11.0) * s0; - const amrex::Real v_2 = -sm1 + amrex::Real(5.0) * s0 + amrex::Real(2.0) * sp1; - const amrex::Real v_3 = amrex::Real(2.0) * s0 + amrex::Real(5.0) * sp1 - sp2; + const amrex::Real v_1 = amrex::Real(2) * sm2 - amrex::Real(7) * sm1 + amrex::Real(11) * s0; + const amrex::Real v_2 = -sm1 + amrex::Real(5) * s0 + amrex::Real(2) * sp1; + const amrex::Real v_3 = amrex::Real(2) * s0 + amrex::Real(5) * sp1 - sp2; - return (omega1 * v_1 + omega2 * v_2 + omega3 * v_3) / (amrex::Real(6.0) * omega); + return (omega1 * v_1 + omega2 * v_2 + omega3 * v_3) / (amrex::Real(6) * omega); } [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE @@ -293,19 +293,19 @@ struct weno_js { constexpr auto eps = amrex::Real(1.0e-6); const amrex::Real beta1 = - amrex::Real(13.0) / amrex::Real(12.0) * (sm2 - amrex::Real(2.0) * sm1 + s0) * - (sm2 - amrex::Real(2.0) * sm1 + s0) + - amrex::Real(0.25) * (sm2 - amrex::Real(4.0) * sm1 + amrex::Real(3.0) * s0) * - (sm2 - amrex::Real(4.0) * sm1 + amrex::Real(3.0) * s0); + amrex::Real(13) / amrex::Real(12) * (sm2 - amrex::Real(2) * sm1 + s0) * + (sm2 - amrex::Real(2) * sm1 + s0) + + amrex::Real(0.25) * (sm2 - amrex::Real(4) * sm1 + amrex::Real(3) * s0) * + (sm2 - amrex::Real(4) * sm1 + amrex::Real(3) * s0); const amrex::Real beta2 = - amrex::Real(13.0) / amrex::Real(12.0) * (sm1 - amrex::Real(2.0) * s0 + sp1) * - (sm1 - amrex::Real(2.0) * s0 + sp1) + + amrex::Real(13) / amrex::Real(12) * (sm1 - amrex::Real(2) * s0 + sp1) * + (sm1 - amrex::Real(2) * s0 + sp1) + amrex::Real(0.25) * (sm1 - sp1) * (sm1 - sp1); const amrex::Real beta3 = - amrex::Real(13.0) / amrex::Real(12.0) * (s0 - amrex::Real(2.0) * sp1 + sp2) * - (s0 - amrex::Real(2.0) * sp1 + sp2) + - amrex::Real(0.25) * (amrex::Real(3.0) * s0 - amrex::Real(4.0) * sp1 + sp2) * - (amrex::Real(3.0) * s0 - amrex::Real(4.0) * sp1 + sp2); + amrex::Real(13) / amrex::Real(12) * (s0 - amrex::Real(2) * sp1 + sp2) * + (s0 - amrex::Real(2) * sp1 + sp2) + + amrex::Real(0.25) * (amrex::Real(3) * s0 - amrex::Real(4) * sp1 + sp2) * + (amrex::Real(3) * s0 - amrex::Real(4) * sp1 + sp2); const amrex::Real omega1 = amrex::Real(0.1) / (eps + beta1); const amrex::Real omega2 = amrex::Real(0.6) / (eps + beta2); @@ -313,11 +313,11 @@ struct weno_js { const amrex::Real omega = omega1 + omega2 + omega3; - const amrex::Real v_1 = amrex::Real(2.0) * sm2 - amrex::Real(7.0) * sm1 + amrex::Real(11.0) * s0; - const amrex::Real v_2 = -sm1 + amrex::Real(5.0) * s0 + amrex::Real(2.0) * sp1; - const amrex::Real v_3 = amrex::Real(2.0) * s0 + amrex::Real(5.0) * sp1 - sp2; + const amrex::Real v_1 = amrex::Real(2) * sm2 - amrex::Real(7) * sm1 + amrex::Real(11) * s0; + const amrex::Real v_2 = -sm1 + amrex::Real(5) * s0 + amrex::Real(2) * sp1; + const amrex::Real v_3 = amrex::Real(2) * s0 + amrex::Real(5) * sp1 - sp2; - return (omega1 * v_1 + omega2 * v_2 + omega3 * v_3) / (amrex::Real(6.0) * omega); + return (omega1 * v_1 + omega2 * v_2 + omega3 * v_3) / (amrex::Real(6) * omega); } [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE @@ -426,7 +426,7 @@ void SetXBCs ( const int i, const int j, const int k, const int n, using namespace amrex; if ( (bclo == BCType::ext_dir) || (bclo == BCType::hoextrap) || - (bclo == BCType::direction_dependent && velp >= 0.0) ) + (bclo == BCType::direction_dependent && velp >= Real(0)) ) { if (i == domlo) { @@ -453,21 +453,21 @@ void SetXBCs ( const int i, const int j, const int k, const int n, sm = sedge1; if constexpr (Limiter::do_limiting) { - if ( (sp - s(domlo+1,j,k,n))*(s(domlo+1,j,k,n) - sm) <= amrex::Real(0.0)) + if ( (sp - s(domlo+1,j,k,n))*(s(domlo+1,j,k,n) - sm) <= amrex::Real(0)) { sp = s(domlo+1,j,k,n); sm = s(domlo+1,j,k,n); - } else if(amrex::Math::abs(sp - s(domlo+1,j,k,n)) >= amrex::Real(2.0)*amrex::Math::abs(sm - s(domlo+1,j,k,n))) { - sp = amrex::Real(3.0)*s(domlo+1,j,k,n) - amrex::Real(2.0)*sm; - } else if(amrex::Math::abs(sm - s(domlo+1,j,k,n)) >= amrex::Real(2.0)*amrex::Math::abs(sp - s(domlo+1,j,k,n))) { - sm = amrex::Real(3.0)*s(domlo+1,j,k,n) - amrex::Real(2.0)*sp; + } else if(amrex::Math::abs(sp - s(domlo+1,j,k,n)) >= amrex::Real(2)*amrex::Math::abs(sm - s(domlo+1,j,k,n))) { + sp = amrex::Real(3)*s(domlo+1,j,k,n) - amrex::Real(2)*sm; + } else if(amrex::Math::abs(sm - s(domlo+1,j,k,n)) >= amrex::Real(2)*amrex::Math::abs(sp - s(domlo+1,j,k,n))) { + sm = amrex::Real(3)*s(domlo+1,j,k,n) - amrex::Real(2)*sp; } } } } if ( (bchi == BCType::ext_dir) || (bchi == BCType::hoextrap) || - (bchi == BCType::direction_dependent && velm <= 0.0) ) + (bchi == BCType::direction_dependent && velm <= Real(0)) ) { if (i == domhi) { @@ -494,14 +494,14 @@ void SetXBCs ( const int i, const int j, const int k, const int n, sm = sedge1; if constexpr (Limiter::do_limiting) { - if( (sp - s(domhi-1,j,k,n))*(s(domhi-1,j,k,n) - sm) <= amrex::Real(0.0)) + if( (sp - s(domhi-1,j,k,n))*(s(domhi-1,j,k,n) - sm) <= amrex::Real(0)) { sp = s(domhi-1,j,k,n); sm = s(domhi-1,j,k,n); } else if(amrex::Math::abs(sp - s(domhi-1,j,k,n)) >= 2.*amrex::Math::abs(sm - s(domhi-1,j,k,n))) { - sp = amrex::Real(3.0)*s(domhi-1,j,k,n) - amrex::Real(2.0)*sm; + sp = amrex::Real(3)*s(domhi-1,j,k,n) - amrex::Real(2)*sm; } else if(amrex::Math::abs(sm - s(domhi-1,j,k,n)) >= 2.*amrex::Math::abs(sp - s(domhi-1,j,k,n))) { - sm = amrex::Real(3.0)*s(domhi-1,j,k,n) - amrex::Real(2.0)*sp; + sm = amrex::Real(3)*s(domhi-1,j,k,n) - amrex::Real(2)*sp; } } } @@ -521,7 +521,7 @@ void SetYBCs ( const int i, const int j, const int k, const int n, using namespace amrex; if ( (bclo == BCType::ext_dir) || (bclo == BCType::hoextrap) || - (bclo == BCType::direction_dependent && velm >= 0.0) ) + (bclo == BCType::direction_dependent && velm >= Real(0)) ) { if (j == domlo) { @@ -548,21 +548,21 @@ void SetYBCs ( const int i, const int j, const int k, const int n, sm = sedge1; if constexpr (Limiter::do_limiting) { - if ( (sp - s(i,domlo+1,k,n))*(s(i,domlo+1,k,n) - sm) <= amrex::Real(0.0)) + if ( (sp - s(i,domlo+1,k,n))*(s(i,domlo+1,k,n) - sm) <= amrex::Real(0)) { sp = s(i,domlo+1,k,n); sm = s(i,domlo+1,k,n); - } else if(amrex::Math::abs(sp - s(i,domlo+1,k,n)) >= amrex::Real(2.0)*amrex::Math::abs(sm - s(i,domlo+1,k,n))) { - sp = amrex::Real(3.0)*s(i,domlo+1,k,n) - amrex::Real(2.0)*sm; - } else if(amrex::Math::abs(sm - s(i,domlo+1,k,n)) >= amrex::Real(2.0)*amrex::Math::abs(sp - s(i,domlo+1,k,n))) { - sm = amrex::Real(3.0)*s(i,domlo+1,k,n) - amrex::Real(2.0)*sp; + } else if(amrex::Math::abs(sp - s(i,domlo+1,k,n)) >= amrex::Real(2)*amrex::Math::abs(sm - s(i,domlo+1,k,n))) { + sp = amrex::Real(3)*s(i,domlo+1,k,n) - amrex::Real(2)*sm; + } else if(amrex::Math::abs(sm - s(i,domlo+1,k,n)) >= amrex::Real(2)*amrex::Math::abs(sp - s(i,domlo+1,k,n))) { + sm = amrex::Real(3)*s(i,domlo+1,k,n) - amrex::Real(2)*sp; } } } } if ( (bchi == BCType::ext_dir) || (bchi == BCType::hoextrap) || - (bchi == BCType::direction_dependent && velp <= 0.0) ) + (bchi == BCType::direction_dependent && velp <= Real(0)) ) { if (j == domhi) { @@ -589,13 +589,13 @@ void SetYBCs ( const int i, const int j, const int k, const int n, sm = sedge1; if constexpr (Limiter::do_limiting) { - if( (sp - s(i,domhi-1,k,n))*(s(i,domhi-1,k,n) - sm) <= amrex::Real(0.0)){ + if( (sp - s(i,domhi-1,k,n))*(s(i,domhi-1,k,n) - sm) <= amrex::Real(0)){ sp = s(i,domhi-1,k,n); sm = s(i,domhi-1,k,n); } else if(amrex::Math::abs(sp - s(i,domhi-1,k,n)) >= 2.*amrex::Math::abs(sm - s(i,domhi-1,k,n))) { - sp = amrex::Real(3.0)*s(i,domhi-1,k,n) - amrex::Real(2.0)*sm; + sp = amrex::Real(3)*s(i,domhi-1,k,n) - amrex::Real(2)*sm; } else if(amrex::Math::abs(sm - s(i,domhi-1,k,n)) >= 2.*amrex::Math::abs(sp - s(i,domhi-1,k,n))) { - sm = amrex::Real(3.0)*s(i,domhi-1,k,n) - amrex::Real(2.0)*sp; + sm = amrex::Real(3)*s(i,domhi-1,k,n) - amrex::Real(2)*sp; } } } @@ -616,7 +616,7 @@ void SetZBCs ( const int i, const int j, const int k, const int n, using namespace amrex; if ( (bclo == BCType::ext_dir) || (bclo == BCType::hoextrap) || - (bclo == BCType::direction_dependent && velm >= 0.0) ) + (bclo == BCType::direction_dependent && velm >= Real(0)) ) { if (k == domlo) { @@ -648,16 +648,16 @@ void SetZBCs ( const int i, const int j, const int k, const int n, sp = s(i,j,domlo+1,n); sm = s(i,j,domlo+1,n); } else if(amrex::Math::abs(sp - s(i,j,domlo+1,n)) >= 2.*amrex::Math::abs(sm - s(i,j,domlo+1,n))) { - sp = amrex::Real(3.0)*s(i,j,domlo+1,n) - amrex::Real(2.0)*sm; + sp = amrex::Real(3)*s(i,j,domlo+1,n) - amrex::Real(2)*sm; } else if(amrex::Math::abs(sm - s(i,j,domlo+1,n)) >= 2.*amrex::Math::abs(sp - s(i,j,domlo+1,n))) { - sm = amrex::Real(3.0)*s(i,j,domlo+1,n) - amrex::Real(2.0)*sp; + sm = amrex::Real(3)*s(i,j,domlo+1,n) - amrex::Real(2)*sp; } } } } if ( (bchi == BCType::ext_dir) || (bchi == BCType::hoextrap) || - (bchi == BCType::direction_dependent && velp <= 0.0) ) + (bchi == BCType::direction_dependent && velp <= Real(0)) ) { if (k == domhi) { @@ -689,9 +689,9 @@ void SetZBCs ( const int i, const int j, const int k, const int n, sp = s(i,j,domhi-1,n); sm = s(i,j,domhi-1,n); } else if(amrex::Math::abs(sp - s(i,j,domhi-1,n)) >= 2.*amrex::Math::abs(sm - s(i,j,domhi-1,n))) { - sp = amrex::Real(3.0)*s(i,j,domhi-1,n) - amrex::Real(2.0)*sm; + sp = amrex::Real(3)*s(i,j,domhi-1,n) - amrex::Real(2)*sm; } else if(amrex::Math::abs(sm - s(i,j,domhi-1,n)) >= 2.*amrex::Math::abs(sp - s(i,j,domhi-1,n))) { - sm = amrex::Real(3.0)*s(i,j,domhi-1,n) - amrex::Real(2.0)*sp; + sm = amrex::Real(3)*s(i,j,domhi-1,n) - amrex::Real(2)*sp; } } @@ -715,8 +715,8 @@ void PredictVelOnXFace ( const int i, const int j, const int k, const int n, const Limiter& /*limiter*/) { constexpr auto half{amrex::Real(0.5)}; - constexpr auto one{amrex::Real(1.0)}; - constexpr auto two3rds{amrex::Real(2.0/3.0)}; + constexpr auto one{amrex::Real(1)}; + constexpr auto two3rds{amrex::Real(2.0)/amrex::Real(3.0)}; const amrex::Real sm2 = S(i-2,j,k,n); const amrex::Real sm1 = S(i-1,j,k,n); @@ -732,7 +732,7 @@ void PredictVelOnXFace ( const int i, const int j, const int k, const int n, SetXBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, v_ad, v_ad, bc.lo(0), bc.hi(0), domlo, domhi); - const amrex::Real s6 = amrex::Real(6.0)*s0 - amrex::Real(3.0)*(sm + sp); + const amrex::Real s6 = amrex::Real(6)*s0 - amrex::Real(3)*(sm + sp); const amrex::Real sigma = amrex::Math::abs(v_ad)*dtdx; @@ -763,8 +763,8 @@ void PredictVelOnYFace ( const int i, const int j, const int k, const int n, const Limiter& /*limiter*/) { constexpr auto half{amrex::Real(0.5)}; - constexpr auto one{amrex::Real(1.0)}; - constexpr auto two3rds{amrex::Real(2.0/3.0)}; + constexpr auto one{amrex::Real(1)}; + constexpr auto two3rds{amrex::Real(2.0)/amrex::Real(3.0)}; const amrex::Real sm2 = S(i,j-2,k,n); const amrex::Real sm1 = S(i,j-1,k,n); @@ -780,7 +780,7 @@ void PredictVelOnYFace ( const int i, const int j, const int k, const int n, SetYBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, v_ad, v_ad, bc.lo(1), bc.hi(1), domlo, domhi); - const amrex::Real s6 = amrex::Real(6.0)*s0 - amrex::Real(3.0)*(sm + sp); + const amrex::Real s6 = amrex::Real(6)*s0 - amrex::Real(3)*(sm + sp); const amrex::Real sigma = amrex::Math::abs(v_ad)*dtdy; @@ -812,8 +812,8 @@ void PredictVelOnZFace ( const int i, const int j, const int k, const int n, const Limiter& /*limiter*/) { constexpr auto half{amrex::Real(0.5)}; - constexpr auto one{amrex::Real(1.0)}; - constexpr auto two3rds{amrex::Real(2.0/3.0)}; + constexpr auto one{amrex::Real(1)}; + constexpr auto two3rds{amrex::Real(2.0)/amrex::Real(3.0)}; const amrex::Real sm2 = S(i,j,k-2,n); const amrex::Real sm1 = S(i,j,k-1,n); @@ -829,7 +829,7 @@ void PredictVelOnZFace ( const int i, const int j, const int k, const int n, SetZBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, v_ad, v_ad, bc.lo(2), bc.hi(2), domlo, domhi); - const amrex::Real s6 = amrex::Real(6.0)*s0 - amrex::Real(3.0)*(sm + sp); + const amrex::Real s6 = amrex::Real(6)*s0 - amrex::Real(3)*(sm + sp); const amrex::Real sigma = amrex::Math::abs(v_ad)*dtdz; @@ -871,12 +871,12 @@ void PredictStateOnXFace ( const int i, const int j, const int k, const int n, const amrex::Real sp1 = S(i+1,j,k,n); constexpr auto half{amrex::Real(0.5)}; - constexpr auto one{amrex::Real(1.0)}; + constexpr auto one{amrex::Real(1)}; const amrex::Real sigmap = amrex::Math::abs(vel_edge(i+1,j,k))*dt/dx; const amrex::Real sigmam = amrex::Math::abs(vel_edge(i ,j,k))*dt/dx; - constexpr auto two3rds{amrex::Real(2.0/3.0)}; + constexpr auto two3rds{amrex::Real(2.0)/amrex::Real(3.0)}; const amrex::Real sm2 = S(i-2,j,k,n); const amrex::Real sp2 = S(i+2,j,k,n); @@ -889,7 +889,7 @@ void PredictStateOnXFace ( const int i, const int j, const int k, const int n, SetXBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k), bc.lo(0), bc.hi(0), domlo, domhi); - const amrex::Real s6 = 6.0*s0 - 3.0*(sm + sp); + const amrex::Real s6 = amrex::Real(6)*s0 - amrex::Real(3)*(sm + sp); if (vel_edge(i+1,j,k) > small_vel) { Ip = sp - (half*sigmap)*((sp - sm) - (one -two3rds*sigmap)*s6); @@ -984,7 +984,7 @@ void PredictStateOnYFace ( const int i, const int j, const int k, const int n, const amrex::Real sigmap = amrex::Math::abs(vel_edge(i,j+1,k))*dt/dx; const amrex::Real sigmam = amrex::Math::abs(vel_edge(i,j ,k))*dt/dx; - constexpr auto two3rds{amrex::Real(2.0/3.0)}; + constexpr auto two3rds{amrex::Real(2.0)/amrex::Real(3.0)}; const amrex::Real sm2 = S(i,j-2,k,n); const amrex::Real sp2 = S(i,j+2,k,n); @@ -997,7 +997,7 @@ void PredictStateOnYFace ( const int i, const int j, const int k, const int n, SetYBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k), bc.lo(1), bc.hi(1), domlo, domhi); - const amrex::Real s6 = 6.0*s0- 3.0*(sm + sp); + const amrex::Real s6 = amrex::Real(6)*s0- amrex::Real(3)*(sm + sp); if (vel_edge(i,j+1,k) > small_vel) { Ip = sp - (half*sigmap)*((sp - sm) - (one -two3rds*sigmap)*s6); @@ -1093,7 +1093,7 @@ void PredictStateOnZFace ( const int i, const int j, const int k, const int n, const amrex::Real sigmap = amrex::Math::abs(vel_edge(i,j,k+1))*dt/dx; const amrex::Real sigmam = amrex::Math::abs(vel_edge(i,j,k ))*dt/dx; - constexpr auto two3rds{amrex::Real(2.0/3.0)}; + constexpr auto two3rds{amrex::Real(2.0)/amrex::Real(3.0)}; const amrex::Real sm2 = S(i,j,k-2,n); const amrex::Real sp2 = S(i,j,k+2,n); @@ -1106,7 +1106,7 @@ void PredictStateOnZFace ( const int i, const int j, const int k, const int n, SetZBCs(i, j, k, n, sm, sp, sedge1, sedge2, S, vel_edge(i,j,k), vel_edge(i,j,k), bc.lo(2), bc.hi(2), domlo, domhi); - const amrex::Real s6 = 6.0*s0- 3.0*(sm + sp); + const amrex::Real s6 = amrex::Real(6)*s0- amrex::Real(3)*(sm + sp); if(vel_edge(i,j,k+1) > small_vel) { Ip = sp - (half*sigmap)*((sp-sm) - (amrex::Real(1.0) -two3rds*sigmap)*s6); diff --git a/MOL/hydro_mol_edge_state_K.H b/MOL/hydro_mol_edge_state_K.H index c96442cef..73f6bfa59 100644 --- a/MOL/hydro_mol_edge_state_K.H +++ b/MOL/hydro_mol_edge_state_K.H @@ -8,8 +8,6 @@ #include #include -using namespace amrex::literals; - namespace MOL { AMREX_GPU_DEVICE AMREX_FORCE_INLINE @@ -40,9 +38,9 @@ amrex::Real hydro_mol_xedge_state_extdir ( int i, int j, int k, int n, } else { - amrex::Real qpls = q(i ,j,k,n) - 0.5 * + amrex::Real qpls = q(i ,j,k,n) - amrex::Real(0.5) * amrex_calc_xslope_extdir( i , j, k, n, order, q, extdir_or_ho_lo, extdir_or_ho_hi, domlo, domhi ); - amrex::Real qmns = q(i-1,j,k,n) + 0.5 * + amrex::Real qmns = q(i-1,j,k,n) + amrex::Real(0.5) * amrex_calc_xslope_extdir( i-1, j, k, n, order, q, extdir_or_ho_lo, extdir_or_ho_hi, domlo, domhi ); HydroBC::SetEdgeBCsLo(0,i,j,k,n,q,qmns,qpls,umac(i,j,k),d_bcrec[n].lo(0),domlo,is_velocity); @@ -50,12 +48,12 @@ amrex::Real hydro_mol_xedge_state_extdir ( int i, int j, int k, int n, if ( (i==domlo) && (d_bcrec[n].lo(0) == amrex::BCType::foextrap || d_bcrec[n].lo(0) == amrex::BCType::hoextrap) ) { - if ( umac(i,j,k) >= 0. && n==XVEL && is_velocity ) qpls = amrex::min(qpls,0.0_rt); + if ( umac(i,j,k) >= 0. && n==XVEL && is_velocity ) qpls = amrex::min(qpls,amrex::Real(0)); qmns = qpls; } if ( (i==domhi+1) && (d_bcrec[n].hi(0) == amrex::BCType::foextrap || d_bcrec[n].hi(0) == amrex::BCType::hoextrap) ) { - if ( umac(i,j,k) <= 0. && n==XVEL && is_velocity ) qmns = amrex::max(qmns,0.0_rt); + if ( umac(i,j,k) <= 0. && n==XVEL && is_velocity ) qmns = amrex::max(qmns,amrex::Real(0)); qpls = qmns; } @@ -69,7 +67,7 @@ amrex::Real hydro_mol_xedge_state_extdir ( int i, int j, int k, int n, } else { - qs = 0.5*(qmns+qpls); + qs = amrex::Real(0.5)*(qmns+qpls); } } @@ -87,20 +85,20 @@ amrex::Real hydro_mol_xedge_state ( int i, int j, int k, int n, int order = 2; amrex::Real qs; - amrex::Real qpls = q(i ,j,k,n) - 0.5 * amrex_calc_xslope( i , j, k, n, order, q ); - amrex::Real qmns = q(i-1,j,k,n) + 0.5 * amrex_calc_xslope( i-1, j, k, n, order, q ); + amrex::Real qpls = q(i ,j,k,n) - amrex::Real(0.5) * amrex_calc_xslope( i , j, k, n, order, q ); + amrex::Real qmns = q(i-1,j,k,n) + amrex::Real(0.5) * amrex_calc_xslope( i-1, j, k, n, order, q ); HydroBC::SetEdgeBCsLo(0,i,j,k,n,q,qmns,qpls,umac(i,j,k),d_bcrec[n].lo(0),domlo,is_velocity); HydroBC::SetEdgeBCsHi(0,i,j,k,n,q,qmns,qpls,umac(i,j,k),d_bcrec[n].hi(0),domhi,is_velocity); if ( (i==domlo) && (d_bcrec[n].lo(0) == amrex::BCType::foextrap || d_bcrec[n].lo(0) == amrex::BCType::hoextrap) ) { - if ( umac(i,j,k) >= 0. && n==XVEL && is_velocity ) qpls = amrex::min(qpls,0.0_rt); + if ( umac(i,j,k) >= 0. && n==XVEL && is_velocity ) qpls = amrex::min(qpls,amrex::Real(0)); qmns = qpls; } if ( (i==domhi+1) && (d_bcrec[n].hi(0) == amrex::BCType::foextrap || d_bcrec[n].hi(0) == amrex::BCType::hoextrap) ) { - if ( umac(i,j,k) <= 0. && n==XVEL && is_velocity ) qmns = amrex::max(qmns,0.0_rt); + if ( umac(i,j,k) <= 0. && n==XVEL && is_velocity ) qmns = amrex::max(qmns,amrex::Real(0)); qpls = qmns; } @@ -114,7 +112,7 @@ amrex::Real hydro_mol_xedge_state ( int i, int j, int k, int n, } else { - qs = 0.5*(qmns+qpls); + qs = amrex::Real(0.5)*(qmns+qpls); } return qs; @@ -149,9 +147,9 @@ amrex::Real hydro_mol_yedge_state_extdir ( int i, int j, int k, int n, } else { - amrex::Real qpls = q(i,j ,k,n) - 0.5 * + amrex::Real qpls = q(i,j ,k,n) - amrex::Real(0.5) * amrex_calc_yslope_extdir( i, j , k, n, order, q, extdir_or_ho_lo, extdir_or_ho_hi, domlo, domhi ); - amrex::Real qmns = q(i,j-1,k,n) + 0.5 * + amrex::Real qmns = q(i,j-1,k,n) + amrex::Real(0.5) * amrex_calc_yslope_extdir( i, j-1, k, n, order, q, extdir_or_ho_lo, extdir_or_ho_hi, domlo, domhi ); HydroBC::SetEdgeBCsLo(1,i,j,k,n,q,qmns,qpls,vmac(i,j,k),d_bcrec[n].lo(1),domlo,is_velocity); @@ -159,12 +157,12 @@ amrex::Real hydro_mol_yedge_state_extdir ( int i, int j, int k, int n, if ( (j==domlo) && (d_bcrec[n].lo(1) == amrex::BCType::foextrap || d_bcrec[n].lo(1) == amrex::BCType::hoextrap) ) { - if ( vmac(i,j,k) >= 0. && n==YVEL && is_velocity ) qpls = amrex::min(qpls,0.0_rt); + if ( vmac(i,j,k) >= 0. && n==YVEL && is_velocity ) qpls = amrex::min(qpls,amrex::Real(0)); qmns = qpls; } if ( (j==domhi+1) && (d_bcrec[n].hi(1) == amrex::BCType::foextrap || d_bcrec[n].hi(1) == amrex::BCType::hoextrap) ) { - if ( vmac(i,j,k) <= 0. && n==YVEL && is_velocity ) qmns = amrex::max(qmns,0.0_rt); + if ( vmac(i,j,k) <= 0. && n==YVEL && is_velocity ) qmns = amrex::max(qmns,amrex::Real(0)); qpls = qmns; } @@ -178,7 +176,7 @@ amrex::Real hydro_mol_yedge_state_extdir ( int i, int j, int k, int n, } else { - qs = 0.5*(qmns+qpls); + qs = amrex::Real(0.5)*(qmns+qpls); } } @@ -198,20 +196,20 @@ amrex::Real hydro_mol_yedge_state ( int i, int j, int k, int n, int order = 2; amrex::Real qs; - amrex::Real qpls = q(i,j ,k,n) - 0.5 * amrex_calc_yslope( i, j , k, n, order, q ); - amrex::Real qmns = q(i,j-1,k,n) + 0.5 * amrex_calc_yslope( i, j-1, k, n, order, q ); + amrex::Real qpls = q(i,j ,k,n) - amrex::Real(0.5) * amrex_calc_yslope( i, j , k, n, order, q ); + amrex::Real qmns = q(i,j-1,k,n) + amrex::Real(0.5) * amrex_calc_yslope( i, j-1, k, n, order, q ); HydroBC::SetEdgeBCsLo(1,i,j,k,n,q,qmns,qpls,vmac(i,j,k), d_bcrec[n].lo(1),domlo,is_velocity); HydroBC::SetEdgeBCsHi(1,i,j,k,n,q,qmns,qpls,vmac(i,j,k), d_bcrec[n].hi(1),domhi,is_velocity); if ( (j==domlo) && (d_bcrec[n].lo(1) == amrex::BCType::foextrap || d_bcrec[n].lo(1) == amrex::BCType::hoextrap) ) { - if ( vmac(i,j,k) >= 0. && n==YVEL && is_velocity ) qpls = amrex::min(qpls,0.0_rt); + if ( vmac(i,j,k) >= 0. && n==YVEL && is_velocity ) qpls = amrex::min(qpls,amrex::Real(0)); qmns = qpls; } if ( (j==domhi+1) && (d_bcrec[n].hi(1) == amrex::BCType::foextrap || d_bcrec[n].hi(1) == amrex::BCType::hoextrap) ) { - if ( vmac(i,j,k) <= 0. && n==YVEL && is_velocity ) qmns = amrex::max(qmns,0.0_rt); + if ( vmac(i,j,k) <= 0. && n==YVEL && is_velocity ) qmns = amrex::max(qmns,amrex::Real(0)); qpls = qmns; } @@ -225,7 +223,7 @@ amrex::Real hydro_mol_yedge_state ( int i, int j, int k, int n, } else { - qs = 0.5*(qmns+qpls); + qs = amrex::Real(0.5)*(qmns+qpls); } return qs; @@ -262,9 +260,9 @@ amrex::Real hydro_mol_zedge_state_extdir ( int i, int j, int k, int n, } else { - amrex::Real qpls = q(i,j,k ,n) - 0.5 * + amrex::Real qpls = q(i,j,k ,n) - amrex::Real(0.5) * amrex_calc_zslope_extdir( i, j, k , n, order, q, extdir_or_ho_lo, extdir_or_ho_hi, domlo, domhi ); - amrex::Real qmns = q(i,j,k-1,n) + 0.5 * + amrex::Real qmns = q(i,j,k-1,n) + amrex::Real(0.5) * amrex_calc_zslope_extdir( i, j, k-1, n, order, q, extdir_or_ho_lo, extdir_or_ho_hi, domlo, domhi ); HydroBC::SetEdgeBCsLo(2,i,j,k,n,q,qmns,qpls,wmac(i,j,k),d_bcrec[n].lo(2),domlo,is_velocity); @@ -272,12 +270,12 @@ amrex::Real hydro_mol_zedge_state_extdir ( int i, int j, int k, int n, if ( (k==domlo) && (d_bcrec[n].lo(2) == amrex::BCType::foextrap || d_bcrec[n].lo(2) == amrex::BCType::hoextrap) ) { - if ( wmac(i,j,k) >= 0. && n==ZVEL && is_velocity ) qpls = amrex::min(qpls,0.0_rt); + if ( wmac(i,j,k) >= 0. && n==ZVEL && is_velocity ) qpls = amrex::min(qpls,amrex::Real(0)); qmns = qpls; } if ( (k==domhi+1) && (d_bcrec[n].hi(2) == amrex::BCType::foextrap || d_bcrec[n].hi(2) == amrex::BCType::hoextrap) ) { - if ( wmac(i,j,k) <= 0. && n==ZVEL && is_velocity ) qmns = amrex::max(qmns,0.0_rt); + if ( wmac(i,j,k) <= 0. && n==ZVEL && is_velocity ) qmns = amrex::max(qmns,amrex::Real(0)); qpls = qmns; } @@ -291,7 +289,7 @@ amrex::Real hydro_mol_zedge_state_extdir ( int i, int j, int k, int n, } else { - qs = 0.5*(qmns+qpls); + qs = amrex::Real(0.5)*(qmns+qpls); } } @@ -311,20 +309,20 @@ amrex::Real hydro_mol_zedge_state ( int i, int j, int k, int n, int order = 2; amrex::Real qs; - amrex::Real qpls = q(i,j,k ,n) - 0.5 * amrex_calc_zslope( i, j, k , n, order, q ); - amrex::Real qmns = q(i,j,k-1,n) + 0.5 * amrex_calc_zslope( i, j, k-1, n, order, q ); + amrex::Real qpls = q(i,j,k ,n) - amrex::Real(0.5) * amrex_calc_zslope( i, j, k , n, order, q ); + amrex::Real qmns = q(i,j,k-1,n) + amrex::Real(0.5) * amrex_calc_zslope( i, j, k-1, n, order, q ); HydroBC::SetEdgeBCsLo(2,i,j,k,n,q,qmns,qpls,wmac(i,j,k),d_bcrec[n].lo(2),domlo,is_velocity); HydroBC::SetEdgeBCsHi(2,i,j,k,n,q,qmns,qpls,wmac(i,j,k),d_bcrec[n].hi(2),domhi,is_velocity); if ( (k==domlo) && (d_bcrec[n].lo(2) == amrex::BCType::foextrap || d_bcrec[n].lo(2) == amrex::BCType::hoextrap) ) { - if ( wmac(i,j,k) >= 0. && n==ZVEL && is_velocity ) qpls = amrex::min(qpls,0.0_rt); + if ( wmac(i,j,k) >= 0. && n==ZVEL && is_velocity ) qpls = amrex::min(qpls,amrex::Real(0)); qmns = qpls; } if ( (k==domhi+1) && (d_bcrec[n].hi(2) == amrex::BCType::foextrap || d_bcrec[n].hi(2) == amrex::BCType::hoextrap) ) { - if ( wmac(i,j,k) <= 0. && n==ZVEL && is_velocity ) qmns = amrex::max(qmns,0.0_rt); + if ( wmac(i,j,k) <= 0. && n==ZVEL && is_velocity ) qmns = amrex::max(qmns,amrex::Real(0)); qpls = qmns; } @@ -338,7 +336,7 @@ amrex::Real hydro_mol_zedge_state ( int i, int j, int k, int n, } else { - qs = 0.5*(qmns+qpls); + qs = amrex::Real(0.5)*(qmns+qpls); } return qs; diff --git a/Projections/hydro_NodalProjector.cpp b/Projections/hydro_NodalProjector.cpp index 49bfbcd86..d62417028 100644 --- a/Projections/hydro_NodalProjector.cpp +++ b/Projections/hydro_NodalProjector.cpp @@ -268,7 +268,7 @@ NodalProjector::project ( Real a_rtol, Real a_atol ) AMREX_ALWAYS_ASSERT(!m_need_bcs); if (m_verbose > 0) - amrex::Print() << "Nodal Projection:" << std::endl; + amrex::Print() << "Nodal Projection:" << '\n'; // // Average fine grid velocity values down to the coarse grid @@ -295,9 +295,9 @@ NodalProjector::project ( Real a_rtol, Real a_atol ) // Print diagnostics if (m_verbose > 0) { - amrex::Print() << " >> Before projection:" << std::endl; + amrex::Print() << " >> Before projection:" << '\n'; printInfo(); - amrex::Print() << std::endl; + amrex::Print() << '\n'; } // Solve @@ -343,9 +343,9 @@ NodalProjector::project ( Real a_rtol, Real a_atol ) if ( (m_verbose > 0) && (!m_has_rhs)) { computeRHS( GetVecOfPtrs(m_rhs), m_vel, m_S_cc, m_S_nd ); - amrex::Print() << " >> After projection:" << std::endl; + amrex::Print() << " >> After projection:" << '\n'; printInfo(); - amrex::Print() << std::endl; + amrex::Print() << '\n'; } } @@ -447,7 +447,7 @@ NodalProjector::printInfo () amrex::Print() << " * On lev " << lev << " max(abs(rhs)) = " << m_rhs[lev].norm0(0,0,false,true) - << std::endl; + << '\n'; } } diff --git a/Redistribution/FOR_PAPER/hydro_redistribution.cpp b/Redistribution/FOR_PAPER/hydro_redistribution.cpp index 1a3873eaf..81beda750 100644 --- a/Redistribution/FOR_PAPER/hydro_redistribution.cpp +++ b/Redistribution/FOR_PAPER/hydro_redistribution.cpp @@ -400,7 +400,7 @@ Redistribution::Make1DProfile ( Box const& bx, int ncomp, for (int i = 0; i < 64; i++) { - amrex::Print() << xloc[i] << " " << prof_ptr[i] << std::endl; + amrex::Print() << xloc[i] << " " << prof_ptr[i] << '\n'; } } diff --git a/Tests/MAC_Projection_EB/main.cpp b/Tests/MAC_Projection_EB/main.cpp index 4990cae06..c41054887 100644 --- a/Tests/MAC_Projection_EB/main.cpp +++ b/Tests/MAC_Projection_EB/main.cpp @@ -19,7 +19,7 @@ void write_plotfile(const Geometry& geom, const MultiFab& plotmf, int regtest) sstream << "plt00000"; std::string plotfile_name = sstream.str(); - amrex::Print() << "Writing " << plotfile_name << std::endl; + amrex::Print() << "Writing " << plotfile_name << '\n'; #if (AMREX_SPACEDIM == 2) EB_WriteSingleLevelPlotfile(plotfile_name, plotmf, @@ -235,12 +235,12 @@ int main (int argc, char* argv[]) macproj.getMLMG().setBottomVerbose(bottom_verbose); - amrex::Print() << " \n********************************************************************" << std::endl; - amrex::Print() << " Let's project the initial velocity to find " << std::endl; - amrex::Print() << " the flow field around the obstacles ... " << std::endl; - amrex::Print() << " The domain has " << n_cell_x << " cells in the x-direction " << std::endl; - amrex::Print() << " The maximum grid size is " << max_grid_size << std::endl; - amrex::Print() << "******************************************************************** \n" << std::endl; + amrex::Print() << " \n********************************************************************" << '\n'; + amrex::Print() << " Let's project the initial velocity to find " << '\n'; + amrex::Print() << " the flow field around the obstacles ... " << '\n'; + amrex::Print() << " The domain has " << n_cell_x << " cells in the x-direction " << '\n'; + amrex::Print() << " The maximum grid size is " << max_grid_size << '\n'; + amrex::Print() << "******************************************************************** \n" << '\n'; // Solve for phi and subtract from the velocity to make it divergence-free // Note that the normal velocities are at face centers (not centroids) @@ -254,9 +254,9 @@ int main (int argc, char* argv[]) // phi_inout.setVal(0.); // macproj.project_center_vels({&phi_inout},reltol,abstol,MLMG::Location::FaceCenter); - amrex::Print() << " \n********************************************************************" << std::endl; - amrex::Print() << " Done!" << std::endl; - amrex::Print() << "******************************************************************** \n" << std::endl; + amrex::Print() << " \n********************************************************************" << '\n'; + amrex::Print() << " Done!" << '\n'; + amrex::Print() << "******************************************************************** \n" << '\n'; for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { vel[idim].FillBoundary(geom.periodicity()); @@ -272,7 +272,7 @@ int main (int argc, char* argv[]) } auto stop_time = amrex::second() - strt_time; - amrex::Print() << "Total run time " << stop_time << std::endl; + amrex::Print() << "Total run time " << stop_time << '\n'; amrex::Finalize(); } diff --git a/Tests/Nodal_Projection_EB/main.cpp b/Tests/Nodal_Projection_EB/main.cpp index 95673682d..b4f6403e7 100644 --- a/Tests/Nodal_Projection_EB/main.cpp +++ b/Tests/Nodal_Projection_EB/main.cpp @@ -16,7 +16,7 @@ void write_plotfile(const Geometry& geom, const MultiFab& plotmf) { std::string plotfile_name("plt00000"); - amrex::Print() << "Writing " << plotfile_name << std::endl; + amrex::Print() << "Writing " << plotfile_name << '\n'; #if (AMREX_SPACEDIM == 2) EB_WriteSingleLevelPlotfile(plotfile_name, plotmf, @@ -204,12 +204,12 @@ int main (int argc, char* argv[]) LinOpBCType::Periodic)}); - amrex::Print() << " \n********************************************************************" << std::endl; - amrex::Print() << " Let's project the initial velocity to find " << std::endl; - amrex::Print() << " the flow field around the obstacles ... " << std::endl; - amrex::Print() << " The domain has " << n_cell_x << " cells in the x-direction " << std::endl; - amrex::Print() << " The maximum grid size is " << max_grid_size << std::endl; - amrex::Print() << "******************************************************************** \n" << std::endl; + amrex::Print() << " \n********************************************************************" << '\n'; + amrex::Print() << " Let's project the initial velocity to find " << '\n'; + amrex::Print() << " the flow field around the obstacles ... " << '\n'; + amrex::Print() << " The domain has " << n_cell_x << " cells in the x-direction " << '\n'; + amrex::Print() << " The maximum grid size is " << max_grid_size << '\n'; + amrex::Print() << "******************************************************************** \n" << '\n'; // // Solve div( sigma * grad(phi) ) = RHS @@ -223,9 +223,9 @@ int main (int argc, char* argv[]) // phi.setVal(0.0); // Must initialize phi; we simply set to 0 for this example. // nodal_proj.project( {&phi}, reltol, abstol); - amrex::Print() << " \n********************************************************************" << std::endl; - amrex::Print() << " Projection complete " << std::endl; - amrex::Print() << "******************************************************************** \n" << std::endl; + amrex::Print() << " \n********************************************************************" << '\n'; + amrex::Print() << " Projection complete " << '\n'; + amrex::Print() << "******************************************************************** \n" << '\n'; // Store plotfile variables; velocity and processor id MultiFab plotfile_mf(grids, dmap, AMREX_SPACEDIM+1, 0, MFInfo(), factory); @@ -241,7 +241,7 @@ int main (int argc, char* argv[]) } auto stop_time = amrex::second() - strt_time; - amrex::Print() << "Total run time " << stop_time << std::endl; + amrex::Print() << "Total run time " << stop_time << '\n'; amrex::Finalize(); }