From 60e21263e9b0769efee7dc326a508e178bf22d78 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sun, 10 Dec 2023 16:29:56 -0800 Subject: [PATCH] =?UTF-8?q?Move=20away=20from=20moisture=20being=20a=20com?= =?UTF-8?q?pile-time=20option=20...=20we=20no=20longer=20=E2=80=A6=20(#133?= =?UTF-8?q?4)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Move away from moisture being a compile-time option ... we no longer need to set USE_MOISTURE=TRUE * fix oops * remove print statement --- .github/workflows/cuda-ci.yml | 43 - .github/workflows/hip.yml | 102 -- Build/cmake_with_moisture.sh | 18 - CMake/BuildERFExe.cmake | 42 +- CMakeLists.txt | 2 +- Docs/sphinx_doc/Inputs.rst | 10 +- Docs/sphinx_doc/RegressionTests.rst | 2 +- Docs/sphinx_doc/building.rst | 4 - Exec/ABL/GNUmakefile | 3 - Exec/ABL/prob.H | 5 - Exec/ABL/prob.cpp | 14 +- Exec/ABL_input_sounding/prob.H | 5 - Exec/ABL_input_sounding/prob.cpp | 14 +- Exec/DevTests/MetGrid/GNUmakefile | 3 - Exec/DevTests/MiguelDev/prob.H | 5 - Exec/DevTests/MiguelDev/prob.cpp | 5 - Exec/DevTests/MovingTerrain/prob.H | 7 +- Exec/DevTests/MovingTerrain/prob.cpp | 14 +- Exec/DevTests/MultiBlock/prob.H | 5 - Exec/DevTests/MultiBlock/prob.cpp | 14 +- Exec/DevTests/ParticlesOverWoA/prob.H | 5 - Exec/DevTests/ParticlesOverWoA/prob.cpp | 14 +- Exec/LLJ/prob.H | 5 - Exec/LLJ/prob.cpp | 5 - Exec/Make.ERF | 48 +- Exec/Radiation/GNUmakefile | 3 - Exec/Radiation/inputs_radiation | 2 +- Exec/Radiation/prob.H | 5 - Exec/Radiation/prob.cpp | 16 +- Exec/RegTests/Bubble/GNUmakefile | 2 - Exec/RegTests/Bubble/prob.H | 5 - Exec/RegTests/Bubble/prob.cpp | 20 +- Exec/RegTests/CouetteFlow/prob.H | 5 - Exec/RegTests/CouetteFlow/prob.cpp | 14 +- Exec/RegTests/DensityCurrent/GNUmakefile | 3 - Exec/RegTests/DensityCurrent/prob.H | 5 - Exec/RegTests/DensityCurrent/prob.cpp | 20 +- Exec/RegTests/DynamicRefinement/prob.H | 5 - Exec/RegTests/DynamicRefinement/prob.cpp | 14 +- Exec/RegTests/EkmanSpiral_custom/prob.H | 5 - Exec/RegTests/EkmanSpiral_custom/prob.cpp | 14 +- Exec/RegTests/EkmanSpiral_ideal/prob.H | 5 - Exec/RegTests/EkmanSpiral_ideal/prob.cpp | 5 - .../EkmanSpiral_input_sounding/prob.H | 5 - .../EkmanSpiral_input_sounding/prob.cpp | 5 - Exec/RegTests/GABLS1/GNUmakefile | 3 - Exec/RegTests/GABLS1/prob.H | 5 - Exec/RegTests/GABLS1/prob.cpp | 14 +- Exec/RegTests/IsentropicVortex/prob.H | 5 - Exec/RegTests/IsentropicVortex/prob.cpp | 14 +- Exec/RegTests/PoiseuilleFlow/prob.H | 5 - Exec/RegTests/PoiseuilleFlow/prob.cpp | 14 +- Exec/RegTests/ScalarAdvDiff/GNUmakefile | 3 - Exec/RegTests/ScalarAdvDiff/prob.H | 5 - Exec/RegTests/ScalarAdvDiff/prob.cpp | 15 +- .../ScalarAdvDiff/prob.cpp.convergence | 218 ++++ Exec/RegTests/StokesSecondProblem/prob.cpp | 14 +- Exec/RegTests/TaylorGreenVortex/prob.H | 5 - Exec/RegTests/TaylorGreenVortex/prob.cpp | 14 +- Exec/RegTests/Terrain2d_Cylinder/prob.H | 5 - Exec/RegTests/Terrain2d_Cylinder/prob.cpp | 14 +- Exec/RegTests/Terrain3d_Hemisphere/prob.H | 5 - Exec/RegTests/Terrain3d_Hemisphere/prob.cpp | 14 +- Exec/RegTests/WPS_Test/GNUmakefile | 2 - Exec/RegTests/WitchOfAgnesi/prob.H | 5 - Exec/RegTests/WitchOfAgnesi/prob.cpp | 14 +- Exec/SquallLine_2D/GNUmakefile | 3 - Exec/SquallLine_2D/prob.H | 5 - Exec/SquallLine_2D/prob.cpp | 12 +- Exec/SuperCell/GNUmakefile | 3 - Exec/SuperCell/prob.H | 5 - Exec/SuperCell/prob.cpp | 16 +- .../BoundaryConditions_metgrid.cpp | 2 +- Source/BoundaryConditions/ERF_FillPatch.cpp | 6 +- Source/DataStructs/AdvStruct.H | 8 - Source/DataStructs/DataStruct.H | 20 +- Source/Derive.H | 56 - Source/Derive.cpp | 129 +- .../Diffusion/ComputeTurbulentViscosity.cpp | 2 +- Source/Diffusion/DiffusionSrcForState_N.cpp | 22 +- Source/Diffusion/DiffusionSrcForState_T.cpp | 22 +- Source/EOS.H | 5 - Source/ERF.H | 37 +- Source/ERF.cpp | 126 +- Source/ERF_Tagging.cpp | 2 - Source/ERF_make_new_level.cpp | 6 - Source/IO/ERF_ReadBndryPlanes.H | 14 +- Source/IO/ERF_ReadBndryPlanes.cpp | 21 +- Source/IO/Plotfile.cpp | 55 +- Source/IO/Plotfile.cpp.convergence | 1071 +++++++++++++++++ Source/IndexDefines.H | 72 +- Source/Initialization/ERF_init1d.cpp | 8 +- Source/Initialization/ERF_init_bcs.cpp | 58 +- Source/Initialization/ERF_init_custom.cpp | 40 +- .../ERF_init_from_input_sounding.cpp | 8 +- .../Initialization/ERF_init_from_metgrid.cpp | 6 +- .../Initialization/ERF_init_from_wrfinput.cpp | 10 +- Source/Microphysics/Kessler/Init_Kessler.cpp | 4 +- .../Microphysics/Kessler/Update_Kessler.cpp | 4 +- .../SAM/{Cloud.cpp => Cloud_SAM.cpp} | 3 - .../SAM/{Diagnose.cpp => Diagnose_SAM.cpp} | 0 .../SAM/{Init.cpp => Init_SAM.cpp} | 4 +- Source/Microphysics/SAM/Make.package | 8 +- .../SAM/{Update.cpp => Update_SAM.cpp} | 4 +- Source/TimeIntegration/ERF_fast_rhs_MT.cpp | 12 +- Source/TimeIntegration/ERF_fast_rhs_N.cpp | 12 +- Source/TimeIntegration/ERF_fast_rhs_T.cpp | 12 +- Source/TimeIntegration/ERF_make_buoyancy.cpp | 18 +- .../TimeIntegration/ERF_make_fast_coeffs.cpp | 8 +- Source/TimeIntegration/ERF_slow_rhs_post.cpp | 4 +- Source/TimeIntegration/ERF_slow_rhs_pre.cpp | 24 +- Source/main.cpp | 3 +- Source/prob_common.H | 5 - 113 files changed, 1678 insertions(+), 1242 deletions(-) delete mode 100755 Build/cmake_with_moisture.sh create mode 100644 Exec/RegTests/ScalarAdvDiff/prob.cpp.convergence create mode 100644 Source/IO/Plotfile.cpp.convergence rename Source/Microphysics/SAM/{Cloud.cpp => Cloud_SAM.cpp} (95%) rename Source/Microphysics/SAM/{Diagnose.cpp => Diagnose_SAM.cpp} (100%) rename Source/Microphysics/SAM/{Init.cpp => Init_SAM.cpp} (98%) rename Source/Microphysics/SAM/{Update.cpp => Update_SAM.cpp} (95%) diff --git a/.github/workflows/cuda-ci.yml b/.github/workflows/cuda-ci.yml index ef7d7d8c6..ea9b0df02 100644 --- a/.github/workflows/cuda-ci.yml +++ b/.github/workflows/cuda-ci.yml @@ -56,46 +56,3 @@ jobs: -DERF_ENABLE_MPI:BOOL=ON \ -DERF_ENABLE_CUDA:BOOL=ON . cmake --build build-${{matrix.cuda_pkg}} -- -j $(nproc) - - cuda-moisture-build: - runs-on: ubuntu-20.04 - name: (MOISTURE ON) CUDA v${{matrix.cuda_ver}} - strategy: - matrix: - cuda_pkg: [11-0] - include: - - cuda_ver: "11.0" - cuda_pkg: 11-0 - cuda_extra: libcurand-dev-11-0 cuda-cupti-dev-11-0 libcusolver-dev-11-0 libcublas-dev-11-0 libcusparse-dev-11-0 - steps: - - name: Cancel previous runs - uses: styfle/cancel-workflow-action@0.6.0 - with: - access_token: ${{github.token}} - - uses: actions/checkout@v3 - with: - submodules: true - - name: Prepare CUDA environment - run: | - export DEBIAN_FRONTEND=noninteractive - sudo apt-key adv --fetch-keys https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2004/x86_64/3bf863cc.pub - sudo apt-key adv --fetch-keys https://developer.download.nvidia.com/compute/machine-learning/repos/ubuntu2004/x86_64/7fa2af80.pub - echo "deb https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2004/x86_64 /" | sudo tee /etc/apt/sources.list.d/cuda.list - echo "deb https://developer.download.nvidia.com/compute/machine-learning/repos/ubuntu2004/x86_64 /" | sudo tee /etc/apt/sources.list.d/nvidia-ml.list - sudo apt-get update - sudo apt-get install -y --no-install-recommends \ - libopenmpi-dev cuda-command-line-tools-${{matrix.cuda_pkg}} cuda-compiler-${{matrix.cuda_pkg}} cuda-minimal-build-${{matrix.cuda_pkg}} cuda-nvml-dev-${{matrix.cuda_pkg}} cuda-nvtx-${{matrix.cuda_pkg}} ${{matrix.cuda_extra}} - - name: Configure and build - run: | - export PATH=/usr/local/nvidia/bin:/usr/local/cuda-${{matrix.cuda_ver}}/bin:${PATH} - export LD_LIBRARY_PATH=/usr/local/nvidia/lib:/usr/local/nvidia/lib64:/usr/local/cuda-${{matrix.cuda_ver}}/lib:${LD_LIBRARY_PATH} - cmake -Bbuild-${{matrix.cuda_pkg}} \ - -DCMAKE_BUILD_TYPE:STRING=RelWithDebInfo \ - -DAMReX_CUDA_ERROR_CROSS_EXECUTION_SPACE_CALL=ON \ - -DAMReX_CUDA_ERROR_CAPTURE_THIS=ON \ - -DERF_DIM:STRING=3 \ - -DERF_ENABLE_MPI:BOOL=OFF \ - -DERF_ENABLE_MOISTURE:BOOL=ON \ - -DERF_ENABLE_CUDA:BOOL=ON . - cmake --build build-${{matrix.cuda_pkg}} -- -j $(nproc) - diff --git a/.github/workflows/hip.yml b/.github/workflows/hip.yml index 35ed8fa2d..8addaf8cb 100644 --- a/.github/workflows/hip.yml +++ b/.github/workflows/hip.yml @@ -119,105 +119,3 @@ jobs: name: build-and-test path: | ${{runner.workspace}}/ERF/regressioncov.xml - - Build-And-Test-HIP-Moisture: - name: (MOISTURE ON) HIP ROCm GFortran@9.3 C++17 [tests] - runs-on: ubuntu-20.04 - # Have to have -Wno-deprecated-declarations due to deprecated atomicAddNoRet - # Have to have -Wno-gnu-zero-variadic-macro-arguments to avoid - # amrex/Src/Base/AMReX_GpuLaunchGlobal.H:15:5: error: must specify at least one argument for '...' parameter of variadic macro [-Werror,-Wgnu-zero-variadic-macro-arguments] - # __launch_bounds__(amrex_launch_bounds_max_threads) - # ^ - # /opt/rocm-4.1.1/hip/include/hip/hcc_detail/hip_runtime.h:178:71: note: expanded from macro '__launch_bounds__' - # select_impl_(__VA_ARGS__, launch_bounds_impl1, launch_bounds_impl0)(__VA_ARGS__) - # ^ - # /opt/rocm-4.1.1/hip/include/hip/hcc_detail/hip_runtime.h:176:9: note: macro 'select_impl_' defined here - # #define select_impl_(_1, _2, impl_, ...) impl_ - # NOTE: -Werror was removed because ERF specifically had a lot of warnings. It will be a separate task to go through and fix them all... - env: {CXXFLAGS: "-Wall -Wextra -Wpedantic -Wnull-dereference -Wfloat-conversion -Wshadow -Woverloaded-virtual -Wextra-semi -Wunreachable-code -Wno-deprecated-declarations -Wno-gnu-zero-variadic-macro-arguments -Wno-pass-failed"} - steps: - - uses: actions/checkout@v3 - with: - submodules: true - - - name: Dependencies - run: .github/workflows/dependencies/dependencies_hip.sh - - - name: Build & Install - run: | - source /etc/profile.d/rocm.sh - hipcc --version - which clang - which clang++ - - # "mpic++ --showme" forgets open-pal in Ubuntu 20.04 + OpenMPI 4.0.3 - # https://bugs.launchpad.net/ubuntu/+source/openmpi/+bug/1941786 - # https://github.com/open-mpi/ompi/issues/9317 - export LDFLAGS="-lopen-pal" - - cmake \ - -B${{runner.workspace}}/ERF/build-${{matrix.os}} \ - -DCMAKE_VERBOSE_MAKEFILE=ON \ - -DCMAKE_INSTALL_PREFIX:PATH=${{runner.workspace}}/ERF/install \ - -DCMAKE_BUILD_TYPE:STRING=RelWithDebInfo \ - -DERF_DIM:STRING=3 \ - -DERF_ENABLE_MPI:BOOL=ON \ - -DERF_ENABLE_HIP:BOOL=ON \ - -DAMReX_AMD_ARCH=gfx908 \ - -DERF_ENABLE_TESTS:BOOL=ON \ - -DERF_ENABLE_ALL_WARNINGS:BOOL=ON \ - -DERF_ENABLE_FCOMPARE:BOOL=ON \ - -DERF_ENABLE_MOISTURE:BOOL=ON \ - -DCODECOV:BOOL=ON \ - -DCMAKE_C_COMPILER=$(which clang) \ - -DCMAKE_CXX_COMPILER=$(which clang++) \ - -DCMAKE_CXX_STANDARD=17 \ - ${{github.workspace}}; - # ${{matrix.mpipreflags}} \ - - # for some reason this cmake command fails to build the code, - # and we need to use the make command that follows instead ... - # cmake --build ${{runner.workspace}}/ERF/build-${{matrix.os}} --parallel ${{env.NPROCS}}; - - pushd ${{runner.workspace}}/ERF/build-${{matrix.os}}; - # make -j ${{env.NPROCS}}; - make -j 2; - - # - name: Regression Tests - # run: | - # ctest -L regression -VV - # working-directory: ${{runner.workspace}}/ERF/build-${{matrix.os}} - - # - name: Generate coverage report - # working-directory: ${{runner.workspace}}/ERF/build-${{matrix.os}} - # run: | - # find . -type f -name '*.gcno' -path "**Source**" -exec gcov -pb {} + - # cd .. - # gcovr -g -k -r . --xml regressioncov.xml # -v - - # - name: Upload coverage to Codecov - # uses: codecov/codecov-action@v3 - # with: - # dry_run: false - # # token: ${{ secrets.CODECOV_TOKEN }} # not required for public repos - # files: ./regressioncov.xml # optional - # flags: regtests # optional - # # name: codecov-umbrella # optional - # fail_ci_if_error: true # optional (default = false) - # verbose: true # optional (default = false) - # directory: ${{runner.workspace}}/ERF - - - name: Success artifacts - uses: actions/upload-artifact@v3 - if: success() - with: - name: build-and-test - path: | - ${{runner.workspace}}/ERF/regressioncov.xml - - name: Failing test artifacts - uses: actions/upload-artifact@v3 - if: failure() - with: - name: build-and-test - path: | - ${{runner.workspace}}/ERF/regressioncov.xml diff --git a/Build/cmake_with_moisture.sh b/Build/cmake_with_moisture.sh deleted file mode 100755 index ccae64360..000000000 --- a/Build/cmake_with_moisture.sh +++ /dev/null @@ -1,18 +0,0 @@ -#!/bin/bash - -# Example CMake config script for an OSX laptop with OpenMPI - -cmake -DCMAKE_INSTALL_PREFIX:PATH=./install \ - -DCMAKE_CXX_COMPILER:STRING=mpicxx \ - -DCMAKE_C_COMPILER:STRING=mpicc \ - -DCMAKE_Fortran_COMPILER:STRING=mpifort \ - -DMPIEXEC_PREFLAGS:STRING=--oversubscribe \ - -DCMAKE_BUILD_TYPE:STRING=Release \ - -DERF_DIM:STRING=3 \ - -DERF_ENABLE_MPI:BOOL=ON \ - -DERF_ENABLE_TESTS:BOOL=ON \ - -DERF_ENABLE_FCOMPARE:BOOL=ON \ - -DERF_ENABLE_DOCUMENTATION:BOOL=OFF \ - -DERF_ENABLE_MOISTURE:BOOL=ON \ - -DCMAKE_EXPORT_COMPILE_COMMANDS:BOOL=ON \ - .. && make -j8 diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index dbf0788ad..358149741 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -15,9 +15,7 @@ function(build_erf_lib erf_lib_name) include(${CMAKE_SOURCE_DIR}/CMake/SetERFCompileFlags.cmake) set_erf_compile_flags(${erf_lib_name}) - if(ERF_ENABLE_MOISTURE) - target_compile_definitions(${erf_lib_name} PUBLIC ERF_USE_MOISTURE) - endif() + target_compile_definitions(${erf_lib_name} PUBLIC ERF_USE_MOISTURE) if(ERF_ENABLE_MULTIBLOCK) target_sources(${erf_lib_name} PRIVATE @@ -57,22 +55,6 @@ function(build_erf_lib erf_lib_name) target_compile_definitions(${erf_lib_name} PUBLIC ERF_USE_NETCDF) endif() - if(ERF_ENABLE_MOISTURE) - target_sources(${erf_lib_name} PRIVATE - ${SRC_DIR}/Microphysics/SAM/Init.cpp - ${SRC_DIR}/Microphysics/SAM/Cloud.cpp - ${SRC_DIR}/Microphysics/SAM/IceFall.cpp - ${SRC_DIR}/Microphysics/SAM/Precip.cpp - ${SRC_DIR}/Microphysics/SAM/PrecipFall.cpp - ${SRC_DIR}/Microphysics/SAM/Diagnose.cpp - ${SRC_DIR}/Microphysics/SAM/Update.cpp - ${SRC_DIR}/Microphysics/Kessler/Init_Kessler.cpp - ${SRC_DIR}/Microphysics/Kessler/Kessler.cpp - ${SRC_DIR}/Microphysics/Kessler/Diagnose_Kessler.cpp - ${SRC_DIR}/Microphysics/Kessler/Update_Kessler.cpp) - target_compile_definitions(${erf_lib_name} PUBLIC ERF_USE_MOISTURE) - endif() - if(ERF_ENABLE_RRTMGP) target_sources(${erf_lib_name} PRIVATE ${SRC_DIR}/Radiation/Init_rrtmgp.cpp @@ -166,6 +148,17 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/Utils/TerrainMetrics.cpp ${SRC_DIR}/Utils/VelocityToMomentum.cpp ${SRC_DIR}/Utils/InteriorGhostCells.cpp + ${SRC_DIR}/Microphysics/SAM/Init_SAM.cpp + ${SRC_DIR}/Microphysics/SAM/Cloud_SAM.cpp + ${SRC_DIR}/Microphysics/SAM/IceFall.cpp + ${SRC_DIR}/Microphysics/SAM/Precip.cpp + ${SRC_DIR}/Microphysics/SAM/PrecipFall.cpp + ${SRC_DIR}/Microphysics/SAM/Diagnose_SAM.cpp + ${SRC_DIR}/Microphysics/SAM/Update_SAM.cpp + ${SRC_DIR}/Microphysics/Kessler/Init_Kessler.cpp + ${SRC_DIR}/Microphysics/Kessler/Kessler.cpp + ${SRC_DIR}/Microphysics/Kessler/Diagnose_Kessler.cpp + ${SRC_DIR}/Microphysics/Kessler/Update_Kessler.cpp ) if(NOT "${erf_exe_name}" STREQUAL "erf_unit_tests") @@ -187,13 +180,6 @@ function(build_erf_lib erf_lib_name) endif() endif() - if(ERF_ENABLE_MOISTURE) - target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics) - target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/Null) - target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/SAM) - target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/Kessler) - endif() - if(ERF_ENABLE_RRTMGP) target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Radiation) endif() @@ -213,6 +199,10 @@ function(build_erf_lib erf_lib_name) target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/TimeIntegration) target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Utils) target_include_directories(${erf_lib_name} PUBLIC ${CMAKE_BINARY_DIR}) + target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics) + target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/Null) + target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/SAM) + target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/Kessler) if(ERF_ENABLE_RRTMGP) target_link_libraries(${erf_lib_name} PUBLIC yakl) diff --git a/CMakeLists.txt b/CMakeLists.txt index c0e4dc7ab..538c10741 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -17,7 +17,7 @@ option(ERF_ENABLE_HDF5 "Enable HDF5 IO" ${ERF_ENABLE_NETCDF}) option(ERF_ENABLE_FCOMPARE "Enable building fcompare when not testing" OFF) set(ERF_PRECISION "DOUBLE" CACHE STRING "Floating point precision SINGLE or DOUBLE") -option(ERF_ENABLE_MOISTURE "Enable Full Moisture" OFF) +option(ERF_ENABLE_MOISTURE "Enable Full Moisture" ON) option(ERF_ENABLE_WARM_NO_PRECIP "Enable Warm Moisture" OFF) option(ERF_ENABLE_RRTMGP "Enable RTE-RRTMGP Radiation" OFF) diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst index 779c7bfb7..ebf539b16 100644 --- a/Docs/sphinx_doc/Inputs.rst +++ b/Docs/sphinx_doc/Inputs.rst @@ -883,14 +883,8 @@ Examples of Usage Moisture ======== -ERF has two different moisture models -- one that includes only water vapor and cloud water, -and a more complete model that includes water vapor, cloud water, cloud ice, rain, snow and graupel. - -If ERF is compiled with ERF_USE_WARM_NO_PRECIP defined, then the first model is used and no -further inputs are required. - -If ERF is compiled with ERF_USE_MOISTURE defined, then the following run-time options control how -the full moisture model is used. +ERF has several different moisture models. +The following run-time options control how the full moisture model is used. List of Parameters ------------------ diff --git a/Docs/sphinx_doc/RegressionTests.rst b/Docs/sphinx_doc/RegressionTests.rst index 56266be2d..ca3e7e111 100644 --- a/Docs/sphinx_doc/RegressionTests.rst +++ b/Docs/sphinx_doc/RegressionTests.rst @@ -24,7 +24,7 @@ The following problems are currently tested in the CI: +-------------------------------+----------+----------+----------+------------+-------+-----------------------+ | Test | nx ny nz | xbc | ybc | zbc | Ext | Other | +===============================+==========+==========+==========+============+=======+=======================+ -| Bubble_Density_Current | 256 4 64 | Symmetry | Periodic | SlipWall | None | USE_MOISTURE=TRUE | +| Bubble_Density_Current | 256 4 64 | Symmetry | Periodic | SlipWall | None | moist bubble | | | | Outflow | | SlipWall | | | +-------------------------------+----------+----------+----------+------------+-------+-----------------------+ | CouetteFlow | 32 4 16 | Periodic | Periodic | SlipWall | None | inhomogeneous | diff --git a/Docs/sphinx_doc/building.rst b/Docs/sphinx_doc/building.rst index ac587667b..3a6cbf5c4 100644 --- a/Docs/sphinx_doc/building.rst +++ b/Docs/sphinx_doc/building.rst @@ -85,8 +85,6 @@ or if using tcsh, +--------------------+------------------------------+------------------+-------------+ | USE_HDF5 | Whether to enable HDF5 | TRUE / FALSE | FALSE | +--------------------+------------------------------+------------------+-------------+ - | USE_MOISTURE | Whether to enable moisture | TRUE / FALSE | FALSE | - +--------------------+------------------------------+------------------+-------------+ | USE_WARM_NO_PRECIP | Whether to use warm moisture | TRUE / FALSE | FALSE | +--------------------+------------------------------+------------------+-------------+ | USE_MULTIBLOCK | Whether to enable multiblock | TRUE / FALSE | FALSE | @@ -176,8 +174,6 @@ Analogous to GNU Make, the list of cmake directives is as follows: +---------------------------+------------------------------+------------------+-------------+ | ERF_ENABLE_HDF5 | Whether to enable HDF5 | TRUE / FALSE | FALSE | +---------------------------+------------------------------+------------------+-------------+ - | ERF_ENABLE_MOISTURE | Whether to enable moisture | TRUE / FALSE | FALSE | - +---------------------------+------------------------------+------------------+-------------+ | ERF_ENABLE_WARM_NO_PRECIP | Whether to use warm moisture | TRUE / FALSE | FALSE | +---------------------------+------------------------------+------------------+-------------+ | ERF_ENABLE_MULTIBLOCK | Whether to enable multiblock | TRUE / FALSE | FALSE | diff --git a/Exec/ABL/GNUmakefile b/Exec/ABL/GNUmakefile index 63dbe29ac..c85179014 100644 --- a/Exec/ABL/GNUmakefile +++ b/Exec/ABL/GNUmakefile @@ -18,9 +18,6 @@ USE_CUDA = FALSE USE_HIP = FALSE USE_SYCL = FALSE -# Physics -USE_MOISTURE = FALSE - # Debugging DEBUG = FALSE diff --git a/Exec/ABL/prob.H b/Exec/ABL/prob.H index 6ceea9f5f..af697ac01 100644 --- a/Exec/ABL/prob.H +++ b/Exec/ABL/prob.H @@ -62,14 +62,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/ABL/prob.cpp b/Exec/ABL/prob.cpp index 46350b794..84724b1e3 100644 --- a/Exec/ABL/prob.cpp +++ b/Exec/ABL/prob.cpp @@ -56,14 +56,9 @@ Problem::init_custom_pert( amrex::Array4 const& /*p_hse*/, amrex::Array4 const& /*z_nd*/, amrex::Array4 const& /*z_cc*/, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& /*qv*/, amrex::Array4 const& /*qc*/, amrex::Array4 const& /*qi*/, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& /*qv*/, - amrex::Array4 const& /*qc*/, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& /*mf_m*/, amrex::Array4 const& /*mf_u*/, @@ -98,13 +93,8 @@ Problem::init_custom_pert( // Set an initial value for QKE state(i, j, k, RhoQKE_comp) = parms.QKE_0; -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = 0.0; - state(i, j, k, RhoQp_comp) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = 0.0; - state(i, j, k, RhoQc_comp) = 0.0; -#endif + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; }); // Set the x-velocity diff --git a/Exec/ABL_input_sounding/prob.H b/Exec/ABL_input_sounding/prob.H index 3c725869a..f3e709ece 100644 --- a/Exec/ABL_input_sounding/prob.H +++ b/Exec/ABL_input_sounding/prob.H @@ -61,14 +61,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/ABL_input_sounding/prob.cpp b/Exec/ABL_input_sounding/prob.cpp index 98deed87a..09c3d0090 100644 --- a/Exec/ABL_input_sounding/prob.cpp +++ b/Exec/ABL_input_sounding/prob.cpp @@ -55,14 +55,9 @@ Problem::init_custom_pert( amrex::Array4 const& /*p_hse*/, amrex::Array4 const& /*z_nd*/, amrex::Array4 const& /*z_cc*/, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& /*qv*/, amrex::Array4 const& /*qc*/, amrex::Array4 const& /*qi*/, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& /*qv*/, - amrex::Array4 const& /*qc*/, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& /*mf_m*/, amrex::Array4 const& /*mf_u*/, @@ -91,13 +86,8 @@ Problem::init_custom_pert( // Set an initial value for QKE state(i, j, k, RhoQKE_comp) = parms.QKE_0; -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = 0.0; - state(i, j, k, RhoQp_comp) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = 0.0; - state(i, j, k, RhoQc_comp) = 0.0; -#endif + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; }); // Set the x-velocity diff --git a/Exec/DevTests/MetGrid/GNUmakefile b/Exec/DevTests/MetGrid/GNUmakefile index 0b412a6bb..ef426569a 100644 --- a/Exec/DevTests/MetGrid/GNUmakefile +++ b/Exec/DevTests/MetGrid/GNUmakefile @@ -27,9 +27,6 @@ USE_ASSERTION = TRUE USE_NETCDF = TRUE -#USE_MOISTURE = FALSE -USE_MOISTURE = TRUE - #USE_WARM_NO_PRECIP = TRUE # GNU Make diff --git a/Exec/DevTests/MiguelDev/prob.H b/Exec/DevTests/MiguelDev/prob.H index a81e8923f..468a61cc3 100644 --- a/Exec/DevTests/MiguelDev/prob.H +++ b/Exec/DevTests/MiguelDev/prob.H @@ -51,14 +51,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/DevTests/MiguelDev/prob.cpp b/Exec/DevTests/MiguelDev/prob.cpp index 00a636cb7..5db590f00 100644 --- a/Exec/DevTests/MiguelDev/prob.cpp +++ b/Exec/DevTests/MiguelDev/prob.cpp @@ -47,14 +47,9 @@ Problem::init_custom_pert( Array4 const&, Array4 const&, Array4 const&, -#if defined(ERF_USE_MOISTURE) Array4 const&, Array4 const&, Array4 const&, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const&, - Array4 const&, -#endif amrex::GeometryData const& geomdata, Array4 const&, Array4 const&, diff --git a/Exec/DevTests/MovingTerrain/prob.H b/Exec/DevTests/MovingTerrain/prob.H index 7c80b7bcd..495a7eea7 100644 --- a/Exec/DevTests/MovingTerrain/prob.H +++ b/Exec/DevTests/MovingTerrain/prob.H @@ -32,14 +32,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, @@ -60,7 +55,7 @@ public: amrex::Geometry const& geom) override; protected: - std::string name() override { return "Supercell"; } + std::string name() override { return "MovingTerrain"; } private: ProbParm parms; diff --git a/Exec/DevTests/MovingTerrain/prob.cpp b/Exec/DevTests/MovingTerrain/prob.cpp index f18760759..fe467cab9 100644 --- a/Exec/DevTests/MovingTerrain/prob.cpp +++ b/Exec/DevTests/MovingTerrain/prob.cpp @@ -36,14 +36,9 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const& z_nd, Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) Array4 const&, Array4 const&, Array4 const&, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const&, - Array4 const&, -#endif GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -88,13 +83,8 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = state(i,j,k,Rho_comp); -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = 0.0; - state(i, j, k, RhoQp_comp) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = 0.0; - state(i, j, k, RhoQc_comp) = 0.0; -#endif + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; }); diff --git a/Exec/DevTests/MultiBlock/prob.H b/Exec/DevTests/MultiBlock/prob.H index 8108c26e1..e44d13361 100644 --- a/Exec/DevTests/MultiBlock/prob.H +++ b/Exec/DevTests/MultiBlock/prob.H @@ -40,14 +40,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/DevTests/MultiBlock/prob.cpp b/Exec/DevTests/MultiBlock/prob.cpp index 4f43e1f37..4607e45bf 100644 --- a/Exec/DevTests/MultiBlock/prob.cpp +++ b/Exec/DevTests/MultiBlock/prob.cpp @@ -40,14 +40,9 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const&, Array4 const&, -#if defined(ERF_USE_MOISTURE) Array4 const&, Array4 const&, Array4 const&, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const&, - Array4 const&, -#endif GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -91,13 +86,8 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = 0.0; -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = 0.0; - state(i, j, k, RhoQp_comp) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = 0.0; - state(i, j, k, RhoQc_comp) = 0.0; -#endif + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; }); // Set the x-velocity diff --git a/Exec/DevTests/ParticlesOverWoA/prob.H b/Exec/DevTests/ParticlesOverWoA/prob.H index 695338d96..b0e9c7c05 100644 --- a/Exec/DevTests/ParticlesOverWoA/prob.H +++ b/Exec/DevTests/ParticlesOverWoA/prob.H @@ -40,14 +40,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/DevTests/ParticlesOverWoA/prob.cpp b/Exec/DevTests/ParticlesOverWoA/prob.cpp index b7c179a91..8d52ee3ca 100644 --- a/Exec/DevTests/ParticlesOverWoA/prob.cpp +++ b/Exec/DevTests/ParticlesOverWoA/prob.cpp @@ -41,14 +41,9 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const& z_nd, Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) Array4 const&, Array4 const&, Array4 const&, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const&, - Array4 const&, -#endif GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -90,13 +85,8 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = 0.0; -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = 0.0; - state(i, j, k, RhoQp_comp) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = 0.0; - state(i, j, k, RhoQc_comp) = 0.0; -#endif + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; }); // Set the x-velocity diff --git a/Exec/LLJ/prob.H b/Exec/LLJ/prob.H index 636de03ba..e3e7ae0b8 100644 --- a/Exec/LLJ/prob.H +++ b/Exec/LLJ/prob.H @@ -33,14 +33,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/LLJ/prob.cpp b/Exec/LLJ/prob.cpp index a7c6c999c..d2d78c631 100644 --- a/Exec/LLJ/prob.cpp +++ b/Exec/LLJ/prob.cpp @@ -34,14 +34,9 @@ Problem::init_custom_pert( Array4 const&, Array4 const&, Array4 const&, -#if defined(ERF_USE_MOISTURE) Array4 const&, Array4 const&, Array4 const&, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const&, - Array4 const&, -#endif GeometryData const&, Array4 const&, Array4 const&, diff --git a/Exec/Make.ERF b/Exec/Make.ERF index 5d38535b7..320ae7927 100644 --- a/Exec/Make.ERF +++ b/Exec/Make.ERF @@ -100,33 +100,27 @@ AMReXpack += $(foreach dir, $(AMReXdirs), $(AMREX_HOME)/Src/$(dir)/M include $(AMReXpack) -ifeq ($(USE_MOISTURE), TRUE) - DEFINES += -DERF_USE_MOISTURE - - ERF_MOISTURE_DIR = $(ERF_SOURCE_DIR)/Microphysics - include $(ERF_MOISTURE_DIR)/Make.package - VPATH_LOCATIONS += $(ERF_MOISTURE_DIR) - INCLUDE_LOCATIONS += $(ERF_MOISTURE_DIR) - - ERF_MOISTURE_NULL_DIR = $(ERF_SOURCE_DIR)/Microphysics/Null - include $(ERF_MOISTURE_NULL_DIR)/Make.package - VPATH_LOCATIONS += $(ERF_MOISTURE_NULL_DIR) - INCLUDE_LOCATIONS += $(ERF_MOISTURE_NULL_DIR) - - ERF_MOISTURE_SAM_DIR = $(ERF_SOURCE_DIR)/Microphysics/SAM - include $(ERF_MOISTURE_SAM_DIR)/Make.package - VPATH_LOCATIONS += $(ERF_MOISTURE_SAM_DIR) - INCLUDE_LOCATIONS += $(ERF_MOISTURE_SAM_DIR) - - ERF_MOISTURE_KESSLER_DIR = $(ERF_SOURCE_DIR)/Microphysics/Kessler - include $(ERF_MOISTURE_KESSLER_DIR)/Make.package - VPATH_LOCATIONS += $(ERF_MOISTURE_KESSLER_DIR) - INCLUDE_LOCATIONS += $(ERF_MOISTURE_KESSLER_DIR) -endif - -ifeq ($(USE_WARM_NO_PRECIP), TRUE) - DEFINES += -DERF_USE_WARM_NO_PRECIP -endif +DEFINES += -DERF_USE_MOISTURE + +ERF_MOISTURE_DIR = $(ERF_SOURCE_DIR)/Microphysics +include $(ERF_MOISTURE_DIR)/Make.package +VPATH_LOCATIONS += $(ERF_MOISTURE_DIR) +INCLUDE_LOCATIONS += $(ERF_MOISTURE_DIR) + +ERF_MOISTURE_NULL_DIR = $(ERF_SOURCE_DIR)/Microphysics/Null +include $(ERF_MOISTURE_NULL_DIR)/Make.package +VPATH_LOCATIONS += $(ERF_MOISTURE_NULL_DIR) +INCLUDE_LOCATIONS += $(ERF_MOISTURE_NULL_DIR) + +ERF_MOISTURE_SAM_DIR = $(ERF_SOURCE_DIR)/Microphysics/SAM +include $(ERF_MOISTURE_SAM_DIR)/Make.package +VPATH_LOCATIONS += $(ERF_MOISTURE_SAM_DIR) +INCLUDE_LOCATIONS += $(ERF_MOISTURE_SAM_DIR) + +ERF_MOISTURE_KESSLER_DIR = $(ERF_SOURCE_DIR)/Microphysics/Kessler +include $(ERF_MOISTURE_KESSLER_DIR)/Make.package +VPATH_LOCATIONS += $(ERF_MOISTURE_KESSLER_DIR) +INCLUDE_LOCATIONS += $(ERF_MOISTURE_KESSLER_DIR) ifeq ($(COMPUTE_ERROR), TRUE) DEFINES += -DERF_COMPUTE_ERROR diff --git a/Exec/Radiation/GNUmakefile b/Exec/Radiation/GNUmakefile index 4472b8d49..968d704e3 100644 --- a/Exec/Radiation/GNUmakefile +++ b/Exec/Radiation/GNUmakefile @@ -24,9 +24,6 @@ DEBUG = FALSE TEST = TRUE USE_ASSERTION = TRUE -USE_MOISTURE = TRUE -#USE_WARM_NO_PRECIP = TRUE - # GNU Make Bpack := ./Make.package Blocs := . diff --git a/Exec/Radiation/inputs_radiation b/Exec/Radiation/inputs_radiation index 63b3f839f..e385a3740 100644 --- a/Exec/Radiation/inputs_radiation +++ b/Exec/Radiation/inputs_radiation @@ -38,7 +38,7 @@ amr.check_int = 10000 # number of timesteps between checkpoints # PLOTFILES erf.plot_file_1 = plt # root name of plotfile erf.plot_int_1 = 1 # number of timesteps between plotfiles -erf.plot_vars_1 = density rhotheta rhoQt rhoQp x_velocity y_velocity z_velocity pressure theta temp qt qp qv qc qi +erf.plot_vars_1 = density rhotheta rhoQ1 rhoQ2 x_velocity y_velocity z_velocity pressure theta temp qt qp qv qc qi # SOLVER CHOICE erf.use_gravity = true diff --git a/Exec/Radiation/prob.H b/Exec/Radiation/prob.H index c7d7b8c72..01cc85b99 100644 --- a/Exec/Radiation/prob.H +++ b/Exec/Radiation/prob.H @@ -40,14 +40,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/Radiation/prob.cpp b/Exec/Radiation/prob.cpp index 4b9269c53..a1e175746 100644 --- a/Exec/Radiation/prob.cpp +++ b/Exec/Radiation/prob.cpp @@ -96,14 +96,9 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const& /*z_nd*/, Array4 const& /*z_cc*/, -#if defined(ERF_USE_MOISTURE) Array4 const& qv, Array4 const& qc, Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const& , - Array4 const& , -#endif GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -187,15 +182,14 @@ Problem::init_custom_pert( state(i, j, k, RhoScalar_comp) = 0.0; // mean states -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = rho*qvapor; - state(i, j, k, RhoQp_comp) = 0.0; + state(i, j, k, RhoQ1_comp) = rho*qvapor; + state(i, j, k, RhoQ2_comp) = 0.0; qv(i, j, k) = qvapor; qc(i, j, k) = 0.0; qi(i, j, k) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = rho*qvapor; - state(i, j, k, RhoQc_comp) = 0.0; +#if defined(ERF_USE_WARM_NO_PRECIP) + state(i, j, k, RhoQ1_comp) = rho*qvapor; + state(i, j, k, RhoQ2_comp) = 0.0; #endif }); diff --git a/Exec/RegTests/Bubble/GNUmakefile b/Exec/RegTests/Bubble/GNUmakefile index 82226520b..fef156b43 100644 --- a/Exec/RegTests/Bubble/GNUmakefile +++ b/Exec/RegTests/Bubble/GNUmakefile @@ -23,8 +23,6 @@ DEBUG = FALSE TEST = TRUE USE_ASSERTION = TRUE -USE_MOISTURE = TRUE - # GNU Make Bpack := ./Make.package Blocs := . diff --git a/Exec/RegTests/Bubble/prob.H b/Exec/RegTests/Bubble/prob.H index ea84220c8..0218aff1d 100644 --- a/Exec/RegTests/Bubble/prob.H +++ b/Exec/RegTests/Bubble/prob.H @@ -56,14 +56,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/Bubble/prob.cpp b/Exec/RegTests/Bubble/prob.cpp index 053c802fa..b01a5e963 100644 --- a/Exec/RegTests/Bubble/prob.cpp +++ b/Exec/RegTests/Bubble/prob.cpp @@ -47,14 +47,9 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const& /*z_nd*/, Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) Array4 const&, Array4 const&, Array4 const&, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const&, - Array4 const&, -#endif GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -136,13 +131,8 @@ Problem::init_custom_pert( state(i, j, k, RhoScalar_comp) = 0.0; -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = 0.0; - state(i, j, k, RhoQp_comp) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = 0.0; - state(i, j, k, RhoQc_comp) = 0.0; -#endif + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; }); } else { @@ -187,10 +177,8 @@ Problem::init_custom_pert( state(i, j, k, RhoScalar_comp) = 0.0; -#ifdef ERF_USE_MOISTURE - state(i, j, k, RhoQt_comp) = 0.0; - state(i, j, k, RhoQp_comp) = 0.0; -#endif + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; }); } diff --git a/Exec/RegTests/CouetteFlow/prob.H b/Exec/RegTests/CouetteFlow/prob.H index 55f01f63d..6bfe6a227 100644 --- a/Exec/RegTests/CouetteFlow/prob.H +++ b/Exec/RegTests/CouetteFlow/prob.H @@ -35,14 +35,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/CouetteFlow/prob.cpp b/Exec/RegTests/CouetteFlow/prob.cpp index 5f9f795bf..bd088184f 100644 --- a/Exec/RegTests/CouetteFlow/prob.cpp +++ b/Exec/RegTests/CouetteFlow/prob.cpp @@ -37,14 +37,9 @@ Problem::init_custom_pert( Array4 const&, Array4 const&, Array4 const&, -#if defined(ERF_USE_MOISTURE) Array4 const&, Array4 const&, Array4 const&, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const&, - Array4 const&, -#endif amrex::GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -56,13 +51,8 @@ Problem::init_custom_pert( // Set scalar to 0 state(i, j, k, RhoScalar_comp) = 0.0; -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = 0.0; - state(i, j, k, RhoQp_comp) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = 0.0; - state(i, j, k, RhoQc_comp) = 0.0; -#endif + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; }); diff --git a/Exec/RegTests/DensityCurrent/GNUmakefile b/Exec/RegTests/DensityCurrent/GNUmakefile index 073213506..75690762b 100644 --- a/Exec/RegTests/DensityCurrent/GNUmakefile +++ b/Exec/RegTests/DensityCurrent/GNUmakefile @@ -24,9 +24,6 @@ DEBUG = FALSE TEST = TRUE USE_ASSERTION = TRUE -#USE_MOISTURE = FALSE -#USE_WARM_NO_PRECIP = TRUE - # GNU Make Bpack := ./Make.package Blocs := . diff --git a/Exec/RegTests/DensityCurrent/prob.H b/Exec/RegTests/DensityCurrent/prob.H index a17015b62..c6ca99772 100644 --- a/Exec/RegTests/DensityCurrent/prob.H +++ b/Exec/RegTests/DensityCurrent/prob.H @@ -38,14 +38,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/DensityCurrent/prob.cpp b/Exec/RegTests/DensityCurrent/prob.cpp index a2adaddbf..87af727b5 100644 --- a/Exec/RegTests/DensityCurrent/prob.cpp +++ b/Exec/RegTests/DensityCurrent/prob.cpp @@ -40,14 +40,9 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const& /*z_nd*/, Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) Array4 const&, Array4 const&, Array4 const&, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const&, - Array4 const&, -#endif GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -93,13 +88,8 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = 0.0; -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = 0.0; - state(i, j, k, RhoQp_comp) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = 0.0; - state(i, j, k, RhoQc_comp) = 0.0; -#endif + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; }); } else { amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept @@ -128,10 +118,8 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = 0.0; -#ifdef ERF_USE_MOISTURE - state(i, j, k, RhoQt_comp) = 0.0; - state(i, j, k, RhoQp_comp) = 0.0; -#endif + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; }); } diff --git a/Exec/RegTests/DynamicRefinement/prob.H b/Exec/RegTests/DynamicRefinement/prob.H index 42ebbd98a..b53f1d829 100644 --- a/Exec/RegTests/DynamicRefinement/prob.H +++ b/Exec/RegTests/DynamicRefinement/prob.H @@ -45,14 +45,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/DynamicRefinement/prob.cpp b/Exec/RegTests/DynamicRefinement/prob.cpp index a7fb04dc1..76cbb4bae 100644 --- a/Exec/RegTests/DynamicRefinement/prob.cpp +++ b/Exec/RegTests/DynamicRefinement/prob.cpp @@ -84,14 +84,9 @@ Problem::init_custom_pert( Array4 const&, Array4 const&, Array4 const&, -#if defined(ERF_USE_MOISTURE) Array4 const&, Array4 const&, Array4 const&, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const&, - Array4 const&, -#endif amrex::GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -131,13 +126,8 @@ Problem::init_custom_pert( // Set scalar = 0 -- unused state(i, j, k, RhoScalar_comp) = 1.0 + Omg; -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = 0.0; - state(i, j, k, RhoQp_comp) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = 0.0; - state(i, j, k, RhoQc_comp) = 0.0; -#endif + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; }); // Set the x-velocity diff --git a/Exec/RegTests/EkmanSpiral_custom/prob.H b/Exec/RegTests/EkmanSpiral_custom/prob.H index fca5887e4..040c0af18 100644 --- a/Exec/RegTests/EkmanSpiral_custom/prob.H +++ b/Exec/RegTests/EkmanSpiral_custom/prob.H @@ -33,14 +33,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/EkmanSpiral_custom/prob.cpp b/Exec/RegTests/EkmanSpiral_custom/prob.cpp index 3a279f0ab..c25554060 100644 --- a/Exec/RegTests/EkmanSpiral_custom/prob.cpp +++ b/Exec/RegTests/EkmanSpiral_custom/prob.cpp @@ -34,14 +34,9 @@ Problem::init_custom_pert( Array4 const&, Array4 const&, Array4 const&, -#if defined(ERF_USE_MOISTURE) Array4 const&, Array4 const&, Array4 const&, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const&, - Array4 const&, -#endif GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -59,13 +54,8 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = 0.0; -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = 0.0; - state(i, j, k, RhoQp_comp) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = 0.0; - state(i, j, k, RhoQc_comp) = 0.0; -#endif + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; }); diff --git a/Exec/RegTests/EkmanSpiral_ideal/prob.H b/Exec/RegTests/EkmanSpiral_ideal/prob.H index d9bbd2ec4..9345870cd 100644 --- a/Exec/RegTests/EkmanSpiral_ideal/prob.H +++ b/Exec/RegTests/EkmanSpiral_ideal/prob.H @@ -33,14 +33,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/EkmanSpiral_ideal/prob.cpp b/Exec/RegTests/EkmanSpiral_ideal/prob.cpp index 1c3322473..7cdc1e5cf 100644 --- a/Exec/RegTests/EkmanSpiral_ideal/prob.cpp +++ b/Exec/RegTests/EkmanSpiral_ideal/prob.cpp @@ -34,14 +34,9 @@ Problem::init_custom_pert( Array4 const&, Array4 const&, Array4 const&, -#if defined(ERF_USE_MOISTURE) Array4 const&, Array4 const&, Array4 const&, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const&, - Array4 const&, -#endif GeometryData const&, Array4 const&, Array4 const&, diff --git a/Exec/RegTests/EkmanSpiral_input_sounding/prob.H b/Exec/RegTests/EkmanSpiral_input_sounding/prob.H index 3047c6f8f..10d9a8836 100644 --- a/Exec/RegTests/EkmanSpiral_input_sounding/prob.H +++ b/Exec/RegTests/EkmanSpiral_input_sounding/prob.H @@ -33,14 +33,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/EkmanSpiral_input_sounding/prob.cpp b/Exec/RegTests/EkmanSpiral_input_sounding/prob.cpp index a9ce69a00..0e691087e 100644 --- a/Exec/RegTests/EkmanSpiral_input_sounding/prob.cpp +++ b/Exec/RegTests/EkmanSpiral_input_sounding/prob.cpp @@ -34,14 +34,9 @@ Problem::init_custom_pert( Array4 const&, Array4 const&, Array4 const&, -#if defined(ERF_USE_MOISTURE) Array4 const&, Array4 const&, Array4 const&, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const&, - Array4 const&, -#endif GeometryData const&, Array4 const&, Array4 const&, diff --git a/Exec/RegTests/GABLS1/GNUmakefile b/Exec/RegTests/GABLS1/GNUmakefile index caf4d648e..ff565d83c 100644 --- a/Exec/RegTests/GABLS1/GNUmakefile +++ b/Exec/RegTests/GABLS1/GNUmakefile @@ -18,9 +18,6 @@ USE_CUDA = FALSE USE_HIP = FALSE USE_SYCL = FALSE -# Physics -USE_MOISTURE = FALSE - # Debugging DEBUG = FALSE diff --git a/Exec/RegTests/GABLS1/prob.H b/Exec/RegTests/GABLS1/prob.H index e4a8ee5be..6017dca83 100644 --- a/Exec/RegTests/GABLS1/prob.H +++ b/Exec/RegTests/GABLS1/prob.H @@ -62,14 +62,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/GABLS1/prob.cpp b/Exec/RegTests/GABLS1/prob.cpp index 314897e03..74da64e9b 100644 --- a/Exec/RegTests/GABLS1/prob.cpp +++ b/Exec/RegTests/GABLS1/prob.cpp @@ -56,14 +56,9 @@ Problem::init_custom_pert( amrex::Array4 const& /*p_hse*/, amrex::Array4 const& /*z_nd*/, amrex::Array4 const& /*z_cc*/, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& /*qv*/, amrex::Array4 const& /*qc*/, amrex::Array4 const& /*qi*/, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& /*qv*/, - amrex::Array4 const& /*qc*/, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& /*mf_m*/, amrex::Array4 const& /*mf_u*/, @@ -127,13 +122,8 @@ Problem::init_custom_pert( } } -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = 0.0; - state(i, j, k, RhoQp_comp) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = 0.0; - state(i, j, k, RhoQc_comp) = 0.0; -#endif + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; }); // Set the x-velocity diff --git a/Exec/RegTests/IsentropicVortex/prob.H b/Exec/RegTests/IsentropicVortex/prob.H index bbc43474b..dfb671e71 100644 --- a/Exec/RegTests/IsentropicVortex/prob.H +++ b/Exec/RegTests/IsentropicVortex/prob.H @@ -44,14 +44,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/IsentropicVortex/prob.cpp b/Exec/RegTests/IsentropicVortex/prob.cpp index 3fce4d0a3..e011626d6 100644 --- a/Exec/RegTests/IsentropicVortex/prob.cpp +++ b/Exec/RegTests/IsentropicVortex/prob.cpp @@ -86,14 +86,9 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const&, Array4 const&, -#if defined(ERF_USE_MOISTURE) Array4 const&, Array4 const&, Array4 const&, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const&, - Array4 const&, -#endif amrex::GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -133,13 +128,8 @@ Problem::init_custom_pert( // Set scalar = 0 -- unused state(i, j, k, RhoScalar_comp) = 0.0; -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = 0.0; - state(i, j, k, RhoQp_comp) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = 0.0; - state(i, j, k, RhoQc_comp) = 0.0; -#endif + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; }); // Set the x-velocity diff --git a/Exec/RegTests/PoiseuilleFlow/prob.H b/Exec/RegTests/PoiseuilleFlow/prob.H index a6b3e7ff4..2ab2bff7b 100644 --- a/Exec/RegTests/PoiseuilleFlow/prob.H +++ b/Exec/RegTests/PoiseuilleFlow/prob.H @@ -34,14 +34,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/PoiseuilleFlow/prob.cpp b/Exec/RegTests/PoiseuilleFlow/prob.cpp index de1d608b8..9302e2f2f 100644 --- a/Exec/RegTests/PoiseuilleFlow/prob.cpp +++ b/Exec/RegTests/PoiseuilleFlow/prob.cpp @@ -36,14 +36,9 @@ Problem::init_custom_pert( amrex::Array4 const&, amrex::Array4 const&, amrex::Array4 const&, -#if defined(ERF_USE_MOISTURE) Array4 const&, Array4 const&, Array4 const&, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const&, - Array4 const&, -#endif amrex::GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -57,13 +52,8 @@ Problem::init_custom_pert( // Set scalar = A_0*exp(-10r^2), where r is distance from center of domain state(i, j, k, RhoScalar_comp) = 0.0; -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = 0.0; - state(i, j, k, RhoQp_comp) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = 0.0; - state(i, j, k, RhoQc_comp) = 0.0; -#endif + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; }); diff --git a/Exec/RegTests/ScalarAdvDiff/GNUmakefile b/Exec/RegTests/ScalarAdvDiff/GNUmakefile index 637a48248..a81273724 100644 --- a/Exec/RegTests/ScalarAdvDiff/GNUmakefile +++ b/Exec/RegTests/ScalarAdvDiff/GNUmakefile @@ -24,9 +24,6 @@ DEBUG = FALSE TEST = TRUE USE_ASSERTION = TRUE -#USE_MOISTURE = TRUE -#USE_WARM_NO_PRECIP = TRUE - # GNU Make Bpack := ./Make.package Blocs := . diff --git a/Exec/RegTests/ScalarAdvDiff/prob.H b/Exec/RegTests/ScalarAdvDiff/prob.H index 3f763fdb1..684050b4f 100644 --- a/Exec/RegTests/ScalarAdvDiff/prob.H +++ b/Exec/RegTests/ScalarAdvDiff/prob.H @@ -46,14 +46,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/ScalarAdvDiff/prob.cpp b/Exec/RegTests/ScalarAdvDiff/prob.cpp index a7d7dd560..b2df00229 100644 --- a/Exec/RegTests/ScalarAdvDiff/prob.cpp +++ b/Exec/RegTests/ScalarAdvDiff/prob.cpp @@ -70,14 +70,9 @@ Problem::init_custom_pert( Array4 const&, Array4 const&, Array4 const&, -#if defined(ERF_USE_MOISTURE) Array4 const&, Array4 const&, Array4 const&, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const&, - Array4 const&, -#endif GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -134,14 +129,8 @@ Problem::init_custom_pert( state(i, j, k, RhoScalar_comp) *= parms.rho_0; -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = 0.0; - state(i, j, k, RhoQp_comp) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = 0.0; - state(i, j, k, RhoQc_comp) = 0.0; -#endif - + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; }); // Set the x-velocity diff --git a/Exec/RegTests/ScalarAdvDiff/prob.cpp.convergence b/Exec/RegTests/ScalarAdvDiff/prob.cpp.convergence new file mode 100644 index 000000000..ecc88cf4c --- /dev/null +++ b/Exec/RegTests/ScalarAdvDiff/prob.cpp.convergence @@ -0,0 +1,218 @@ +#include "prob.H" +#include "prob_common.H" + +#include "AMReX_ParmParse.H" +#include "AMReX_MultiFab.H" +#include "ERF_Constants.H" +#include "IndexDefines.H" + +using namespace amrex; + +ProbParm parms; + +void +erf_init_rayleigh(amrex::Vector& tau, + amrex::Vector& ubar, + amrex::Vector& vbar, + amrex::Vector& wbar, + amrex::Vector& thetabar, + amrex::Geometry const& geom) +{ + const int khi = geom.Domain().bigEnd()[2]; + + // We just use these values to test the Rayleigh damping + for (int k = 0; k <= khi; k++) + { + tau[k] = 1.0; + ubar[k] = 2.0; + vbar[k] = 1.0; + wbar[k] = 0.0; + thetabar[k] = parms.Theta_0; + } +} + +void +erf_init_dens_hse(MultiFab& rho_hse, + std::unique_ptr&, + std::unique_ptr&, + amrex::Geometry const&) +{ + Real rho_0 = parms.rho_0; +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(rho_hse,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.growntilebox(1); + const Array4 rho_hse_arr = rho_hse[mfi].array(); + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + rho_hse_arr(i,j,k) = rho_0; + }); + } +} + +void +init_custom_prob( + const Box& bx, + const Box& xbx, + const Box& ybx, + const Box& zbx, + Array4 const& state, + Array4 const& x_vel, + Array4 const& y_vel, + Array4 const& z_vel, + Array4 const&, + Array4 const&, + Array4 const&, + Array4 const&, + Array4 const&, + Array4 const&, + Array4 const&, + GeometryData const& geomdata, + Array4 const& /*mf_m*/, + Array4 const& /*mf_u*/, + Array4 const& /*mf_v*/, + const SolverChoice&) +{ + amrex::ParallelFor(bx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Geometry + const Real* prob_lo = geomdata.ProbLo(); + const Real* prob_hi = geomdata.ProbHi(); + const Real* dx = geomdata.CellSize(); + const Real x = prob_lo[0] + (i + 0.5) * dx[0]; + const Real y = prob_lo[1] + (j + 0.5) * dx[1]; + const Real z = prob_lo[2] + (k + 0.5) * dx[2]; + + // Define a point (xc,yc,zc) at the center of the domain + const Real xc = parms.xc_frac * (prob_lo[0] + prob_hi[0]); + const Real yc = parms.yc_frac * (prob_lo[1] + prob_hi[1]); + const Real zc = parms.zc_frac * (prob_lo[2] + prob_hi[2]); + + // Define ellipse parameters + const Real r0 = parms.rad_0 * (prob_hi[0] - prob_lo[0]); + + const Real r3d = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc) + (z-zc)*(z-zc)); + const Real r2d_xy = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc)); + const Real r2d_xz = std::sqrt((x-xc)*(x-xc) + (z-zc)*(z-zc)); + + // Set the density + state(i, j, k, Rho_comp) = parms.rho_0; + + // Initial potential temperature + state(i, j, k, RhoTheta_comp) = parms.rho_0 * parms.Theta_0; + + if (parms.prob_type == 10) + { + // Set scalar = A_0*exp(-10r^2), where r is distance from center of domain, + // + B_0*sin(x) + state(i, j, k, RhoScalar_comp) = parms.A_0 * exp(-10.*r3d*r3d) + parms.B_0*sin(x); + + } else if (parms.prob_type == 11) { + // state(i, j, k, RhoScalar_comp) = parms.A_0 * 0.25 * (1.0 + std::cos(PI * std::min(r2d_xy, r0) / r0)); + state(i, j, k, RhoScalar_comp) = std::cos(PI * x); + } else if (parms.prob_type == 12) { + state(i, j, k, RhoScalar_comp) = parms.A_0 * 0.25 * (1.0 + std::cos(PI * std::min(r2d_xz, r0) / r0)); + } else if (parms.prob_type == 13) { + const Real r0_z = parms.rad_0 * (prob_hi[2] - prob_lo[2]); + const Real r2d_xz = std::sqrt((x-xc)*(x-xc)/(r0*r0) + (z-zc)*(z-zc)/(r0_z*r0_z)); //ellipse for mapfac shear validation + state(i, j, k, RhoScalar_comp) = parms.A_0 * 0.25 * (1.0 + std::cos(PI * std::min(r2d_xz, r0_z) / r0_z)); + } else { + // Set scalar = A_0 in a ball of radius r0 and 0 elsewhere + if (r3d < r0) { + state(i, j, k, RhoScalar_comp) = parms.A_0; + } else { + state(i, j, k, RhoScalar_comp) = 0.0; + } + } + + state(i, j, k, RhoScalar_comp) *= parms.rho_0; + + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; + + }); + + // Set the x-velocity + amrex::ParallelFor(xbx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + x_vel(i, j, k) = parms.u_0; + + const Real* prob_lo = geomdata.ProbLo(); + const Real* dx = geomdata.CellSize(); + const Real z = prob_lo[2] + (k + 0.5) * dx[2]; + + // Set the x-velocity + x_vel(i, j, k) = parms.u_0 + parms.uRef * + std::log((z + parms.z0)/parms.z0)/ + std::log((parms.zRef +parms.z0)/parms.z0); + }); + + // Set the y-velocity + amrex::ParallelFor(ybx, [=, parms=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + y_vel(i, j, k) = parms.v_0; + }); + + // Set the z-velocity + amrex::ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + z_vel(i, j, k) = 0.0; + }); +} + +void +init_custom_terrain(const Geometry& geom, MultiFab& z_phys_nd, + const Real& /*time*/) +{ + auto dx = geom.CellSizeArray(); + + for ( MFIter mfi(z_phys_nd, TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + const Box& gbx = mfi.growntilebox(1); + Array4 z_arr = z_phys_nd.array(mfi); + ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + + Real z = k * dx[2]; + + // Flat terrain with z = 0 at k = 0 + z_arr(i,j,k) = z; + }); + } + z_phys_nd.FillBoundary(geom.periodicity()); +} + +void +amrex_probinit( + const amrex_real* /*problo*/, + const amrex_real* /*probhi*/) +{ + // Parse params + amrex::ParmParse pp("prob"); + pp.query("rho_0", parms.rho_0); + pp.query("T_0", parms.Theta_0); + pp.query("A_0", parms.A_0); + pp.query("B_0", parms.B_0); + pp.query("u_0", parms.u_0); + pp.query("v_0", parms.v_0); + pp.query("rad_0", parms.rad_0); + pp.query("z0", parms.z0); + pp.query("zRef", parms.zRef); + pp.query("uRef", parms.uRef); + + pp.query("xc_frac", parms.xc_frac); + pp.query("yc_frac", parms.yc_frac); + pp.query("zc_frac", parms.zc_frac); + + pp.query("prob_type", parms.prob_type); +} + +AMREX_GPU_DEVICE +Real +dhdt(int /*i*/, int /*j*/, + const GpuArray /*dx*/, + const Real /*time_mt*/, const Real /*delta_t*/) +{ + return 0.; +} diff --git a/Exec/RegTests/StokesSecondProblem/prob.cpp b/Exec/RegTests/StokesSecondProblem/prob.cpp index cfa8418fd..80ca94bd9 100644 --- a/Exec/RegTests/StokesSecondProblem/prob.cpp +++ b/Exec/RegTests/StokesSecondProblem/prob.cpp @@ -32,14 +32,9 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const&, Array4 const&, -#if defined(ERF_USE_MOISTURE) Array4 const&, Array4 const&, Array4 const&, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const&, - Array4 const&, -#endif GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -59,13 +54,8 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = 0.0; -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = 0.0; - state(i, j, k, RhoQp_comp) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = 0.0; - state(i, j, k, RhoQc_comp) = 0.0; -#endif + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; }); // Set the x-velocity diff --git a/Exec/RegTests/TaylorGreenVortex/prob.H b/Exec/RegTests/TaylorGreenVortex/prob.H index f663c1591..5fca0d226 100644 --- a/Exec/RegTests/TaylorGreenVortex/prob.H +++ b/Exec/RegTests/TaylorGreenVortex/prob.H @@ -33,14 +33,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/TaylorGreenVortex/prob.cpp b/Exec/RegTests/TaylorGreenVortex/prob.cpp index 994a33b44..4b55a5a23 100644 --- a/Exec/RegTests/TaylorGreenVortex/prob.cpp +++ b/Exec/RegTests/TaylorGreenVortex/prob.cpp @@ -36,14 +36,9 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const&, Array4 const&, -#if defined(ERF_USE_MOISTURE) Array4 const&, Array4 const&, Array4 const&, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const&, - Array4 const&, -#endif GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -69,13 +64,8 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = 1.0 * parms.rho_0; -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = 0.0; - state(i, j, k, RhoQp_comp) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = 0.0; - state(i, j, k, RhoQc_comp) = 0.0; -#endif + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; }); // Set the x-velocity diff --git a/Exec/RegTests/Terrain2d_Cylinder/prob.H b/Exec/RegTests/Terrain2d_Cylinder/prob.H index acf61fc1a..46cf99f0c 100644 --- a/Exec/RegTests/Terrain2d_Cylinder/prob.H +++ b/Exec/RegTests/Terrain2d_Cylinder/prob.H @@ -34,14 +34,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/Terrain2d_Cylinder/prob.cpp b/Exec/RegTests/Terrain2d_Cylinder/prob.cpp index fd4c2d18e..5bb9cb8ac 100644 --- a/Exec/RegTests/Terrain2d_Cylinder/prob.cpp +++ b/Exec/RegTests/Terrain2d_Cylinder/prob.cpp @@ -36,14 +36,9 @@ Problem::init_custom_pert( Array4 const&, Array4 const& z_nd, Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) Array4 const&, Array4 const&, Array4 const&, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const&, - Array4 const&, -#endif GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -65,13 +60,8 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = 0.0; -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = 0.0; - state(i, j, k, RhoQp_comp) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = 0.0; - state(i, j, k, RhoQc_comp) = 0.0; -#endif + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; }); // Set the x-velocity diff --git a/Exec/RegTests/Terrain3d_Hemisphere/prob.H b/Exec/RegTests/Terrain3d_Hemisphere/prob.H index 60e622470..a36f54712 100644 --- a/Exec/RegTests/Terrain3d_Hemisphere/prob.H +++ b/Exec/RegTests/Terrain3d_Hemisphere/prob.H @@ -34,14 +34,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/Terrain3d_Hemisphere/prob.cpp b/Exec/RegTests/Terrain3d_Hemisphere/prob.cpp index 338206958..3283b9af7 100644 --- a/Exec/RegTests/Terrain3d_Hemisphere/prob.cpp +++ b/Exec/RegTests/Terrain3d_Hemisphere/prob.cpp @@ -35,14 +35,9 @@ Problem::init_custom_pert( Array4 const&, Array4 const& z_nd, Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) Array4 const&, Array4 const&, Array4 const&, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const&, - Array4 const&, -#endif GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -64,13 +59,8 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = 0.0; -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = 0.0; - state(i, j, k, RhoQp_comp) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = 0.0; - state(i, j, k, RhoQc_comp) = 0.0; -#endif + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; }); // Set the x-velocity diff --git a/Exec/RegTests/WPS_Test/GNUmakefile b/Exec/RegTests/WPS_Test/GNUmakefile index f83bb0ff5..3e501d0a0 100644 --- a/Exec/RegTests/WPS_Test/GNUmakefile +++ b/Exec/RegTests/WPS_Test/GNUmakefile @@ -24,8 +24,6 @@ DEBUG = FALSE TEST = TRUE USE_ASSERTION = TRUE -USE_MOISTURE = TRUE - USE_NETCDF = TRUE #USE_HDF5 = TRUE diff --git a/Exec/RegTests/WitchOfAgnesi/prob.H b/Exec/RegTests/WitchOfAgnesi/prob.H index 0bac29daa..1e4322ecf 100644 --- a/Exec/RegTests/WitchOfAgnesi/prob.H +++ b/Exec/RegTests/WitchOfAgnesi/prob.H @@ -31,14 +31,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/RegTests/WitchOfAgnesi/prob.cpp b/Exec/RegTests/WitchOfAgnesi/prob.cpp index 65fa868c5..04a2e8f6b 100644 --- a/Exec/RegTests/WitchOfAgnesi/prob.cpp +++ b/Exec/RegTests/WitchOfAgnesi/prob.cpp @@ -34,14 +34,9 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const& z_nd, Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) Array4 const&, Array4 const&, Array4 const&, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const&, - Array4 const&, -#endif GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -63,13 +58,8 @@ Problem::init_custom_pert( // Set scalar = 0 everywhere state(i, j, k, RhoScalar_comp) = 0.0; -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = 0.0; - state(i, j, k, RhoQp_comp) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = 0.0; - state(i, j, k, RhoQc_comp) = 0.0; -#endif + state(i, j, k, RhoQ1_comp) = 0.0; + state(i, j, k, RhoQ2_comp) = 0.0; }); // Set the x-velocity diff --git a/Exec/SquallLine_2D/GNUmakefile b/Exec/SquallLine_2D/GNUmakefile index 40f807cd1..5705fc951 100644 --- a/Exec/SquallLine_2D/GNUmakefile +++ b/Exec/SquallLine_2D/GNUmakefile @@ -24,9 +24,6 @@ DEBUG = FALSE TEST = TRUE USE_ASSERTION = TRUE -USE_MOISTURE = TRUE -#USE_WARM_NO_PRECIP = TRUE - # GNU Make Bpack := ./Make.package Blocs := . diff --git a/Exec/SquallLine_2D/prob.H b/Exec/SquallLine_2D/prob.H index 90f8f73f5..f9f1e8617 100644 --- a/Exec/SquallLine_2D/prob.H +++ b/Exec/SquallLine_2D/prob.H @@ -50,14 +50,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/SquallLine_2D/prob.cpp b/Exec/SquallLine_2D/prob.cpp index 973c924f4..e9e90c08d 100644 --- a/Exec/SquallLine_2D/prob.cpp +++ b/Exec/SquallLine_2D/prob.cpp @@ -282,14 +282,9 @@ Problem::init_custom_pert ( Array4 const& /*p_hse*/, Array4 const& /*z_nd*/, Array4 const& /*z_cc*/, -#if defined(ERF_USE_MOISTURE) Array4 const& qv, Array4 const& qc, Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const& , - Array4 const& , -#endif GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -380,13 +375,12 @@ Problem::init_custom_pert ( state(i, j, k, RhoScalar_comp) = rho*scalar; // mean states -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = rho*q_v_hot; - state(i, j, k, RhoQp_comp) = 0.0; + state(i, j, k, RhoQ1_comp) = rho*q_v_hot; + state(i, j, k, RhoQ2_comp) = 0.0; qv(i, j, k) = q_v_hot; qc(i, j, k) = 0.0; qi(i, j, k) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) +#if defined(ERF_USE_WARM_NO_PRECIP) state(i, j, k, RhoQv_comp) = 0.0;//rho*qvapor; state(i, j, k, RhoQc_comp) = 0.0; #endif diff --git a/Exec/SuperCell/GNUmakefile b/Exec/SuperCell/GNUmakefile index 4472b8d49..968d704e3 100644 --- a/Exec/SuperCell/GNUmakefile +++ b/Exec/SuperCell/GNUmakefile @@ -24,9 +24,6 @@ DEBUG = FALSE TEST = TRUE USE_ASSERTION = TRUE -USE_MOISTURE = TRUE -#USE_WARM_NO_PRECIP = TRUE - # GNU Make Bpack := ./Make.package Blocs := . diff --git a/Exec/SuperCell/prob.H b/Exec/SuperCell/prob.H index c7d7b8c72..01cc85b99 100644 --- a/Exec/SuperCell/prob.H +++ b/Exec/SuperCell/prob.H @@ -40,14 +40,9 @@ public: amrex::Array4 const& p_hse, amrex::Array4 const& z_nd, amrex::Array4 const& z_cc, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& qv, amrex::Array4 const& qc, amrex::Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& qv, - amrex::Array4 const& qc, -#endif amrex::GeometryData const& geomdata, amrex::Array4 const& mf_m, amrex::Array4 const& mf_u, diff --git a/Exec/SuperCell/prob.cpp b/Exec/SuperCell/prob.cpp index 4b9269c53..a1e175746 100644 --- a/Exec/SuperCell/prob.cpp +++ b/Exec/SuperCell/prob.cpp @@ -96,14 +96,9 @@ Problem::init_custom_pert( Array4 const& p_hse, Array4 const& /*z_nd*/, Array4 const& /*z_cc*/, -#if defined(ERF_USE_MOISTURE) Array4 const& qv, Array4 const& qc, Array4 const& qi, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Array4 const& , - Array4 const& , -#endif GeometryData const& geomdata, Array4 const& /*mf_m*/, Array4 const& /*mf_u*/, @@ -187,15 +182,14 @@ Problem::init_custom_pert( state(i, j, k, RhoScalar_comp) = 0.0; // mean states -#if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = rho*qvapor; - state(i, j, k, RhoQp_comp) = 0.0; + state(i, j, k, RhoQ1_comp) = rho*qvapor; + state(i, j, k, RhoQ2_comp) = 0.0; qv(i, j, k) = qvapor; qc(i, j, k) = 0.0; qi(i, j, k) = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - state(i, j, k, RhoQv_comp) = rho*qvapor; - state(i, j, k, RhoQc_comp) = 0.0; +#if defined(ERF_USE_WARM_NO_PRECIP) + state(i, j, k, RhoQ1_comp) = rho*qvapor; + state(i, j, k, RhoQ2_comp) = 0.0; #endif }); diff --git a/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp b/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp index 37f4dc04f..1654da092 100644 --- a/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp +++ b/Source/BoundaryConditions/BoundaryConditions_metgrid.cpp @@ -31,7 +31,7 @@ ERF::fill_from_metgrid (const Vector& mfs, // Flags for read vars and index mapping // See IndexDefines.H #if defined(ERF_USE_MOISTURE) - // Cons includes [Rho RhoTheta RhoKE RhoQKE RhoScalar RhoQt RhoQp NumVars] + // Cons includes [Rho RhoTheta RhoKE RhoQKE RhoScalar RhoQ1 RhoQ2 NumVars] Vector cons_read = {1, 1, 0, 0, 0, 1, 0}; Vector cons_map = {MetGridBdyVars::R, MetGridBdyVars::T, 0, 0, 0, MetGridBdyVars::QV, 0}; #elif defined(ERF_USE_WARM_NO_PRECIP) diff --git a/Source/BoundaryConditions/ERF_FillPatch.cpp b/Source/BoundaryConditions/ERF_FillPatch.cpp index 9b1b6dd37..b12c33fdf 100644 --- a/Source/BoundaryConditions/ERF_FillPatch.cpp +++ b/Source/BoundaryConditions/ERF_FillPatch.cpp @@ -110,7 +110,6 @@ ERF::FillPatch (int lev, Real time, const Vector& mfs, bool fillset) * @param[out] mf MultiFab to be filled (qmoist[lev]) */ -#ifdef ERF_USE_MOISTURE void ERF::FillPatchMoistVars (int lev, MultiFab& mf) { @@ -123,8 +122,8 @@ ERF::FillPatchMoistVars (int lev, MultiFab& mf) int ncomp_cons = 1; // We only fill qv, the first component // Note that we are filling qv, stored in qmoist[lev], with the input data (if there is any), stored - // in RhoQt_comp. - int bccomp_cons = BCVars::RhoQt_bc_comp; + // in RhoQ1_comp. + int bccomp_cons = BCVars::RhoQ1_bc_comp; IntVect ngvect_cons = mf.nGrowVect(); IntVect ngvect_vels = {0,0,0}; @@ -135,7 +134,6 @@ ERF::FillPatchMoistVars (int lev, MultiFab& mf) mf.FillBoundary(geom[lev].periodicity()); } -#endif /* * Fill valid and ghost data diff --git a/Source/DataStructs/AdvStruct.H b/Source/DataStructs/AdvStruct.H index 8aa297708..614801b3c 100644 --- a/Source/DataStructs/AdvStruct.H +++ b/Source/DataStructs/AdvStruct.H @@ -31,11 +31,9 @@ struct AdvChoice { pp.query("dryscal_horiz_adv_type" , dryscal_horiz_adv_string); pp.query("dryscal_vert_adv_type" , dryscal_vert_adv_string); -#if defined(ERF_USE_MOISTURE) or defined(ERF_USE_WARM_NO_PRECIP) std::string moistscal_horiz_adv_string = ""; std::string moistscal_vert_adv_string = ""; pp.query("moistscal_horiz_adv_type", moistscal_horiz_adv_string); pp.query("moistscal_vert_adv_type" , moistscal_vert_adv_string); -#endif if (use_efficient_advection){ amrex::Print() << "Using efficient advection scheme" << std::endl;; @@ -99,7 +97,6 @@ struct AdvChoice { amrex::Print() << "Using default dryscal_vert_adv_type" << std::endl;; } -#if defined(ERF_USE_MOISTURE) or defined(ERF_USE_WARM_NO_PRECIP) if ( (moistscal_horiz_adv_string == "Centered_2nd") || (moistscal_horiz_adv_string == "Upwind_3rd" ) || (moistscal_horiz_adv_string == "Centered_4th") || @@ -133,7 +130,6 @@ struct AdvChoice { } else { amrex::Print() << "Using default moistscal_vert_adv_type" << std::endl;; } -#endif } void display() @@ -143,10 +139,8 @@ struct AdvChoice { amrex::Print() << "dycore_vert_adv_type : " << adv_type_convert_int_to_string(dycore_vert_adv_type) << std::endl; amrex::Print() << "dryscal_horiz_adv_type : " << adv_type_convert_int_to_string(dryscal_horiz_adv_type) << std::endl; amrex::Print() << "dryscal_vert_adv_type : " << adv_type_convert_int_to_string(dryscal_vert_adv_type) << std::endl; -#if defined(ERF_USE_MOISTURE) or defined(ERF_USE_WARM_NO_PRECIP) amrex::Print() << "moistscal_horiz_adv_type : " << adv_type_convert_int_to_string(moistscal_horiz_adv_type) << std::endl; amrex::Print() << "moistscal_vert_adv_type : " << adv_type_convert_int_to_string(moistscal_vert_adv_type) << std::endl; -#endif } std::string @@ -214,9 +208,7 @@ struct AdvChoice { AdvType dycore_vert_adv_type = AdvType::Centered_2nd; AdvType dryscal_horiz_adv_type = AdvType::Centered_2nd; AdvType dryscal_vert_adv_type = AdvType::Centered_2nd; -#if defined(ERF_USE_MOISTURE) or defined(ERF_USE_WARM_NO_PRECIP) AdvType moistscal_horiz_adv_type = AdvType::Weno_3; AdvType moistscal_vert_adv_type = AdvType::Weno_3; -#endif }; #endif diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index 970691474..d24d138b7 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -28,6 +28,10 @@ enum struct TerrainType { Static, Moving }; +enum struct MoistureType { + Kessler, SAM, None +}; + enum struct Coord { x, y, z }; @@ -79,6 +83,17 @@ struct SolverChoice { amrex::Abort("terrain_type can be either Moving/moving or Static/static"); } + // What type of moisture model to use + static std::string moisture_model_string = "None"; + pp.query("moisture_model", moisture_model_string); + if (moisture_model_string == "SAM") { + moisture_type = MoistureType::SAM; + } else if (moisture_model_string == "Kessler") { + moisture_type = MoistureType::Kessler; + } else { + moisture_type = MoistureType::None; + } + // Use lagged_delta_rt in the fast integrator? pp.query("use_lagged_delta_rt", use_lagged_delta_rt); @@ -159,11 +174,9 @@ struct SolverChoice { pp.query("Ave_Plane", ave_plane); -#ifdef ERF_USE_MOISTURE pp.query("mp_clouds", do_cloud); pp.query("mp_precip", do_precip); pp.query("use_moist_background", use_moist_background); -#endif // Use numerical diffusion? pp.query("use_NumDiff",use_NumDiff); @@ -350,17 +363,16 @@ struct SolverChoice { CouplingType coupling_type; TerrainType terrain_type; + MoistureType moisture_type; ABLDriverType abl_driver_type; amrex::GpuArray abl_pressure_grad; amrex::GpuArray abl_geo_forcing; int ave_plane {2}; -#ifdef ERF_USE_MOISTURE // Microphysics params bool do_cloud {true}; bool do_precip {true}; bool use_moist_background {false}; -#endif }; #endif diff --git a/Source/Derive.H b/Source/Derive.H index f405ed4f1..f5693bec8 100644 --- a/Source/Derive.H +++ b/Source/Derive.H @@ -23,17 +23,6 @@ void erf_dernull ( const int* bcrec, const int level); -void erf_derpres ( - const amrex::Box& bx, - amrex::FArrayBox& derfab, - int dcomp, - int ncomp, - const amrex::FArrayBox& datfab, - const amrex::Geometry& geomdata, - amrex::Real time, - const int* bcrec, - const int level); - void erf_dersoundspeed ( const amrex::Box& bx, amrex::FArrayBox& derfab, @@ -100,50 +89,5 @@ void erf_derQKE ( const int* bcrec, const int level); -#if defined(ERF_USE_MOISTURE) -void erf_derQt ( - const amrex::Box& bx, - amrex::FArrayBox& derfab, - int dcomp, - int ncomp, - const amrex::FArrayBox& datfab, - const amrex::Geometry& geomdata, - amrex::Real time, - const int* bcrec, - const int level); - -void erf_derQp ( - const amrex::Box& bx, - amrex::FArrayBox& derfab, - int dcomp, - int ncomp, - const amrex::FArrayBox& datfab, - const amrex::Geometry& geomdata, - amrex::Real time, - const int* bcrec, - const int level); -#elif defined(ERF_USE_WARM_NO_PRECIP) -void erf_derQv ( - const amrex::Box& bx, - amrex::FArrayBox& derfab, - int dcomp, - int ncomp, - const amrex::FArrayBox& datfab, - const amrex::Geometry& geomdata, - amrex::Real time, - const int* bcrec, - const int level); -void erf_derQc ( - const amrex::Box& bx, - amrex::FArrayBox& derfab, - int dcomp, - int ncomp, - const amrex::FArrayBox& datfab, - const amrex::Geometry& geomdata, - amrex::Real time, - const int* bcrec, - const int level); -#endif - } #endif diff --git a/Source/Derive.cpp b/Source/Derive.cpp index 9ed1f85c8..dfc232078 100644 --- a/Source/Derive.cpp +++ b/Source/Derive.cpp @@ -46,44 +46,6 @@ erf_dernull ( const int /*level*/) { } -/** - * Function to define pressure by calling an EOS routine - * - * @params[in] bx box on which to divide by density - * @params[out] derfab array of derived quantity -- here it holds pressure - * @params[in] datfab array of data used to construct derived quantity -*/ -void -erf_derpres ( - const amrex::Box& bx, - amrex::FArrayBox& derfab, - int /*dcomp*/, - int /*ncomp*/, - const amrex::FArrayBox& datfab, - const amrex::Geometry& /*geomdata*/, - amrex::Real /*time*/, - const int* /*bcrec*/, - const int /*level*/) -{ - auto const dat = datfab.array(); - auto pfab = derfab.array(); - - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - const amrex::Real rhotheta = dat(i, j, k, RhoTheta_comp); - AMREX_ALWAYS_ASSERT(rhotheta > 0.); -#if defined(ERF_USE_MOISTURE) - // HACK HACK HACK -- FOR NOW THIS IS ONLY THE PARTIAL PRESSURE OF THE DRY AIR - amrex::Real qv = 0.; - pfab(i,j,k) = getPgivenRTh(rhotheta,qv); -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Real qv = dat(i,j,k,RhoQv_comp) / dat(i,j,k,Rho_comp); - pfab(i,j,k) = getPgivenRTh(rhotheta,qv); -#else - pfab(i,j,k) = getPgivenRTh(rhotheta); -#endif - }); -} - /** * Function to define the sound speed by calling an EOS routine * @@ -235,93 +197,4 @@ erf_derQKE ( erf_derrhodivide(bx, derfab, datfab, RhoQKE_comp); } -#if defined(ERF_USE_MOISTURE) -/** - * Function to define total water Qt by dividing (rho Qt) by rho - * - * @params[in] bx box on which to divide by density - * @params[out] derfab array of derived quantity -- here it holds Qt - * @params[in] datfab array of data used to construct derived quantity -*/ -void -erf_derQt ( - const amrex::Box& bx, - amrex::FArrayBox& derfab, - int /*dcomp*/, - int /*ncomp*/, - const amrex::FArrayBox& datfab, - const amrex::Geometry& /*geomdata*/, - amrex::Real /*time*/, - const int* /*bcrec*/, - const int /*level*/) -{ - erf_derrhodivide(bx, derfab, datfab, RhoQt_comp); -} - -/** - * Function to define precipitating water Qp by dividing (rho Qp) by rho - * - * @params[in] bx box on which to divide by density - * @params[out] derfab array of derived quantity -- here it holds Qp - * @params[in] datfab array of data used to construct derived quantity -*/ -void -erf_derQp ( - const amrex::Box& bx, - amrex::FArrayBox& derfab, - int /*dcomp*/, - int /*ncomp*/, - const amrex::FArrayBox& datfab, - const amrex::Geometry& /*geomdata*/, - amrex::Real /*time*/, - const int* /*bcrec*/, - const int /*level*/) -{ - erf_derrhodivide(bx, derfab, datfab, RhoQp_comp); -} -#elif defined(ERF_USE_WARM_NO_PRECIP) -/** - * Function to define water vapor Qv by dividing (rho Qv) by rho - * - * @params[in] bx box on which to divide by density - * @params[out] derfab array of derived quantity -- here it holds Qv - * @params[in] datfab array of data used to construct derived quantity -*/ -void -erf_derQv ( - const amrex::Box& bx, - amrex::FArrayBox& derfab, - int /*dcomp*/, - int /*ncomp*/, - const amrex::FArrayBox& datfab, - const amrex::Geometry& /*geomdata*/, - amrex::Real /*time*/, - const int* /*bcrec*/, - const int /*level*/) -{ - erf_derrhodivide(bx, derfab, datfab, RhoQv_comp); -} -/** - * Function to define cloud water Qc by dividing (rho Qc) by rho - * - * @params[in] bx box on which to divide by density - * @params[out] derfab array of derived quantity -- here it holds Qc - * @params[in] datfab array of data used to construct derived quantity -*/ -void -erf_derQc ( - const amrex::Box& bx, - amrex::FArrayBox& derfab, - int /*dcomp*/, - int /*ncomp*/, - const amrex::FArrayBox& datfab, - const amrex::Geometry& /*geomdata*/, - amrex::Real /*time*/, - const int* /*bcrec*/, - const int /*level*/) -{ - erf_derrhodivide(bx, derfab, datfab, RhoQc_comp); -} -#endif - -} +} // namespace diff --git a/Source/Diffusion/ComputeTurbulentViscosity.cpp b/Source/Diffusion/ComputeTurbulentViscosity.cpp index 17ad13032..029b36c71 100644 --- a/Source/Diffusion/ComputeTurbulentViscosity.cpp +++ b/Source/Diffusion/ComputeTurbulentViscosity.cpp @@ -170,7 +170,7 @@ void ComputeTurbulentViscosityLES (const amrex::MultiFab& Tau11, const amrex::Mu Real inv_Sc_t = turbChoice.Sc_t_inv; Real inv_sigma_k = 1.0 / turbChoice.sigma_k; #if defined(ERF_USE_MOISTURE) - // EddyDiff mapping : Theta_h Scalar_h KE_h QKE_h Qt_h Qp_h + // EddyDiff mapping : Theta_h Scalar_h KE_h QKE_h Q1_h Q2_h Vector Factors = {inv_Pr_t, inv_Sc_t, inv_sigma_k, inv_sigma_k, inv_Sc_t, inv_Sc_t}; // alpha = mu/Pr #elif defined(ERF_USE_WARM_NO_PRECIP) // EddyDiff mapping : Theta_h Scalar_h KE_h QKE_h Qv_h Qc_h diff --git a/Source/Diffusion/DiffusionSrcForState_N.cpp b/Source/Diffusion/DiffusionSrcForState_N.cpp index 0df9470f0..eb097f984 100644 --- a/Source/Diffusion/DiffusionSrcForState_N.cpp +++ b/Source/Diffusion/DiffusionSrcForState_N.cpp @@ -95,11 +95,11 @@ DiffusionSrcForState_N (const amrex::Box& bx, const amrex::Box& domain, alpha_eff[PrimScalar_comp] = diffChoice.alpha_C; break; #if defined(ERF_USE_MOISTURE) - case PrimQt_comp: - alpha_eff[PrimQt_comp] = diffChoice.alpha_C; + case PrimQ1_comp: + alpha_eff[PrimQ1_comp] = diffChoice.alpha_C; break; - case PrimQp_comp: - alpha_eff[PrimQp_comp] = diffChoice.alpha_C; + case PrimQ2_comp: + alpha_eff[PrimQ2_comp] = diffChoice.alpha_C; break; #elif defined(ERF_USE_WARM_NO_PRECIP) case PrimQv_comp: @@ -124,11 +124,11 @@ DiffusionSrcForState_N (const amrex::Box& bx, const amrex::Box& domain, alpha_eff[PrimScalar_comp] = diffChoice.rhoAlpha_C; break; #if defined(ERF_USE_MOISTURE) - case PrimQt_comp: - alpha_eff[PrimQt_comp] = diffChoice.rhoAlpha_C; + case PrimQ1_comp: + alpha_eff[PrimQ1_comp] = diffChoice.rhoAlpha_C; break; - case PrimQp_comp: - alpha_eff[PrimQp_comp] = diffChoice.rhoAlpha_C; + case PrimQ2_comp: + alpha_eff[PrimQ2_comp] = diffChoice.rhoAlpha_C; break; #elif defined(ERF_USE_WARM_NO_PRECIP) case PrimQv_comp: @@ -145,9 +145,9 @@ DiffusionSrcForState_N (const amrex::Box& bx, const amrex::Box& domain, } } #if defined(ERF_USE_MOISTURE) - Vector eddy_diff_idx{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Qt_h, EddyDiff::Qp_h}; - Vector eddy_diff_idy{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Qt_h, EddyDiff::Qp_h}; - Vector eddy_diff_idz{EddyDiff::Theta_v, EddyDiff::KE_v, EddyDiff::QKE_v, EddyDiff::Scalar_v, EddyDiff::Qt_v, EddyDiff::Qp_v}; + Vector eddy_diff_idx{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Q1_h, EddyDiff::Q2_h}; + Vector eddy_diff_idy{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Q1_h, EddyDiff::Q2_h}; + Vector eddy_diff_idz{EddyDiff::Theta_v, EddyDiff::KE_v, EddyDiff::QKE_v, EddyDiff::Scalar_v, EddyDiff::Q1_v, EddyDiff::Q2_v}; #elif defined(ERF_USE_WARM_NO_PRECIP) Vector eddy_diff_idx{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Qv_h, EddyDiff::Qc_h}; Vector eddy_diff_idy{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Qv_h, EddyDiff::Qc_h}; diff --git a/Source/Diffusion/DiffusionSrcForState_T.cpp b/Source/Diffusion/DiffusionSrcForState_T.cpp index edf0a8d21..884937916 100644 --- a/Source/Diffusion/DiffusionSrcForState_T.cpp +++ b/Source/Diffusion/DiffusionSrcForState_T.cpp @@ -102,11 +102,11 @@ DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain, alpha_eff[PrimScalar_comp] = diffChoice.alpha_C; break; #if defined(ERF_USE_MOISTURE) - case PrimQt_comp: - alpha_eff[PrimQt_comp] = diffChoice.alpha_C; + case PrimQ1_comp: + alpha_eff[PrimQ1_comp] = diffChoice.alpha_C; break; - case PrimQp_comp: - alpha_eff[PrimQp_comp] = diffChoice.alpha_C; + case PrimQ2_comp: + alpha_eff[PrimQ2_comp] = diffChoice.alpha_C; break; #elif defined(ERF_USE_WARM_NO_PRECIP) case PrimQv_comp: @@ -131,11 +131,11 @@ DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain, alpha_eff[PrimScalar_comp] = diffChoice.rhoAlpha_C; break; #if defined(ERF_USE_MOISTURE) - case PrimQt_comp: - alpha_eff[PrimQt_comp] = diffChoice.rhoAlpha_C; + case PrimQ1_comp: + alpha_eff[PrimQ1_comp] = diffChoice.rhoAlpha_C; break; - case PrimQp_comp: - alpha_eff[PrimQp_comp] = diffChoice.rhoAlpha_C; + case PrimQ2_comp: + alpha_eff[PrimQ2_comp] = diffChoice.rhoAlpha_C; break; #elif defined(ERF_USE_WARM_NO_PRECIP) case PrimQv_comp: @@ -153,9 +153,9 @@ DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain, } #if defined(ERF_USE_MOISTURE) - Vector eddy_diff_idx{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Qt_h, EddyDiff::Qp_h}; - Vector eddy_diff_idy{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Qt_h, EddyDiff::Qp_h}; - Vector eddy_diff_idz{EddyDiff::Theta_v, EddyDiff::KE_v, EddyDiff::QKE_v, EddyDiff::Scalar_v, EddyDiff::Qt_v, EddyDiff::Qp_v}; + Vector eddy_diff_idx{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Q1_h, EddyDiff::Q2_h}; + Vector eddy_diff_idy{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Q1_h, EddyDiff::Q2_h}; + Vector eddy_diff_idz{EddyDiff::Theta_v, EddyDiff::KE_v, EddyDiff::QKE_v, EddyDiff::Scalar_v, EddyDiff::Q1_v, EddyDiff::Q2_v}; #elif defined(ERF_USE_WARM_NO_PRECIP) Vector eddy_diff_idx{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Qv_h, EddyDiff::Qc_h}; Vector eddy_diff_idy{EddyDiff::Theta_h, EddyDiff::KE_h, EddyDiff::QKE_h, EddyDiff::Scalar_h, EddyDiff::Qv_h, EddyDiff::Qc_h}; diff --git a/Source/EOS.H b/Source/EOS.H index 2ce071025..c2e8c6828 100644 --- a/Source/EOS.H +++ b/Source/EOS.H @@ -61,12 +61,7 @@ amrex::Real getThgivenRandT(const amrex::Real rho, const amrex::Real T, const am AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv = 0.) { -#if defined(ERF_USE_WARM_NO_PRECIP) || defined(ERF_USE_MOISTURE) return p_0 * std::pow(R_d * rhotheta * (1.0+(R_v/R_d)*qv) * ip_0, Gamma); -#else - amrex::ignore_unused(qv); - return p_0 * std::pow(R_d * rhotheta * ip_0, Gamma); -#endif } /** diff --git a/Source/ERF.H b/Source/ERF.H index 10020f20a..e02d106d7 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -44,9 +44,7 @@ #include #include -#ifdef ERF_USE_MOISTURE #include "Microphysics.H" -#endif #ifdef ERF_USE_RRTMGP #include "Radiation.H" @@ -215,10 +213,8 @@ public: amrex::Real dt, amrex::Real time, amrex::InterpFaceRegister* ifr); -#if defined(ERF_USE_MOISTURE) void advance_microphysics (int lev, amrex::MultiFab& cons_in, const amrex::Real& dt_advance); -#endif #if defined(ERF_USE_RRTMGP) void advance_radiation (int lev, @@ -362,11 +358,9 @@ private: // works for single level and 2-level cases (fill fine grid ghost by interpolating from coarse) void FillPatch (int lev, amrex::Real time, const amrex::Vector& mf, bool fillset=true); -#if defined(ERF_USE_MOISTURE) // Compute a new MultiFab by copying from valid region and filling ghost cells - // works for single level and 2-level cases (fill fine grid ghost by interpolating from coarse) void FillPatchMoistVars (int lev, amrex::MultiFab& mf); -#endif // compute new multifabs by copying in data from valid region and filling ghost cells // works for single level and 2-level cases (fill fine grid ghost by interpolating from coarse) @@ -535,11 +529,8 @@ private: amrex::Vector rW_old; amrex::Vector rW_new; -#if defined(ERF_USE_MOISTURE) - std::string moisture_model = "NullMoist"; Microphysics micro; amrex::Vector qmoist; // This has 6 components: qv, qc, qi, qr, qs, qg -#endif #if defined(ERF_USE_RRTMGP) Radiation rad; @@ -659,28 +650,19 @@ private: amrex::Vector plot_var_names_1; amrex::Vector plot_var_names_2; - const amrex::Vector cons_names {"density", "rhotheta", "rhoKE", "rhoQKE", "rhoadv_0" -#if defined(ERF_USE_MOISTURE) - ,"rhoQt", "rhoQp" -#elif defined(ERF_USE_WARM_NO_PRECIP) - ,"rhoQv", "rhoQc" -#endif - }; + const amrex::Vector cons_names {"density", "rhotheta", "rhoKE", "rhoQKE", "rhoadv_0", + "rhoQt", "rhoQp"}; // Note that the order of variable names here must match the order in Derive.cpp - const amrex::Vector derived_names {"pressure", "soundspeed", "temp", "theta", "KE", "QKE", "scalar", - "pres_hse", "dens_hse", "pert_pres", "pert_dens", + const amrex::Vector derived_names {"soundspeed", "temp", "theta", "KE", "QKE", "scalar", + "pres_hse", "dens_hse", "pressure", "pert_pres", "pert_dens", "dpdx", "dpdy", "pres_hse_x", "pres_hse_y", "z_phys", "detJ" , "mapfac", // eddy viscosity "Kmv","Kmh", // eddy diffusivity of heat "Khv","Khh" -#if defined(ERF_USE_MOISTURE) ,"qt", "qp", "qv", "qc", "qi", "qrain", "qsnow", "qgraup" -#elif defined(ERF_USE_WARM_NO_PRECIP) - ,"qv", "qc" -#endif #ifdef ERF_COMPUTE_ERROR ,"xvel_err", "yvel_err", "zvel_err", "pp_err" #endif @@ -763,17 +745,14 @@ private: amrex::Vector h_havg_density; amrex::Vector h_havg_temperature; amrex::Vector h_havg_pressure; -#ifdef ERF_USE_MOISTURE amrex::Vector h_havg_qv; amrex::Vector h_havg_qc; -#endif + amrex::Gpu::DeviceVector d_havg_density; amrex::Gpu::DeviceVector d_havg_temperature; amrex::Gpu::DeviceVector d_havg_pressure; -#ifdef ERF_USE_MOISTURE amrex::Gpu::DeviceVector d_havg_qv; amrex::Gpu::DeviceVector d_havg_qc; -#endif void refinement_criteria_setup (); @@ -805,32 +784,26 @@ private: if ( (advChoice.dycore_horiz_adv_type == AdvType::Centered_6th) || (advChoice.dycore_vert_adv_type == AdvType::Centered_6th) -#if defined(ERF_USE_MOISTURE) or defined(ERF_USE_WARM_NO_PRECIP) || (advChoice.moistscal_horiz_adv_type == AdvType::Centered_6th) || (advChoice.moistscal_horiz_adv_type == AdvType::Centered_6th) -#endif || (advChoice.dryscal_horiz_adv_type == AdvType::Centered_6th) || (advChoice.dryscal_vert_adv_type == AdvType::Centered_6th) ) { return 3; } else if ( (advChoice.dycore_horiz_adv_type == AdvType::Upwind_5th) || (advChoice.dycore_vert_adv_type == AdvType::Upwind_5th) -#if defined(ERF_USE_MOISTURE) or defined(ERF_USE_WARM_NO_PRECIP) || (advChoice.moistscal_horiz_adv_type == AdvType::Upwind_5th) || (advChoice.moistscal_horiz_adv_type == AdvType::Upwind_5th) -#endif || (advChoice.dryscal_horiz_adv_type == AdvType::Upwind_5th) || (advChoice.dryscal_vert_adv_type == AdvType::Upwind_5th) ) { return 3; } else if ( (advChoice.dryscal_horiz_adv_type == AdvType::Weno_5) || (advChoice.dryscal_vert_adv_type == AdvType::Weno_5) -#if defined(ERF_USE_MOISTURE) or defined(ERF_USE_WARM_NO_PRECIP) || (advChoice.moistscal_horiz_adv_type == AdvType::Weno_5) || (advChoice.moistscal_vert_adv_type == AdvType::Weno_5) || (advChoice.moistscal_horiz_adv_type == AdvType::Weno_5Z) || (advChoice.moistscal_vert_adv_type == AdvType::Weno_5Z) -#endif || (advChoice.dryscal_horiz_adv_type == AdvType::Weno_5Z) || (advChoice.dryscal_vert_adv_type == AdvType::Weno_5Z) ) { return 3; } diff --git a/Source/ERF.cpp b/Source/ERF.cpp index ab5040d55..55f7229fb 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -157,9 +157,7 @@ ERF::ERF () vars_old[lev].resize(Vars::NumTypes); } -#if defined(ERF_USE_MOISTURE) qmoist.resize(nlevs_max); -#endif mri_integrator_mem.resize(nlevs_max); physbcs.resize(nlevs_max); @@ -546,12 +544,13 @@ ERF::InitData () restart(); -#ifdef ERF_USE_MOISTURE // Need to fill ghost cells here since we will use this qmoist in advance - for (int lev = 0; lev <= finest_level; lev++) { - FillPatchMoistVars(lev, qmoist[lev]); + if (solverChoice.moisture_type != MoistureType::None) + { + for (int lev = 0; lev <= finest_level; lev++) { + FillPatchMoistVars(lev, qmoist[lev]); + } } -#endif } if (input_bndry_planes) { @@ -578,21 +577,22 @@ ERF::InitData () } } -#ifdef ERF_USE_MOISTURE // Initialize microphysics here micro.define(solverChoice); // Call Init which will call Diagnose to fill qmoist - for (int lev = 0; lev <= finest_level; ++lev) + if (solverChoice.moisture_type != MoistureType::None) { - // If not restarting we need to fill qmoist given qt and qp. - if (restart_chkfile.empty()) { - micro.Init(vars_new[lev][Vars::cons], qmoist[lev], - grids[lev], Geom(lev), 0.0); // dummy value, not needed just to diagnose - micro.Update(vars_new[lev][Vars::cons], qmoist[lev]); + for (int lev = 0; lev <= finest_level; ++lev) + { + // If not restarting we need to fill qmoist given qt and qp. + if (restart_chkfile.empty()) { + micro.Init(vars_new[lev][Vars::cons], qmoist[lev], + grids[lev], Geom(lev), 0.0); // dummy value, not needed just to diagnose + micro.Update(vars_new[lev][Vars::cons], qmoist[lev]); + } } } -#endif // Configure ABLMost params if used MostWall boundary condition // NOTE: we must set up the MOST routine before calling WritePlotFile because @@ -889,9 +889,7 @@ ERF::init_only (int lev, Real time) init_from_hse(lev); } -#if defined(ERF_USE_MOISTURE) qmoist[lev].setVal(0.); -#endif // Add problem-specific flow features // @@ -907,6 +905,8 @@ ERF::init_only (int lev, Real time) lev_new[Vars::xvel].OverrideSync(geom[lev].periodicity()); lev_new[Vars::yvel].OverrideSync(geom[lev].periodicity()); lev_new[Vars::zvel].OverrideSync(geom[lev].periodicity()); + + print_state(vars_new[0][Vars::cons], IntVect(30,30,2)); } // read in some parameters from inputs file @@ -960,18 +960,17 @@ ERF::ReadParameters () particleData.init_particle_params(); #endif -#ifdef ERF_USE_MOISTURE // What type of moisture model to use - pp.query("moisture_model", moisture_model); - if (moisture_model == "SAM") { + if (solverChoice.moisture_type == MoistureType::SAM) { micro.SetModel(); - } else if (moisture_model == "Kessler") { + } else if (solverChoice.moisture_type == MoistureType::Kessler) { micro.SetModel(); - } else { + } else if (solverChoice.moisture_type == MoistureType::None) { micro.SetModel(); amrex::Print() << "WARNING: Compiled with moisture but using NullMoist model!\n"; + } else { + amrex::Abort("Dont know this moisture_type!") ; } -#endif // If this is set, it must be even if (fixed_mri_dt_ratio > 0 && (fixed_mri_dt_ratio%2 != 0) ) @@ -1152,76 +1151,85 @@ ERF::MakeHorizontalAverages () MultiFab mf(grids[lev], dmap[lev], 5, 0); -#if defined(ERF_USE_MOISTURE) - MultiFab qv(qmoist[lev], make_alias, 0, 1); -#endif + int zdir = 2; + auto domain = geom[0].Domain(); + + bool use_moisture = (solverChoice.moisture_type != MoistureType::None); for (MFIter mfi(mf); mfi.isValid(); ++mfi) { const Box& bx = mfi.validbox(); auto fab_arr = mf.array(mfi); - auto cons_arr = vars_new[lev][Vars::cons].array(mfi); -#if defined(ERF_USE_MOISTURE) - auto qv_arr = qv.array(mfi); -#endif + auto const cons_arr = vars_new[lev][Vars::cons].const_array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Real dens = cons_arr(i, j, k, Cons::Rho); fab_arr(i, j, k, 0) = dens; fab_arr(i, j, k, 1) = cons_arr(i, j, k, Cons::RhoTheta) / dens; -#if defined(ERF_USE_MOISTURE) - fab_arr(i, j, k, 2) = getPgivenRTh(cons_arr(i, j, k, Cons::RhoTheta), qv_arr(i,j,k)); -#else - fab_arr(i, j, k, 2) = getPgivenRTh(cons_arr(i, j, k, Cons::RhoTheta)); -#endif -#if defined(ERF_USE_MOISTURE) - fab_arr(i, j, k, 3) = cons_arr(i, j, k, Cons::RhoQt) / dens; - fab_arr(i, j, k, 4) = cons_arr(i, j, k, Cons::RhoQp) / dens; -#endif + if (!use_moisture) { + fab_arr(i, j, k, 2) = getPgivenRTh(cons_arr(i, j, k, Cons::RhoTheta)); + } }); } - int zdir = 2; - auto domain = geom[0].Domain(); + if (use_moisture) + { + MultiFab qv(qmoist[lev], make_alias, 0, 1); + + for (MFIter mfi(mf); mfi.isValid(); ++mfi) { + const Box& bx = mfi.validbox(); + auto fab_arr = mf.array(mfi); + auto const cons_arr = vars_new[lev][Vars::cons].const_array(mfi); + auto const qv_arr = qv.const_array(mfi); + + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + Real dens = cons_arr(i, j, k, Cons::Rho); + fab_arr(i, j, k, 2) = getPgivenRTh(cons_arr(i, j, k, Cons::RhoTheta), qv_arr(i,j,k)); + fab_arr(i, j, k, 3) = cons_arr(i, j, k, Cons::RhoQ1) / dens; + fab_arr(i, j, k, 4) = cons_arr(i, j, k, Cons::RhoQ2) / dens; + }); + } + + Gpu::HostVector h_avg_qv = sumToLine(mf,3,1,domain,zdir); + Gpu::HostVector h_avg_qc = sumToLine(mf,4,1,domain,zdir); + } // Sum in the horizontal plane Gpu::HostVector h_avg_density = sumToLine(mf,0,1,domain,zdir); Gpu::HostVector h_avg_temperature = sumToLine(mf,1,1,domain,zdir); Gpu::HostVector h_avg_pressure = sumToLine(mf,2,1,domain,zdir); -#ifdef ERF_USE_MOISTURE - Gpu::HostVector h_avg_qv = sumToLine(mf,3,1,domain,zdir); - Gpu::HostVector h_avg_qc = sumToLine(mf,4,1,domain,zdir); -#endif // Divide by the total number of cells we are averaging over - int size_z = domain.length(zdir); + int size_z = domain.length(zdir); Real area_z = static_cast(domain.length(0)*domain.length(1)); int klen = static_cast(h_avg_density.size()); + for (int k = 0; k < klen; ++k) { h_havg_density[k] /= area_z; h_havg_temperature[k] /= area_z; h_havg_pressure[k] /= area_z; -#if defined(ERF_USE_MOISTURE) - h_havg_qc[k] /= area_z; - h_havg_qv[k] /= area_z; -#endif - } + if (solverChoice.moisture_type != MoistureType::None) + { + h_havg_qc[k] /= area_z; + h_havg_qv[k] /= area_z; + } + } // k // resize device vectors d_havg_density.resize(size_z, 0.0_rt); d_havg_temperature.resize(size_z, 0.0_rt); d_havg_pressure.resize(size_z, 0.0_rt); -#if defined(ERF_USE_MOISTURE) - d_havg_qv.resize(size_z, 0.0_rt); - d_havg_qc.resize(size_z, 0.0_rt); -#endif // copy host vectors to device vectors amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_havg_density.begin(), h_havg_density.end(), d_havg_density.begin()); amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_havg_temperature.begin(), h_havg_temperature.end(), d_havg_temperature.begin()); amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_havg_pressure.begin(), h_havg_pressure.end(), d_havg_pressure.begin()); -#if defined(ERF_USE_MOISTURE) - amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_havg_qv.begin(), h_havg_qv.end(), d_havg_qv.begin()); - amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_havg_qc.begin(), h_havg_qc.end(), d_havg_qc.begin()); -#endif + + if (solverChoice.moisture_type != MoistureType::None) + { + d_havg_qv.resize(size_z, 0.0_rt); + d_havg_qc.resize(size_z, 0.0_rt); + amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_havg_qv.begin(), h_havg_qv.end(), d_havg_qv.begin()); + amrex::Gpu::copy(amrex::Gpu::hostToDevice, h_havg_qc.begin(), h_havg_qc.end(), d_havg_qc.begin()); + } } // Create horizontal average quantities for the MultiFab passed in @@ -1491,9 +1499,7 @@ ERF::ERF (const amrex::RealBox& rb, int max_level_in, rV_old.resize(nlevs_max); rW_old.resize(nlevs_max); -#if defined(ERF_USE_MOISTURE) qmoist.resize(nlevs_max); -#endif mri_integrator_mem.resize(nlevs_max); physbcs.resize(nlevs_max); diff --git a/Source/ERF_Tagging.cpp b/Source/ERF_Tagging.cpp index b065afd79..7682744fa 100644 --- a/Source/ERF_Tagging.cpp +++ b/Source/ERF_Tagging.cpp @@ -39,8 +39,6 @@ ERF::ErrorEst (int level, TagBoxArray& tags, Real time, int /*ngrow*/) auto& sfab = vars_new[level][Vars::cons][mfi]; if (ref_tags[j].Field() == "scalar") { derived::erf_derscalar(bx, dfab, 0, 1, sfab, Geom(level), time, nullptr, level); - } else if (ref_tags[j].Field() == "pressure") { - derived::erf_derpres(bx, dfab, 0, 1, sfab, Geom(level), time, nullptr, level); } else if (ref_tags[j].Field() == "theta") { derived::erf_dertheta(bx, dfab, 0, 1, sfab, Geom(level), time, nullptr, level); } diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index 5ba0f8dc2..ff62dd201 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -59,9 +59,7 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba, //******************************************************************************************** // Microphysics // ******************************************************************************************* -#if defined(ERF_USE_MOISTURE) qmoist[lev].define(ba, dm, 6, ngrow_state); // qv, qc, qi, qr, qs, qg -#endif // ******************************************************************************************** // Build the data structures for calculating diffusive/turbulent terms @@ -195,10 +193,8 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, //******************************************************************************************** // Microphysics // ******************************************************************************************* -#if defined(ERF_USE_MOISTURE) int ngrow_state = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_NumDiff) + 1; qmoist[lev].define(ba, dm, 6, ngrow_state); // qv, qc, qi, qr, qs, qg -#endif init_stuff(lev, ba, dm); @@ -300,9 +296,7 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp //******************************************************************************************** // Microphysics // ******************************************************************************************* -#if defined(ERF_USE_MOISTURE) qmoist[lev].define(ba, dm, 6, ngrow_state); // qv, qc, qi, qr, qs, qg -#endif init_stuff(lev,ba,dm); diff --git a/Source/IO/ERF_ReadBndryPlanes.H b/Source/IO/ERF_ReadBndryPlanes.H index 62aae4676..f2da53fa1 100644 --- a/Source/IO/ERF_ReadBndryPlanes.H +++ b/Source/IO/ERF_ReadBndryPlanes.H @@ -40,10 +40,9 @@ public: [[nodiscard]] int ingested_theta() const {return (is_temperature_read || is_theta_read);} [[nodiscard]] int ingested_density() const {return is_density_read;} [[nodiscard]] int ingested_scalar() const {return is_scalar_read;} -#if defined(ERF_USE_MOISTURE) - [[nodiscard]] int ingested_qt() const {return is_qt_read;} - [[nodiscard]] int ingested_qp() const {return is_qp_read;} -#elif defined(ERF_USE_WARM_NO_PRECIP) + [[nodiscard]] int ingested_q1() const {return is_q1_read;} + [[nodiscard]] int ingested_q2() const {return is_q2_read;} +#if defined(ERF_USE_WARM_NO_PRECIP) int ingested_qv() const {return is_qv_read;} int ingested_qc() const {return is_qc_read;} #endif @@ -101,10 +100,9 @@ private: int is_temperature_read; int is_theta_read; int is_scalar_read; -#if defined(ERF_USE_MOISTURE) - int is_qt_read; - int is_qp_read; -#elif defined(ERF_USE_WARM_NO_PRECIP) + int is_q1_read; + int is_q2_read; +#if defined(ERF_USE_WARM_NO_PRECIP) int is_qv_read; int is_qc_read; #endif diff --git a/Source/IO/ERF_ReadBndryPlanes.cpp b/Source/IO/ERF_ReadBndryPlanes.cpp index 4a2fcd8bc..f495e2ba6 100644 --- a/Source/IO/ERF_ReadBndryPlanes.cpp +++ b/Source/IO/ERF_ReadBndryPlanes.cpp @@ -155,10 +155,9 @@ ReadBndryPlanes::ReadBndryPlanes(const Geometry& geom, const Real& rdOcp_in) is_temperature_read = 0; is_theta_read = 0; is_scalar_read = 0; -#if defined(ERF_USE_MOISTURE) - is_qt_read = 0; - is_qp_read = 0; -#elif defined(ERF_USE_WARM_NO_PRECIP) + is_q1_read = 0; + is_q2_read = 0; +#if defined(ERF_USE_WARM_NO_PRECIP) is_qv_read = 0; is_qc_read = 0; #endif @@ -176,10 +175,9 @@ ReadBndryPlanes::ReadBndryPlanes(const Geometry& geom, const Real& rdOcp_in) if (m_var_names[i] == "temperature") is_temperature_read = 1; if (m_var_names[i] == "theta") is_theta_read = 1; if (m_var_names[i] == "scalar") is_scalar_read = 1; -#if defined(ERF_USE_MOISTURE) - if (m_var_names[i] == "qt") is_qt_read = 1; - if (m_var_names[i] == "qp") is_qp_read = 1; -#elif defined(ERF_USE_WARM_NO_PRECIP) + if (m_var_names[i] == "qt") is_q1_read = 1; + if (m_var_names[i] == "qp") is_q2_read = 1; +#if defined(ERF_USE_WARM_NO_PRECIP) if (m_var_names[i] == "qv") is_qv_read = 1; if (m_var_names[i] == "qc") is_qc_read = 1; #endif @@ -413,10 +411,9 @@ void ReadBndryPlanes::read_file(const int idx, Vector plot_var_names) }; // Note: All derived variables must be computed in order of "derived_names" defined in ERF.H - calculate_derived("pressure", derived::erf_derpres); calculate_derived("soundspeed", derived::erf_dersoundspeed); calculate_derived("temp", derived::erf_dertemp); calculate_derived("theta", derived::erf_dertheta); @@ -195,6 +194,28 @@ ERF::WritePlotFile (int which, Vector plot_var_names) MultiFab::Copy(mf[lev],r_hse,0,mf_comp,1,0); mf_comp += 1; } + + if (containerHasElement(plot_var_names, "pressure")) + { +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4& S_arr = vars_new[lev][Vars::cons].const_array(mfi); + const Array4& qv_arr = qmoist[0].const_array(mfi); + + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + Real qv_for_p = qv_arr(i,j,k); + const Real rhotheta = S_arr(i,j,k,RhoTheta_comp); + derdat(i, j, k, mf_comp) = getPgivenRTh(rhotheta,qv_for_p); + }); + } + mf_comp += 1; + } if (containerHasElement(plot_var_names, "pert_pres")) { #ifdef _OPENMP @@ -206,19 +227,11 @@ ERF::WritePlotFile (int which, Vector plot_var_names) const Array4& derdat = mf[lev].array(mfi); const Array4& p0_arr = p_hse.const_array(mfi); const Array4& S_arr = vars_new[lev][Vars::cons].const_array(mfi); -#if defined(ERF_USE_MOISTURE) const Array4 & qv_arr = qmoist[0].const_array(mfi); -#endif ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { -#if defined(ERF_USE_MOISTURE) Real qv_for_p = qv_arr(i,j,k); -#elif defined(ERF_USE_WARM_NO_PRECIP) - Real qv_for_p = S_arr(i,j,k,RhoQv_comp) / S_arr(i,j,k,Rho_comp); -#else - Real qv_for_p = 0.; -#endif const Real rhotheta = S_arr(i,j,k,RhoTheta_comp); derdat(i, j, k, mf_comp) = getPgivenRTh(rhotheta,qv_for_p) - p0_arr(i,j,k); }); @@ -580,10 +593,6 @@ ERF::WritePlotFile (int which, Vector plot_var_names) mf_comp ++; } -#if defined(ERF_USE_MOISTURE) - calculate_derived("qt", derived::erf_derQt); - calculate_derived("qp", derived::erf_derQp); - MultiFab qv_mf(qmoist[lev], make_alias, 0, 1); MultiFab qc_mf(qmoist[lev], make_alias, 1, 1); MultiFab qi_mf(qmoist[lev], make_alias, 2, 1); @@ -591,6 +600,22 @@ ERF::WritePlotFile (int which, Vector plot_var_names) MultiFab qs_mf(qmoist[lev], make_alias, 4, 1); MultiFab qg_mf(qmoist[lev], make_alias, 5, 1); + if (containerHasElement(plot_var_names, "qt")) + { + MultiFab::Copy(mf[lev],qv_mf,0,mf_comp,1,0); + MultiFab::Add (mf[lev],qc_mf,0,mf_comp,1,0); + MultiFab::Add (mf[lev],qi_mf,0,mf_comp,1,0); + mf_comp += 1; + } + + if (containerHasElement(plot_var_names, "qp")) + { + MultiFab::Copy(mf[lev],qr_mf,0,mf_comp,1,0); + MultiFab::Add (mf[lev],qs_mf,0,mf_comp,1,0); + MultiFab::Add (mf[lev],qg_mf,0,mf_comp,1,0); + mf_comp += 1; + } + if (containerHasElement(plot_var_names, "qv")) { MultiFab::Copy(mf[lev],qv_mf,0,mf_comp,1,0); @@ -626,10 +651,6 @@ ERF::WritePlotFile (int which, Vector plot_var_names) MultiFab::Copy(mf[lev],qg_mf,0,mf_comp,1,0); mf_comp += 1; } -#elif defined(ERF_USE_WARM_NO_PRECIP) - calculate_derived("qv", derived::erf_derQv); - calculate_derived("qc", derived::erf_derQc); -#endif #ifdef ERF_USE_PARTICLES if (containerHasElement(plot_var_names, "tracer_particle_count")) diff --git a/Source/IO/Plotfile.cpp.convergence b/Source/IO/Plotfile.cpp.convergence new file mode 100644 index 000000000..8396e370d --- /dev/null +++ b/Source/IO/Plotfile.cpp.convergence @@ -0,0 +1,1071 @@ +#include +#include +#include "AMReX_Interp_3D_C.H" +#include "AMReX_PlotFileUtil.H" +#include "TerrainMetrics.H" +#include "ERF_Constants.H" + +using namespace amrex; + +template +bool containerHasElement(const V& iterable, const T& query) { + return std::find(iterable.begin(), iterable.end(), query) != iterable.end(); +} + +void +ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector& plot_var_names) +{ + ParmParse pp(pp_prefix); + + if (pp.contains(pp_plot_var_names.c_str())) + { + std::string nm; + + int nPltVars = pp.countval(pp_plot_var_names.c_str()); + + for (int i = 0; i < nPltVars; i++) + { + pp.get(pp_plot_var_names.c_str(), nm, i); + + // Add the named variable to our list of plot variables + // if it is not already in the list + if (!containerHasElement(plot_var_names, nm)) { + plot_var_names.push_back(nm); + } + } + } else { + // + // The default is to add none of the variables to the list + // + plot_var_names.clear(); + } + + // Get state variables in the same order as we define them, + // since they may be in any order in the input list + Vector tmp_plot_names; + for (int i = 0; i < Cons::NumVars; ++i) { + if ( containerHasElement(plot_var_names, cons_names[i]) ) { + tmp_plot_names.push_back(cons_names[i]); + } + } + // check for velocity since it's not in cons_names + // if we are asked for any velocity component, we will need them all + if (containerHasElement(plot_var_names, "x_velocity") || + containerHasElement(plot_var_names, "y_velocity") || + containerHasElement(plot_var_names, "z_velocity")) { + tmp_plot_names.push_back("x_velocity"); + tmp_plot_names.push_back("y_velocity"); + tmp_plot_names.push_back("z_velocity"); + } + for (int i = 0; i < derived_names.size(); ++i) { + if ( containerHasElement(plot_var_names, derived_names[i]) ) { + if (solverChoice.use_terrain || (derived_names[i] != "z_phys" && derived_names[i] != "detJ") ) { + tmp_plot_names.push_back(derived_names[i]); + } + } + } + + // Check to see if we found all the requested variables + for (const auto& plot_name : plot_var_names) { + if (!containerHasElement(tmp_plot_names, plot_name)) { + if (amrex::ParallelDescriptor::IOProcessor()) { + Warning("\nWARNING: Requested to plot variable '" + plot_name + "' but it is not available"); + } + } + } + plot_var_names = tmp_plot_names; +} + +// set plotfile variable names +Vector +ERF::PlotFileVarNames ( Vector plot_var_names ) +{ + Vector names; + + names.insert(names.end(), plot_var_names.begin(), plot_var_names.end()); + + return names; + +} + +// Write plotfile to disk +void +ERF::WritePlotFile (int which, Vector plot_var_names) +{ + const Vector varnames = PlotFileVarNames(plot_var_names); + const int ncomp_mf = varnames.size(); + + // We fillpatch here because some of the derived quantities require derivatives + // which require ghost cells to be filled + for (int lev = 0; lev <= finest_level; ++lev) { + FillPatch(lev, t_new[lev], {&vars_new[lev][Vars::cons], &vars_new[lev][Vars::xvel], + &vars_new[lev][Vars::yvel], &vars_new[lev][Vars::zvel]}); + } + + if (ncomp_mf == 0) + return; + + Vector mf(finest_level+1); + for (int lev = 0; lev <= finest_level; ++lev) { + mf[lev].define(grids[lev], dmap[lev], ncomp_mf, 0); + } + + Vector mf_nd(finest_level+1); + if (solverChoice.use_terrain) { + for (int lev = 0; lev <= finest_level; ++lev) { + BoxArray nodal_grids(grids[lev]); nodal_grids.surroundingNodes(); + mf_nd[lev].define(nodal_grids, dmap[lev], 3, 0); + mf_nd[lev].setVal(0.); + } + } + + for (int lev = 0; lev <= finest_level; ++lev) { + int mf_comp = 0; + + // First, copy any of the conserved state variables into the output plotfile + AMREX_ALWAYS_ASSERT(cons_names.size() == Cons::NumVars); + for (int i = 0; i < Cons::NumVars; ++i) { + if (containerHasElement(plot_var_names, cons_names[i])) { + MultiFab::Copy(mf[lev],vars_new[lev][Vars::cons],i,mf_comp,1,0); + mf_comp++; + } + } + + // Next, check for velocities and if desired, output them -- note we output none or all, not just some + if (containerHasElement(plot_var_names, "x_velocity") || + containerHasElement(plot_var_names, "y_velocity") || + containerHasElement(plot_var_names, "z_velocity")) { + + average_face_to_cellcenter(mf[lev],mf_comp, + Array{&vars_new[lev][Vars::xvel],&vars_new[lev][Vars::yvel],&vars_new[lev][Vars::zvel]}); + mf_comp += AMREX_SPACEDIM; + } + + // Finally, check for any derived quantities and compute them, inserting + // them into our output multifab + auto calculate_derived = [&](const std::string& der_name, + decltype(derived::erf_dernull)& der_function) + { + if (containerHasElement(plot_var_names, der_name)) { + MultiFab dmf(mf[lev], make_alias, mf_comp, 1); +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(dmf, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + auto& dfab = dmf[mfi]; + auto& sfab = vars_new[lev][Vars::cons][mfi]; + der_function(bx, dfab, 0, 1, sfab, Geom(lev), t_new[0], nullptr, lev); + } + + mf_comp++; + } + }; + + // Note: All derived variables must be computed in order of "derived_names" defined in ERF.H + calculate_derived("pressure", derived::erf_derpres); + calculate_derived("soundspeed", derived::erf_dersoundspeed); + calculate_derived("temp", derived::erf_dertemp); + calculate_derived("theta", derived::erf_dertheta); + calculate_derived("KE", derived::erf_derKE); + calculate_derived("QKE", derived::erf_derQKE); + calculate_derived("scalar", derived::erf_derscalar); + + MultiFab r_hse(base_state[lev], make_alias, 0, 1); // r_0 is first component + MultiFab p_hse(base_state[lev], make_alias, 1, 1); // p_0 is second component + if (containerHasElement(plot_var_names, "pres_hse")) + { + // p_0 is second component of base_state + MultiFab::Copy(mf[lev],p_hse,0,mf_comp,1,0); + mf_comp += 1; + } + if (containerHasElement(plot_var_names, "dens_hse")) + { + // r_0 is first component of base_state + MultiFab::Copy(mf[lev],r_hse,0,mf_comp,1,0); + mf_comp += 1; + } + if (containerHasElement(plot_var_names, "pert_pres")) + { +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4& p0_arr = p_hse.const_array(mfi); + const Array4& S_arr = vars_new[lev][Vars::cons].const_array(mfi); + + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + + const Real rhotheta = S_arr(i,j,k,RhoTheta_comp); + derdat(i, j, k, mf_comp) = getPgivenRTh(rhotheta) - p0_arr(i,j,k); + }); + } + mf_comp += 1; + } + if (containerHasElement(plot_var_names, "pert_dens")) + { +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4& S_arr = vars_new[lev][Vars::cons].const_array(mfi); + const Array4& r0_arr = r_hse.const_array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + derdat(i, j, k, mf_comp) = S_arr(i,j,k,Rho_comp) - r0_arr(i,j,k); + }); + } + mf_comp ++; + } + + int klo = geom[lev].Domain().smallEnd(2); + int khi = geom[lev].Domain().bigEnd(2); + + if (containerHasElement(plot_var_names, "dpdx")) + { + auto dxInv = geom[lev].InvCellSizeArray(); + MultiFab pres(vars_new[lev][Vars::cons].boxArray(), vars_new[lev][Vars::cons].DistributionMap(), 1, 1); +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + // First define pressure on grown box + const Box& gbx = mfi.growntilebox(1); + const Array4 & p_arr = pres.array(mfi); + const Array4& S_arr = vars_new[lev][Vars::cons].const_array(mfi); + amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + p_arr(i,j,k) = getPgivenRTh(S_arr(i,j,k,RhoTheta_comp)); + }); + } + pres.FillBoundary(geom[lev].periodicity()); + +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + // Now compute pressure gradient on valid box + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4 & p_arr = pres.array(mfi); + + if (solverChoice.use_terrain) { + const Array4& z_nd = z_phys_nd[lev]->const_array(mfi); + + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + + // Pgrad at lower I face + Real met_h_xi_lo = Compute_h_xi_AtIface (i, j, k, dxInv, z_nd); + Real met_h_zeta_lo = Compute_h_zeta_AtIface(i, j, k, dxInv, z_nd); + Real gp_xi_lo = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k)); + Real gp_zeta_on_iface_lo; + if(k == klo) { + gp_zeta_on_iface_lo = 0.5 * dxInv[2] * ( + p_arr(i-1,j,k+1) + p_arr(i,j,k+1) + - p_arr(i-1,j,k ) - p_arr(i,j,k ) ); + } else if (k == khi) { + gp_zeta_on_iface_lo = 0.5 * dxInv[2] * ( + p_arr(i-1,j,k ) + p_arr(i,j,k ) + - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) ); + } else { + gp_zeta_on_iface_lo = 0.25 * dxInv[2] * ( + p_arr(i-1,j,k+1) + p_arr(i,j,k+1) + - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) ); + } + amrex::Real gpx_lo = gp_xi_lo - (met_h_xi_lo/ met_h_zeta_lo) * gp_zeta_on_iface_lo; + + // Pgrad at higher I face + Real met_h_xi_hi = Compute_h_xi_AtIface (i+1, j, k, dxInv, z_nd); + Real met_h_zeta_hi = Compute_h_zeta_AtIface(i+1, j, k, dxInv, z_nd); + Real gp_xi_hi = dxInv[0] * (p_arr(i+1,j,k) - p_arr(i,j,k)); + Real gp_zeta_on_iface_hi; + if(k == klo) { + gp_zeta_on_iface_hi = 0.5 * dxInv[2] * ( + p_arr(i+1,j,k+1) + p_arr(i,j,k+1) + - p_arr(i+1,j,k ) - p_arr(i,j,k ) ); + } else if (k == khi) { + gp_zeta_on_iface_hi = 0.5 * dxInv[2] * ( + p_arr(i+1,j,k ) + p_arr(i,j,k ) + - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) ); + } else { + gp_zeta_on_iface_hi = 0.25 * dxInv[2] * ( + p_arr(i+1,j,k+1) + p_arr(i,j,k+1) + - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) ); + } + amrex::Real gpx_hi = gp_xi_hi - (met_h_xi_hi/ met_h_zeta_hi) * gp_zeta_on_iface_hi; + + // Average P grad to CC + derdat(i ,j ,k, mf_comp) = 0.5 * (gpx_lo + gpx_hi); + }); + } else { + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + derdat(i ,j ,k, mf_comp) = 0.5 * (p_arr(i+1,j,k) - p_arr(i-1,j,k)) * dxInv[0]; + }); + } + } // mfi + mf_comp ++; + } // dpdx + + if (containerHasElement(plot_var_names, "dpdy")) + { + auto dxInv = geom[lev].InvCellSizeArray(); + + MultiFab pres(vars_new[lev][Vars::cons].boxArray(), vars_new[lev][Vars::cons].DistributionMap(), 1, 1); +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + // First define pressure on grown box + const Box& gbx = mfi.growntilebox(1); + const Array4 & p_arr = pres.array(mfi); + const Array4& S_arr = vars_new[lev][Vars::cons].const_array(mfi); + amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + p_arr(i,j,k) = getPgivenRTh(S_arr(i,j,k,RhoTheta_comp)); + }); + } + pres.FillBoundary(geom[lev].periodicity()); + +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + // Now compute pressure gradient on valid box + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4 & p_arr = pres.array(mfi); + + if (solverChoice.use_terrain) { + const Array4& z_nd = z_phys_nd[lev]->const_array(mfi); + + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + + Real met_h_eta_lo = Compute_h_eta_AtJface (i, j, k, dxInv, z_nd); + Real met_h_zeta_lo = Compute_h_zeta_AtJface(i, j, k, dxInv, z_nd); + Real gp_eta_lo = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k)); + Real gp_zeta_on_jface_lo; + if (k == klo) { + gp_zeta_on_jface_lo = 0.5 * dxInv[2] * ( + p_arr(i,j,k+1) + p_arr(i,j-1,k+1) + - p_arr(i,j,k ) - p_arr(i,j-1,k ) ); + } else if (k == khi) { + gp_zeta_on_jface_lo = 0.5 * dxInv[2] * ( + p_arr(i,j,k ) + p_arr(i,j-1,k ) + - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) ); + } else { + gp_zeta_on_jface_lo = 0.25 * dxInv[2] * ( + p_arr(i,j,k+1) + p_arr(i,j-1,k+1) + - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) ); + } + amrex::Real gpy_lo = gp_eta_lo - (met_h_eta_lo / met_h_zeta_lo) * gp_zeta_on_jface_lo; + + Real met_h_eta_hi = Compute_h_eta_AtJface (i, j+1, k, dxInv, z_nd); + Real met_h_zeta_hi = Compute_h_zeta_AtJface(i, j+1, k, dxInv, z_nd); + Real gp_eta_hi = dxInv[1] * (p_arr(i,j+1,k) - p_arr(i,j,k)); + Real gp_zeta_on_jface_hi; + if (k == klo) { + gp_zeta_on_jface_hi = 0.5 * dxInv[2] * ( + p_arr(i,j+1,k+1) + p_arr(i,j,k+1) + - p_arr(i,j+1,k ) - p_arr(i,j,k ) ); + } else if (k == khi) { + gp_zeta_on_jface_hi = 0.5 * dxInv[2] * ( + p_arr(i,j+1,k ) + p_arr(i,j,k ) + - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) ); + } else { + gp_zeta_on_jface_hi = 0.25 * dxInv[2] * ( + p_arr(i,j+1,k+1) + p_arr(i,j,k+1) + - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) ); + } + amrex::Real gpy_hi = gp_eta_hi - (met_h_eta_hi / met_h_zeta_hi) * gp_zeta_on_jface_hi; + + derdat(i ,j ,k, mf_comp) = 0.5 * (gpy_lo + gpy_hi); + }); + } else { + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + derdat(i ,j ,k, mf_comp) = 0.5 * (p_arr(i,j+1,k) - p_arr(i,j-1,k)) * dxInv[1]; + }); + } + } // mf + mf_comp ++; + } // dpdy + + if (containerHasElement(plot_var_names, "pres_hse_x")) + { + auto dxInv = geom[lev].InvCellSizeArray(); +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4& p_arr = p_hse.const_array(mfi); + + //USE_TERRAIN POSSIBLE ISSUE HERE + const Array4& z_nd = z_phys_nd[lev]->const_array(mfi); + + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + Real met_h_xi_lo = Compute_h_xi_AtIface (i, j, k, dxInv, z_nd); + Real met_h_zeta_lo = Compute_h_zeta_AtIface(i, j, k, dxInv, z_nd); + Real gp_xi_lo = dxInv[0] * (p_arr(i,j,k) - p_arr(i-1,j,k)); + Real gp_zeta_on_iface_lo; + if (k == klo) { + gp_zeta_on_iface_lo = 0.5 * dxInv[2] * ( + p_arr(i-1,j,k+1) + p_arr(i,j,k+1) + - p_arr(i-1,j,k ) - p_arr(i,j,k ) ); + } else if (k == khi) { + gp_zeta_on_iface_lo = 0.5 * dxInv[2] * ( + p_arr(i-1,j,k ) + p_arr(i,j,k ) + - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) ); + } else { + gp_zeta_on_iface_lo = 0.25 * dxInv[2] * ( + p_arr(i-1,j,k+1) + p_arr(i,j,k+1) + - p_arr(i-1,j,k-1) - p_arr(i,j,k-1) ); + } + amrex::Real gpx_lo = gp_xi_lo - (met_h_xi_lo/ met_h_zeta_lo) * gp_zeta_on_iface_lo; + + Real met_h_xi_hi = Compute_h_xi_AtIface (i+1, j, k, dxInv, z_nd); + Real met_h_zeta_hi = Compute_h_zeta_AtIface(i+1, j, k, dxInv, z_nd); + Real gp_xi_hi = dxInv[0] * (p_arr(i+1,j,k) - p_arr(i,j,k)); + Real gp_zeta_on_iface_hi; + if (k == klo) { + gp_zeta_on_iface_hi = 0.5 * dxInv[2] * ( + p_arr(i+1,j,k+1) + p_arr(i,j,k+1) + - p_arr(i+1,j,k ) - p_arr(i,j,k ) ); + } else if (k == khi) { + gp_zeta_on_iface_hi = 0.5 * dxInv[2] * ( + p_arr(i+1,j,k ) + p_arr(i,j,k ) + - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) ); + } else { + gp_zeta_on_iface_hi = 0.25 * dxInv[2] * ( + p_arr(i+1,j,k+1) + p_arr(i,j,k+1) + - p_arr(i+1,j,k-1) - p_arr(i,j,k-1) ); + } + amrex::Real gpx_hi = gp_xi_hi - (met_h_xi_hi/ met_h_zeta_hi) * gp_zeta_on_iface_hi; + + derdat(i ,j ,k, mf_comp) = 0.5 * (gpx_lo + gpx_hi); + }); + } + mf_comp += 1; + } // pres_hse_x + + if (containerHasElement(plot_var_names, "pres_hse_y")) + { + auto dxInv = geom[lev].InvCellSizeArray(); +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4& p_arr = p_hse.const_array(mfi); + const Array4& z_nd = z_phys_nd[lev]->const_array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + Real met_h_eta_lo = Compute_h_eta_AtJface (i, j, k, dxInv, z_nd); + Real met_h_zeta_lo = Compute_h_zeta_AtJface(i, j, k, dxInv, z_nd); + Real gp_eta_lo = dxInv[1] * (p_arr(i,j,k) - p_arr(i,j-1,k)); + Real gp_zeta_on_jface_lo; + if (k == klo) { + gp_zeta_on_jface_lo = 0.5 * dxInv[2] * ( + p_arr(i,j,k+1) + p_arr(i,j-1,k+1) + - p_arr(i,j,k ) - p_arr(i,j-1,k ) ); + } else if (k == khi) { + gp_zeta_on_jface_lo = 0.5 * dxInv[2] * ( + p_arr(i,j,k ) + p_arr(i,j-1,k ) + - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) ); + } else { + gp_zeta_on_jface_lo = 0.25 * dxInv[2] * ( + p_arr(i,j,k+1) + p_arr(i,j-1,k+1) + - p_arr(i,j,k-1) - p_arr(i,j-1,k-1) ); + } + amrex::Real gpy_lo = gp_eta_lo - (met_h_eta_lo / met_h_zeta_lo) * gp_zeta_on_jface_lo; + + Real met_h_eta_hi = Compute_h_eta_AtJface (i, j+1, k, dxInv, z_nd); + Real met_h_zeta_hi = Compute_h_zeta_AtJface(i, j+1, k, dxInv, z_nd); + Real gp_eta_hi = dxInv[1] * (p_arr(i,j+1,k) - p_arr(i,j,k)); + Real gp_zeta_on_jface_hi; + if (k == klo) { + gp_zeta_on_jface_hi = 0.5 * dxInv[2] * ( + p_arr(i,j+1,k+1) + p_arr(i,j,k+1) + - p_arr(i,j+1,k ) - p_arr(i,j,k ) ); + } else if (k == khi) { + gp_zeta_on_jface_hi = 0.5 * dxInv[2] * ( + p_arr(i,j+1,k ) + p_arr(i,j,k ) + - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) ); + } else { + gp_zeta_on_jface_hi = 0.25 * dxInv[2] * ( + p_arr(i,j+1,k+1) + p_arr(i,j,k+1) + - p_arr(i,j+1,k-1) - p_arr(i,j,k-1) ); + } + amrex::Real gpy_hi = gp_eta_hi - (met_h_eta_hi / met_h_zeta_hi) * gp_zeta_on_jface_hi; + + derdat(i ,j ,k, mf_comp) = 0.5 * (gpy_lo + gpy_hi); + }); + } + mf_comp += 1; + } // pres_hse_y + + if (solverChoice.use_terrain) { + if (containerHasElement(plot_var_names, "z_phys")) + { + MultiFab::Copy(mf[lev],*z_phys_cc[lev],0,mf_comp,1,0); + mf_comp ++; + } + + if (containerHasElement(plot_var_names, "detJ")) + { + MultiFab::Copy(mf[lev],*detJ_cc[lev],0,mf_comp,1,0); + mf_comp ++; + } + } // use_terrain + + if (containerHasElement(plot_var_names, "mapfac")) { +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4& mf_m = mapfac_m[lev]->array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + derdat(i ,j ,k, mf_comp) = mf_m(i,j,0); + }); + } + mf_comp ++; + } + +#if defined(ERF_USE_MOISTURE) + calculate_derived("qt", derived::erf_derQt); + calculate_derived("qp", derived::erf_derQp); + + MultiFab qv_mf(qmoist[lev], make_alias, 0, 1); + MultiFab qc_mf(qmoist[lev], make_alias, 1, 1); + MultiFab qi_mf(qmoist[lev], make_alias, 2, 1); + MultiFab qr_mf(qmoist[lev], make_alias, 3, 1); + MultiFab qs_mf(qmoist[lev], make_alias, 4, 1); + MultiFab qg_mf(qmoist[lev], make_alias, 5, 1); + + if (containerHasElement(plot_var_names, "qv")) + { + MultiFab::Copy(mf[lev],qv_mf,0,mf_comp,1,0); + mf_comp += 1; + } + + if (containerHasElement(plot_var_names, "qc")) + { + MultiFab::Copy(mf[lev],qc_mf,0,mf_comp,1,0); + mf_comp += 1; + } + + if (containerHasElement(plot_var_names, "qi")) + { + MultiFab::Copy(mf[lev],qi_mf,0,mf_comp,1,0); + mf_comp += 1; + } + + if (containerHasElement(plot_var_names, "qrain")) + { + MultiFab::Copy(mf[lev],qr_mf,0,mf_comp,1,0); + mf_comp += 1; + } + + if (containerHasElement(plot_var_names, "qsnow")) + { + MultiFab::Copy(mf[lev],qs_mf,0,mf_comp,1,0); + mf_comp += 1; + } + + if (containerHasElement(plot_var_names, "qgraup")) + { + MultiFab::Copy(mf[lev],qg_mf,0,mf_comp,1,0); + mf_comp += 1; + } +#elif defined(ERF_USE_WARM_NO_PRECIP) + calculate_derived("qv", derived::erf_derQv); + calculate_derived("qc", derived::erf_derQc); +#endif + +#ifdef ERF_COMPUTE_ERROR + // Next, check for error in velocities and if desired, output them -- note we output none or all, not just some + if (containerHasElement(plot_var_names, "xvel_err") || + containerHasElement(plot_var_names, "yvel_err") || + containerHasElement(plot_var_names, "zvel_err")) + { + // + // Moving terrain ANALYTICAL + // + Real H = geom[lev].ProbHi()[2]; + Real Ampl = 0.16; + Real wavelength = 100.; + Real kp = 2. * PI / wavelength; + Real g = CONST_GRAV; + Real omega = std::sqrt(g * kp); + Real omega_t = omega * t_new[lev]; + + const auto dx = geom[lev].CellSizeArray(); + +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.validbox(); + Box xbx(bx); xbx.surroundingNodes(0); + const Array4 xvel_arr = vars_new[lev][Vars::xvel].array(mfi); + const Array4 zvel_arr = vars_new[lev][Vars::zvel].array(mfi); + + const Array4& z_nd = z_phys_nd[lev]->const_array(mfi); + + ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + Real x = i * dx[0]; + Real z = 0.25 * (z_nd(i,j,k) + z_nd(i,j+1,k) + z_nd(i,j,k+1) + z_nd(i,j+1,k+1)); + + Real z_base = Ampl * std::sin(kp * x - omega_t); + z -= z_base; + + Real fac = std::cosh( kp * (z - H) ) / std::sinh(kp * H); + + xvel_arr(i,j,k) -= -Ampl * omega * fac * std::sin(kp * x - omega_t); + }); + + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + Real x = (i + 0.5) * dx[0]; + Real z = 0.25 * ( z_nd(i,j,k) + z_nd(i+1,j,k) + z_nd(i,j+1,k) + z_nd(i+1,j+1,k)); + + Real z_base = Ampl * std::sin(kp * x - omega_t); + z -= z_base; + + Real fac = std::sinh( kp * (z - H) ) / std::sinh(kp * H); + + zvel_arr(i,j,k) -= Ampl * omega * fac * std::cos(kp * x - omega_t); + }); + } + + MultiFab temp_mf(mf[lev].boxArray(), mf[lev].DistributionMap(), AMREX_SPACEDIM, 0); + average_face_to_cellcenter(temp_mf,0, + Array{&vars_new[lev][Vars::xvel],&vars_new[lev][Vars::yvel],&vars_new[lev][Vars::zvel]}); + + if (containerHasElement(plot_var_names, "xvel_err")) { + MultiFab::Copy(mf[lev],temp_mf,0,mf_comp,1,0); + mf_comp += 1; + } + if (containerHasElement(plot_var_names, "yvel_err")) { + MultiFab::Copy(mf[lev],temp_mf,1,mf_comp,1,0); + mf_comp += 1; + } + if (containerHasElement(plot_var_names, "zvel_err")) { + MultiFab::Copy(mf[lev],temp_mf,2,mf_comp,1,0); + mf_comp += 1; + } + + // Now restore the velocities to what they were +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(mf[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.validbox(); + Box xbx(bx); xbx.surroundingNodes(0); + + const Array4 xvel_arr = vars_new[lev][Vars::xvel].array(mfi); + const Array4 zvel_arr = vars_new[lev][Vars::zvel].array(mfi); + + const Array4& z_nd = z_phys_nd[lev]->const_array(mfi); + + ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + Real x = i * dx[0]; + Real z = 0.25 * (z_nd(i,j,k) + z_nd(i,j+1,k) + z_nd(i,j,k+1) + z_nd(i,j+1,k+1)); + Real z_base = Ampl * std::sin(kp * x - omega_t); + + z -= z_base; + + Real fac = std::cosh( kp * (z - H) ) / std::sinh(kp * H); + xvel_arr(i,j,k) += -Ampl * omega * fac * std::sin(kp * x - omega_t); + }); + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + Real x = (i + 0.5) * dx[0]; + Real z = 0.25 * ( z_nd(i,j,k) + z_nd(i+1,j,k) + z_nd(i,j+1,k) + z_nd(i+1,j+1,k)); + Real z_base = Ampl * std::sin(kp * x - omega_t); + + z -= z_base; + Real fac = std::sinh( kp * (z - H) ) / std::sinh(kp * H); + + zvel_arr(i,j,k) += Ampl * omega * fac * std::cos(kp * x - omega_t); + }); + } + } // end xvel_err, yvel_err, zvel_err + + if (containerHasElement(plot_var_names, "pp_err")) + { + // Moving terrain ANALYTICAL +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for ( MFIter mfi(mf[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const Array4& derdat = mf[lev].array(mfi); + const Array4& S_arr = vars_new[lev][Vars::cons].const_array(mfi); + + const auto dx = geom[lev].CellSizeArray(); +#if 0 + const Array4& z_nd = z_phys_nd[lev]->const_array(mfi); + const Array4& p0_arr = p_hse.const_array(mfi); + const Array4& r0_arr = r_hse.const_array(mfi); + + Real H = geom[lev].ProbHi()[2]; + Real Ampl = 0.16; + Real wavelength = 100.; + Real kp = 2. * PI / wavelength; + Real g = CONST_GRAV; + Real omega = std::sqrt(g * kp); + Real omega_t = omega * t_new[lev]; + + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + const Real rhotheta = S_arr(i,j,k,RhoTheta_comp); + derdat(i, j, k, mf_comp) = getPgivenRTh(rhotheta) - p0_arr(i,j,k); + + Real rho_hse = r0_arr(i,j,k); + + Real x = (i + 0.5) * dx[0]; + Real z = 0.125 * ( z_nd(i,j,k ) + z_nd(i+1,j,k ) + z_nd(i,j+1,k ) + z_nd(i+1,j+1,k ) + +z_nd(i,j,k+1) + z_nd(i+1,j,k+1) + z_nd(i,j+1,k+1) + z_nd(i+1,j+1,k+1) ); + Real z_base = Ampl * std::sin(kp * x - omega_t); + + z -= z_base; + Real fac = std::cosh( kp * (z - H) ) / std::sinh(kp * H); + Real pprime_exact = -(Ampl * omega * omega / kp) * fac * + std::sin(kp * x - omega_t) * r0_arr(i,j,k); + + derdat(i,j,k,mf_comp) -= pprime_exact; + }); +#endif + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + Real x = (i + 0.5) * dx[0]; + Real y = (j + 0.5) * dx[1]; + Real S_exact = cos(PI * x); + + derdat(i,j,k,mf_comp) = S_arr(i,j,k,RhoScalar_comp) - S_exact; + }); + } // mfi +#endif + mf_comp += 1; + + int nx = geom[lev].Domain().length(0); + int nv = std::sqrt(static_cast(geom[lev].Domain().numPts())); + Real dx = 2.0 / static_cast(nx); + Real norm_0 = mf[0].norm0(mf_comp-1,0); + Real norm_2 = mf[0].norm2(mf_comp-1) / nv; + amrex::Print() << " \n ERR: dx = " << dx << " " << norm_2 << std::endl; +// << " order = " << solverChoice.horiz_spatial_order +// << " 2-norm = " << norm_2 +// << " 0-norm = " << norm_0 << "\n" << std::endl; + } // pp_err + } // lev + + std::string plotfilename; + if (which == 1) + plotfilename = Concatenate(plot_file_1, istep[0], 5); + else if (which == 2) + plotfilename = Concatenate(plot_file_2, istep[0], 5); + + if (finest_level == 0) + { + if (plotfile_type == "amrex") { + amrex::Print() << "Writing plotfile " << plotfilename << "\n"; + if (solverChoice.use_terrain) { + // We started with mf_nd holding 0 in every component; here we fill only the offset in z + int lev = 0; + MultiFab::Copy(mf_nd[lev],*z_phys_nd[lev],0,2,1,0); + Real dz = Geom()[lev].CellSizeArray()[2]; + for (MFIter mfi(mf_nd[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const Box& bx = mfi.tilebox(); + Array4< Real> mf_arr = mf_nd[lev].array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + mf_arr(i,j,k,2) -= k * dz; + }); + } + WriteMultiLevelPlotfileWithTerrain(plotfilename, finest_level+1, + GetVecOfConstPtrs(mf), + GetVecOfConstPtrs(mf_nd), + varnames, + t_new[0], istep); + } else { + WriteMultiLevelPlotfile(plotfilename, finest_level+1, + GetVecOfConstPtrs(mf), + varnames, + Geom(), t_new[0], istep, refRatio()); + } + writeJobInfo(plotfilename); +#ifdef ERF_USE_HDF5 + } else if (plotfile_type == "hdf5" || plotfile_type == "HDF5") { + amrex::Print() << "Writing plotfile " << plotfilename+"d01.h5" << "\n"; + WriteMultiLevelPlotfileHDF5(plotfilename, finest_level+1, + GetVecOfConstPtrs(mf), + varnames, + Geom(), t_new[0], istep, refRatio()); +#endif +#ifdef ERF_USE_NETCDF + } else if (plotfile_type == "netcdf" || plotfile_type == "NetCDF") { + int lev = 0; + int l_which = 0; + writeNCPlotFile(lev, l_which, plotfilename, GetVecOfConstPtrs(mf), varnames, istep, t_new[0]); +#endif + } else { + amrex::Print() << "User specified plot_filetype = " << plotfile_type << std::endl; + amrex::Abort("Dont know this plot_filetype"); + } + + } else { // multilevel + + Vector r2(finest_level); + Vector g2(finest_level+1); + Vector mf2(finest_level+1); + + mf2[0].define(grids[0], dmap[0], ncomp_mf, 0); + + // Copy level 0 as is + MultiFab::Copy(mf2[0],mf[0],0,0,mf[0].nComp(),0); + + // Define a new multi-level array of Geometry's so that we pass the new "domain" at lev > 0 + Array periodicity = + {Geom()[0].isPeriodic(0),Geom()[0].isPeriodic(1),Geom()[0].isPeriodic(2)}; + g2[0].define(Geom()[0].Domain(),&(Geom()[0].ProbDomain()),0,periodicity.data()); + + if (plotfile_type == "amrex") { + r2[0] = IntVect(1,1,ref_ratio[0][0]); + for (int lev = 1; lev <= finest_level; ++lev) { + if (lev > 1) { + r2[lev-1][0] = 1; + r2[lev-1][1] = 1; + r2[lev-1][2] = r2[lev-2][2] * ref_ratio[lev-1][0]; + } + + mf2[lev].define(refine(grids[lev],r2[lev-1]), dmap[lev], ncomp_mf, 0); + + // Set the new problem domain + Box d2(Geom()[lev].Domain()); + d2.refine(r2[lev-1]); + + g2[lev].define(d2,&(Geom()[lev].ProbDomain()),0,periodicity.data()); + } + + // Do piecewise interpolation of mf into mf2 + for (int lev = 1; lev <= finest_level; ++lev) { +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(mf2[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const Box& bx = mfi.tilebox(); + pcinterp_interp(bx,mf2[lev].array(mfi), 0, mf[lev].nComp(), mf[lev].const_array(mfi),0,r2[lev-1]); + } + } + + // Define an effective ref_ratio which is isotropic to be passed into WriteMultiLevelPlotfile + Vector rr(finest_level); + for (int lev = 0; lev < finest_level; ++lev) { + rr[lev] = IntVect(ref_ratio[lev][0],ref_ratio[lev][1],ref_ratio[lev][0]); + } + + WriteMultiLevelPlotfile(plotfilename, finest_level+1, GetVecOfConstPtrs(mf2), varnames, + g2, t_new[0], istep, rr); + writeJobInfo(plotfilename); +#ifdef ERF_USE_NETCDF + } else if (plotfile_type == "netcdf" || plotfile_type == "NetCDF") { + for (int lev = 0; lev <= finest_level; ++lev) { + for (int which_box = 0; which_box < num_boxes_at_level[lev]; which_box++) { + writeNCPlotFile(lev, which_box, plotfilename, GetVecOfConstPtrs(mf), varnames, istep, t_new[0]); + } + } +#endif + } + } // end multi-level +} + +void +ERF::WriteMultiLevelPlotfileWithTerrain (const std::string& plotfilename, int nlevels, + const Vector& mf, + const Vector& mf_nd, + const Vector& varnames, + Real time, + const Vector& level_steps, + const std::string &versionName, + const std::string &levelPrefix, + const std::string &mfPrefix, + const Vector& extra_dirs) const +{ + BL_PROFILE("WriteMultiLevelPlotfileWithTerrain()"); + + BL_ASSERT(nlevels <= mf.size()); + BL_ASSERT(nlevels <= ref_ratio.size()+1); + BL_ASSERT(nlevels <= level_steps.size()); + BL_ASSERT(mf[0]->nComp() == varnames.size()); + + bool callBarrier(false); + PreBuildDirectorHierarchy(plotfilename, levelPrefix, nlevels, callBarrier); + if (!extra_dirs.empty()) { + for (const auto& d : extra_dirs) { + const std::string ed = plotfilename+"/"+d; + PreBuildDirectorHierarchy(ed, levelPrefix, nlevels, callBarrier); + } + } + ParallelDescriptor::Barrier(); + + if (ParallelDescriptor::MyProc() == ParallelDescriptor::NProcs()-1) { + Vector boxArrays(nlevels); + for(int level(0); level < boxArrays.size(); ++level) { + boxArrays[level] = mf[level]->boxArray(); + } + + auto f = [=]() { + VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size); + std::string HeaderFileName(plotfilename + "/Header"); + std::ofstream HeaderFile; + HeaderFile.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size()); + HeaderFile.open(HeaderFileName.c_str(), std::ofstream::out | + std::ofstream::trunc | + std::ofstream::binary); + if( ! HeaderFile.good()) FileOpenFailed(HeaderFileName); + WriteGenericPlotfileHeaderWithTerrain(HeaderFile, nlevels, boxArrays, varnames, + time, level_steps, versionName, + levelPrefix, mfPrefix); + }; + + if (AsyncOut::UseAsyncOut()) { + AsyncOut::Submit(std::move(f)); + } else { + f(); + } + } + + std::string mf_nodal_prefix = "Nu_nd"; + for (int level = 0; level <= finest_level; ++level) + { + if (AsyncOut::UseAsyncOut()) { + VisMF::AsyncWrite(*mf[level], + MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mfPrefix), + true); + VisMF::AsyncWrite(*mf_nd[level], + MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_nodal_prefix), + true); + } else { + const MultiFab* data; + std::unique_ptr mf_tmp; + if (mf[level]->nGrowVect() != 0) { + mf_tmp = std::make_unique(mf[level]->boxArray(), + mf[level]->DistributionMap(), + mf[level]->nComp(), 0, MFInfo(), + mf[level]->Factory()); + MultiFab::Copy(*mf_tmp, *mf[level], 0, 0, mf[level]->nComp(), 0); + data = mf_tmp.get(); + } else { + data = mf[level]; + } + VisMF::Write(*data , MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mfPrefix)); + VisMF::Write(*mf_nd[level], MultiFabFileFullPrefix(level, plotfilename, levelPrefix, mf_nodal_prefix)); + } + } +} + +void +ERF::WriteGenericPlotfileHeaderWithTerrain (std::ostream &HeaderFile, + int nlevels, + const Vector &bArray, + const Vector &varnames, + Real time, + const Vector &level_steps, + const std::string &versionName, + const std::string &levelPrefix, + const std::string &mfPrefix) const +{ + BL_ASSERT(nlevels <= bArray.size()); + BL_ASSERT(nlevels <= ref_ratio.size()+1); + BL_ASSERT(nlevels <= level_steps.size()); + + HeaderFile.precision(17); + + // ---- this is the generic plot file type name + HeaderFile << versionName << '\n'; + + HeaderFile << varnames.size() << '\n'; + + for (int ivar = 0; ivar < varnames.size(); ++ivar) { + HeaderFile << varnames[ivar] << "\n"; + } + HeaderFile << AMREX_SPACEDIM << '\n'; + HeaderFile << time << '\n'; + HeaderFile << finest_level << '\n'; + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + HeaderFile << geom[0].ProbLo(i) << ' '; + } + HeaderFile << '\n'; + for (int i = 0; i < AMREX_SPACEDIM; ++i) { + HeaderFile << geom[0].ProbHi(i) << ' '; + } + HeaderFile << '\n'; + for (int i = 0; i < finest_level; ++i) { + HeaderFile << ref_ratio[i][0] << ' '; + } + HeaderFile << '\n'; + for (int i = 0; i <= finest_level; ++i) { + HeaderFile << geom[i].Domain() << ' '; + } + HeaderFile << '\n'; + for (int i = 0; i <= finest_level; ++i) { + HeaderFile << level_steps[i] << ' '; + } + HeaderFile << '\n'; + for (int i = 0; i <= finest_level; ++i) { + for (int k = 0; k < AMREX_SPACEDIM; ++k) { + HeaderFile << geom[i].CellSize()[k] << ' '; + } + HeaderFile << '\n'; + } + HeaderFile << (int) geom[0].Coord() << '\n'; + HeaderFile << "0\n"; + + for (int level = 0; level <= finest_level; ++level) { + HeaderFile << level << ' ' << bArray[level].size() << ' ' << time << '\n'; + HeaderFile << level_steps[level] << '\n'; + + const IntVect& domain_lo = geom[level].Domain().smallEnd(); + for (int i = 0; i < bArray[level].size(); ++i) + { + // Need to shift because the RealBox ctor we call takes the + // physical location of index (0,0,0). This does not affect + // the usual cases where the domain index starts with 0. + const Box& b = shift(bArray[level][i], -domain_lo); + RealBox loc = RealBox(b, geom[level].CellSize(), geom[level].ProbLo()); + for (int n = 0; n < AMREX_SPACEDIM; ++n) { + HeaderFile << loc.lo(n) << ' ' << loc.hi(n) << '\n'; + } + } + + HeaderFile << MultiFabHeaderPath(level, levelPrefix, mfPrefix) << '\n'; + } + HeaderFile << "1" << "\n"; + HeaderFile << "3" << "\n"; + HeaderFile << "amrexvec_nu_x" << "\n"; + HeaderFile << "amrexvec_nu_y" << "\n"; + HeaderFile << "amrexvec_nu_z" << "\n"; + std::string mf_nodal_prefix = "Nu_nd"; + for (int level = 0; level <= finest_level; ++level) { + HeaderFile << MultiFabHeaderPath(level, levelPrefix, mf_nodal_prefix) << '\n'; + } +} diff --git a/Source/IndexDefines.H b/Source/IndexDefines.H index 9c90436a1..2bccf05dd 100644 --- a/Source/IndexDefines.H +++ b/Source/IndexDefines.H @@ -14,30 +14,19 @@ #define RhoKE_comp 2 // for Deardorff LES Model #define RhoQKE_comp 3 // for MYNN or YSU PBL Model #define RhoScalar_comp 4 -#if defined(ERF_USE_MOISTURE) - #define RhoQt_comp 5 - #define RhoQp_comp 6 - #define NVAR 7 -#elif defined(ERF_USE_WARM_NO_PRECIP) - #define RhoQv_comp 5 - #define RhoQc_comp 6 - #define NVAR 7 -#else - #define NVAR 5 -#endif +#define RhoQ1_comp 5 +#define RhoQ2_comp 6 + +#define NVAR 7 // Cell-centered primitive variables #define PrimTheta_comp (RhoTheta_comp -1) #define PrimKE_comp (RhoKE_comp -1) #define PrimQKE_comp (RhoQKE_comp -1) #define PrimScalar_comp (RhoScalar_comp-1) -#if defined(ERF_USE_MOISTURE) - #define PrimQt_comp (RhoQt_comp-1) - #define PrimQp_comp (RhoQp_comp-1) -#elif defined(ERF_USE_WARM_NO_PRECIP) - #define PrimQv_comp (RhoQv_comp-1) - #define PrimQc_comp (RhoQc_comp-1) -#endif +#define PrimQ1_comp (RhoQ1_comp-1) +#define PrimQ2_comp (RhoQ2_comp-1) + #define NUM_PRIM (NVAR-1) namespace BCVars { @@ -48,13 +37,8 @@ namespace BCVars { RhoKE_bc_comp, RhoQKE_bc_comp, RhoScalar_bc_comp, -#if defined(ERF_USE_MOISTURE) - RhoQt_bc_comp, - RhoQp_bc_comp, -#elif defined(ERF_USE_WARM_NO_PRECIP) - RhoQv_bc_comp, - RhoQc_bc_comp, -#endif + RhoQ1_bc_comp, + RhoQ2_bc_comp, xvel_bc = NVAR, yvel_bc, zvel_bc, @@ -106,33 +90,25 @@ namespace Cons { RhoKE, RhoQKE, RhoScalar, -#if defined(ERF_USE_MOISTURE) - RhoQt, - RhoQp, -#elif defined(ERF_USE_WARM_NO_PRECIP) - RhoQv, - RhoQc, -#endif + RhoQ1, + RhoQ2, NumVars }; } +#if 0 namespace Prim { enum { Theta = 0, KE, QKE, Scalar, -#if defined(ERF_USE_MOISTURE) - Qt, - Qp, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Qv, - Qc, -#endif + Q1, + Q2, NumVars }; } +#endif // We separate out horizontal and vertical turbulent diffusivities // These are the same for LES, but different for PBL models @@ -143,25 +119,15 @@ namespace EddyDiff { Scalar_h, KE_h, QKE_h, -#if defined(ERF_USE_MOISTURE) - Qt_h, - Qp_h, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Qv_h, - Qc_h, -#endif + Q1_h, + Q2_h, Mom_v, Theta_v, Scalar_v, KE_v, QKE_v, -#if defined(ERF_USE_MOISTURE) - Qt_v, - Qp_v, -#elif defined(ERF_USE_WARM_NO_PRECIP) - Qv_v, - Qc_v, -#endif + Q1_v, + Q2_v, PBL_lengthscale, NumDiffs }; diff --git a/Source/Initialization/ERF_init1d.cpp b/Source/Initialization/ERF_init1d.cpp index 2046e3954..35d343cb6 100644 --- a/Source/Initialization/ERF_init1d.cpp +++ b/Source/Initialization/ERF_init1d.cpp @@ -129,16 +129,12 @@ ERF::initHSE (int lev) MultiFab p_hse (base_state[lev], make_alias, 1, 1); // p_0 is second component MultiFab pi_hse(base_state[lev], make_alias, 2, 1); // pi_0 is third component -#ifdef ERF_USE_MOISTURE // Initial r_hse may or may not be in HSE -- defined in prob.cpp - if(solverChoice.use_moist_background){ + if (solverChoice.use_moist_background){ prob->erf_init_dens_hse_moist(r_hse, z_phys_nd[lev], geom[lev]); - }else{ + } else { prob->erf_init_dens_hse(r_hse, z_phys_nd[lev], z_phys_cc[lev], geom[lev]); } -#else - prob->erf_init_dens_hse(r_hse, z_phys_nd[lev], z_phys_cc[lev], geom[lev]); -#endif // This integrates up through column to update p_hse, pi_hse; // r_hse is not const b/c FillBoundary is called at the end for r_hse and p_hse diff --git a/Source/Initialization/ERF_init_bcs.cpp b/Source/Initialization/ERF_init_bcs.cpp index de6712bfd..7ae16d317 100644 --- a/Source/Initialization/ERF_init_bcs.cpp +++ b/Source/Initialization/ERF_init_bcs.cpp @@ -29,13 +29,8 @@ void ERF::init_bcs () m_bc_extdir_vals[BCVars::RhoKE_bc_comp][ori] = 0.0; m_bc_extdir_vals[BCVars::RhoQKE_bc_comp][ori] = 0.0; m_bc_extdir_vals[BCVars::RhoScalar_bc_comp][ori] = 0.0; -#if defined(ERF_USE_MOISTURE) - m_bc_extdir_vals[BCVars::RhoQt_bc_comp][ori] = 0.0; - m_bc_extdir_vals[BCVars::RhoQp_bc_comp][ori] = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - m_bc_extdir_vals[BCVars::RhoQv_bc_comp][ori] = 0.0; - m_bc_extdir_vals[BCVars::RhoQc_bc_comp][ori] = 0.0; -#endif + m_bc_extdir_vals[BCVars::RhoQ1_bc_comp][ori] = 0.0; + m_bc_extdir_vals[BCVars::RhoQ2_bc_comp][ori] = 0.0; m_bc_extdir_vals[BCVars::xvel_bc][ori] = 0.0; // default m_bc_extdir_vals[BCVars::yvel_bc][ori] = 0.0; @@ -47,13 +42,8 @@ void ERF::init_bcs () m_bc_neumann_vals[BCVars::RhoKE_bc_comp][ori] = 0.0; m_bc_neumann_vals[BCVars::RhoQKE_bc_comp][ori] = 0.0; m_bc_neumann_vals[BCVars::RhoScalar_bc_comp][ori] = 0.0; -#if defined(ERF_USE_MOISTURE) - m_bc_neumann_vals[BCVars::RhoQt_bc_comp][ori] = 0.0; - m_bc_neumann_vals[BCVars::RhoQp_bc_comp][ori] = 0.0; -#elif defined(ERF_USE_WARM_NO_PRECIP) - m_bc_neumann_vals[BCVars::RhoQv_bc_comp][ori] = 0.0; - m_bc_neumann_vals[BCVars::RhoQc_bc_comp][ori] = 0.0; -#endif + m_bc_neumann_vals[BCVars::RhoQ1_bc_comp][ori] = 0.0; + m_bc_neumann_vals[BCVars::RhoQ2_bc_comp][ori] = 0.0; m_bc_neumann_vals[BCVars::xvel_bc][ori] = 0.0; m_bc_neumann_vals[BCVars::yvel_bc][ori] = 0.0; m_bc_neumann_vals[BCVars::zvel_bc][ori] = 0.0; @@ -120,22 +110,23 @@ void ERF::init_bcs () if (pp.query("scalar", scalar_in)) m_bc_extdir_vals[BCVars::RhoScalar_bc_comp][ori] = rho_in*scalar_in; } -#if defined(ERF_USE_MOISTURE) + Real qt_in = 0.; - if (input_bndry_planes && m_r2d->ingested_qt()) { - m_bc_extdir_vals[BCVars::RhoQt_bc_comp][ori] = 0.; + if (input_bndry_planes && m_r2d->ingested_q1()) { + m_bc_extdir_vals[BCVars::RhoQ1_bc_comp][ori] = 0.; } else { if (pp.query("qt", qt_in)) - m_bc_extdir_vals[BCVars::RhoQt_bc_comp][ori] = rho_in*qt_in; + m_bc_extdir_vals[BCVars::RhoQ1_bc_comp][ori] = rho_in*qt_in; } Real qp_in = 0.; - if (input_bndry_planes && m_r2d->ingested_qp()) { - m_bc_extdir_vals[BCVars::RhoQp_bc_comp][ori] = 0.; + if (input_bndry_planes && m_r2d->ingested_q2()) { + m_bc_extdir_vals[BCVars::RhoQ2_bc_comp][ori] = 0.; } else { if (pp.query("qp", qp_in)) - m_bc_extdir_vals[BCVars::RhoQp_bc_comp][ori] = rho_in*qp_in; + m_bc_extdir_vals[BCVars::RhoQ2_bc_comp][ori] = rho_in*qp_in; } -#elif defined(ERF_USE_WARM_NO_PRECIP) + +#if defined(ERF_USE_WARM_NO_PRECIP) Real qv_in = 0.; if (input_bndry_planes && m_r2d->ingested_qv()) { m_bc_extdir_vals[BCVars::RhoQv_bc_comp][ori] = 0.; @@ -451,15 +442,9 @@ void ERF::init_bcs () ( (BCVars::cons_bc+i == BCVars::RhoTheta_bc_comp) && m_r2d->ingested_theta() ) || ( (BCVars::cons_bc+i == BCVars::RhoKE_bc_comp) && m_r2d->ingested_KE() ) || ( (BCVars::cons_bc+i == BCVars::RhoQKE_bc_comp) && m_r2d->ingested_QKE() ) || - ( (BCVars::cons_bc+i == BCVars::RhoScalar_bc_comp) && m_r2d->ingested_scalar() ) -#if defined(ERF_USE_MOISTURE) - || ( (BCVars::cons_bc+i == BCVars::RhoQt_bc_comp) && m_r2d->ingested_qt() ) - || ( (BCVars::cons_bc+i == BCVars::RhoQp_bc_comp) && m_r2d->ingested_qp() ) -#elif defined(ERF_USE_WARM_NO_PRECIP) - || ( (BCVars::cons_bc+i == BCVars::RhoQv_bc_comp) && m_r2d->ingested_qv() ) - || ( (BCVars::cons_bc+i == BCVars::RhoQc_bc_comp) && m_r2d->ingested_qc() ) -#endif - ) ) { + ( (BCVars::cons_bc+i == BCVars::RhoScalar_bc_comp) && m_r2d->ingested_scalar() ) || + ( (BCVars::cons_bc+i == BCVars::RhoQ1_bc_comp) && m_r2d->ingested_q1() ) || + ( (BCVars::cons_bc+i == BCVars::RhoQ2_bc_comp) && m_r2d->ingested_q2() )) ) { domain_bcs_type[BCVars::cons_bc+i].setLo(dir, ERFBCType::ext_dir_ingested); } } @@ -471,14 +456,9 @@ void ERF::init_bcs () ( (BCVars::cons_bc+i == BCVars::RhoTheta_bc_comp) && m_r2d->ingested_theta() ) || ( (BCVars::cons_bc+i == BCVars::RhoKE_bc_comp) && m_r2d->ingested_KE() ) || ( (BCVars::cons_bc+i == BCVars::RhoQKE_bc_comp) && m_r2d->ingested_QKE() ) || - ( (BCVars::cons_bc+i == BCVars::RhoScalar_bc_comp) && m_r2d->ingested_scalar() ) -#if defined(ERF_USE_MOISTURE) - || ( (BCVars::cons_bc+i == BCVars::RhoQt_bc_comp) && m_r2d->ingested_qt() ) - || ( (BCVars::cons_bc+i == BCVars::RhoQp_bc_comp) && m_r2d->ingested_qp() ) -#elif defined(ERF_USE_WARM_NO_PRECIP) - || ( (BCVars::cons_bc+i == BCVars::RhoQv_bc_comp) && m_r2d->ingested_qv() ) - || ( (BCVars::cons_bc+i == BCVars::RhoQc_bc_comp) && m_r2d->ingested_qc() ) -#endif + ( (BCVars::cons_bc+i == BCVars::RhoScalar_bc_comp) && m_r2d->ingested_scalar() ) || + ( (BCVars::cons_bc+i == BCVars::RhoQ1_bc_comp) && m_r2d->ingested_q1() ) || + ( (BCVars::cons_bc+i == BCVars::RhoQ2_bc_comp) && m_r2d->ingested_q2() ) ) ) { domain_bcs_type[BCVars::cons_bc+i].setHi(dir, ERFBCType::ext_dir_ingested); } diff --git a/Source/Initialization/ERF_init_custom.cpp b/Source/Initialization/ERF_init_custom.cpp index ce9a3a036..fcc158133 100644 --- a/Source/Initialization/ERF_init_custom.cpp +++ b/Source/Initialization/ERF_init_custom.cpp @@ -24,10 +24,9 @@ using namespace amrex; void ERF::init_custom (int lev) { - auto& lev_new = vars_new[lev]; -#if defined(ERF_USE_MOISTURE) - auto& qmoist_new = qmoist[lev]; -#endif + auto& lev_new = vars_new[lev]; + auto& qmoist_new = qmoist[lev]; + MultiFab r_hse(base_state[lev], make_alias, 0, 1); // r_0 is first component MultiFab p_hse(base_state[lev], make_alias, 1, 1); // p_0 is second component @@ -43,18 +42,12 @@ ERF::init_custom (int lev) yvel_pert.setVal(0.); zvel_pert.setVal(0.); -#if defined(ERF_USE_MOISTURE) MultiFab qmoist_pert(qmoist[lev].boxArray(), qmoist[lev].DistributionMap(), 3, qmoist[lev].nGrow()); qmoist_pert.setVal(0.); + MultiFab qv_pert(qmoist_pert, amrex::make_alias, 0, 1); MultiFab qc_pert(qmoist_pert, amrex::make_alias, 1, 1); MultiFab qi_pert(qmoist_pert, amrex::make_alias, 2, 1); -#elif defined(ERF_USE_WARM_NO_PRECIP) - MultiFab qmoist_pert(cons_pert.boxArray(), cons_pert.DistributionMap(), 2, cons_pert.nGrow()); - qmoist_pert.setVal(0.); - MultiFab qv_pert(qmoist_pert, amrex::make_alias, 0, 1); - MultiFab qc_pert(qmoist_pert, amrex::make_alias, 1, 1); -#endif int fix_random_seed = 0; ParmParse pp("erf"); pp.query("fix_random_seed", fix_random_seed); @@ -90,21 +83,12 @@ ERF::init_custom (int lev) Array4 r_hse_arr = r_hse.array(mfi); Array4 p_hse_arr = p_hse.array(mfi); -#if defined(ERF_USE_MOISTURE) const auto &qv_pert_arr = qv_pert.array(mfi); const auto &qc_pert_arr = qc_pert.array(mfi); const auto &qi_pert_arr = qi_pert.array(mfi); -#elif defined(ERF_USE_WARM_NO_PRECIP) - const auto &qv_pert_arr = qv_pert.array(mfi); - const auto &qc_pert_arr = qc_pert.array(mfi); -#endif prob->init_custom_pert(bx, xbx, ybx, zbx, cons_pert_arr, xvel_pert_arr, yvel_pert_arr, zvel_pert_arr, r_hse_arr, p_hse_arr, z_nd_arr, z_cc_arr, -#if defined(ERF_USE_MOISTURE) qv_pert_arr, qc_pert_arr, qi_pert_arr, -#elif defined(ERF_USE_WARM_NO_PRECIP) - qv_pert_arr, qc_pert_arr, -#endif geom[lev].data(), mf_m, mf_u, mf_v, solverChoice); } //mfi @@ -121,14 +105,16 @@ ERF::init_custom (int lev) MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQKE_comp, RhoQKE_comp, 1, cons_pert.nGrow()); } -#if defined(ERF_USE_MOISTURE) - MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQt_comp, RhoQt_comp, 1, cons_pert.nGrow()); - MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQp_comp, RhoQp_comp, 1, cons_pert.nGrow()); - MultiFab::Add( qmoist_new, qmoist_pert, 0, 0, 3, qmoist_pert.nGrow()); // qv, qc, qi -#elif defined(ERF_USE_WARM_NO_PRECIP) - MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQv_comp, RhoQv_comp, 1, cons_pert.nGrow()); - MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQc_comp, RhoQc_comp, 1, cons_pert.nGrow()); + if (solverChoice.moisture_type != MoistureType::None) { + MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQ1_comp, RhoQ1_comp, 1, cons_pert.nGrow()); + MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQ2_comp, RhoQ2_comp, 1, cons_pert.nGrow()); + MultiFab::Add( qmoist_new, qmoist_pert, 0, 0, 3, qmoist_pert.nGrow()); // qv, qc, qi +#if defined(ERF_USE_WARM_NO_PRECIP) + MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQv_comp, RhoQv_comp, 1, cons_pert.nGrow()); + MultiFab::Add(lev_new[Vars::cons], cons_pert, RhoQc_comp, RhoQc_comp, 1, cons_pert.nGrow()); #endif + } + MultiFab::Add(lev_new[Vars::xvel], xvel_pert, 0, 0, 1, xvel_pert.nGrowVect()); MultiFab::Add(lev_new[Vars::yvel], yvel_pert, 0, 0, 1, yvel_pert.nGrowVect()); MultiFab::Add(lev_new[Vars::zvel], zvel_pert, 0, 0, 1, zvel_pert.nGrowVect()); diff --git a/Source/Initialization/ERF_init_from_input_sounding.cpp b/Source/Initialization/ERF_init_from_input_sounding.cpp index 568d87e01..593f77478 100644 --- a/Source/Initialization/ERF_init_from_input_sounding.cpp +++ b/Source/Initialization/ERF_init_from_input_sounding.cpp @@ -143,9 +143,9 @@ init_bx_scalars_from_input_sounding (const amrex::Box &bx, // Set scalar = A_0*exp(-10r^2), where r is distance from center of domain state(i, j, k, RhoScalar_comp) = 0; - // total nonprecipitating water (Qt) == water vapor (Qv), i.e., there is no cloud water or cloud ice + // total nonprecipitating water (Q1) == water vapor (Qv), i.e., there is no cloud water or cloud ice #if defined(ERF_USE_MOISTURE) - state(i, j, k, RhoQt_comp) = rho_0 * interpolate_1d(z_inp_sound, qv_inp_sound, z, inp_sound_size); + state(i, j, k, RhoQ1_comp) = rho_0 * interpolate_1d(z_inp_sound, qv_inp_sound, z, inp_sound_size); #elif defined(ERF_USE_WARM_NO_PRECIP) state(i, j, k, RhoQv_comp) = rho_0 * interpolate_1d(z_inp_sound, qv_inp_sound, z, inp_sound_size); #endif @@ -235,9 +235,9 @@ init_bx_scalars_from_input_sounding_hse (const amrex::Box &bx, } #if defined(ERF_USE_MOISTURE) - // total nonprecipitating water (Qt) == water vapor (Qv), i.e., there + // total nonprecipitating water (Q1) == water vapor (Qv), i.e., there // is no cloud water or cloud ice - state(i, j, k, RhoQt_comp) = rho_k * interpolate_1d(z_inp_sound, qv_inp_sound, z, inp_sound_size); + state(i, j, k, RhoQ1_comp) = rho_k * interpolate_1d(z_inp_sound, qv_inp_sound, z, inp_sound_size); #elif defined(ERF_USE_WARM_NO_PRECIP) state(i, j, k, RhoQv_comp) = rho_k * interpolate_1d(z_inp_sound, qv_inp_sound, z, inp_sound_size); #endif diff --git a/Source/Initialization/ERF_init_from_metgrid.cpp b/Source/Initialization/ERF_init_from_metgrid.cpp index 6037c0de7..10be9dc82 100644 --- a/Source/Initialization/ERF_init_from_metgrid.cpp +++ b/Source/Initialization/ERF_init_from_metgrid.cpp @@ -772,7 +772,7 @@ init_state_from_metgrid (const Real l_rdOcp, int kmax = amrex::ubound(tbxc).z; #if defined(ERF_USE_MOISTURE) - int state_indx = RhoQt_comp; + int state_indx = RhoQ1_comp; #elif defined(ERF_USE_WARM_NO_PRECIP) int state_indx = RhoQv_comp; #endif @@ -829,7 +829,7 @@ init_base_state_from_metgrid (const Real l_rdOcp, const amrex::Array4& mask_c_arr) { #if defined(ERF_USE_MOISTURE) - int RhoQ_comp = RhoQt_comp; + int RhoQ_comp = RhoQ1_comp; #elif defined(ERF_USE_WARM_NO_PRECIP) int RhoQ_comp = RhoQv_comp; #endif @@ -894,7 +894,7 @@ init_base_state_from_metgrid (const Real l_rdOcp, { new_data(i,j,k,Rho_comp) = r_hse_arr(i,j,k); new_data(i,j,k,RhoScalar_comp) = 0.0; - // RhoTheta and RhoQt or RhoQv currently hold Theta and Qt or Qv. Multiply by Rho. + // RhoTheta and RhoQ1 or RhoQv currently hold Theta and Q1 or Qv. Multiply by Rho. Real RhoTheta = r_hse_arr(i,j,k)*new_data(i,j,k,RhoTheta_comp); new_data(i,j,k,RhoTheta_comp) = RhoTheta; #if defined(ERF_USE_MOISTURE) || defined(ERF_USE_WARM_NO_PRECIP) diff --git a/Source/Initialization/ERF_init_from_wrfinput.cpp b/Source/Initialization/ERF_init_from_wrfinput.cpp index d29fe922d..b2af0162b 100644 --- a/Source/Initialization/ERF_init_from_wrfinput.cpp +++ b/Source/Initialization/ERF_init_from_wrfinput.cpp @@ -318,12 +318,12 @@ init_state_from_wrfinput (int lev, state_fab.template copy(NC_rhotheta_fab[idx], 0, RhoTheta_comp, 1); #if defined(ERF_USE_MOISTURE) - state_fab.template copy(NC_QVAPOR_fab[idx], 0, RhoQt_comp, 1); - state_fab.template plus(NC_QCLOUD_fab[idx], 0, RhoQt_comp, 1); - state_fab.template mult(NC_rho_fab[idx] , 0, RhoQt_comp, 1); + state_fab.template copy(NC_QVAPOR_fab[idx], 0, RhoQ1_comp, 1); + state_fab.template plus(NC_QCLOUD_fab[idx], 0, RhoQ1_comp, 1); + state_fab.template mult(NC_rho_fab[idx] , 0, RhoQ1_comp, 1); - state_fab.template copy(NC_QRAIN_fab[idx], 0, RhoQp_comp, 1); - state_fab.template mult(NC_rho_fab[idx] , 0, RhoQp_comp, 1); + state_fab.template copy(NC_QRAIN_fab[idx], 0, RhoQ2_comp, 1); + state_fab.template mult(NC_rho_fab[idx] , 0, RhoQ2_comp, 1); # elif defined(ERF_USE_WARM_NO_PRECIP) #endif } // idx diff --git a/Source/Microphysics/Kessler/Init_Kessler.cpp b/Source/Microphysics/Kessler/Init_Kessler.cpp index 0b97652cf..f8a3d835f 100644 --- a/Source/Microphysics/Kessler/Init_Kessler.cpp +++ b/Source/Microphysics/Kessler/Init_Kessler.cpp @@ -137,8 +137,8 @@ void Kessler::Init (const MultiFab& cons_in, MultiFab& qmoist, amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { rho_array(i,j,k) = states_array(i,j,k,Rho_comp); theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp); - qt_array(i,j,k) = states_array(i,j,k,RhoQt_comp)/states_array(i,j,k,Rho_comp); - qp_array(i,j,k) = states_array(i,j,k,RhoQp_comp)/states_array(i,j,k,Rho_comp); + qt_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp); + qp_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp); qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k); qv_array(i,j,k) = qv_array_from_moist(i,j,k); temp_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)); diff --git a/Source/Microphysics/Kessler/Update_Kessler.cpp b/Source/Microphysics/Kessler/Update_Kessler.cpp index a32d44452..304c812a3 100644 --- a/Source/Microphysics/Kessler/Update_Kessler.cpp +++ b/Source/Microphysics/Kessler/Update_Kessler.cpp @@ -42,8 +42,8 @@ void Kessler::Update (amrex::MultiFab& cons, amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { states_arr(i,j,k,Rho_comp) = rho_arr(i,j,k); states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k); - states_arr(i,j,k,RhoQt_comp) = rho_arr(i,j,k)*qt_arr(i,j,k); - states_arr(i,j,k,RhoQp_comp) = rho_arr(i,j,k)*qp_arr(i,j,k); + states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qt_arr(i,j,k); + states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qp_arr(i,j,k); // Graupel == precip total - rain - snow (but must be >= 0) qgraup_arr(i,j,k) = 0.0;// diff --git a/Source/Microphysics/SAM/Cloud.cpp b/Source/Microphysics/SAM/Cloud_SAM.cpp similarity index 95% rename from Source/Microphysics/SAM/Cloud.cpp rename to Source/Microphysics/SAM/Cloud_SAM.cpp index 6ef1cc15b..9001a62b8 100644 --- a/Source/Microphysics/SAM/Cloud.cpp +++ b/Source/Microphysics/SAM/Cloud_SAM.cpp @@ -62,9 +62,6 @@ void SAM::Cloud () { erf_qsati(tabs1, pres1d_t(k), qsatt2); qsatt = om*qsatt1 + (1.-om)*qsatt2; } -//if(i==2 && j==2) -// printf("cloud: %d, %d, %d, %13.6f, %13.6f, %13.6f, %13.6f, %13.6f, %13.6f\n", -// i,j,k,tabs1,pres1d_t(k),qt_array(i,j,k),qsatt,qn_array(i,j,k),qp_array(i,j,k)); int niter; Real dtabs, lstarn, dlstarn, omp, lstarp, dlstarp, fff, dfff, dqsat; diff --git a/Source/Microphysics/SAM/Diagnose.cpp b/Source/Microphysics/SAM/Diagnose_SAM.cpp similarity index 100% rename from Source/Microphysics/SAM/Diagnose.cpp rename to Source/Microphysics/SAM/Diagnose_SAM.cpp diff --git a/Source/Microphysics/SAM/Init.cpp b/Source/Microphysics/SAM/Init_SAM.cpp similarity index 98% rename from Source/Microphysics/SAM/Init.cpp rename to Source/Microphysics/SAM/Init_SAM.cpp index 385eff7f8..ad5f0c484 100644 --- a/Source/Microphysics/SAM/Init.cpp +++ b/Source/Microphysics/SAM/Init_SAM.cpp @@ -134,8 +134,8 @@ void SAM::Init (const MultiFab& cons_in, MultiFab& qmoist, amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { rho_array(i,j,k) = states_array(i,j,k,Rho_comp); theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp); - qt_array(i,j,k) = states_array(i,j,k,RhoQt_comp)/states_array(i,j,k,Rho_comp); - qp_array(i,j,k) = states_array(i,j,k,RhoQp_comp)/states_array(i,j,k,Rho_comp); + qt_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp); + qp_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp); qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k); temp_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp)); pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp))/100.; diff --git a/Source/Microphysics/SAM/Make.package b/Source/Microphysics/SAM/Make.package index f55996515..11c88cf81 100644 --- a/Source/Microphysics/SAM/Make.package +++ b/Source/Microphysics/SAM/Make.package @@ -1,7 +1,7 @@ -CEXE_sources += Cloud.cpp -CEXE_sources += Init.cpp -CEXE_sources += Diagnose.cpp -CEXE_sources += Update.cpp +CEXE_sources += Cloud_SAM.cpp +CEXE_sources += Init_SAM.cpp +CEXE_sources += Diagnose_SAM.cpp +CEXE_sources += Update_SAM.cpp CEXE_sources += IceFall.cpp CEXE_sources += Precip.cpp CEXE_sources += PrecipFall.cpp diff --git a/Source/Microphysics/SAM/Update.cpp b/Source/Microphysics/SAM/Update_SAM.cpp similarity index 95% rename from Source/Microphysics/SAM/Update.cpp rename to Source/Microphysics/SAM/Update_SAM.cpp index 3c99d2802..7e391d1dd 100644 --- a/Source/Microphysics/SAM/Update.cpp +++ b/Source/Microphysics/SAM/Update_SAM.cpp @@ -44,8 +44,8 @@ void SAM::Update (amrex::MultiFab& cons, amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { states_arr(i,j,k,Rho_comp) = rho_arr(i,j,k); states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k); - states_arr(i,j,k,RhoQt_comp) = rho_arr(i,j,k)*qt_arr(i,j,k); - states_arr(i,j,k,RhoQp_comp) = rho_arr(i,j,k)*qp_arr(i,j,k); + states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qt_arr(i,j,k); + states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qp_arr(i,j,k); // Graupel == precip total - rain - snow (but must be >= 0) qgraup_arr(i,j,k) = std::max(0.0, qp_arr(i,j,k)-qpl_arr(i,j,k)-qpi_arr(i,j,k)); diff --git a/Source/TimeIntegration/ERF_fast_rhs_MT.cpp b/Source/TimeIntegration/ERF_fast_rhs_MT.cpp index 9b41d7426..28675be13 100644 --- a/Source/TimeIntegration/ERF_fast_rhs_MT.cpp +++ b/Source/TimeIntegration/ERF_fast_rhs_MT.cpp @@ -259,8 +259,8 @@ void erf_fast_rhs_MT (int step, int nrk, gpx *= mf_u(i,j,0); #if defined(ERF_USE_MOISTURE) - Real q = 0.5 * ( prim(i,j,k,PrimQt_comp) + prim(i-1,j,k,PrimQt_comp) - +prim(i,j,k,PrimQp_comp) + prim(i-1,j,k,PrimQp_comp) ); + Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i-1,j,k,PrimQ1_comp) + +prim(i,j,k,PrimQ2_comp) + prim(i-1,j,k,PrimQ2_comp) ); gpx /= (1.0 + q); #elif defined(ERF_USE_WARM_NO_PRECIP) Real q = 0.5 * ( prim(i,j,k,PrimQv_comp) + prim(i-1,j,k,PrimQv_comp) @@ -289,8 +289,8 @@ void erf_fast_rhs_MT (int step, int nrk, gpy *= mf_v(i,j,0); #if defined(ERF_USE_MOISTURE) - Real q = 0.5 * ( prim(i,j,k,PrimQt_comp) + prim(i,j-1,k,PrimQt_comp) - +prim(i,j,k,PrimQp_comp) + prim(i,j-1,k,PrimQp_comp) ); + Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i,j-1,k,PrimQ1_comp) + +prim(i,j,k,PrimQ2_comp) + prim(i,j-1,k,PrimQ2_comp) ); gpy /= (1.0 + q); #elif defined(ERF_USE_WARM_NO_PRECIP) Real q = 0.5 * ( prim(i,j,k,PrimQv_comp) + prim(i,j-1,k,PrimQv_comp) @@ -402,8 +402,8 @@ void erf_fast_rhs_MT (int step, int nrk, Real coeff_Q = coeffQ_a(i,j,k); #if defined(ERF_USE_MOISTURE) - Real q = 0.5 * ( prim(i,j,k,PrimQt_comp) + prim(i,j,k-1,PrimQt_comp) - +prim(i,j,k,PrimQp_comp) + prim(i,j,k-1,PrimQp_comp) ); + Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i,j,k-1,PrimQ1_comp) + +prim(i,j,k,PrimQ2_comp) + prim(i,j,k-1,PrimQ2_comp) ); coeff_P /= (1.0 + q); coeff_Q /= (1.0 + q); #elif defined(ERF_USE_WARM_NO_PRECIP) diff --git a/Source/TimeIntegration/ERF_fast_rhs_N.cpp b/Source/TimeIntegration/ERF_fast_rhs_N.cpp index 4f27f325b..729bf0d33 100644 --- a/Source/TimeIntegration/ERF_fast_rhs_N.cpp +++ b/Source/TimeIntegration/ERF_fast_rhs_N.cpp @@ -223,8 +223,8 @@ void erf_fast_rhs_N (int step, int nrk, gpx *= mf_u(i,j,0); #if defined(ERF_USE_MOISTURE) - Real q = 0.5 * ( prim(i,j,k,PrimQt_comp) + prim(i-1,j,k,PrimQt_comp) - +prim(i,j,k,PrimQp_comp) + prim(i-1,j,k,PrimQp_comp) ); + Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i-1,j,k,PrimQ1_comp) + +prim(i,j,k,PrimQ2_comp) + prim(i-1,j,k,PrimQ2_comp) ); gpx /= (1.0 + q); #elif defined(ERF_USE_WARM_NO_PRECIP) Real q = 0.5 * ( prim(i,j,k,PrimQv_comp) + prim(i-1,j,k,PrimQv_comp) @@ -249,8 +249,8 @@ void erf_fast_rhs_N (int step, int nrk, gpy *= mf_v(i,j,0); #if defined(ERF_USE_MOISTURE) - Real q = 0.5 * ( prim(i,j,k,PrimQt_comp) + prim(i,j-1,k,PrimQt_comp) - +prim(i,j,k,PrimQp_comp) + prim(i,j-1,k,PrimQp_comp) ); + Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i,j-1,k,PrimQ1_comp) + +prim(i,j,k,PrimQ2_comp) + prim(i,j-1,k,PrimQ2_comp) ); gpy /= (1.0 + q); #elif defined(ERF_USE_WARM_NO_PRECIP) Real q = 0.5 * ( prim(i,j,k,PrimQv_comp) + prim(i,j-1,k,PrimQv_comp) @@ -383,8 +383,8 @@ void erf_fast_rhs_N (int step, int nrk, Real coeff_Q = coeffQ_a(i,j,k); #if defined(ERF_USE_MOISTURE) - Real q = 0.5 * ( prim(i,j,k,PrimQt_comp) + prim(i,j,k-1,PrimQt_comp) - +prim(i,j,k,PrimQp_comp) + prim(i,j,k-1,PrimQp_comp) ); + Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i,j,k-1,PrimQ1_comp) + +prim(i,j,k,PrimQ2_comp) + prim(i,j,k-1,PrimQ2_comp) ); coeff_P /= (1.0 + q); coeff_Q /= (1.0 + q); #elif defined(ERF_USE_WARM_NO_PRECIP) diff --git a/Source/TimeIntegration/ERF_fast_rhs_T.cpp b/Source/TimeIntegration/ERF_fast_rhs_T.cpp index 40aea842a..b1ac58d6b 100644 --- a/Source/TimeIntegration/ERF_fast_rhs_T.cpp +++ b/Source/TimeIntegration/ERF_fast_rhs_T.cpp @@ -258,8 +258,8 @@ void erf_fast_rhs_T (int step, int nrk, gpx *= mf_u(i,j,0); #if defined(ERF_USE_MOISTURE) - Real q = 0.5 * ( prim(i,j,k,PrimQt_comp) + prim(i-1,j,k,PrimQt_comp) - +prim(i,j,k,PrimQp_comp) + prim(i-1,j,k,PrimQp_comp) ); + Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i-1,j,k,PrimQ1_comp) + +prim(i,j,k,PrimQ2_comp) + prim(i-1,j,k,PrimQ2_comp) ); gpx /= (1.0 + q); #elif defined(ERF_USE_WARM_NO_PRECIP) Real q = 0.5 * ( prim(i,j,k,PrimQv_comp) + prim(i-1,j,k,PrimQv_comp) @@ -291,8 +291,8 @@ void erf_fast_rhs_T (int step, int nrk, gpy *= mf_v(i,j,0); #if defined(ERF_USE_MOISTURE) - Real q = 0.5 * ( prim(i,j,k,PrimQt_comp) + prim(i,j-1,k,PrimQt_comp) - +prim(i,j,k,PrimQp_comp) + prim(i,j-1,k,PrimQp_comp) ); + Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i,j-1,k,PrimQ1_comp) + +prim(i,j,k,PrimQ2_comp) + prim(i,j-1,k,PrimQ2_comp) ); gpy /= (1.0 + q); #elif defined(ERF_USE_WARM_NO_PRECIP) Real q = 0.5 * ( prim(i,j,k,PrimQv_comp) + prim(i,j-1,k,PrimQv_comp) @@ -469,8 +469,8 @@ void erf_fast_rhs_T (int step, int nrk, Real coeff_Q = coeffQ_a(i,j,k); #if defined(ERF_USE_MOISTURE) - Real q = 0.5 * ( prim(i,j,k,PrimQt_comp) + prim(i,j,k-1,PrimQt_comp) - +prim(i,j,k,PrimQp_comp) + prim(i,j,k-1,PrimQp_comp) ); + Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i,j,k-1,PrimQ1_comp) + +prim(i,j,k,PrimQ2_comp) + prim(i,j,k-1,PrimQ2_comp) ); coeff_P /= (1.0 + q); coeff_Q /= (1.0 + q); #elif defined(ERF_USE_WARM_NO_PRECIP) diff --git a/Source/TimeIntegration/ERF_make_buoyancy.cpp b/Source/TimeIntegration/ERF_make_buoyancy.cpp index b5c320621..e0e8aec42 100644 --- a/Source/TimeIntegration/ERF_make_buoyancy.cpp +++ b/Source/TimeIntegration/ERF_make_buoyancy.cpp @@ -173,9 +173,9 @@ void make_buoyancy (Vector& S_data, amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Real rhop_hi = cell_data(i,j,k ,Rho_comp) * (1.0 + qv_data(i,j,k ) + qc_data(i,j,k ) - + qi_data(i,j,k )) + cell_data(i,j,k,RhoQp_comp) - r0_arr(i,j,k ); + + qi_data(i,j,k )) + cell_data(i,j,k,RhoQ2_comp) - r0_arr(i,j,k ); Real rhop_lo = cell_data(i,j,k-1,Rho_comp) * (1.0 + qv_data(i,j,k-1) + qc_data(i,j,k-1) - + qi_data(i,j,k-1)) + cell_data(i,j,k-1,RhoQp_comp) - r0_arr(i,j,k-1); + + qi_data(i,j,k-1)) + cell_data(i,j,k-1,RhoQ2_comp) - r0_arr(i,j,k-1); buoyancy_fab(i, j, k) = grav_gpu[2] * 0.5 * ( rhop_hi + rhop_lo ); }); } // mfi @@ -205,7 +205,7 @@ void make_buoyancy (Vector& S_data, if (solverChoice.buoyancy_type == 2 || solverChoice.buoyancy_type == 4 ) { - prim_ave.line_average(PrimQp_comp, qp_h); + prim_ave.line_average(PrimQ2_comp, qp_h); Gpu::DeviceVector qp_d(ncell); @@ -251,13 +251,13 @@ void make_buoyancy (Vector& S_data, qplus = 0.61* ( qv_data(i,j,k)-qv_d_ptr[k]) - (qc_data(i,j,k)-qc_d_ptr[k]+ qi_data(i,j,k)-qi_d_ptr[k]+ - cell_prim(i,j,k,PrimQp_comp)-qp_d_ptr[k]) + cell_prim(i,j,k,PrimQ2_comp)-qp_d_ptr[k]) + (tempp3d-tempp1d)/tempp1d*(Real(1.0) + Real(0.61)*qv_d_ptr[k]-qc_d_ptr[k]-qi_d_ptr[k]-qp_d_ptr[k]); qminus = 0.61 *( qv_data(i,j,k-1)-qv_d_ptr[k-1]) - (qc_data(i,j,k-1)-qc_d_ptr[k-1]+ qi_data(i,j,k-1)-qi_d_ptr[k-1]+ - cell_prim(i,j,k-1,PrimQp_comp)-qp_d_ptr[k-1]) + cell_prim(i,j,k-1,PrimQ2_comp)-qp_d_ptr[k-1]) + (tempm3d-tempm1d)/tempm1d*(Real(1.0) + Real(0.61)*qv_d_ptr[k-1]-qi_d_ptr[k-1]-qc_d_ptr[k-1]-qp_d_ptr[k-1]); } @@ -265,13 +265,13 @@ void make_buoyancy (Vector& S_data, qplus = 0.61* ( qv_data(i,j,k)-qv_d_ptr[k]) - (qc_data(i,j,k)-qc_d_ptr[k]+ qi_data(i,j,k)-qi_d_ptr[k]+ - cell_prim(i,j,k,PrimQp_comp)-qp_d_ptr[k]) + cell_prim(i,j,k,PrimQ2_comp)-qp_d_ptr[k]) + (cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp) - theta_d_ptr[k ])/theta_d_ptr[k ]; qminus = 0.61 *( qv_data(i,j,k-1)-qv_d_ptr[k-1]) - (qc_data(i,j,k-1)-qc_d_ptr[k-1]+ qi_data(i,j,k-1)-qi_d_ptr[k-1]+ - cell_prim(i,j,k-1,PrimQp_comp)-qp_d_ptr[k-1]) + cell_prim(i,j,k-1,PrimQ2_comp)-qp_d_ptr[k-1]) + (cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) - theta_d_ptr[k-1])/theta_d_ptr[k-1]; } @@ -313,10 +313,10 @@ void make_buoyancy (Vector& S_data, Real tempp3d = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), cell_data(i,j,k ,RhoTheta_comp)); Real tempm3d = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp)); - Real qplus = 0.61 * qv_data(i,j,k) - (qc_data(i,j,k)+ qi_data(i,j,k)+ cell_prim(i,j,k,PrimQp_comp)) + Real qplus = 0.61 * qv_data(i,j,k) - (qc_data(i,j,k)+ qi_data(i,j,k)+ cell_prim(i,j,k,PrimQ2_comp)) + (tempp3d-tempp1d)/tempp1d; - Real qminus = 0.61 *qv_data(i,j,k-1) - (qc_data(i,j,k-1)+ qi_data(i,j,k-1)+ cell_prim(i,j,k-1,PrimQp_comp)) + Real qminus = 0.61 *qv_data(i,j,k-1) - (qc_data(i,j,k-1)+ qi_data(i,j,k-1)+ cell_prim(i,j,k-1,PrimQ2_comp)) + (tempm3d-tempm1d)/tempm1d; Real qavg = Real(0.5) * (qplus + qminus); diff --git a/Source/TimeIntegration/ERF_make_fast_coeffs.cpp b/Source/TimeIntegration/ERF_make_fast_coeffs.cpp index 311868e38..0871b4fb0 100644 --- a/Source/TimeIntegration/ERF_make_fast_coeffs.cpp +++ b/Source/TimeIntegration/ERF_make_fast_coeffs.cpp @@ -136,8 +136,8 @@ void make_fast_coeffs (int /*level*/, coeffQ_a(i,j,k) = coeff_Q; #if defined(ERF_USE_MOISTURE) - Real q = 0.5 * ( prim(i,j,k,PrimQt_comp) + prim(i,j,k-1,PrimQt_comp) - +prim(i,j,k,PrimQp_comp) + prim(i,j,k-1,PrimQp_comp) ); + Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i,j,k-1,PrimQ1_comp) + +prim(i,j,k,PrimQ2_comp) + prim(i,j,k-1,PrimQ2_comp) ); coeff_P /= (1.0 + q); coeff_Q /= (1.0 + q); #elif defined(ERF_USE_WARM_NO_PRECIP) @@ -185,8 +185,8 @@ void make_fast_coeffs (int /*level*/, coeffQ_a(i,j,k) = coeff_Q; #if defined(ERF_USE_MOISTURE) - Real q = 0.5 * ( prim(i,j,k,PrimQt_comp) + prim(i,j,k-1,PrimQt_comp) - +prim(i,j,k,PrimQp_comp) + prim(i,j,k-1,PrimQp_comp) ); + Real q = 0.5 * ( prim(i,j,k,PrimQ1_comp) + prim(i,j,k-1,PrimQ1_comp) + +prim(i,j,k,PrimQ2_comp) + prim(i,j,k-1,PrimQ2_comp) ); coeff_P /= (1.0 + q); coeff_Q /= (1.0 + q); #elif defined(ERF_USE_WARM_NO_PRECIP) diff --git a/Source/TimeIntegration/ERF_slow_rhs_post.cpp b/Source/TimeIntegration/ERF_slow_rhs_post.cpp index ff1b961ae..b8128f6d9 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_post.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_post.cpp @@ -295,7 +295,7 @@ void erf_slow_rhs_post (int level, int finest_level, horiz_adv_type, vert_adv_type, l_use_terrain, flx_arr); #ifdef ERF_USE_MOISTURE - start_comp = RhoQt_comp; + start_comp = RhoQ1_comp; num_comp = 2; AdvType moist_horiz_adv_type = ac.moistscal_horiz_adv_type; @@ -413,7 +413,7 @@ void erf_slow_rhs_post (int level, int finest_level, bx_ylo, bx_yhi); int icomp; #if defined(ERF_USE_MOISTURE) - icomp = RhoQt_comp; + icomp = RhoQ1_comp; #elif defined(ERF_USE_WARM_NO_PRECIP) icomp = RhoQv_comp; #endif diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index dd80373af..b406a7028 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -804,8 +804,8 @@ void erf_slow_rhs_pre (int level, int finest_level, Real q = 0.0; #if defined(ERF_USE_MOISTURE) - q = 0.5 * ( cell_prim(i,j,k,PrimQt_comp) + cell_prim(i-1,j,k,PrimQt_comp) - +cell_prim(i,j,k,PrimQp_comp) + cell_prim(i-1,j,k,PrimQp_comp) ); + q = 0.5 * ( cell_prim(i,j,k,PrimQ1_comp) + cell_prim(i-1,j,k,PrimQ1_comp) + +cell_prim(i,j,k,PrimQ2_comp) + cell_prim(i-1,j,k,PrimQ2_comp) ); #elif defined(ERF_USE_WARM_NO_PRECIP) q = 0.5 * ( cell_prim(i,j,k,PrimQv_comp) + cell_prim(i-1,j,k,PrimQv_comp) +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i-1,j,k,PrimQc_comp) ); @@ -848,8 +848,8 @@ void erf_slow_rhs_pre (int level, int finest_level, Real q = 0.0; #if defined(ERF_USE_MOISTURE) - q = 0.5 * ( cell_prim(i,j,k,PrimQt_comp) + cell_prim(i-1,j,k,PrimQt_comp) - +cell_prim(i,j,k,PrimQp_comp) + cell_prim(i-1,j,k,PrimQp_comp) ); + q = 0.5 * ( cell_prim(i,j,k,PrimQ1_comp) + cell_prim(i-1,j,k,PrimQ1_comp) + +cell_prim(i,j,k,PrimQ2_comp) + cell_prim(i-1,j,k,PrimQ2_comp) ); #elif defined(ERF_USE_WARM_NO_PRECIP) q = 0.5 * ( cell_prim(i,j,k,PrimQv_comp) + cell_prim(i-1,j,k,PrimQv_comp) +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i-1,j,k,PrimQc_comp) ); @@ -912,8 +912,8 @@ void erf_slow_rhs_pre (int level, int finest_level, Real q = 0.0; #if defined(ERF_USE_MOISTURE) - q = 0.5 * ( cell_prim(i,j,k,PrimQt_comp) + cell_prim(i,j-1,k,PrimQt_comp) - +cell_prim(i,j,k,PrimQp_comp) + cell_prim(i,j-1,k,PrimQp_comp) ); + q = 0.5 * ( cell_prim(i,j,k,PrimQ1_comp) + cell_prim(i,j-1,k,PrimQ1_comp) + +cell_prim(i,j,k,PrimQ2_comp) + cell_prim(i,j-1,k,PrimQ2_comp) ); #elif defined(ERF_USE_WARM_NO_PRECIP) q = 0.5 * ( cell_prim(i,j,k,PrimQv_comp) + cell_prim(i,j-1,k,PrimQv_comp) +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i,j-1,k,PrimQc_comp) ); @@ -954,8 +954,8 @@ void erf_slow_rhs_pre (int level, int finest_level, Real q = 0.0; #if defined(ERF_USE_MOISTURE) - q = 0.5 * ( cell_prim(i,j,k,PrimQt_comp) + cell_prim(i,j-1,k,PrimQt_comp) - +cell_prim(i,j,k,PrimQp_comp) + cell_prim(i,j-1,k,PrimQp_comp) ); + q = 0.5 * ( cell_prim(i,j,k,PrimQ1_comp) + cell_prim(i,j-1,k,PrimQ1_comp) + +cell_prim(i,j,k,PrimQ2_comp) + cell_prim(i,j-1,k,PrimQ2_comp) ); #elif defined(ERF_USE_WARM_NO_PRECIP) q = 0.5 * ( cell_prim(i,j,k,PrimQv_comp) + cell_prim(i,j-1,k,PrimQv_comp) +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i,j-1,k,PrimQc_comp) ); @@ -1008,8 +1008,8 @@ void erf_slow_rhs_pre (int level, int finest_level, Real q = 0.0; #if defined(ERF_USE_MOISTURE) - q = 0.5 * ( cell_prim(i,j,k,PrimQt_comp) + cell_prim(i,j,k-1,PrimQt_comp) - +cell_prim(i,j,k,PrimQp_comp) + cell_prim(i,j,k-1,PrimQp_comp) ); + q = 0.5 * ( cell_prim(i,j,k,PrimQ1_comp) + cell_prim(i,j,k-1,PrimQ1_comp) + +cell_prim(i,j,k,PrimQ2_comp) + cell_prim(i,j,k-1,PrimQ2_comp) ); #elif defined(ERF_USE_WARM_NO_PRECIP) q = 0.5 * ( cell_prim(i,j,k,PrimQv_comp) + cell_prim(i,j,k-1,PrimQv_comp) +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i,j,k-1,PrimQc_comp) ); @@ -1049,8 +1049,8 @@ void erf_slow_rhs_pre (int level, int finest_level, Real q = 0.0; #if defined(ERF_USE_MOISTURE) - q = 0.5 * ( cell_prim(i,j,k,PrimQt_comp) + cell_prim(i,j,k-1,PrimQt_comp) - +cell_prim(i,j,k,PrimQp_comp) + cell_prim(i,j,k-1,PrimQp_comp) ); + q = 0.5 * ( cell_prim(i,j,k,PrimQ1_comp) + cell_prim(i,j,k-1,PrimQ1_comp) + +cell_prim(i,j,k,PrimQ2_comp) + cell_prim(i,j,k-1,PrimQ2_comp) ); #elif defined(ERF_USE_WARM_NO_PRECIP) q = 0.5 * ( cell_prim(i,j,k,PrimQv_comp) + cell_prim(i,j,k-1,PrimQv_comp) +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i,j,k-1,PrimQc_comp) ); diff --git a/Source/main.cpp b/Source/main.cpp index fd96fcc9c..9c00f7ff0 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -36,7 +36,8 @@ void add_par () { // (rather than the amrex default of 1). pp.add("n_proper",2); - pp.add("max_grid_size",2048); + int max_grid_size = 2048; + pp.queryAdd("max_grid_size",max_grid_size); // This will set the default value of blocking_factor to be 1, but will allow // the user to override it in the inputs file or on command line diff --git a/Source/prob_common.H b/Source/prob_common.H index ca365830b..a90aebee6 100644 --- a/Source/prob_common.H +++ b/Source/prob_common.H @@ -93,14 +93,9 @@ public: amrex::Array4 const& /*p_hse*/, amrex::Array4 const& /*z_nd*/, amrex::Array4 const& /*z_cc*/, -#if defined(ERF_USE_MOISTURE) amrex::Array4 const& /*qv*/, amrex::Array4 const& /*qc*/, amrex::Array4 const& /*qi*/, -#elif defined(ERF_USE_WARM_NO_PRECIP) - amrex::Array4 const& /*qv*/, - amrex::Array4 const& /*qc*/, -#endif amrex::GeometryData const& /*geomdata*/, amrex::Array4 const& /*mf_m*/, amrex::Array4 const& /*mf_u*/,