Kessler Microphysics (#1324)

* Changes to EOS.H with moisture * Change all functions in EOS.H to work with moisture * Adding qp for buoyancy type 1 * Updating docs for buoyancy term * Updating docs for buoyancy term * Update docs for buoyancy * Avoiding extra divide * Update buoyancy docs * Some corrections to buoyancy docs * Some corrections to Docs in buoyancy * Correcting docs issue * Correcting docs * adding Kessler microphysics docs * adding Kessler microphysics docs * adding Kessler microphysics docs * adding Kessler microphysics docs * WIP: Adding Kessler microphysics * Guards for microphysics * WIP: Adding Kessler microphysics * WIP: kessler microphysics * Compiled and verified initial condition for Kessler * theta pert to 0 for hydrostatic case * Calling moist initial conditions * Call moist for moisture cases * Pert to density and theta * Avoiding name clashes * Complete Kessler microphysics draft * Updating squall line test case * Moving towards Newton Raphson for IC * Updated Kessler microphysics for push * Updating inputs * Adding inputs with outflow bcs * Updating periodic inputs * Removing old prob.cpp * Minor changes and removing whitespaces and tabs * Updating periodic inputs * Finalizing squall line Kessler * Updating inputs for Kessler * Updating squall line test case * enforce qv >= 0.0 * Comment out unused vars in problemBase. * Initi1D an only access use_moist_bkgd when compiled with moisture support. * Fix HIP issues. * Add Kessler to CMake compile. * Temporary patch, doesn't fix device capture. * Not ideal but compiled with nvcc. * Use PI instead of M_PI for windows. --------- Co-authored-by: Mahesh Natarajan <[email protected]> Co-authored-by: Ann Almgren <[email protected]> Co-authored-by: Aaron Lattanzi <[email protected]>
erf-model · Nov 29, 2023 · 0e5e440 · 0e5e440
1 parent 7ecb265
commit 0e5e440
Show file tree

Hide file tree

Showing 25 changed files with 1,612 additions and 18 deletions.
diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake
@@ -65,7 +65,11 @@ function(build_erf_lib erf_lib_name)
        ${SRC_DIR}/Microphysics/SAM/Precip.cpp
        ${SRC_DIR}/Microphysics/SAM/PrecipFall.cpp
        ${SRC_DIR}/Microphysics/SAM/Diagnose.cpp
-       ${SRC_DIR}/Microphysics/SAM/Update.cpp)
+       ${SRC_DIR}/Microphysics/SAM/Update.cpp
+       ${SRC_DIR}/Microphysics/Kessler/Init_Kessler.cpp
+       ${SRC_DIR}/Microphysics/Kessler/Kessler.cpp
+       ${SRC_DIR}/Microphysics/Kessler/Diagnose_Kessler.cpp
+       ${SRC_DIR}/Microphysics/Kessler/Update_Kessler.cpp)
     target_compile_definitions(${erf_lib_name} PUBLIC ERF_USE_MOISTURE)
   endif()
 
@@ -187,6 +191,7 @@ function(build_erf_lib erf_lib_name)
     target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics)
     target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/Null)
     target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/SAM)
+    target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Microphysics/Kessler)
   endif()
 
   if(ERF_ENABLE_RRTMGP)

diff --git a/Exec/CMakeLists.txt b/Exec/CMakeLists.txt
@@ -8,6 +8,7 @@ else ()
   add_subdirectory(ABL)
   add_subdirectory(SuperCell)
   add_subdirectory(Radiation)
+  add_subdirectory(SquallLine_2D)
   add_subdirectory(RegTests/Bubble)
   add_subdirectory(RegTests/CouetteFlow)
   add_subdirectory(RegTests/DensityCurrent)

diff --git a/Exec/Make.ERF b/Exec/Make.ERF
@@ -117,6 +117,11 @@ ifeq ($(USE_MOISTURE), TRUE)
   include $(ERF_MOISTURE_SAM_DIR)/Make.package
   VPATH_LOCATIONS   += $(ERF_MOISTURE_SAM_DIR)
   INCLUDE_LOCATIONS += $(ERF_MOISTURE_SAM_DIR) 
