diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake index 0a926b1bb..dc4888ff7 100644 --- a/CMake/BuildERFExe.cmake +++ b/CMake/BuildERFExe.cmake @@ -170,6 +170,7 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/SourceTerms/ERF_MoistSetRhs.cpp ${SRC_DIR}/SourceTerms/ERF_NumericalDiffusion.cpp ${SRC_DIR}/SourceTerms/ERF_ForestDrag.cpp + ${SRC_DIR}/SourceTerms/ERF_TerrainDrag.cpp ${SRC_DIR}/TimeIntegration/ERF_ComputeTimestep.cpp ${SRC_DIR}/TimeIntegration/ERF_Advance.cpp ${SRC_DIR}/TimeIntegration/ERF_TimeStep.cpp @@ -206,6 +207,9 @@ function(build_erf_lib erf_lib_name) ${SRC_DIR}/Microphysics/Kessler/ERF_InitKessler.cpp ${SRC_DIR}/Microphysics/Kessler/ERF_Kessler.cpp ${SRC_DIR}/Microphysics/Kessler/ERF_UpdateKessler.cpp + ${SRC_DIR}/Microphysics/SatAdj/ERF_InitSatAdj.cpp + ${SRC_DIR}/Microphysics/SatAdj/ERF_SatAdj.cpp + ${SRC_DIR}/Microphysics/SatAdj/ERF_UpdateSatAdj.cpp ${SRC_DIR}/WindFarmParametrization/Fitch/ERF_AdvanceFitch.cpp ${SRC_DIR}/WindFarmParametrization/EWP/ERF_AdvanceEWP.cpp ${SRC_DIR}/WindFarmParametrization/SimpleActuatorDisk/ERF_AdvanceSimpleAD.cpp @@ -254,6 +258,7 @@ endif() target_include_directories(${erf_lib_name} PUBLIC $) target_include_directories(${erf_lib_name} PUBLIC $) target_include_directories(${erf_lib_name} PUBLIC $) + target_include_directories(${erf_lib_name} PUBLIC $) target_include_directories(${erf_lib_name} PUBLIC $) target_include_directories(${erf_lib_name} PUBLIC $) target_include_directories(${erf_lib_name} PUBLIC $) diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst index 46ec4148a..abf362ad9 100644 --- a/Docs/sphinx_doc/Inputs.rst +++ b/Docs/sphinx_doc/Inputs.rst @@ -432,10 +432,12 @@ List of Parameters | **erf.no_substepping** | Should we turn off | int (0 or 1) | 0 | | | substepping in time? | | | +----------------------------+----------------------+----------------+-------------------+ -| **erf.cfl** | CFL number for | Real > 0 and | 0.8 | -| | hydro | <= 1 | | -| | | | | -| | | | | +| **erf.cfl** | CFL number used to | Real > 0 and | 0.8 | +| | compute level 0 dt | <= 1 | | ++----------------------------+----------------------+----------------+-------------------+ +| **erf.substepping_cfl** | CFL number used to | Real > 0 and | 1.0 | +| | compute the number | <= 1 | | +| | of substeps | | | +----------------------------+----------------------+----------------+-------------------+ | **erf.fixed_dt** | set level 0 dt | Real > 0 | unused if not | | | as this value | | set | @@ -1443,7 +1445,7 @@ List of Parameters | | | Values | | +=============================+==========================+====================+============+ | **erf.moisture_model** | Name of moisture model | "SAM", "Kessler", | "Null" | -| | | "FastEddy" | | +| | | "SatAdj" | | +-----------------------------+--------------------------+--------------------+------------+ | **erf.do_cloud** | use basic moisture model | true / false | true | +-----------------------------+--------------------------+--------------------+------------+ diff --git a/Docs/sphinx_doc/theory/Microphysics.rst b/Docs/sphinx_doc/theory/Microphysics.rst index b62d6704a..afe563eff 100644 --- a/Docs/sphinx_doc/theory/Microphysics.rst +++ b/Docs/sphinx_doc/theory/Microphysics.rst @@ -170,3 +170,8 @@ The evaporation rate of rain is .. math:: Q_{revp} = 2\pi(S-1)n_{0R}[0.78\lambda_{R}^{-2}+0.31S_{c}^{1/3}\Gamma[(b+5)/2]a^{1/2}\mu^{-1/2}(\frac{\rho_{0}}{\rho})^{1/4}\lambda_{R}^{(b+5)/2}](\frac{1}{\rho})(\frac{L_{v}^{2}}{K_{0}R_{w}T^{2}}+\frac{1}{\rho r_{s}\psi})^{-1} + + Saturation Adjustment Microphysics Model +---------------------------------- +The saturation adjustment microphysics model is the simplest possible moisture model and only transports the +water vapor mixing ratio, :math:`q_v`, and the cloud water mixing ration, :math:`q_c`. Evaporation, :math:`q_v \longrightarrow q_c`, and condensation, :math:`q_c \longrightarrow q_v`, are the only relevant mechanisms. The final saturation state, :math:`q_v = q_{vs}(T)` is obtained from Newton-Raphson iterations on the thermal temperature :math:`T`. diff --git a/Exec/ABL/erf_terrain_def b/Exec/ABL/erf_terrain_def new file mode 100644 index 000000000..1a649774c --- /dev/null +++ b/Exec/ABL/erf_terrain_def @@ -0,0 +1,1024 @@ +768 768 108 +784 768 108 +800 768 108 +816 768 108 +832 768 108 +848 768 108 +864 768 108 +880 768 108 +896 768 108 +912 768 108 +928 768 108 +944 768 108 +960 768 108 +976 768 108 +992 768 108 +1008 768 108 +1024 768 108 +1040 768 108 +1056 768 108 +1072 768 108 +1088 768 108 +1104 768 108 +1120 768 108 +1136 768 108 +1152 768 108 +1168 768 108 +1184 768 108 +1200 768 108 +1216 768 108 +1232 768 108 +1248 768 108 +1264 768 108 +768 784 108 +784 784 108 +800 784 108 +816 784 108 +832 784 108 +848 784 108 +864 784 108 +880 784 108 +896 784 108 +912 784 108 +928 784 108 +944 784 108 +960 784 108 +976 784 108 +992 784 108 +1008 784 108 +1024 784 108 +1040 784 108 +1056 784 108 +1072 784 108 +1088 784 108 +1104 784 108 +1120 784 108 +1136 784 108 +1152 784 108 +1168 784 108 +1184 784 108 +1200 784 108 +1216 784 108 +1232 784 108 +1248 784 108 +1264 784 108 +768 800 108 +784 800 108 +800 800 108 +816 800 108 +832 800 108 +848 800 108 +864 800 108 +880 800 108 +896 800 108 +912 800 108 +928 800 108 +944 800 108 +960 800 108 +976 800 108 +992 800 108 +1008 800 108 +1024 800 108 +1040 800 108 +1056 800 108 +1072 800 108 +1088 800 108 +1104 800 108 +1120 800 108 +1136 800 108 +1152 800 108 +1168 800 108 +1184 800 108 +1200 800 108 +1216 800 108 +1232 800 108 +1248 800 108 +1264 800 108 +768 816 108 +784 816 108 +800 816 108 +816 816 108 +832 816 108 +848 816 108 +864 816 108 +880 816 108 +896 816 108 +912 816 108 +928 816 108 +944 816 108 +960 816 108 +976 816 108 +992 816 108 +1008 816 108 +1024 816 108 +1040 816 108 +1056 816 108 +1072 816 108 +1088 816 108 +1104 816 108 +1120 816 108 +1136 816 108 +1152 816 108 +1168 816 108 +1184 816 108 +1200 816 108 +1216 816 108 +1232 816 108 +1248 816 108 +1264 816 108 +768 832 108 +784 832 108 +800 832 108 +816 832 108 +832 832 108 +848 832 108 +864 832 108 +880 832 108 +896 832 108 +912 832 108 +928 832 108 +944 832 108 +960 832 108 +976 832 108 +992 832 108 +1008 832 108 +1024 832 108 +1040 832 108 +1056 832 108 +1072 832 108 +1088 832 108 +1104 832 108 +1120 832 108 +1136 832 108 +1152 832 108 +1168 832 108 +1184 832 108 +1200 832 108 +1216 832 108 +1232 832 108 +1248 832 108 +1264 832 108 +768 848 108 +784 848 108 +800 848 108 +816 848 108 +832 848 108 +848 848 108 +864 848 108 +880 848 108 +896 848 108 +912 848 108 +928 848 108 +944 848 108 +960 848 108 +976 848 108 +992 848 108 +1008 848 108 +1024 848 108 +1040 848 108 +1056 848 108 +1072 848 108 +1088 848 108 +1104 848 108 +1120 848 108 +1136 848 108 +1152 848 108 +1168 848 108 +1184 848 108 +1200 848 108 +1216 848 108 +1232 848 108 +1248 848 108 +1264 848 108 +768 864 108 +784 864 108 +800 864 108 +816 864 108 +832 864 108 +848 864 108 +864 864 108 +880 864 108 +896 864 108 +912 864 108 +928 864 108 +944 864 108 +960 864 108 +976 864 108 +992 864 108 +1008 864 108 +1024 864 108 +1040 864 108 +1056 864 108 +1072 864 108 +1088 864 108 +1104 864 108 +1120 864 108 +1136 864 108 +1152 864 108 +1168 864 108 +1184 864 108 +1200 864 108 +1216 864 108 +1232 864 108 +1248 864 108 +1264 864 108 +768 880 108 +784 880 108 +800 880 108 +816 880 108 +832 880 108 +848 880 108 +864 880 108 +880 880 108 +896 880 108 +912 880 108 +928 880 108 +944 880 108 +960 880 108 +976 880 108 +992 880 108 +1008 880 108 +1024 880 108 +1040 880 108 +1056 880 108 +1072 880 108 +1088 880 108 +1104 880 108 +1120 880 108 +1136 880 108 +1152 880 108 +1168 880 108 +1184 880 108 +1200 880 108 +1216 880 108 +1232 880 108 +1248 880 108 +1264 880 108 +768 896 108 +784 896 108 +800 896 108 +816 896 108 +832 896 108 +848 896 108 +864 896 108 +880 896 108 +896 896 108 +912 896 108 +928 896 108 +944 896 108 +960 896 108 +976 896 108 +992 896 108 +1008 896 108 +1024 896 108 +1040 896 108 +1056 896 108 +1072 896 108 +1088 896 108 +1104 896 108 +1120 896 108 +1136 896 108 +1152 896 108 +1168 896 108 +1184 896 108 +1200 896 108 +1216 896 108 +1232 896 108 +1248 896 108 +1264 896 108 +768 912 108 +784 912 108 +800 912 108 +816 912 108 +832 912 108 +848 912 108 +864 912 108 +880 912 108 +896 912 108 +912 912 108 +928 912 108 +944 912 108 +960 912 108 +976 912 108 +992 912 108 +1008 912 108 +1024 912 108 +1040 912 108 +1056 912 108 +1072 912 108 +1088 912 108 +1104 912 108 +1120 912 108 +1136 912 108 +1152 912 108 +1168 912 108 +1184 912 108 +1200 912 108 +1216 912 108 +1232 912 108 +1248 912 108 +1264 912 108 +768 928 108 +784 928 108 +800 928 108 +816 928 108 +832 928 108 +848 928 108 +864 928 108 +880 928 108 +896 928 108 +912 928 108 +928 928 108 +944 928 108 +960 928 108 +976 928 108 +992 928 108 +1008 928 108 +1024 928 108 +1040 928 108 +1056 928 108 +1072 928 108 +1088 928 108 +1104 928 108 +1120 928 108 +1136 928 108 +1152 928 108 +1168 928 108 +1184 928 108 +1200 928 108 +1216 928 108 +1232 928 108 +1248 928 108 +1264 928 108 +768 944 108 +784 944 108 +800 944 108 +816 944 108 +832 944 108 +848 944 108 +864 944 108 +880 944 108 +896 944 108 +912 944 108 +928 944 108 +944 944 108 +960 944 108 +976 944 108 +992 944 108 +1008 944 108 +1024 944 108 +1040 944 108 +1056 944 108 +1072 944 108 +1088 944 108 +1104 944 108 +1120 944 108 +1136 944 108 +1152 944 108 +1168 944 108 +1184 944 108 +1200 944 108 +1216 944 108 +1232 944 108 +1248 944 108 +1264 944 108 +768 960 108 +784 960 108 +800 960 108 +816 960 108 +832 960 108 +848 960 108 +864 960 108 +880 960 108 +896 960 108 +912 960 108 +928 960 108 +944 960 108 +960 960 108 +976 960 108 +992 960 108 +1008 960 108 +1024 960 108 +1040 960 108 +1056 960 108 +1072 960 108 +1088 960 108 +1104 960 108 +1120 960 108 +1136 960 108 +1152 960 108 +1168 960 108 +1184 960 108 +1200 960 108 +1216 960 108 +1232 960 108 +1248 960 108 +1264 960 108 +768 976 108 +784 976 108 +800 976 108 +816 976 108 +832 976 108 +848 976 108 +864 976 108 +880 976 108 +896 976 108 +912 976 108 +928 976 108 +944 976 108 +960 976 108 +976 976 108 +992 976 108 +1008 976 108 +1024 976 108 +1040 976 108 +1056 976 108 +1072 976 108 +1088 976 108 +1104 976 108 +1120 976 108 +1136 976 108 +1152 976 108 +1168 976 108 +1184 976 108 +1200 976 108 +1216 976 108 +1232 976 108 +1248 976 108 +1264 976 108 +768 992 108 +784 992 108 +800 992 108 +816 992 108 +832 992 108 +848 992 108 +864 992 108 +880 992 108 +896 992 108 +912 992 108 +928 992 108 +944 992 108 +960 992 108 +976 992 108 +992 992 108 +1008 992 108 +1024 992 108 +1040 992 108 +1056 992 108 +1072 992 108 +1088 992 108 +1104 992 108 +1120 992 108 +1136 992 108 +1152 992 108 +1168 992 108 +1184 992 108 +1200 992 108 +1216 992 108 +1232 992 108 +1248 992 108 +1264 992 108 +768 1008 108 +784 1008 108 +800 1008 108 +816 1008 108 +832 1008 108 +848 1008 108 +864 1008 108 +880 1008 108 +896 1008 108 +912 1008 108 +928 1008 108 +944 1008 108 +960 1008 108 +976 1008 108 +992 1008 108 +1008 1008 108 +1024 1008 108 +1040 1008 108 +1056 1008 108 +1072 1008 108 +1088 1008 108 +1104 1008 108 +1120 1008 108 +1136 1008 108 +1152 1008 108 +1168 1008 108 +1184 1008 108 +1200 1008 108 +1216 1008 108 +1232 1008 108 +1248 1008 108 +1264 1008 108 +768 1024 108 +784 1024 108 +800 1024 108 +816 1024 108 +832 1024 108 +848 1024 108 +864 1024 108 +880 1024 108 +896 1024 108 +912 1024 108 +928 1024 108 +944 1024 108 +960 1024 108 +976 1024 108 +992 1024 108 +1008 1024 108 +1024 1024 108 +1040 1024 108 +1056 1024 108 +1072 1024 108 +1088 1024 108 +1104 1024 108 +1120 1024 108 +1136 1024 108 +1152 1024 108 +1168 1024 108 +1184 1024 108 +1200 1024 108 +1216 1024 108 +1232 1024 108 +1248 1024 108 +1264 1024 108 +768 1040 108 +784 1040 108 +800 1040 108 +816 1040 108 +832 1040 108 +848 1040 108 +864 1040 108 +880 1040 108 +896 1040 108 +912 1040 108 +928 1040 108 +944 1040 108 +960 1040 108 +976 1040 108 +992 1040 108 +1008 1040 108 +1024 1040 108 +1040 1040 108 +1056 1040 108 +1072 1040 108 +1088 1040 108 +1104 1040 108 +1120 1040 108 +1136 1040 108 +1152 1040 108 +1168 1040 108 +1184 1040 108 +1200 1040 108 +1216 1040 108 +1232 1040 108 +1248 1040 108 +1264 1040 108 +768 1056 108 +784 1056 108 +800 1056 108 +816 1056 108 +832 1056 108 +848 1056 108 +864 1056 108 +880 1056 108 +896 1056 108 +912 1056 108 +928 1056 108 +944 1056 108 +960 1056 108 +976 1056 108 +992 1056 108 +1008 1056 108 +1024 1056 108 +1040 1056 108 +1056 1056 108 +1072 1056 108 +1088 1056 108 +1104 1056 108 +1120 1056 108 +1136 1056 108 +1152 1056 108 +1168 1056 108 +1184 1056 108 +1200 1056 108 +1216 1056 108 +1232 1056 108 +1248 1056 108 +1264 1056 108 +768 1072 108 +784 1072 108 +800 1072 108 +816 1072 108 +832 1072 108 +848 1072 108 +864 1072 108 +880 1072 108 +896 1072 108 +912 1072 108 +928 1072 108 +944 1072 108 +960 1072 108 +976 1072 108 +992 1072 108 +1008 1072 108 +1024 1072 108 +1040 1072 108 +1056 1072 108 +1072 1072 108 +1088 1072 108 +1104 1072 108 +1120 1072 108 +1136 1072 108 +1152 1072 108 +1168 1072 108 +1184 1072 108 +1200 1072 108 +1216 1072 108 +1232 1072 108 +1248 1072 108 +1264 1072 108 +768 1088 108 +784 1088 108 +800 1088 108 +816 1088 108 +832 1088 108 +848 1088 108 +864 1088 108 +880 1088 108 +896 1088 108 +912 1088 108 +928 1088 108 +944 1088 108 +960 1088 108 +976 1088 108 +992 1088 108 +1008 1088 108 +1024 1088 108 +1040 1088 108 +1056 1088 108 +1072 1088 108 +1088 1088 108 +1104 1088 108 +1120 1088 108 +1136 1088 108 +1152 1088 108 +1168 1088 108 +1184 1088 108 +1200 1088 108 +1216 1088 108 +1232 1088 108 +1248 1088 108 +1264 1088 108 +768 1104 108 +784 1104 108 +800 1104 108 +816 1104 108 +832 1104 108 +848 1104 108 +864 1104 108 +880 1104 108 +896 1104 108 +912 1104 108 +928 1104 108 +944 1104 108 +960 1104 108 +976 1104 108 +992 1104 108 +1008 1104 108 +1024 1104 108 +1040 1104 108 +1056 1104 108 +1072 1104 108 +1088 1104 108 +1104 1104 108 +1120 1104 108 +1136 1104 108 +1152 1104 108 +1168 1104 108 +1184 1104 108 +1200 1104 108 +1216 1104 108 +1232 1104 108 +1248 1104 108 +1264 1104 108 +768 1120 108 +784 1120 108 +800 1120 108 +816 1120 108 +832 1120 108 +848 1120 108 +864 1120 108 +880 1120 108 +896 1120 108 +912 1120 108 +928 1120 108 +944 1120 108 +960 1120 108 +976 1120 108 +992 1120 108 +1008 1120 108 +1024 1120 108 +1040 1120 108 +1056 1120 108 +1072 1120 108 +1088 1120 108 +1104 1120 108 +1120 1120 108 +1136 1120 108 +1152 1120 108 +1168 1120 108 +1184 1120 108 +1200 1120 108 +1216 1120 108 +1232 1120 108 +1248 1120 108 +1264 1120 108 +768 1136 108 +784 1136 108 +800 1136 108 +816 1136 108 +832 1136 108 +848 1136 108 +864 1136 108 +880 1136 108 +896 1136 108 +912 1136 108 +928 1136 108 +944 1136 108 +960 1136 108 +976 1136 108 +992 1136 108 +1008 1136 108 +1024 1136 108 +1040 1136 108 +1056 1136 108 +1072 1136 108 +1088 1136 108 +1104 1136 108 +1120 1136 108 +1136 1136 108 +1152 1136 108 +1168 1136 108 +1184 1136 108 +1200 1136 108 +1216 1136 108 +1232 1136 108 +1248 1136 108 +1264 1136 108 +768 1152 108 +784 1152 108 +800 1152 108 +816 1152 108 +832 1152 108 +848 1152 108 +864 1152 108 +880 1152 108 +896 1152 108 +912 1152 108 +928 1152 108 +944 1152 108 +960 1152 108 +976 1152 108 +992 1152 108 +1008 1152 108 +1024 1152 108 +1040 1152 108 +1056 1152 108 +1072 1152 108 +1088 1152 108 +1104 1152 108 +1120 1152 108 +1136 1152 108 +1152 1152 108 +1168 1152 108 +1184 1152 108 +1200 1152 108 +1216 1152 108 +1232 1152 108 +1248 1152 108 +1264 1152 108 +768 1168 108 +784 1168 108 +800 1168 108 +816 1168 108 +832 1168 108 +848 1168 108 +864 1168 108 +880 1168 108 +896 1168 108 +912 1168 108 +928 1168 108 +944 1168 108 +960 1168 108 +976 1168 108 +992 1168 108 +1008 1168 108 +1024 1168 108 +1040 1168 108 +1056 1168 108 +1072 1168 108 +1088 1168 108 +1104 1168 108 +1120 1168 108 +1136 1168 108 +1152 1168 108 +1168 1168 108 +1184 1168 108 +1200 1168 108 +1216 1168 108 +1232 1168 108 +1248 1168 108 +1264 1168 108 +768 1184 108 +784 1184 108 +800 1184 108 +816 1184 108 +832 1184 108 +848 1184 108 +864 1184 108 +880 1184 108 +896 1184 108 +912 1184 108 +928 1184 108 +944 1184 108 +960 1184 108 +976 1184 108 +992 1184 108 +1008 1184 108 +1024 1184 108 +1040 1184 108 +1056 1184 108 +1072 1184 108 +1088 1184 108 +1104 1184 108 +1120 1184 108 +1136 1184 108 +1152 1184 108 +1168 1184 108 +1184 1184 108 +1200 1184 108 +1216 1184 108 +1232 1184 108 +1248 1184 108 +1264 1184 108 +768 1200 108 +784 1200 108 +800 1200 108 +816 1200 108 +832 1200 108 +848 1200 108 +864 1200 108 +880 1200 108 +896 1200 108 +912 1200 108 +928 1200 108 +944 1200 108 +960 1200 108 +976 1200 108 +992 1200 108 +1008 1200 108 +1024 1200 108 +1040 1200 108 +1056 1200 108 +1072 1200 108 +1088 1200 108 +1104 1200 108 +1120 1200 108 +1136 1200 108 +1152 1200 108 +1168 1200 108 +1184 1200 108 +1200 1200 108 +1216 1200 108 +1232 1200 108 +1248 1200 108 +1264 1200 108 +768 1216 108 +784 1216 108 +800 1216 108 +816 1216 108 +832 1216 108 +848 1216 108 +864 1216 108 +880 1216 108 +896 1216 108 +912 1216 108 +928 1216 108 +944 1216 108 +960 1216 108 +976 1216 108 +992 1216 108 +1008 1216 108 +1024 1216 108 +1040 1216 108 +1056 1216 108 +1072 1216 108 +1088 1216 108 +1104 1216 108 +1120 1216 108 +1136 1216 108 +1152 1216 108 +1168 1216 108 +1184 1216 108 +1200 1216 108 +1216 1216 108 +1232 1216 108 +1248 1216 108 +1264 1216 108 +768 1232 108 +784 1232 108 +800 1232 108 +816 1232 108 +832 1232 108 +848 1232 108 +864 1232 108 +880 1232 108 +896 1232 108 +912 1232 108 +928 1232 108 +944 1232 108 +960 1232 108 +976 1232 108 +992 1232 108 +1008 1232 108 +1024 1232 108 +1040 1232 108 +1056 1232 108 +1072 1232 108 +1088 1232 108 +1104 1232 108 +1120 1232 108 +1136 1232 108 +1152 1232 108 +1168 1232 108 +1184 1232 108 +1200 1232 108 +1216 1232 108 +1232 1232 108 +1248 1232 108 +1264 1232 108 +768 1248 108 +784 1248 108 +800 1248 108 +816 1248 108 +832 1248 108 +848 1248 108 +864 1248 108 +880 1248 108 +896 1248 108 +912 1248 108 +928 1248 108 +944 1248 108 +960 1248 108 +976 1248 108 +992 1248 108 +1008 1248 108 +1024 1248 108 +1040 1248 108 +1056 1248 108 +1072 1248 108 +1088 1248 108 +1104 1248 108 +1120 1248 108 +1136 1248 108 +1152 1248 108 +1168 1248 108 +1184 1248 108 +1200 1248 108 +1216 1248 108 +1232 1248 108 +1248 1248 108 +1264 1248 108 +768 1264 108 +784 1264 108 +800 1264 108 +816 1264 108 +832 1264 108 +848 1264 108 +864 1264 108 +880 1264 108 +896 1264 108 +912 1264 108 +928 1264 108 +944 1264 108 +960 1264 108 +976 1264 108 +992 1264 108 +1008 1264 108 +1024 1264 108 +1040 1264 108 +1056 1264 108 +1072 1264 108 +1088 1264 108 +1104 1264 108 +1120 1264 108 +1136 1264 108 +1152 1264 108 +1168 1264 108 +1184 1264 108 +1200 1264 108 +1216 1264 108 +1232 1264 108 +1248 1264 108 +1264 1264 108 diff --git a/Exec/ABL/inputs_terrain b/Exec/ABL/inputs_terrain new file mode 100644 index 000000000..e91f8c3f1 --- /dev/null +++ b/Exec/ABL/inputs_terrain @@ -0,0 +1,63 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 1 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 1024 1024 1024 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_extent = 2048 2048 1024 +amr.n_cell = 64 64 32 + +geometry.is_periodic = 1 1 0 + +zlo.type = "SlipWall" +zhi.type = "SlipWall" + +# TIME STEP CONTROL +erf.cfl = 0.5 + +# DIAGNOSTICS & VERBOSITY +erf.sum_interval = 1 # timesteps between computing mass +erf.v = 1 # verbosity in ERF.cpp +amr.v = 1 # verbosity in Amr.cpp + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed + +# CHECKPOINT FILES +erf.check_file = chk # root name of checkpoint file +erf.check_int = 512 # number of timesteps between checkpoints + +# PLOTFILES +erf.plot_file_1 = plt # prefix of plotfile name +erf.plot_int_1 = 512 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhoadv_0 x_velocity y_velocity z_velocity pressure temp theta terrain_IB_mask + +# SOLVER CHOICE +erf.alpha_T = 0.0 +erf.alpha_C = 1.0 +erf.use_gravity = false + +erf.molec_diff_type = "None" +erf.les_type = "Smagorinsky" +erf.Cs = 0.1 + +erf.init_type = "uniform" + +erf.terrain_file = erf_terrain_def + +# PROBLEM PARAMETERS +prob.rho_0 = 1.0 +prob.A_0 = 1.0 + +prob.U_0 = 10.0 +prob.V_0 = 0.0 +prob.W_0 = 0.0 +prob.T_0 = 300.0 + +# Higher values of perturbations lead to instability +# Instability seems to be coming from BC +prob.U_0_Pert_Mag = 0.0 +prob.V_0_Pert_Mag = 0.0 # +prob.W_0_Pert_Mag = 0.0 diff --git a/Exec/ABL/write_terrain.py b/Exec/ABL/write_terrain.py new file mode 100644 index 000000000..e85e89ff7 --- /dev/null +++ b/Exec/ABL/write_terrain.py @@ -0,0 +1,15 @@ +import numpy as np + +target=open("erf_terrain_def","w") + +xc=1024 +yc=1024 +x=np.arange(xc-256,xc+256,16) +y=np.arange(yc-256,yc+256,16) + +X,Y=np.meshgrid(x,y) + +for i in range(0,X.shape[0]): + for j in range(0,X.shape[1]): + target.write("%g %g %g\n"%(X[i,j],Y[i,j],108)) +target.close() diff --git a/Exec/DevTests/EB_Test/ERF_Prob.H b/Exec/DevTests/EB_Test/ERF_Prob.H index 8d670810a..9dc736590 100644 --- a/Exec/DevTests/EB_Test/ERF_Prob.H +++ b/Exec/DevTests/EB_Test/ERF_Prob.H @@ -23,6 +23,9 @@ struct ProbParm : ProbParmDefaults { amrex::Real yc_frac = 0.5; // Location of "center" of scalar (multiplies domain length) amrex::Real zc_frac = 0.5; // Location of "center" of scalar (multiplies domain length) + amrex::Real xradius = 10.0; // x-radius of scalar bubble + amrex::Real zradius = 10.0; // z-radius of scalar bubble + int prob_type = -1; }; // namespace ProbParm diff --git a/Exec/DevTests/EB_Test/ERF_Prob.cpp b/Exec/DevTests/EB_Test/ERF_Prob.cpp index 1c4ec28cc..e15bd4222 100644 --- a/Exec/DevTests/EB_Test/ERF_Prob.cpp +++ b/Exec/DevTests/EB_Test/ERF_Prob.cpp @@ -32,6 +32,9 @@ Problem::Problem() pp.query("prob_type", parms.prob_type); + pp.query("xradius", parms.xradius); + pp.query("zradius", parms.zradius); + init_base_parms(parms.rho_0, parms.T_0); } @@ -81,16 +84,71 @@ Problem::init_custom_pert( AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1); - ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + // ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + // { + // // Set scalar = 0 everywhere + // state_pert(i, j, k, RhoScalar_comp) = 0.0; + + // if (use_moisture) { + // state_pert(i, j, k, RhoQ1_comp) = 0.0; + // state_pert(i, j, k, RhoQ2_comp) = 0.0; + // } + // }); + + // Set the state_pert + ParallelFor(bx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - // Set scalar = 0 everywhere - state_pert(i, j, k, RhoScalar_comp) = 0.0; + // 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_d.xc_frac * (prob_lo[0] + prob_hi[0]); + const Real yc = parms_d.yc_frac * (prob_lo[1] + prob_hi[1]); + const Real zc = parms_d.zc_frac * (prob_lo[2] + prob_hi[2]); + + // Define ellipse parameters + const Real r0 = parms_d.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)); + const Real r2d_xz_nd = std::sqrt((x-xc)*(x-xc)/parms_d.xradius/parms_d.xradius + + (z-zc)*(z-zc)/parms_d.zradius/parms_d.zradius); + + if (parms_d.prob_type == 10) + { + // Set scalar = A_0*exp(-10r^2), where r is distance from center of domain, + // + B_0*sin(x) + // state_pert(i, j, k, RhoScalar_comp) = parms_d.A_0 * exp(-10.*r3d*r3d) + parms_d.B_0*sin(x); + state_pert(i, j, k, RhoScalar_comp) = parms_d.A_0 * exp(-0.1*r2d_xz*r2d_xz) + parms_d.B_0*sin(x); + + } else if (parms_d.prob_type == 11) { + if (r2d_xz_nd < 1.0) + { + state_pert(i, j, k, RhoScalar_comp) = 0.5 * parms_d.A_0 * (1.0 + std::cos(PI*r2d_xz_nd)); + } else { + state_pert(i, j, k, RhoScalar_comp) = 0.0; + } + } else { + // Set scalar = A_0 in a ball of radius r0 and 0 elsewhere + if (r3d < r0) { + state_pert(i, j, k, RhoScalar_comp) = parms_d.A_0; + } else { + state_pert(i, j, k, RhoScalar_comp) = 0.0; + } + } + + state_pert(i, j, k, RhoScalar_comp) *= parms_d.rho_0; if (use_moisture) { state_pert(i, j, k, RhoQ1_comp) = 0.0; state_pert(i, j, k, RhoQ2_comp) = 0.0; } - }); + }); // Set the x-velocity ParallelFor(xbx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept diff --git a/Exec/DevTests/EB_Test/inputs b/Exec/DevTests/EB_Test/inputs index 6cf5daf48..75e3e23f7 100644 --- a/Exec/DevTests/EB_Test/inputs +++ b/Exec/DevTests/EB_Test/inputs @@ -1,6 +1,5 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -max_step = 100 -max_step = 10 +max_step = 800 amr.max_grid_size = 256 256 256 @@ -17,16 +16,16 @@ fabarray.mfiter_tile_size = 1024 1024 1024 # PROBLEM SIZE & GEOMETRY geometry.prob_lo = 0. 0. 0. -geometry.prob_hi = 16. 1. 16. +geometry.prob_hi = 16. 4. 8. -amr.n_cell = 64 4 64 +amr.n_cell = 64 4 32 geometry.is_periodic = 0 1 0 xlo.type = Inflow xhi.type = Outflow -xlo.velocity = 1. 0. 0. +xlo.velocity = 10. 0. 0. xlo.density = 1.0 xlo.theta = 1.0 xlo.scalar = 0. @@ -36,7 +35,7 @@ zhi.type = SlipWall # TIME STEP CONTROL erf.substepping_type = None -erf.fixed_dt = 1.e-5 +erf.fixed_dt = 1.e-3 # DIAGNOSTICS & VERBOSITY erf.sum_interval = 1 # timesteps between computing mass @@ -72,4 +71,10 @@ erf.init_type = "uniform" # PROBLEM PARAMETERS prob.rho_0 = 1.0 prob.T_0 = 1.0 -prob.u_0 = 1.0 +prob.u_0 = 10.0 +prob.A_0 = 1.0 +prob.prob_type = 11 +prob.xradius = 3.0 +prob.zradius = 1.5 +prob.xc_frac = 0.25 +prob.zc_frac = 0.15 diff --git a/Exec/Make.ERF b/Exec/Make.ERF index d90ccc172..d25135f68 100644 --- a/Exec/Make.ERF +++ b/Exec/Make.ERF @@ -149,6 +149,11 @@ include $(ERF_MOISTURE_KESSLER_DIR)/Make.package VPATH_LOCATIONS += $(ERF_MOISTURE_KESSLER_DIR) INCLUDE_LOCATIONS += $(ERF_MOISTURE_KESSLER_DIR) +ERF_MOISTURE_SATADJ_DIR = $(ERF_SOURCE_DIR)/Microphysics/SatAdj +include $(ERF_MOISTURE_SATADJ_DIR)/Make.package +VPATH_LOCATIONS += $(ERF_MOISTURE_SATADJ_DIR) +INCLUDE_LOCATIONS += $(ERF_MOISTURE_SATADJ_DIR) + # If using windfarm parametrization, then compile all models and choose # at runtime from the inputs ifeq ($(USE_WINDFARM), TRUE) diff --git a/Exec/MoistRegTests/Bubble/ERF_Prob.cpp b/Exec/MoistRegTests/Bubble/ERF_Prob.cpp index a28cdfdaa..f079ba5d9 100644 --- a/Exec/MoistRegTests/Bubble/ERF_Prob.cpp +++ b/Exec/MoistRegTests/Bubble/ERF_Prob.cpp @@ -303,7 +303,6 @@ Problem::init_custom_pert( if (use_moisture) { state_pert(i, j, k, RhoQ1_comp) = 0.0; state_pert(i, j, k, RhoQ2_comp) = 0.0; - state_pert(i, j, k, RhoQ3_comp) = 0.0; } }); } // do_moist_bubble @@ -389,7 +388,6 @@ Problem::init_custom_pert( // mean states state_pert(i, j, k, RhoQ1_comp) = rho*q_v_hot; state_pert(i, j, k, RhoQ2_comp) = rho*(parms_d.qt_init - q_v_hot); - state_pert(i, j, k, RhoQ3_comp) = 0.0; // Cold microphysics are present int nstate = state_pert.ncomp; @@ -429,7 +427,6 @@ Problem::init_custom_pert( if (use_moisture) { state_pert(i, j, k, RhoQ1_comp) = 0.0; state_pert(i, j, k, RhoQ2_comp) = 0.0; - state_pert(i, j, k, RhoQ3_comp) = 0.0; } }); } // do_moist_bubble diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index 3e7349b88..c7e9e6b32 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -29,6 +29,10 @@ AMREX_ENUM(SubsteppingType, None, Explicit, Implicit ); +AMREX_ENUM(MeshType, + ConstantDz, StretchedDz, VariableDz +); + AMREX_ENUM(TerrainType, None, Static, Moving ); @@ -38,7 +42,7 @@ AMREX_ENUM(MoistureModelType, ); AMREX_ENUM(MoistureType, - Kessler, SAM, SAM_NoIce, SAM_NoPrecip_NoIce, Kessler_NoRain, None + SAM, SAM_NoIce, SAM_NoPrecip_NoIce, Kessler, Kessler_NoRain, SatAdj, None ); AMREX_ENUM(WindFarmType, @@ -87,6 +91,9 @@ struct SolverChoice { "The grid stretching ratio must be greater than 1"); } if (grid_stretching_ratio >= 1) { + if (mesh_type == MeshType::ConstantDz) { + mesh_type = MeshType::StretchedDz; + } if (terrain_type != TerrainType::Static) { amrex::Print() << "Turning terrain on to enable grid stretching" << std::endl; terrain_type = TerrainType::Static; @@ -123,15 +130,19 @@ struct SolverChoice { } else if (moisture_type == MoistureType::Kessler_NoRain) { RhoQv_comp = RhoQ1_comp; RhoQc_comp = RhoQ2_comp; + } else if (moisture_type == MoistureType::SatAdj) { + RhoQv_comp = RhoQ1_comp; + RhoQc_comp = RhoQ2_comp; } // TODO: should we set default for dry?? // Set a different default for moist vs dry if (moisture_type != MoistureType::None) { - if (moisture_type == MoistureType::Kessler_NoRain || - moisture_type == MoistureType::SAM || - moisture_type == MoistureType::SAM_NoIce || - moisture_type == MoistureType::SAM_NoPrecip_NoIce) + if (moisture_type == MoistureType::Kessler_NoRain || + moisture_type == MoistureType::SAM || + moisture_type == MoistureType::SAM_NoIce || + moisture_type == MoistureType::SAM_NoPrecip_NoIce || + moisture_type == MoistureType::SatAdj) { buoyancy_type = 1; // uses Rhoprime } else { @@ -148,10 +159,20 @@ struct SolverChoice { // Is the terrain none, static or moving? pp.query_enum_case_insensitive("terrain_type",terrain_type); + + if (terrain_type != TerrainType::None) { + mesh_type = MeshType::VariableDz; + } + int n_zlevels = pp.countval("terrain_z_levels"); - if (n_zlevels > 0 and terrain_type == TerrainType::None) + if (n_zlevels > 0) { - terrain_type = TerrainType::Static; + if (terrain_type == TerrainType::None) { + terrain_type = TerrainType::Static; + } + if (mesh_type == MeshType::ConstantDz) { + mesh_type = MeshType::StretchedDz; + } } // Use lagged_delta_rt in the fast integrator? @@ -434,6 +455,17 @@ struct SolverChoice { amrex::Print() << " None" << std::endl; } + amrex::Print() << " Mesh Type: " << std::endl; + if (mesh_type == MeshType::ConstantDz) { + amrex::Print() << " ConstantDz" << std::endl; + } else if (mesh_type == MeshType::StretchedDz) { + amrex::Print() << " StretchedDz" << std::endl; + } else if (mesh_type == MeshType::VariableDz) { + amrex::Print() << " VariableDz" << std::endl; + } else { + amrex::Abort("No mesh_type set!"); + } + amrex::Print() << "ABL Driver Type: " << std::endl; if (abl_driver_type == ABLDriverType::None) { amrex::Print() << " None" << std::endl; @@ -533,6 +565,9 @@ struct SolverChoice { inline static TerrainType terrain_type = TerrainType::None; + inline static + MeshType mesh_type = MeshType::ConstantDz; + static void set_flat_terrain_flag () { @@ -657,7 +692,10 @@ struct SolverChoice { amrex::Real windfarm_x_shift = -1.0; amrex::Real windfarm_y_shift = -1.0; - // Flag for valid canopy model - bool do_forest {false}; + // Use forest canopy model? + bool do_forest_drag {false}; + + // Use immersed forcing representation of terrain? + bool do_terrain_drag {false}; }; #endif diff --git a/Source/EB/ERF_redistribute.cpp b/Source/EB/ERF_Redistribute.cpp similarity index 77% rename from Source/EB/ERF_redistribute.cpp rename to Source/EB/ERF_Redistribute.cpp index 184d73ad1..aa98941c1 100644 --- a/Source/EB/ERF_redistribute.cpp +++ b/Source/EB/ERF_Redistribute.cpp @@ -9,10 +9,11 @@ using namespace amrex; void ERF::redistribute_term ( int lev, - MultiFab& result, - MultiFab& result_tmp, // Saves doing a MF::copy. does this matter??? - MultiFab const& state, - BCRec const* bc) // this is bc for the state (needed for SRD slopes) + MultiFab& result, + MultiFab& result_tmp, // Saves doing a MF::copy. does this matter??? + MultiFab const& state, + BCRec const* bc, // this is bc for the state (needed for SRD slopes) + Real const dt) { // ************************************************************************ // Redistribute result_tmp and pass out result @@ -26,16 +27,17 @@ ERF::redistribute_term ( int lev, #endif for (MFIter mfi(state,TilingIfNotGPU()); mfi.isValid(); ++mfi) { - redistribute_term(mfi, result, result_tmp, state, bc, lev); + redistribute_term(mfi, lev, result, result_tmp, state, bc, dt); } } void ERF::redistribute_term ( MFIter const& mfi, int lev, - MultiFab& result, - MultiFab& result_tmp, - MultiFab const& state, - BCRec const* bc) // this is bc for the state (needed for SRD slopes) + MultiFab& result, + MultiFab& result_tmp, + MultiFab const& state, + BCRec const* bc, // this is bc for the state (needed for SRD slopes) + Real const dt) { AMREX_ASSERT(result.nComp() == state.nComp()); @@ -57,13 +59,13 @@ ERF::redistribute_term ( MFIter const& mfi, int lev, auto const& vfrac = ebfact.getVolFrac().const_array(mfi); auto const& ccc = ebfact.getCentroid().const_array(mfi); - auto const& apx = ebfact.getAreaFrac()[0]->const_array(mfi);, - auto const& apy = ebfact.getAreaFrac()[1]->const_array(mfi);, - auto const& apz = ebfact.getAreaFrac()[2]->const_array(mfi);); + auto const& apx = ebfact.getAreaFrac()[0]->const_array(mfi); + auto const& apy = ebfact.getAreaFrac()[1]->const_array(mfi); + auto const& apz = ebfact.getAreaFrac()[2]->const_array(mfi); - auto const& fcx = ebfact.getFaceCent()[0]->const_array(mfi);, - auto const& fcy = ebfact.getFaceCent()[1]->const_array(mfi);, - auto const& fcz = ebfact.getFaceCent()[2]->const_array(mfi);); + auto const& fcx = ebfact.getFaceCent()[0]->const_array(mfi); + auto const& fcy = ebfact.getFaceCent()[1]->const_array(mfi); + auto const& fcz = ebfact.getFaceCent()[2]->const_array(mfi); Box gbx = bx; gbx.grow(3); @@ -77,7 +79,7 @@ ERF::redistribute_term ( MFIter const& mfi, int lev, // scratch(i,j,k) = 1.; // }); - std::string redistribution_type = "StateRedistribution"; + std::string redistribution_type = "StateRedist"; // State redist acts on the state. Array4 state_arr = state.const_array(mfi); @@ -85,7 +87,7 @@ ERF::redistribute_term ( MFIter const& mfi, int lev, scratch, flag, apx, apy, apz, vfrac, fcx, fcy, fcz, ccc, - bc, geom[lev], m_dt, edistribution_type); + bc, geom[lev], dt, redistribution_type); } else { diff --git a/Source/EB/Make.package b/Source/EB/Make.package index a9e8662af..8bae36093 100644 --- a/Source/EB/Make.package +++ b/Source/EB/Make.package @@ -6,6 +6,5 @@ CEXE_sources += ERF_EBCylinder.cpp CEXE_sources += ERF_EBRegular.cpp CEXE_sources += ERF_Redistribute.cpp CEXE_sources += ERF_WriteEBSurface.cpp - CEXE_headers += ERF_EBIF.H CEXE_headers += ERF_TerrainIF.H diff --git a/Source/ERF.H b/Source/ERF.H index d49fffe32..27dd95b14 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -46,6 +46,7 @@ #include #include #include +#include #ifdef ERF_USE_PARTICLES #include "ERF_ParticleData.H" @@ -475,6 +476,18 @@ public: void make_eb_box (); void make_eb_cylinder (); void make_eb_regular (); + void redistribute_term ( int lev, + amrex::MultiFab& result, + amrex::MultiFab& result_tmp, + amrex::MultiFab const& state, + amrex::BCRec const* bc, + amrex::Real const dt); + void redistribute_term ( amrex::MFIter const& mfi, int lev, + amrex::MultiFab& result, + amrex::MultiFab& result_tmp, + amrex::MultiFab const& state, + amrex::BCRec const* bc, + amrex::Real const dt); #endif // more flexible version of AverageDown() that lets you average down across multiple levels @@ -814,8 +827,10 @@ private: amrex::Vector> SFS_q1fx1_lev, SFS_q1fx2_lev, SFS_q1fx3_lev; amrex::Vector> SFS_q2fx3_lev; - // Terrain / grid stretching + // Grid stretching amrex::Vector> zlevels_stag; // nominal height levels + + // Terrain data structures amrex::Vector> z_phys_nd; amrex::Vector> z_phys_cc; @@ -838,6 +853,7 @@ private: amrex::Vector> z_t_rk; + // Map scale factors amrex::Vector> mapfac_m; amrex::Vector> mapfac_u; amrex::Vector> mapfac_v; @@ -906,6 +922,7 @@ private: // Time step controls static amrex::Real cfl; + static amrex::Real sub_cfl; static amrex::Real init_shrink; static amrex::Real change_max; @@ -965,7 +982,9 @@ private: "Lturb", // moisture vars "qt", "qv", "qc", "qi", "qp", "qrain", "qsnow", "qgraup", "qsat", - "rain_accum", "snow_accum", "graup_accum" + "rain_accum", "snow_accum", "graup_accum", + // Terrain IB mask + "terrain_IB_mask" #ifdef ERF_USE_EB // EB variables ,"volfrac", @@ -1139,7 +1158,8 @@ private: std::unique_ptr m_w2d = nullptr; std::unique_ptr m_r2d = nullptr; std::unique_ptr m_most = nullptr; - amrex::Vector> m_forest; + amrex::Vector> m_forest_drag; + amrex::Vector> m_terrain_drag; // // Holds info for dynamically generated tagging criteria diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 3c759adcf..3759ede46 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -25,6 +26,7 @@ SolverChoice ERF::solverChoice; // Time step control Real ERF::cfl = 0.8; +Real ERF::sub_cfl = 1.0; Real ERF::init_shrink = 1.0; Real ERF::change_max = 1.1; int ERF::fixed_mri_dt_ratio = 0; @@ -84,6 +86,15 @@ Vector BCNames = {"xlo", "ylo", "zlo", "xhi", "yhi", "zhi"}; // - initializes BCRec boundary condition object ERF::ERF () { + int fix_random_seed = 0; + ParmParse pp("erf"); pp.query("fix_random_seed", fix_random_seed); + // Note that the value of 1024UL is not significant -- the point here is just to set the + // same seed for all MPI processes for the purpose of regression testing + if (fix_random_seed) { + Print() << "Fixing the random seed" << std::endl; + InitRandom(1024UL); + } + ERF_shared(); } @@ -130,8 +141,13 @@ ERF::ERF_shared () lsm_flux.resize(nlevs_max); // NOTE: size canopy model before readparams (if file exists, we construct) - m_forest.resize(nlevs_max); - for (int lev = 0; lev < max_level; ++lev) { m_forest[lev] = nullptr; } + m_forest_drag.resize(nlevs_max); + for (int lev = 0; lev < max_level; ++lev) { m_forest_drag[lev] = nullptr;} + + // Immersed Forcing Representation of Terrain + m_terrain_drag.resize(nlevs_max); + for (int lev = 0; lev < max_level; ++lev) { m_terrain_drag[lev] = nullptr;} + ReadParameters(); initializeMicrophysics(nlevs_max); @@ -143,7 +159,7 @@ ERF::ERF_shared () const std::string& pv1 = "plot_vars_1"; setPlotVariables(pv1,plot_var_names_1); const std::string& pv2 = "plot_vars_2"; setPlotVariables(pv2,plot_var_names_2); - // This is only used when we have flat terrain and stretched z levels. + // This is only used when we have mesh_type == MeshType::StretchedDz stretched_dz_h.resize(nlevs_max); stretched_dz_d.resize(nlevs_max); @@ -159,7 +175,8 @@ ERF::ERF_shared () solverChoice.zsurf, solverChoice.dz0); - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type == MeshType::StretchedDz || + SolverChoice::mesh_type == MeshType::VariableDz) { int nz = geom[0].Domain().length(2) + 1; // staggered if (std::fabs(zlevels_stag[0][nz-1]-geom[0].ProbHi(2)) > 1.0e-4) { Print() << "Note: prob_hi[2]=" << geom[0].ProbHi(2) @@ -448,7 +465,6 @@ ERF::post_timestep (int nstep, Real time, Real dt_lev0) if (solverChoice.coupling_type == CouplingType::TwoWay) { - bool use_terrain = (SolverChoice::terrain_type != TerrainType::None); int ncomp = vars_new[0][Vars::cons].nComp(); for (int lev = finest_level-1; lev >= 0; lev--) { @@ -459,16 +475,16 @@ ERF::post_timestep (int nstep, Real time, Real dt_lev0) const Box& bx = mfi.tilebox(); const Array4< Real> cons_arr = vars_new[lev][Vars::cons].array(mfi); const Array4 mapfac_arr = mapfac_m[lev]->const_array(mfi); - if (use_terrain) { - const Array4 detJ_arr = detJ_cc[lev]->const_array(mfi); + if (solverChoice.mesh_type == MeshType::ConstantDz) { ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - cons_arr(i,j,k,n) *= detJ_arr(i,j,k) / (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)); + cons_arr(i,j,k,n) /= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)); }); } else { + const Array4 detJ_arr = detJ_cc[lev]->const_array(mfi); ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - cons_arr(i,j,k,n) /= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)); + cons_arr(i,j,k,n) *= detJ_arr(i,j,k) / (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)); }); } } // mfi @@ -482,16 +498,16 @@ ERF::post_timestep (int nstep, Real time, Real dt_lev0) const Box& bx = mfi.tilebox(); const Array4< Real> cons_arr = vars_new[lev][Vars::cons].array(mfi); const Array4 mapfac_arr = mapfac_m[lev]->const_array(mfi); - if (use_terrain) { - const Array4 detJ_arr = detJ_cc[lev]->const_array(mfi); + if (solverChoice.mesh_type == MeshType::ConstantDz) { ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - cons_arr(i,j,k,n) *= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)) / detJ_arr(i,j,k); + cons_arr(i,j,k,n) *= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)); }); } else { + const Array4 detJ_arr = detJ_cc[lev]->const_array(mfi); ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept { - cons_arr(i,j,k,n) *= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)); + cons_arr(i,j,k,n) *= (mapfac_arr(i,j,0)*mapfac_arr(i,j,0)) / detJ_arr(i,j,k); }); } } // mfi @@ -633,16 +649,16 @@ void ERF::InitData_post () { if (restart_chkfile.empty()) { - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { if (init_type == InitType::Ideal) { - Abort("We do not currently support init_type = ideal with terrain"); + Abort("We do not currently support init_type = ideal with non-constant dz"); } } // // Make sure that detJ and z_phys_cc are the average of the data on a finer level if there is one // - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { for (int crse_lev = finest_level-1; crse_lev >= 0; crse_lev--) { average_down( *detJ_cc[crse_lev+1], *detJ_cc[crse_lev], 0, 1, refRatio(crse_lev)); average_down(*z_phys_cc[crse_lev+1], *z_phys_cc[crse_lev], 0, 1, refRatio(crse_lev)); @@ -704,8 +720,8 @@ ERF::InitData_post () { AMREX_ALWAYS_ASSERT_WITH_MESSAGE(finest_level == 0, "Thin immersed body with refinement not currently supported."); - if (SolverChoice::terrain_type != TerrainType::None) { - amrex::Print() << "NOTE: Thin immersed body with terrain has not been tested." << std::endl; + if (SolverChoice::mesh_type != MeshType::ConstantDz) { + amrex::Print() << "NOTE: Thin immersed body with non-constant dz has not been tested." << std::endl; } } @@ -782,8 +798,8 @@ ERF::InitData_post () h_v_geos[lev], d_v_geos[lev], geom[lev], z_phys_cc[lev]); } else { - if (SolverChoice::terrain_type != TerrainType::None && !SolverChoice::terrain_is_flat) { - amrex::Print() << "Note: 1-D geostrophic wind profile input is only defined for grid stretching, not real terrain" << std::endl; + if (SolverChoice::mesh_type == MeshType::VariableDz) { + amrex::Print() << "Note: 1-D geostrophic wind profile input is not defined for real terrain" << std::endl; } init_geo_wind_profile(solverChoice.abl_geo_wind_table, h_u_geos[lev], d_u_geos[lev], @@ -875,7 +891,6 @@ ERF::InitData_post () // if (restart_chkfile == "") { - // Note -- this projection is only defined for no terrain if (solverChoice.project_initial_velocity) { Real dummy_dt = 1.0; for (int lev = 0; lev <= finest_level; ++lev) @@ -917,7 +932,7 @@ ERF::InitData_post () for (int lev = 0; lev <= finest_level; ++lev) { dz_min[lev] = geom[lev].CellSize(2); - if ( SolverChoice::terrain_type != TerrainType::None ) { + if ( SolverChoice::mesh_type != MeshType::ConstantDz ) { dz_min[lev] *= (*detJ_cc[lev]).min(0); } } @@ -1441,6 +1456,7 @@ ERF::ReadParameters () // Time step controls pp.query("cfl", cfl); + pp.query("substepping_cfl", sub_cfl); pp.query("init_shrink", init_shrink); pp.query("change_max", change_max); @@ -1614,10 +1630,19 @@ ERF::ReadParameters () // Query the canopy model file name std::string forestfile; - solverChoice.do_forest = pp.query("forest_file", forestfile); - if (solverChoice.do_forest) { + solverChoice.do_forest_drag = pp.query("forest_file", forestfile); + if (solverChoice.do_forest_drag) { + for (int lev = 0; lev <= max_level; ++lev) { + m_forest_drag[lev] = std::make_unique(forestfile); + } + } + + //Query the terrain file name + std::string terrainfile; + solverChoice.do_terrain_drag = pp.query("terrain_file", terrainfile); + if (solverChoice.do_terrain_drag) { for (int lev = 0; lev <= max_level; ++lev) { - m_forest[lev] = std::make_unique(forestfile); + m_terrain_drag[lev] = std::make_unique(terrainfile); } } diff --git a/Source/ERF_MakeNewArrays.cpp b/Source/ERF_MakeNewArrays.cpp index c0874e153..40a67bbb5 100644 --- a/Source/ERF_MakeNewArrays.cpp +++ b/Source/ERF_MakeNewArrays.cpp @@ -42,7 +42,8 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, // ******************************************************************************************** // Allocate terrain arrays // ******************************************************************************************** - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type == MeshType::StretchedDz || + SolverChoice::mesh_type == MeshType::VariableDz) { z_phys_cc[lev] = std::make_unique(ba,dm,1,1); if (solverChoice.terrain_type == TerrainType::Moving) @@ -452,7 +453,8 @@ ERF::update_diffusive_arrays (int lev, const BoxArray& ba, const DistributionMap void ERF::init_zphys (int lev, Real time) { - if (SolverChoice::terrain_type != TerrainType::None) + if (SolverChoice::mesh_type == MeshType::StretchedDz || + SolverChoice::mesh_type == MeshType::VariableDz) { if (init_type != InitType::Real && init_type != InitType::Metgrid) { @@ -490,7 +492,7 @@ ERF::init_zphys (int lev, Real time) void ERF::remake_zphys (int lev, std::unique_ptr& temp_zphys_nd) { - if (lev > 0 && SolverChoice::terrain_type != TerrainType::None) + if (lev > 0 && SolverChoice::mesh_type == MeshType::VariableDz) { // // First interpolate from coarser level @@ -517,7 +519,8 @@ ERF::remake_zphys (int lev, std::unique_ptr& temp_zphys_nd) void ERF::update_terrain_arrays (int lev) { - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type == MeshType::StretchedDz || + SolverChoice::mesh_type == MeshType::VariableDz) { make_J(geom[lev],*z_phys_nd[lev],*detJ_cc[lev]); make_areas(geom[lev],*z_phys_nd[lev],*ax[lev],*ay[lev],*az[lev]); make_zcc(geom[lev],*z_phys_nd[lev],*z_phys_cc[lev]); @@ -549,7 +552,7 @@ ERF::initialize_integrator (int lev, MultiFab& cons_mf, MultiFab& vel_mf) void ERF::make_physbcs (int lev) { - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type == MeshType::VariableDz) { AMREX_ALWAYS_ASSERT(z_phys_nd[lev] != nullptr); } diff --git a/Source/ERF_MakeNewLevel.cpp b/Source/ERF_MakeNewLevel.cpp index c524b2434..51ceac884 100644 --- a/Source/ERF_MakeNewLevel.cpp +++ b/Source/ERF_MakeNewLevel.cpp @@ -141,8 +141,9 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba_in, // ******************************************************************************************** // Build the data structures for canopy model (depends upon z_phys) // ******************************************************************************************** - if (solverChoice.do_forest) { m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } + if (solverChoice.do_forest_drag) { m_forest_drag[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } + if (solverChoice.do_terrain_drag) { m_terrain_drag[lev]->define_terrain_blank_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } //******************************************************************************************** // Microphysics // ******************************************************************************************* @@ -210,7 +211,7 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, // // Make sure that detJ and z_phys_cc are the average of the data on a finer level if there is one // - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { for (int crse_lev = lev-1; crse_lev >= 0; crse_lev--) { average_down( *detJ_cc[crse_lev+1], *detJ_cc[crse_lev], 0, 1, refRatio(crse_lev)); average_down(*z_phys_cc[crse_lev+1], *z_phys_cc[crse_lev], 0, 1, refRatio(crse_lev)); @@ -220,8 +221,9 @@ ERF::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba, // ******************************************************************************************** // Build the data structures for canopy model (depends upon z_phys) // ******************************************************************************************** - if (solverChoice.do_forest) { m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } + if (solverChoice.do_forest_drag) { m_forest_drag[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } + if (solverChoice.do_terrain_drag) { m_terrain_drag[lev]->define_terrain_blank_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } //******************************************************************************************** // Microphysics // ******************************************************************************************* @@ -339,7 +341,7 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp // // Make sure that detJ and z_phys_cc are the average of the data on a finer level if there is one // - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { for (int crse_lev = lev-1; crse_lev >= 0; crse_lev--) { average_down( *detJ_cc[crse_lev+1], *detJ_cc[crse_lev], 0, 1, refRatio(crse_lev)); average_down(*z_phys_cc[crse_lev+1], *z_phys_cc[crse_lev], 0, 1, refRatio(crse_lev)); @@ -349,8 +351,9 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp // ******************************************************************************************** // Build the data structures for canopy model (depends upon z_phys) // ******************************************************************************************** - if (solverChoice.do_forest) { m_forest[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } + if (solverChoice.do_forest_drag) { m_forest_drag[lev]->define_drag_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } + if (solverChoice.do_terrain_drag) { m_terrain_drag[lev]->define_terrain_blank_field(ba, dm, geom[lev], z_phys_nd[lev].get()); } // ***************************************************************************************************** // Create the physbcs objects (after initializing the terrain but before calling FillCoarsePatch // ***************************************************************************************************** diff --git a/Source/ERF_Tagging.cpp b/Source/ERF_Tagging.cpp index e6d839911..db63f88eb 100644 --- a/Source/ERF_Tagging.cpp +++ b/Source/ERF_Tagging.cpp @@ -161,7 +161,7 @@ ERF::refinement_criteria_setup () jlo = static_cast((rbox_lo[1] - plo[1])/dx[1]); ihi = static_cast((rbox_hi[0] - plo[0])/dx[0]-1); jhi = static_cast((rbox_hi[1] - plo[1])/dx[1]-1); - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { // Search for k indices corresponding to nominal grid // AGL heights const Box& domain = geom[lev_for_box].Domain(); diff --git a/Source/IO/ERF_Checkpoint.cpp b/Source/IO/ERF_Checkpoint.cpp index 72c3a9da1..d2294277f 100644 --- a/Source/IO/ERF_Checkpoint.cpp +++ b/Source/IO/ERF_Checkpoint.cpp @@ -145,7 +145,7 @@ ERF::WriteCheckpointFile () const MultiFab::Copy(base,base_state[lev],0,0,ncomp_base,ng_base); VisMF::Write(base, MultiFabFileFullPrefix(lev, checkpointname, "Level_", "BaseState")); - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { // Note that we also write the ghost cells of z_phys_nd IntVect ng = z_phys_nd[lev]->nGrowVect(); MultiFab z_height(convert(grids[lev],IntVect(1,1,1)),dmap[lev],1,ng); @@ -443,7 +443,7 @@ ERF::ReadCheckpointFile () } base_state[lev].FillBoundary(geom[lev].periodicity()); - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { // Note that we also read the ghost cells of z_phys_nd IntVect ng = z_phys_nd[lev]->nGrowVect(); MultiFab z_height(convert(grids[lev],IntVect(1,1,1)),dmap[lev],1,ng); diff --git a/Source/IO/ERF_Plotfile.cpp b/Source/IO/ERF_Plotfile.cpp index 4b9ab1084..26ed76e16 100644 --- a/Source/IO/ERF_Plotfile.cpp +++ b/Source/IO/ERF_Plotfile.cpp @@ -72,7 +72,8 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector // for (int i = 0; i < derived_names.size(); ++i) { if ( containerHasElement(plot_var_names, derived_names[i]) ) { - if (SolverChoice::terrain_type != TerrainType::None || (derived_names[i] != "z_phys" && derived_names[i] != "detJ") ) { + if ( SolverChoice::mesh_type != MeshType::ConstantDz || + (derived_names[i] != "z_phys" && derived_names[i] != "detJ") ) { if ( (solverChoice.moisture_type == MoistureType::SAM || solverChoice.moisture_type == MoistureType::SAM_NoIce) || (derived_names[i] != "qi" && @@ -225,7 +226,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p // Vector of MultiFabs for nodal data Vector mf_nd(finest_level+1); - if (SolverChoice::terrain_type != TerrainType::None) { + if ( SolverChoice::mesh_type != MeshType::ConstantDz) { 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); @@ -488,6 +489,24 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p mf_comp ++; } + if (containerHasElement(plot_var_names, "terrain_IB_mask")) + { + MultiFab* terrain_blank = m_terrain_drag[lev]->get_terrain_blank_field(); +#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& src = terrain_blank->const_array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + derdat(i, j, k, mf_comp) = src(i,j,k); + }); + } + mf_comp ++; + } + #ifdef ERF_USE_WINDFARM if (containerHasElement(plot_var_names, "num_turb") and (solverChoice.windfarm_type == WindFarmType::Fitch or solverChoice.windfarm_type == WindFarmType::EWP or @@ -578,7 +597,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p const Array4& derdat = mf[lev].array(mfi); const Array4 & p_arr = pres.array(mfi); - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { const Array4& z_nd = z_phys_nd[lev]->const_array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { @@ -672,7 +691,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p const Array4& derdat = mf[lev].array(mfi); const Array4 & p_arr = pres.array(mfi); - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { const Array4& z_nd = z_phys_nd[lev]->const_array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { @@ -843,7 +862,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p mf_comp += 1; } // pres_hse_y - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { if (containerHasElement(plot_var_names, "z_phys")) { MultiFab::Copy(mf[lev],*z_phys_cc[lev],0,mf_comp,1,0); @@ -1053,7 +1072,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p // Precipitating components //-------------------------------------------------------------------------- - if(containerHasElement(plot_var_names, "qp")) + if(containerHasElement(plot_var_names, "qp") && (n_qstate >= 3)) { int n_start = (n_qstate > 3) ? RhoQ4_comp : RhoQ3_comp; int n_end = ncomp_cons - 1; @@ -1350,7 +1369,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p #endif // Fill terrain distortion MF - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { for (int lev(0); lev <= finest_level; ++lev) { MultiFab::Copy(mf_nd[lev],*z_phys_nd[lev],0,2,1,0); Real dz = Geom()[lev].CellSizeArray()[2]; @@ -1399,7 +1418,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p if (plotfile_type == PlotFileType::Amrex) { Print() << "Writing native plotfile " << plotfilename << "\n"; - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { WriteMultiLevelPlotfileWithTerrain(plotfilename, finest_level+1, GetVecOfConstPtrs(mf), GetVecOfConstPtrs(mf_nd), @@ -1510,7 +1529,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p } Print() << "Writing plotfile " << plotfilename << "\n"; - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { WriteMultiLevelPlotfileWithTerrain(plotfilename, finest_level+1, GetVecOfConstPtrs(mf2), GetVecOfConstPtrs(mf_nd), @@ -1523,7 +1542,7 @@ ERF::WritePlotFile (int which, PlotFileType plotfile_type, Vector p } } else { - if (SolverChoice::terrain_type != TerrainType::None) { + if (SolverChoice::mesh_type != MeshType::ConstantDz) { WriteMultiLevelPlotfileWithTerrain(plotfilename, finest_level+1, GetVecOfConstPtrs(mf), GetVecOfConstPtrs(mf_nd), diff --git a/Source/IO/ERF_WriteScalarProfiles.cpp b/Source/IO/ERF_WriteScalarProfiles.cpp index d839640ab..a8ef59498 100644 --- a/Source/IO/ERF_WriteScalarProfiles.cpp +++ b/Source/IO/ERF_WriteScalarProfiles.cpp @@ -406,7 +406,7 @@ ERF::volWgtSumMF (int lev, auto const& dx = geom[lev].CellSizeArray(); Real cell_vol = dx[0]*dx[1]*dx[2]; volume.setVal(cell_vol); - if (solverChoice.terrain_type != TerrainType::None) { + if (solverChoice.mesh_type != MeshType::ConstantDz) { MultiFab::Multiply(volume, *detJ_cc[lev], 0, 0, 1, 0); } sum = MultiFab::Dot(tmp, 0, volume, 0, 1, 0, local); diff --git a/Source/Initialization/ERF_Init1D.cpp b/Source/Initialization/ERF_Init1D.cpp index 6af99e5c7..c1ad1e081 100644 --- a/Source/Initialization/ERF_Init1D.cpp +++ b/Source/Initialization/ERF_Init1D.cpp @@ -95,7 +95,7 @@ ERF::initHSE (int lev) std::unique_ptr new_z_phys_cc; std::unique_ptr new_z_phys_nd; - if (solverChoice.terrain_type != TerrainType::None) { + if (solverChoice.mesh_type != MeshType::ConstantDz) { new_z_phys_cc = std::make_unique(ba_new,dm_new,1,1); new_z_phys_cc->ParallelCopy(*z_phys_cc[lev],0,0,1,1,1); @@ -151,7 +151,7 @@ ERF::erf_enforce_hse (int lev, std::unique_ptr& z_cc) { Real l_gravity = solverChoice.gravity; - bool l_use_terrain = (solverChoice.terrain_type != TerrainType::None); + bool l_use_terrain = (solverChoice.mesh_type != MeshType::ConstantDz); const auto geomdata = geom[lev].data(); const Real dz = geomdata.CellSize(2); diff --git a/Source/Initialization/ERF_InitCustom.cpp b/Source/Initialization/ERF_InitCustom.cpp index f7b7f1eb7..4d7d4ed44 100644 --- a/Source/Initialization/ERF_InitCustom.cpp +++ b/Source/Initialization/ERF_InitCustom.cpp @@ -42,15 +42,6 @@ ERF::init_custom (int lev) yvel_pert.setVal(0.); zvel_pert.setVal(0.); - int fix_random_seed = 0; - ParmParse pp("erf"); pp.query("fix_random_seed", fix_random_seed); - // Note that the value of 1024UL is not significant -- the point here is just to set the - // same seed for all MPI processes for the purpose of regression testing - if (fix_random_seed) { - Print() << "Fixing the random seed" << std::endl; - InitRandom(1024UL); - } - #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif diff --git a/Source/LinearSolvers/ERF_ComputeDivergence.cpp b/Source/LinearSolvers/ERF_ComputeDivergence.cpp index 292b0a62d..fbbeb133e 100644 --- a/Source/LinearSolvers/ERF_ComputeDivergence.cpp +++ b/Source/LinearSolvers/ERF_ComputeDivergence.cpp @@ -11,8 +11,6 @@ void ERF::compute_divergence (int lev, MultiFab& rhs, Array& rho0w_arr = rho0_u_const[2]->const_array(mfi); const Array4& rhs_arr = rhs.array(mfi); - if (SolverChoice::terrain_is_flat) { // flat terrain + if (SolverChoice::mesh_type == MeshType::StretchedDz) { Real* stretched_dz_d_ptr = stretched_dz_d[lev].data(); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real dz = stretched_dz_d_ptr[k]; @@ -41,7 +43,7 @@ void ERF::compute_divergence (int lev, MultiFab& rhs, Array& mom_mf, Mult auto const dom_hi = ubound(geom[lev].Domain()); #endif - bool l_use_terrain = SolverChoice::terrain_type != TerrainType::None; - // Make sure the solver only sees the levels over which we are solving Vector ba_tmp; ba_tmp.push_back(mom_mf[Vars::cons].boxArray()); Vector dm_tmp; dm_tmp.push_back(mom_mf[Vars::cons].DistributionMap()); @@ -44,7 +42,8 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& mom_mf, Mult // **************************************************************************** // #ifndef ERF_USE_EB - if (l_use_terrain && !SolverChoice::terrain_is_flat) { + if (solverChoice.mesh_type == MeshType::VariableDz) + { for ( MFIter mfi(rhs[0],TilingIfNotGPU()); mfi.isValid(); ++mfi) { const Array4& rho0u_arr = mom_mf[IntVars::xmom].const_array(mfi); @@ -133,7 +132,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& mom_mf, Mult // **************************************************************************** // No terrain or grid stretching // **************************************************************************** - if (!l_use_terrain) { + if (solverChoice.mesh_type == MeshType::ConstantDz) { #ifdef ERF_USE_FFT if (use_fft) { if (boxes_make_rectangle) { @@ -156,7 +155,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& mom_mf, Mult // **************************************************************************** // Grid stretching (flat terrain) // **************************************************************************** - else if (l_use_terrain && SolverChoice::terrain_is_flat) { + else if (solverChoice.mesh_type == MeshType::StretchedDz) { #ifndef ERF_USE_FFT amrex::Abort("Rebuild with USE_FFT = TRUE so you can use the FFT solver"); #else @@ -174,7 +173,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& mom_mf, Mult // **************************************************************************** // General terrain // **************************************************************************** - else if (l_use_terrain && !SolverChoice::terrain_is_flat) { + else if (solverChoice.mesh_type == MeshType::VariableDz) { #ifdef ERF_USE_FFT if (use_fft) { @@ -243,7 +242,7 @@ void ERF::project_velocities (int lev, Real l_dt, Vector& mom_mf, Mult // Now convert the rho0w MultiFab back to holding (rho0w) rather than Omega // **************************************************************************** // - if (l_use_terrain && !solverChoice.terrain_is_flat) + if (solverChoice.mesh_type == MeshType::VariableDz) { for (MFIter mfi(mom_mf[Vars::cons],TilingIfNotGPU()); mfi.isValid(); ++mfi) { diff --git a/Source/LinearSolvers/ERF_PoissonSolve_tb.cpp b/Source/LinearSolvers/ERF_PoissonSolve_tb.cpp index 100fdf75b..107f6c0b7 100644 --- a/Source/LinearSolvers/ERF_PoissonSolve_tb.cpp +++ b/Source/LinearSolvers/ERF_PoissonSolve_tb.cpp @@ -20,7 +20,7 @@ bool ERF::projection_has_dirichlet (Array bcs) const void ERF::project_velocities_tb (int lev, Real l_dt, Vector& vmf, MultiFab& pmf) { BL_PROFILE("ERF::project_velocities_tb()"); - AMREX_ALWAYS_ASSERT(solverChoice.terrain_type == TerrainType::None); + AMREX_ALWAYS_ASSERT(solverChoice.mesh_type == MeshType::ConstantDz); // Make sure the solver only sees the levels over which we are solving Vector ba_tmp; ba_tmp.push_back(vmf[Vars::cons].boxArray()); diff --git a/Source/LinearSolvers/ERF_SolveWithFFT.cpp b/Source/LinearSolvers/ERF_SolveWithFFT.cpp index e499a7f07..22228ea40 100644 --- a/Source/LinearSolvers/ERF_SolveWithFFT.cpp +++ b/Source/LinearSolvers/ERF_SolveWithFFT.cpp @@ -12,8 +12,6 @@ void ERF::solve_with_fft (int lev, MultiFab& rhs, MultiFab& phi, Array 0) { amrex::Print() << "Using the 3D FFT solver..." << std::endl; @@ -49,7 +46,8 @@ void ERF::solve_with_fft (int lev, MultiFab& rhs, MultiFab& phi, Array 0) { amrex::Print() << "Using the hybrid FFT solver..." << std::endl; @@ -89,7 +87,7 @@ void ERF::solve_with_fft (int lev, MultiFab& rhs, MultiFab& phi, Array const& fz_arr = fluxes[2].array(mfi); - if (l_use_terrain && SolverChoice::terrain_is_flat) { + if (solverChoice.mesh_type == MeshType::StretchedDz) { Real* stretched_dz_d_ptr = stretched_dz_d[lev].data(); ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { diff --git a/Source/LinearSolvers/ERF_TerrainPoisson.cpp b/Source/LinearSolvers/ERF_TerrainPoisson.cpp index 7eea971f7..f7bf18849 100644 --- a/Source/LinearSolvers/ERF_TerrainPoisson.cpp +++ b/Source/LinearSolvers/ERF_TerrainPoisson.cpp @@ -60,48 +60,70 @@ void TerrainPoisson::apply_bcs (MultiFab& phi) { Box bx = mfi.tilebox(); const Array4& phi_arr = phi.array(mfi); - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - if (i == domlo.x) { - if (bc_fft[0].first == FFT::Boundary::even) { + if (bx.smallEnd(0) <= domlo.x) { + if (bc_fft[0].first == FFT::Boundary::even) { + ParallelFor(makeSlab(bx,0,domlo.x), [=] AMREX_GPU_DEVICE (int i, int j, int k) + { phi_arr(i-1,j,k) = phi_arr(i,j,k); - } else if (bc_fft[0].first == FFT::Boundary::odd) { - phi_arr(i-1,j,k) = -phi_arr(i,j,k); - } - } else if (i == domhi.x) { - if (bc_fft[0].second == FFT::Boundary::even) { + }); + } else if (bc_fft[0].first == FFT::Boundary::odd) { + ParallelFor(makeSlab(bx,0,domlo.x), [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + phi_arr(i-1,j,k) = -phi_arr(i,j,k); + }); + } + } // lo x + if (bx.bigEnd(0) >= domhi.x) { + if (bc_fft[0].second == FFT::Boundary::even) { + ParallelFor(makeSlab(bx,0,domhi.x), [=] AMREX_GPU_DEVICE (int i, int j, int k) + { phi_arr(i+1,j,k) = phi_arr(i,j,k); - } else if (bc_fft[0].second == FFT::Boundary::odd) { - phi_arr(i+1,j,k) = -phi_arr(i,j,k); - } + }); + } else if (bc_fft[0].second == FFT::Boundary::odd) { + ParallelFor(makeSlab(bx,0,domhi.x), [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + phi_arr(i+1,j,k) = -phi_arr(i,j,k); + }); } - }); + } // lo x } // mfi - } + } // not periodic in x + if (!m_geom.isPeriodic(1)) { for (MFIter mfi(phi,true); mfi.isValid(); ++mfi) { Box bx = mfi.tilebox(); Box bx2(bx); bx2.grow(0,1); const Array4& phi_arr = phi.array(mfi); - ParallelFor(bx2, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - if (j == domlo.y) { - if (bc_fft[1].first == FFT::Boundary::even) { + if (bx.smallEnd(1) <= domlo.y) { + if (bc_fft[1].first == FFT::Boundary::even) { + ParallelFor(makeSlab(bx2,1,domlo.y), [=] AMREX_GPU_DEVICE (int i, int j, int k) + { phi_arr(i,j-1,k) = phi_arr(i,j,k); - } else if (bc_fft[1].first == FFT::Boundary::odd) { - phi_arr(i,j-1,k) = -phi_arr(i,j,k); - } - } else if (j == domhi.y) { - if (bc_fft[1].second == FFT::Boundary::even) { + }); + } else if (bc_fft[1].first == FFT::Boundary::odd) { + ParallelFor(makeSlab(bx2,1,domlo.y), [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + phi_arr(i,j-1,k) = -phi_arr(i,j,k); + }); + } + } // lo x + if (bx.bigEnd(1) >= domhi.y) { + if (bc_fft[1].second == FFT::Boundary::even) { + ParallelFor(makeSlab(bx2,1,domhi.y), [=] AMREX_GPU_DEVICE (int i, int j, int k) + { phi_arr(i,j+1,k) = phi_arr(i,j,k); - } else if (bc_fft[1].second == FFT::Boundary::odd) { - phi_arr(i,j+1,k) = -phi_arr(i,j,k); - } + }); + } else if (bc_fft[1].second == FFT::Boundary::odd) { + ParallelFor(makeSlab(bx2,1,domhi.y), [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + phi_arr(i,j+1,k) = -phi_arr(i,j,k); + }); } - }); + } // lo x + } // mfi - } + } // not periodic in y for (MFIter mfi(phi,true); mfi.isValid(); ++mfi) { diff --git a/Source/Microphysics/ERF_EulerianMicrophysics.H b/Source/Microphysics/ERF_EulerianMicrophysics.H index cb91854b4..edf9c34cd 100644 --- a/Source/Microphysics/ERF_EulerianMicrophysics.H +++ b/Source/Microphysics/ERF_EulerianMicrophysics.H @@ -8,6 +8,7 @@ #include "ERF_NullMoist.H" #include "ERF_SAM.H" #include "ERF_Kessler.H" +#include "ERF_SatAdj.H" #include "ERF_Microphysics.H" /*! \brief Eulerian microphysics interface */ @@ -34,6 +35,8 @@ public: } else if (a_model_type == MoistureType::Kessler || a_model_type == MoistureType::Kessler_NoRain) { SetModel(); + } else if (a_model_type == MoistureType::SatAdj) { + SetModel(); } else if (a_model_type == MoistureType::None) { SetModel(); } else { @@ -108,12 +111,12 @@ public: /*! \brief get the indices and names of moisture model variables for restart at a given level */ - void Get_Qmoist_Restart_Vars ( const int a_lev, /*!< level */ - const SolverChoice& a_sc, /*!< Solver choice object */ - std::vector& a_idx, /*!< indices */ - std::vector& a_names /*!< names */ ) const override + void Get_Qmoist_Restart_Vars (const int a_lev, /*!< level */ + const SolverChoice& a_sc, /*!< Solver choice object */ + std::vector& a_idx, /*!< indices */ + std::vector& a_names /*!< names */ ) const override { - m_moist_model[a_lev]->Qmoist_Restart_Vars( a_sc, a_idx, a_names ); + m_moist_model[a_lev]->Qmoist_Restart_Vars(a_sc, a_idx, a_names); } protected: diff --git a/Source/Microphysics/ERF_Microphysics.H b/Source/Microphysics/ERF_Microphysics.H index 29eab027d..084572e0f 100644 --- a/Source/Microphysics/ERF_Microphysics.H +++ b/Source/Microphysics/ERF_Microphysics.H @@ -58,7 +58,10 @@ public: /*! \brief get the indices and names of moisture model variables for restart at a given level */ - virtual void Get_Qmoist_Restart_Vars ( int, const SolverChoice&, std::vector&, std::vector& ) const = 0; + virtual void Get_Qmoist_Restart_Vars (int, + const SolverChoice&, + std::vector&, + std::vector& ) const = 0; /*! \brief query if a specified moisture model is Eulerian or Lagrangian */ static MoistureModelType modelType (const MoistureType a_moisture_type) @@ -68,6 +71,7 @@ public: || (a_moisture_type == MoistureType::SAM_NoPrecip_NoIce) || (a_moisture_type == MoistureType::Kessler) || (a_moisture_type == MoistureType::Kessler_NoRain) + || (a_moisture_type == MoistureType::SatAdj) || (a_moisture_type == MoistureType::None) ) { return MoistureModelType::Eulerian; } else { diff --git a/Source/Microphysics/Kessler/ERF_Kessler.H b/Source/Microphysics/Kessler/ERF_Kessler.H index f361f013a..f890ef35b 100644 --- a/Source/Microphysics/Kessler/ERF_Kessler.H +++ b/Source/Microphysics/Kessler/ERF_Kessler.H @@ -125,9 +125,9 @@ public: Qstate_Size () override { return Kessler::m_qstate_size; } void - Qmoist_Restart_Vars ( const SolverChoice& a_sc, - std::vector& a_idx, - std::vector& a_names) const override + Qmoist_Restart_Vars (const SolverChoice& a_sc, + std::vector& a_idx, + std::vector& a_names) const override { a_idx.clear(); a_names.clear(); diff --git a/Source/Microphysics/SAM/ERF_SAM.H b/Source/Microphysics/SAM/ERF_SAM.H index e514b509a..7652b3093 100644 --- a/Source/Microphysics/SAM/ERF_SAM.H +++ b/Source/Microphysics/SAM/ERF_SAM.H @@ -264,9 +264,9 @@ public: } void - Qmoist_Restart_Vars ( const SolverChoice& a_sc, - std::vector& a_idx, - std::vector& a_names) const override + Qmoist_Restart_Vars (const SolverChoice& a_sc, + std::vector& a_idx, + std::vector& a_names) const override { a_idx.clear(); a_names.clear(); @@ -283,7 +283,7 @@ private: // Number of qmoist variables (qt, qv, qcl, qci, qp, qpr, qps, qpg) int m_qmoist_size = 11; - // Number of qmoist variables + // Number of qstate variables int m_qstate_size = 6; // CFL MAX for vertical advection diff --git a/Source/Microphysics/SatAdj/ERF_InitSatAdj.cpp b/Source/Microphysics/SatAdj/ERF_InitSatAdj.cpp new file mode 100644 index 000000000..5ce0b67c9 --- /dev/null +++ b/Source/Microphysics/SatAdj/ERF_InitSatAdj.cpp @@ -0,0 +1,76 @@ +#include "ERF_SatAdj.H" + +using namespace amrex; + + +/** + * Initializes the Microphysics module. + * + * @param[in] cons_in Conserved variables input + * @param[in] qc_in Cloud variables input + * @param[in,out] qv_in Vapor variables input + * @param[in] qi_in Ice variables input + * @param[in] grids The boxes on which we will evolve the solution + * @param[in] geom Geometry associated with these MultiFabs and grids + * @param[in] dt_advance Timestep for the advance + */ +void SatAdj::Init (const MultiFab& cons_in, + const BoxArray& /*grids*/, + const Geometry& geom, + const Real& dt_advance, + std::unique_ptr& /*z_phys_nd*/, + std::unique_ptr& /*detJ_cc*/) +{ + dt = dt_advance; + m_geom = geom; + + MicVarMap.resize(m_qmoist_size); + MicVarMap = {MicVar_SatAdj::qv, MicVar_SatAdj::qc}; + + // initialize microphysics variables + for (auto ivar = 0; ivar < MicVar_SatAdj::NumVars; ++ivar) { + mic_fab_vars[ivar] = std::make_shared(cons_in.boxArray(), cons_in.DistributionMap(), + 1, cons_in.nGrowVect()); + mic_fab_vars[ivar]->setVal(0.); + } +} + +/** + * Initializes the Microphysics module. + * + * @param[in] cons_in Conserved variables input + */ +void SatAdj::Copy_State_to_Micro (const MultiFab& cons_in) +{ + // Get the temperature, density, theta, qt and qp from input + for (MFIter mfi(cons_in); mfi.isValid(); ++mfi) { + const auto& tbx = mfi.tilebox(); + + auto states_array = cons_in.array(mfi); + + auto qv_array = mic_fab_vars[MicVar_SatAdj::qv]->array(mfi); + auto qc_array = mic_fab_vars[MicVar_SatAdj::qc]->array(mfi); + + auto rho_array = mic_fab_vars[MicVar_SatAdj::rho]->array(mfi); + auto theta_array = mic_fab_vars[MicVar_SatAdj::theta]->array(mfi); + auto tabs_array = mic_fab_vars[MicVar_SatAdj::tabs]->array(mfi); + auto pres_array = mic_fab_vars[MicVar_SatAdj::pres]->array(mfi); + + // Get pressure, theta, temperature, density, and qt, qp + ParallelFor(tbx, [=] 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); + qv_array(i,j,k) = states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp); + qc_array(i,j,k) = states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp); + + tabs_array(i,j,k) = getTgivenRandRTh(states_array(i,j,k,Rho_comp), + states_array(i,j,k,RhoTheta_comp), + qv_array(i,j,k)); + + // Pressure in [mbar] for qsat evaluation + pres_array(i,j,k) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)) * 0.01; + }); + } +} + diff --git a/Source/Microphysics/SatAdj/ERF_SatAdj.H b/Source/Microphysics/SatAdj/ERF_SatAdj.H new file mode 100644 index 000000000..8b42743b5 --- /dev/null +++ b/Source/Microphysics/SatAdj/ERF_SatAdj.H @@ -0,0 +1,214 @@ +#ifndef ERF_SATADJ_H +#define ERF_SATADJ_H + +/* + * Implementation saturation adjustment microphysics model + * The model transports qv and qc and does Newton iterations + * to complete condensation/evaporation. + */ + +#include +#include +#include + +#include +#include +#include +#include + +#include "ERF_EOS.H" +#include "ERF_Constants.H" +#include "ERF_MicrophysicsUtils.H" +#include "ERF_IndexDefines.H" +#include "ERF_DataStruct.H" +#include "ERF_NullMoist.H" +#include "ERF_TileNoZ.H" + +namespace MicVar_SatAdj { + enum { + // independent variables + rho=0, // density + theta, // liquid/ice water potential temperature + tabs, // temperature + pres, // pressure + // non-precipitating vars + qv, // cloud vapor + qc, // cloud water + NumVars + }; +} + +class SatAdj : public NullMoist { + + using FabPtr = std::shared_ptr; + +public: + // constructor + SatAdj () {} + + // destructor + virtual ~SatAdj () = default; + + // cloud physics + void AdvanceSatAdj (const SolverChoice& /*solverChoice*/); + + // Set up for first time + void + Define (SolverChoice& sc) override + { + m_fac_cond = lcond / sc.c_p; + m_rdOcp = sc.rdOcp; + } + + // init + void + Init (const amrex::MultiFab& cons_in, + const amrex::BoxArray& /*grids*/, + const amrex::Geometry& geom, + const amrex::Real& dt_advance, + std::unique_ptr& /*z_phys_nd*/, + std::unique_ptr& /*detJ_cc*/) override; + + // Copy state into micro vars + void + Copy_State_to_Micro (const amrex::MultiFab& cons_in) override; + + // Copy state into micro vars + void + Copy_Micro_to_State (amrex::MultiFab& cons_in) override; + + // update micro vars + void + Update_Micro_Vars (amrex::MultiFab& cons_in) override + { + this->Copy_State_to_Micro(cons_in); + } + + // update state vars + void + Update_State_Vars (amrex::MultiFab& cons_in) override + { + this->Copy_Micro_to_State(cons_in); + } + + // wrapper to advance micro vars + void + Advance (const amrex::Real& dt_advance, + const SolverChoice& solverChoice) override + { + dt = dt_advance; + + this->AdvanceSatAdj(solverChoice); + } + + amrex::MultiFab* + Qmoist_Ptr (const int& varIdx) override + { + AMREX_ALWAYS_ASSERT(varIdx < m_qmoist_size); + return mic_fab_vars[MicVarMap[varIdx]].get(); + } + + int + Qmoist_Size () override { return SatAdj::m_qmoist_size; } + + int + Qstate_Size () override { return SatAdj::m_qstate_size; } + + void + Qmoist_Restart_Vars (const SolverChoice& /*a_sc*/, + std::vector& a_idx, + std::vector& a_names) const override + { + a_idx.clear(); + a_names.clear(); + } + + AMREX_GPU_HOST_DEVICE + AMREX_FORCE_INLINE + static amrex::Real + NewtonIterSat (int& i, int& j, int& k, + const amrex::Real& fac_cond, + const amrex::Array4& tabs_array, + const amrex::Array4& pres_array, + const amrex::Array4& qv_array, + const amrex::Array4& qc_array) + { + // Solution tolerance + amrex::Real tol = 1.0e-8; + + // Saturation moisture fractions + amrex::Real qsat, dqsat; + + // Newton iteration vars + int niter; + amrex::Real fff, dfff; + amrex::Real dtabs, delta_qv; + + // Initial guess for temperature & pressure + amrex::Real tabs = tabs_array(i,j,k); + amrex::Real pres = pres_array(i,j,k); + + niter = 0; + dtabs = 1; + + //================================================== + // Newton iteration to qv=qsat (cloud phase only) + //================================================== + do { + // Saturation moisture fractions + erf_qsatw(tabs, pres, qsat); + erf_dtqsatw(tabs, pres, dqsat); + + // Function for root finding: + // 0 = -T_new + T_old + L_eff/C_p * (qv - qsat) + fff = -tabs + tabs_array(i,j,k) + fac_cond*(qv_array(i,j,k) - qsat); + + // Derivative of function (T_new iterated on) + dfff = -1.0 - fac_cond*dqsat; + + // Update the temperature + dtabs = -fff/dfff; + tabs += dtabs; + + // Update iteration + niter = niter+1; + } while(std::abs(dtabs) > tol && niter < 20); + + // Update qsat from last iteration (dq = dq/dt * dt) + qsat += dqsat*dtabs; + + // Changes in each component + delta_qv = qv_array(i,j,k) - qsat; + + // Partition the change in non-precipitating q + qv_array(i,j,k) = qsat; + qc_array(i,j,k) += delta_qv; + + // Return to temperature + return tabs; + } + +private: + // Number of qmoist variables (qv, qc) + int m_qmoist_size = 2; + + // Number of qstate variables + int m_qstate_size = 2; + + // MicVar map (Qmoist indices -> MicVar enum) + amrex::Vector MicVarMap; + + // geometry + amrex::Geometry m_geom; + + // timestep + amrex::Real dt; + + // constants + amrex::Real m_fac_cond; + amrex::Real m_rdOcp ; + + // independent variables + amrex::Array mic_fab_vars; +}; +#endif diff --git a/Source/Microphysics/SatAdj/ERF_SatAdj.cpp b/Source/Microphysics/SatAdj/ERF_SatAdj.cpp new file mode 100644 index 000000000..0956e4156 --- /dev/null +++ b/Source/Microphysics/SatAdj/ERF_SatAdj.cpp @@ -0,0 +1,83 @@ +#include "ERF_SatAdj.H" + +using namespace amrex; + +/** + * Compute Precipitation-related Microphysics quantities. + */ +void SatAdj::AdvanceSatAdj (const SolverChoice& /*solverChoice*/) +{ + auto tabs = mic_fab_vars[MicVar_SatAdj::tabs]; + + // Expose for GPU + Real d_fac_cond = m_fac_cond; + Real rdOcp = m_rdOcp; + + // get the temperature, dentisy, theta, qt and qc from input + for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + + const auto& tbx = mfi.tilebox(); + + auto qv_array = mic_fab_vars[MicVar_SatAdj::qv]->array(mfi); + auto qc_array = mic_fab_vars[MicVar_SatAdj::qc]->array(mfi); + auto tabs_array = mic_fab_vars[MicVar_SatAdj::tabs]->array(mfi); + auto theta_array = mic_fab_vars[MicVar_SatAdj::theta]->array(mfi); + auto pres_array = mic_fab_vars[MicVar_SatAdj::pres]->array(mfi); + + ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k)); + + //------- Evaporation/condensation + Real qsat; + erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsat); + + // There is enough moisutre to drive to equilibrium + if ((qv_array(i,j,k)+qc_array(i,j,k)) > qsat) { + + // Update temperature + tabs_array(i,j,k) = NewtonIterSat(i, j, k , + d_fac_cond, tabs_array, pres_array, + qv_array , qc_array ); + + // Update theta (constant pressure) + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); + + // + // We cannot blindly relax to qsat, but we can convert qc/qi -> qv. + // The concept here is that if we put all the moisture into qv and modify + // the temperature, we can then check if qv > qsat occurs (for final T/P/qv). + // If the reduction in T/qsat and increase in qv does trigger the + // aforementioned condition, we can do Newton iteration to drive qv = qsat. + // + } else { + // Changes in each component + Real delta_qc = qc_array(i,j,k); + + // Partition the change in non-precipitating q + qv_array(i,j,k) += qc_array(i,j,k); + qc_array(i,j,k) = 0.0; + + // Update temperature (endothermic since we evap/sublime) + tabs_array(i,j,k) -= d_fac_cond * delta_qc; + + // Update theta + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); + + // Verify assumption that qv > qsat does not occur + erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsat); + if (qv_array(i,j,k) > qsat) { + + // Update temperature + tabs_array(i,j,k) = NewtonIterSat(i, j, k , + d_fac_cond , tabs_array, pres_array, + qv_array , qc_array ); + + // Update theta + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); + + } + } + }); + } +} diff --git a/Source/Microphysics/SatAdj/ERF_UpdateSatAdj.cpp b/Source/Microphysics/SatAdj/ERF_UpdateSatAdj.cpp new file mode 100644 index 000000000..b45d619ef --- /dev/null +++ b/Source/Microphysics/SatAdj/ERF_UpdateSatAdj.cpp @@ -0,0 +1,37 @@ +#include "ERF_SatAdj.H" + +using namespace amrex; + +/** + * Updates conserved and microphysics variables in the provided MultiFabs from + * the internal MultiFabs that store Microphysics module data. + * + * @param[out] cons Conserved variables + * @param[out] qmoist: qv, qc + */ +void SatAdj::Copy_Micro_to_State (MultiFab& cons) +{ + // Get the temperature, density, theta, qt and qp from input + for (MFIter mfi(cons,TilingIfNotGPU()); mfi.isValid(); ++mfi) { + const auto& tbx = mfi.tilebox(); + + auto states_arr = cons.array(mfi); + + auto rho_arr = mic_fab_vars[MicVar_SatAdj::rho]->array(mfi); + auto theta_arr = mic_fab_vars[MicVar_SatAdj::theta]->array(mfi); + auto qv_arr = mic_fab_vars[MicVar_SatAdj::qv]->array(mfi); + auto qc_arr = mic_fab_vars[MicVar_SatAdj::qc]->array(mfi); + + // get potential total density, temperature, qt, qp + ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k); + states_arr(i,j,k,RhoQ1_comp) = rho_arr(i,j,k)*qv_arr(i,j,k); + states_arr(i,j,k,RhoQ2_comp) = rho_arr(i,j,k)*qc_arr(i,j,k); + }); + } + + // Fill interior ghost cells and periodic boundaries + cons.FillBoundary(m_geom.periodicity()); +} + diff --git a/Source/Microphysics/SatAdj/Make.package b/Source/Microphysics/SatAdj/Make.package new file mode 100644 index 000000000..5b9ff43e4 --- /dev/null +++ b/Source/Microphysics/SatAdj/Make.package @@ -0,0 +1,4 @@ +CEXE_sources += ERF_InitSatAdj.cpp +CEXE_sources += ERF_UpdateSatAdj.cpp +CEXE_sources += ERF_SatAdj.cpp +CEXE_headers += ERF_SatAdj.H diff --git a/Source/Particles/ERFPC.H b/Source/Particles/ERFPC.H index 986e04ac0..dd3ef7cd9 100644 --- a/Source/Particles/ERFPC.H +++ b/Source/Particles/ERFPC.H @@ -223,6 +223,8 @@ class ERFPC : public amrex::ParticleContainer< ERFParticlesRealIdxAoS::ncomps, std::string m_initialization_type; /*!< initial particle distribution type */ int m_ppc_init; /*!< initial number of particles per cell */ + bool m_stable_redistribute; /*!< use stable redistribute for deterministic simulations */ + /*! read inputs from file */ virtual void readInputs (); diff --git a/Source/Particles/ERFPCInitializations.cpp b/Source/Particles/ERFPCInitializations.cpp index feca8c476..3d4abcc91 100644 --- a/Source/Particles/ERFPCInitializations.cpp +++ b/Source/Particles/ERFPCInitializations.cpp @@ -49,6 +49,14 @@ void ERFPC::readInputs () m_advect_w_gravity = (m_name == ERFParticleNames::hydro ? true : false); pp.query("advect_with_gravity", m_advect_w_gravity); + m_stable_redistribute = false; + pp.query("stable_redistribute", m_stable_redistribute); + + if (m_stable_redistribute) { + Print() << "Note: using stable Redistribute() for particle container " << m_name << ".\n"; + } + setStableRedistribute(m_stable_redistribute); + return; } @@ -65,6 +73,7 @@ void ERFPC::InitializeParticles (const std::unique_ptr& a_height_ptr) << m_name << " particle species.\n"; Error("See error message!"); } + return; } diff --git a/Source/SourceTerms/ERF_MakeMomSources.cpp b/Source/SourceTerms/ERF_MakeMomSources.cpp index 8133f0482..a6c05c654 100644 --- a/Source/SourceTerms/ERF_MakeMomSources.cpp +++ b/Source/SourceTerms/ERF_MakeMomSources.cpp @@ -53,6 +53,7 @@ void make_mom_sources (int level, MultiFab & zmom_src, const MultiFab& base_state, MultiFab* forest_drag, + MultiFab* terrain_blank, const Geometry geom, const SolverChoice& solverChoice, std::unique_ptr& /*mapfac_m*/, @@ -86,10 +87,21 @@ void make_mom_sources (int level, // 6. nudging towards input sounding data // 7. numerical diffusion for (xmom,ymom,zmom) // 8. sponge + // 9. Forest canopy + // 10. Immersed Forcing // ***************************************************************************** - const bool l_use_ndiff = solverChoice.use_NumDiff; - const bool use_terrain = solverChoice.terrain_type != TerrainType::None; - const bool l_do_forest = solverChoice.do_forest; + const bool l_use_ndiff = solverChoice.use_NumDiff; + const bool l_use_zphys = (solverChoice.mesh_type != MeshType::ConstantDz); + const bool l_do_forest_drag = solverChoice.do_forest_drag; + const bool l_do_terrain_drag = solverChoice.do_terrain_drag; + + // Check if terrain and immersed terrain clash + if(l_use_zphys && l_do_terrain_drag){ + amrex::Error(" Cannot use immersed forcing for terrain with terrain-fitted coordinates"); + } + if(l_do_forest_drag && l_do_terrain_drag){ + amrex::Error(" Currently forest canopy cannot be used with immersed forcing"); + } // ***************************************************************************** // Data for Coriolis forcing @@ -228,10 +240,12 @@ void make_mom_sources (int level, const Array4& f_drag_arr = (forest_drag) ? forest_drag->const_array(mfi) : Array4{}; + const Array4& t_blank_arr = (terrain_blank) ? terrain_blank->const_array(mfi) : + Array4{}; - const Array4& z_nd_arr = (use_terrain) ? z_phys_nd->const_array(mfi) : + const Array4& z_nd_arr = (l_use_zphys) ? z_phys_nd->const_array(mfi) : Array4{}; - const Array4& z_cc_arr = (use_terrain) ? z_phys_cc->const_array(mfi) : + const Array4& z_cc_arr = (l_use_zphys) ? z_phys_cc->const_array(mfi) : Array4{}; // ***************************************************************************** @@ -491,7 +505,7 @@ void make_mom_sources (int level, // ***************************************************************************** // 9. Add CANOPY source terms // ***************************************************************************** - if (l_do_forest) { + if (l_do_forest_drag) { ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { const Real ux = u(i, j, k); @@ -526,5 +540,48 @@ void make_mom_sources (int level, zmom_src_arr(i, j, k) -= f_drag * uz * windspeed; }); } + // ***************************************************************************** + // 10. Add Immersed source terms + // ***************************************************************************** + if (l_do_terrain_drag) { + const Real drag_coefficient=10.0/dz; + const Real tiny = std::numeric_limits::epsilon(); + ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + const Real ux = u(i, j, k); + const Real uy = 0.25 * ( v(i, j , k ) + v(i-1, j , k ) + + v(i, j+1, k ) + v(i-1, j+1, k ) ); + const Real uz = 0.25 * ( w(i, j , k ) + w(i-1, j , k ) + + w(i, j , k+1) + w(i-1, j , k+1) ); + const Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz); + const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i-1, j, k)); + const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0); + xmom_src_arr(i, j, k) -= t_blank * CdM * ux * windspeed; + }); + ParallelFor(tby, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + const Real ux = 0.25 * ( u(i , j , k ) + u(i , j-1, k ) + + u(i+1, j , k ) + u(i+1, j-1, k ) ); + const Real uy = v(i, j, k); + const Real uz = 0.25 * ( w(i , j , k ) + w(i , j-1, k ) + + w(i , j , k+1) + w(i , j-1, k+1) ); + const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz); + const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j-1, k)); + const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0); + ymom_src_arr(i, j, k) -= t_blank * CdM * uy * windspeed; + }); + ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + const amrex::Real ux = 0.25 * ( u(i , j , k ) + u(i+1, j , k ) + + u(i , j , k-1) + u(i+1, j , k-1) ); + const amrex::Real uy = 0.25 * ( v(i , j , k ) + v(i , j+1, k ) + + v(i , j , k-1) + v(i , j+1, k-1) ); + const amrex::Real uz = w(i, j, k); + const amrex::Real windspeed = std::sqrt(ux * ux + uy * uy + uz * uz); + const Real t_blank = 0.5 * (t_blank_arr(i, j, k) + t_blank_arr(i, j, k-1)); + const Real CdM = std::min(drag_coefficient / (windspeed + tiny), 1000.0); + zmom_src_arr(i, j, k) -= t_blank * CdM * uz * windspeed; + }); + } } // mfi } diff --git a/Source/SourceTerms/ERF_SrcHeaders.H b/Source/SourceTerms/ERF_SrcHeaders.H index 5ea93913d..4992e80cd 100644 --- a/Source/SourceTerms/ERF_SrcHeaders.H +++ b/Source/SourceTerms/ERF_SrcHeaders.H @@ -59,6 +59,7 @@ void make_mom_sources (int level, int nrk, amrex::MultiFab& zmom_source, const amrex::MultiFab& base_state, amrex::MultiFab* forest_drag, + amrex::MultiFab* terrain_blank, const amrex::Geometry geom, const SolverChoice& solverChoice, std::unique_ptr& mapfac_m, diff --git a/Source/SourceTerms/ERF_TerrainDrag.H b/Source/SourceTerms/ERF_TerrainDrag.H new file mode 100644 index 000000000..289db4a55 --- /dev/null +++ b/Source/SourceTerms/ERF_TerrainDrag.H @@ -0,0 +1,33 @@ +#ifndef ERF_TERRAINDRAG_H_ +#define ERF_TERRAINDRAG_H_ + +#include +#include + +/* + Adapted from: https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2020MS002141 + */ +class TerrainDrag +{ +public: + + explicit TerrainDrag (std::string terrainfile); + + ~TerrainDrag () = default; + + void + define_terrain_blank_field (const amrex::BoxArray& ba, + const amrex::DistributionMapping& dm, + amrex::Geometry& geom, + amrex::MultiFab* z_phys_nd); + + amrex::MultiFab* + get_terrain_blank_field () { return m_terrain_blank.get(); } + +private: + amrex::Vector m_x_terrain; + amrex::Vector m_y_terrain; + amrex::Vector m_height_terrain; + std::unique_ptr m_terrain_blank; +}; +#endif diff --git a/Source/SourceTerms/ERF_TerrainDrag.cpp b/Source/SourceTerms/ERF_TerrainDrag.cpp new file mode 100644 index 000000000..36bc04b07 --- /dev/null +++ b/Source/SourceTerms/ERF_TerrainDrag.cpp @@ -0,0 +1,97 @@ +#include + +using namespace amrex; + +/* + Constructor to get the terrain parameters: + TreeType xc, yc, height +*/ +TerrainDrag::TerrainDrag(std::string terrainfile) +{ + std::ifstream file(terrainfile, std::ios::in); + if (!file.good()) { + Abort("Cannot find terrain file: " + terrainfile); + } + // xc yc height + Real value1, value2, value3; + while (file >> value1 >> value2 >> value3) { + m_x_terrain.push_back(value1); + m_y_terrain.push_back(value2); + m_height_terrain.push_back(value3); + } + file.close(); +} + +void +TerrainDrag::define_terrain_blank_field( + const BoxArray& ba, + const DistributionMapping& dm, + Geometry& geom, + MultiFab* z_phys_nd) +{ + // Geometry params + const auto& dx = geom.CellSizeArray(); + const auto& prob_lo = geom.ProbLoArray(); + + // Allocate the terrain blank MF + // NOTE: 1 ghost cell for averaging to faces + m_terrain_blank.reset(); + m_terrain_blank = std::make_unique(ba, dm, 1, 1); + m_terrain_blank->setVal(0.); + const auto xterrain_size = m_x_terrain.size(); + amrex::Gpu::DeviceVector d_xterrain(xterrain_size); + amrex::Gpu::DeviceVector d_yterrain(xterrain_size); + amrex::Gpu::DeviceVector d_zterrain(xterrain_size); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_x_terrain.begin(), m_x_terrain.end(), + d_xterrain.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_y_terrain.begin(), m_y_terrain.end(), + d_yterrain.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_height_terrain.begin(), m_height_terrain.end(), + d_zterrain.begin()); + const auto* xterrain_ptr = d_xterrain.data(); + const auto* yterrain_ptr = d_yterrain.data(); + const auto* zterrain_ptr = d_zterrain.data(); + // Set the terrain blank data + for (MFIter mfi(*m_terrain_blank); mfi.isValid(); ++mfi) { + Box gtbx = mfi.growntilebox(); + const Array4& levelBlank = m_terrain_blank->array(mfi); + const Array4& z_nd = + (z_phys_nd) ? z_phys_nd->const_array(mfi) : Array4{}; + ParallelFor(gtbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + // Loop over terrain points + for (unsigned ii = 0; ii < xterrain_size; ++ii) { + Real ht = zterrain_ptr[ii]; + Real xt = xterrain_ptr[ii]; + Real yt = yterrain_ptr[ii]; + // Physical positions of cell-centers + const Real x = prob_lo[0] + (i + 0.5) * dx[0]; + const Real y = prob_lo[1] + (j + 0.5) * dx[1]; + Real z = prob_lo[2] + (k + 0.5) * dx[2]; + if (z_nd) { + 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)); + } + z = std::max(z, 0.0); + const Real cell_radius = std::sqrt(dx[0] * dx[0] + dx[1] * dx[1]); + // Proximity to the terrain + const Real radius = + std::sqrt((x - xt) * (x - xt) + (y - yt) * (y - yt)); + const Real terrain_point = + (radius <= cell_radius && z <= ht) ? 1.0 : 0.0; + levelBlank(i, j, k) = terrain_point; + if(terrain_point==1){ + break; + } + } + }); + } // mfi + +// Fillboundary for periodic ghost cell copy +m_terrain_blank->FillBoundary(geom.periodicity()); + +} // init_terrain_blank_field diff --git a/Source/SourceTerms/Make.package b/Source/SourceTerms/Make.package index d70c5c9e0..92b757433 100644 --- a/Source/SourceTerms/Make.package +++ b/Source/SourceTerms/Make.package @@ -6,6 +6,7 @@ CEXE_sources += ERF_ApplySpongeZoneBCs.cpp CEXE_sources += ERF_ApplySpongeZoneBCs_ReadFromFile.cpp CEXE_sources += ERF_NumericalDiffusion.cpp CEXE_sources += ERF_ForestDrag.cpp +CEXE_sources += ERF_TerrainDrag.cpp ifeq ($(USE_NETCDF),TRUE) CEXE_sources += ERF_MoistSetRhs.cpp @@ -15,3 +16,4 @@ CEXE_headers += ERF_NumericalDiffusion.H CEXE_headers += ERF_SrcHeaders.H CEXE_headers += ERF_BuoyancyUtils.H CEXE_headers += ERF_ForestDrag.H +CEXE_headers += ERF_TerrainDrag.H \ No newline at end of file diff --git a/Source/TimeIntegration/ERF_ComputeTimestep.cpp b/Source/TimeIntegration/ERF_ComputeTimestep.cpp index 96f37c1e4..5c79fb85a 100644 --- a/Source/TimeIntegration/ERF_ComputeTimestep.cpp +++ b/Source/TimeIntegration/ERF_ComputeTimestep.cpp @@ -115,15 +115,16 @@ ERF::estTimeStep (int level, long& dt_fast_ratio) const // If we are doing implicit acoustic substepping, then the z-direction does not contribute // to the computation of the time step if (l_implicit_substepping) { - if (nxc > 1 && nyc > 1) { - new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0))+c)*dxinv[0]), - ((amrex::Math::abs(u(i,j,k,1))+c)*dxinv[1]), new_comp_dt); - } else if (nxc > 1) { + if ((nxc > 1) && (nyc==1)) { + // 2-D in x-z new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0))+c)*dxinv[0]), new_comp_dt); - } else if (nyc > 1) { + } else if ((nyc > 1) && (nxc==1)) { + // 2-D in y-z new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,1))+c)*dxinv[1]), new_comp_dt); } else { - amrex::Abort("Not sure how to compute dt for this case"); + // 3-D or SCM + new_comp_dt = amrex::max(((amrex::Math::abs(u(i,j,k,0))+c)*dxinv[0]), + ((amrex::Math::abs(u(i,j,k,1))+c)*dxinv[1]), new_comp_dt); } // If we are not doing implicit acoustic substepping, then the z-direction contributes @@ -201,8 +202,9 @@ ERF::estTimeStep (int level, long& dt_fast_ratio) const if (fixed_dt[level] > 0. && fixed_fast_dt[level] > 0.) { dt_fast_ratio = static_cast( fixed_dt[level] / fixed_fast_dt[level] ); } else if (fixed_dt[level] > 0.) { - // Max CFL_c = 1.0 for substeps, but we enforce a min of 4 substeps - auto dt_sub_max = (estdt_comp/cfl); + // Max CFL_c = 1.0 for substeps by default, but we enforce a min of 4 substeps + auto dt_sub_max = (estdt_comp/cfl * sub_cfl); + Print() << "fixed_dt="<