+
+  ERF_MOISTURE_KESSLER_DIR = $(ERF_SOURCE_DIR)/Microphysics/Kessler
+  include $(ERF_MOISTURE_KESSLER_DIR)/Make.package
+  VPATH_LOCATIONS   += $(ERF_MOISTURE_KESSLER_DIR)
+  INCLUDE_LOCATIONS += $(ERF_MOISTURE_KESSLER_DIR)
 endif
 
 ifeq ($(USE_WARM_NO_PRECIP), TRUE)

diff --git a/Exec/SquallLine_2D/CMakeLists.txt b/Exec/SquallLine_2D/CMakeLists.txt
@@ -0,0 +1,12 @@
+set(erf_exe_name squallline_2d)
+
+add_executable(${erf_exe_name} "")
+target_sources(${erf_exe_name}
+   PRIVATE
+     prob.cpp
+)
+
+target_include_directories(${erf_exe_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
+
+include(${CMAKE_SOURCE_DIR}/CMake/BuildERFExe.cmake)
+build_erf_exe(${erf_exe_name})
diff --git a/Exec/SquallLine_2D/GNUmakefile b/Exec/SquallLine_2D/GNUmakefile
@@ -0,0 +1,35 @@
+# AMReX
+COMP = gnu
+PRECISION = DOUBLE
+
+# Profiling
+PROFILE       = FALSE
+TINY_PROFILE  = FALSE
+COMM_PROFILE  = FALSE
+TRACE_PROFILE = FALSE
+MEM_PROFILE   = FALSE
+USE_GPROF     = FALSE
+
+# Performance
+USE_MPI  = TRUE
+USE_OMP  = FALSE
+
+USE_CUDA = FALSE
+USE_HIP  = FALSE
+USE_SYCL = FALSE
+
+# Debugging
+DEBUG = FALSE
+
+TEST = TRUE
+USE_ASSERTION = TRUE
+
+USE_MOISTURE = TRUE
+#USE_WARM_NO_PRECIP = TRUE
+
+# GNU Make
+Bpack := ./Make.package
+Blocs := .
+ERF_HOME := ../..
+ERF_PROBLEM_DIR = $(ERF_HOME)/Exec/SquallLine_2D
+include $(ERF_HOME)/Exec/Make.ERF
diff --git a/Exec/SquallLine_2D/Make.package b/Exec/SquallLine_2D/Make.package
@@ -0,0 +1,2 @@
+CEXE_headers += prob.H
+CEXE_sources += prob.cpp
diff --git a/Exec/SquallLine_2D/README b/Exec/SquallLine_2D/README
@@ -0,0 +1,2 @@
+This problem setup is the evolution of a supercell, which primarily tests the ability
+of ERF to model moisture physics.
diff --git a/Exec/SquallLine_2D/inputs_moisture b/Exec/SquallLine_2D/inputs_moisture
@@ -0,0 +1,84 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 9000
+stop_time = 90000.0
+
+amrex.fpe_trap_invalid = 1
+
+fabarray.mfiter_tile_size = 2048 1024 2048
+
+# PROBLEM SIZE & GEOMETRY
+geometry.prob_lo     = -25000.   0.    0.
+geometry.prob_hi     =  25000. 400. 20000.
+amr.n_cell           =  192    4    81    # dx=dy=dz=100 m
+
+# periodic in x to match WRF setup
+# - as an alternative, could use symmetry at x=0 and outflow at x=25600
+geometry.is_periodic = 1 1 0
+#xlo.type = "Outflow"
+#xhi.type = "Outflow"
+zlo.type = "SlipWall"
+zhi.type = "Outflow"
+
+erf.sponge_strength = 2.0
+#erf.use_zhi_sponge_damping = true
+erf.zhi_sponge_start = 12000.0
+
+erf.sponge_density = 1.2
+erf.sponge_x_velocity = 0.0
+erf.sponge_y_velocity = 0.0
+erf.sponge_z_velocity = 0.0
+
+# TIME STEP CONTROL
+erf.use_native_mri = 1
+erf.fixed_dt       = 1.0      # fixed time step [s] -- Straka et al 1993
+erf.fixed_fast_dt  = 0.5     # fixed time step [s] -- Straka et al 1993
+#erf.no_substepping  = 1
+
+# 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
+amr.check_file      = chk        # root name of checkpoint file
+amr.check_int       = 1000       # number of timesteps between checkpoints
+#amr.restart         = chk12000
+
+# PLOTFILES
+erf.plot_file_1         = plt        # root name of plotfile
+erf.plot_int_1          = 100         # number of timesteps between plotfiles
+erf.plot_vars_1         = density rhotheta rhoQt rhoQp x_velocity y_velocity z_velocity pressure theta temp qt qp qv qc qi scalar pert_dens
+
+# SOLVER CHOICE
+erf.use_gravity = true
+erf.buoyancy_type = 4
+erf.use_coriolis = false
+erf.use_rayleigh_damping = false
+
+#erf.les_type = "Smagorinsky"
+erf.Cs              = 0.25
+erf.les_type = "None"
+
+#
+# diffusion coefficient from Straka, K = 75 m^2/s
+#
+erf.molec_diff_type = "ConstantAlpha"
+#erf.molec_diff_type = "Constant"
+erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities
+erf.dynamicViscosity = 200.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s
+erf.alpha_T = 00.0 # [m^2/s]
+erf.alpha_C = 50.0
+
+erf.moisture_model = "Kessler"
+erf.use_moist_background = true
+
+erf.moistscal_horiz_adv_string = "Centered_2nd"
+erf.moistscal_vert_adv_string = "Centered_2nd"
+
+# PROBLEM PARAMETERS (optional)
+prob.T_0 = 300.0
+prob.U_0 = 0
+prob.T_pert = 3
diff --git a/Exec/SquallLine_2D/inputs_moisture_outflow b/Exec/SquallLine_2D/inputs_moisture_outflow
@@ -0,0 +1,86 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 9000
+stop_time = 90000.0
+
+amrex.fpe_trap_invalid = 1
+
+fabarray.mfiter_tile_size = 2048 1024 2048
+
+# PROBLEM SIZE & GEOMETRY
+geometry.prob_lo     = -25000.   0.    0.
+geometry.prob_hi     =  25000. 400. 20000.
+amr.n_cell           =  192    4    81    # dx=dy=dz=100 m
+
+# periodic in x to match WRF setup
+# - as an alternative, could use symmetry at x=0 and outflow at x=25600
+geometry.is_periodic = 0 1 0
+xlo.type = "Outflow"
+xhi.type = "Outflow"
+zlo.type = "SlipWall"
+zhi.type = "Outflow"
+
+erf.sponge_strength = 2.0
+#erf.use_zhi_sponge_damping = true
+erf.zhi_sponge_start = 12000.0
+
+erf.sponge_density = 1.2
+erf.sponge_x_velocity = 0.0
+erf.sponge_y_velocity = 0.0
+erf.sponge_z_velocity = 0.0
+
+
+
+# TIME STEP CONTROL
+erf.use_native_mri = 1
+erf.fixed_dt       = 1.0      # fixed time step [s] -- Straka et al 1993
+erf.fixed_fast_dt  = 0.5     # fixed time step [s] -- Straka et al 1993
+#erf.no_substepping  = 1
+
+# 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
+amr.check_file      = chk        # root name of checkpoint file
+amr.check_int       = 1000       # number of timesteps between checkpoints
+#amr.restart         = chk09000
+
+# PLOTFILES
+erf.plot_file_1         = plt        # root name of plotfile
+erf.plot_int_1          = 100         # number of timesteps between plotfiles
+erf.plot_vars_1         = density rhotheta rhoQt rhoQp x_velocity y_velocity z_velocity pressure theta temp qt qp qv qc qi scalar pert_dens
+
+# SOLVER CHOICE
+erf.use_gravity = true
+erf.buoyancy_type = 4
+erf.use_coriolis = false
+erf.use_rayleigh_damping = false
+
+#erf.les_type = "Smagorinsky"
+erf.Cs              = 0.25
+erf.les_type = "None"
+
+#
+# diffusion coefficient from Straka, K = 75 m^2/s
+#
+erf.molec_diff_type = "ConstantAlpha"
+#erf.molec_diff_type = "Constant"
+erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities
+erf.dynamicViscosity = 200.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s
+erf.alpha_T = 00.0 # [m^2/s]
+erf.alpha_C = 0.0
+
+erf.moisture_model = "Kessler"
+erf.use_moist_background = true
+
+erf.moistscal_horiz_adv_string = "Centered_2nd"
+erf.moistscal_vert_adv_string = "Centered_2nd"
+
+# PROBLEM PARAMETERS (optional)
+prob.T_0 = 300.0
+prob.U_0 = 0
+prob.T_pert = 3
diff --git a/Exec/SquallLine_2D/prob.H b/Exec/SquallLine_2D/prob.H
@@ -0,0 +1,89 @@
+#ifndef _PROB_H_
+#define _PROB_H_
+
+#include <string>
+
+#include "AMReX_REAL.H"
+#include "AMReX_ParmParse.H"
+#include "AMReX_MultiFab.H"
+
+#include "prob_common.H"
+#include "EOS.H"
+#include "IndexDefines.H"
+#include "TileNoZ.H"
+
+struct ProbParm : ProbParmDefaults {
+  amrex::Real T_0 = 300.0; // surface temperature == mean potential temperature
+  amrex::Real U_0 = 10.0;
+  amrex::Real V_0 = 0.0;
+  amrex::Real x_c = 0.0; // center of thermal perturbation
+  amrex::Real z_c = 3200.0;
+  amrex::Real x_r = 1000.0;
+  amrex::Real z_r = 1000.0;
+  amrex::Real T_pert = 5.0; // perturbation temperature
+  // overridden physical constants
+  amrex::Real C_p = 1004.0;
+}; // namespace ProbParm
+
+
+class Problem : public ProblemBase
+{
+public:
+    Problem();
+
+#include "Prob/init_density_hse_dry.H"
+
+    void init_custom_pert (
+        const amrex::Box&  bx,
+        const amrex::Box& xbx,
+        const amrex::Box& ybx,
+        const amrex::Box& zbx,
+        amrex::Array4<amrex::Real      > const& state,
+        amrex::Array4<amrex::Real      > const& x_vel,
+        amrex::Array4<amrex::Real      > const& y_vel,
+        amrex::Array4<amrex::Real      > const& z_vel,
+        amrex::Array4<amrex::Real      > const& r_hse,
+        amrex::Array4<amrex::Real      > const& p_hse,
+        amrex::Array4<amrex::Real const> const& z_nd,
+        amrex::Array4<amrex::Real const> const& z_cc,
+#if defined(ERF_USE_MOISTURE)
+        amrex::Array4<amrex::Real      > const& qv,
+        amrex::Array4<amrex::Real      > const& qc,
+        amrex::Array4<amrex::Real      > const& qi,
+#elif defined(ERF_USE_WARM_NO_PRECIP)
+        amrex::Array4<amrex::Real      > const& qv,
+        amrex::Array4<amrex::Real      > const& qc,
+#endif
+        amrex::GeometryData const& geomdata,
+        amrex::Array4<amrex::Real const> const& mf_m,
+        amrex::Array4<amrex::Real const> const& mf_u,
+        amrex::Array4<amrex::Real const> const& mf_v,
+        const SolverChoice& sc) override;
+
+
+    void erf_init_dens_hse_moist (amrex::MultiFab& rho_hse,
+                                  std::unique_ptr<amrex::MultiFab>& z_phys_nd,
+                                  std::unique_ptr<amrex::MultiFab>& z_phys_cc,
+                                  amrex::Geometry const& geom) override;
+
+    void init_custom_terrain (
+        const amrex::Geometry& geom,
+        amrex::MultiFab& z_phys_nd,
+        const amrex::Real& time) override;
+
+    void erf_init_rayleigh (
+        amrex::Vector<amrex::Real>& tau,
+        amrex::Vector<amrex::Real>& ubar,
+        amrex::Vector<amrex::Real>& vbar,
+        amrex::Vector<amrex::Real>& wbar,
+        amrex::Vector<amrex::Real>& thetabar,
+        amrex::Geometry      const& geom) override;
+
+protected:
+    std::string name () override { return "Supercell"; }
+
+private:
+    ProbParm parms;
+};
+
+#endif
Original file line number	Diff line number	Diff line change
		@@ -0,0 +1,2 @@
		CEXE_headers += prob.H
		CEXE_sources += prob.cpp
Original file line number	Diff line number	Diff line change
		@@ -0,0 +1,2 @@
		This problem setup is the evolution of a supercell, which primarily tests the ability
		of ERF to model moisture physics.