From 7c03e4b5ddc4a79e135434cd8f944435d4330feb Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Tue, 10 Oct 2023 12:58:18 -0700 Subject: [PATCH 01/37] Changes to EOS.H with moisture --- Source/EOS.H | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/Source/EOS.H b/Source/EOS.H index ad53ddb30..189303e7e 100644 --- a/Source/EOS.H +++ b/Source/EOS.H @@ -13,10 +13,11 @@ * @params[in] rhotheta density times potential temperature theta */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta) +amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta, const amrex::Real qv=0.0) { - amrex::Real p_loc = p_0 * std::pow(R_d * rhotheta * ip_0, Gamma); - return p_loc / (R_d * rho); + // rho and rhotheta are dry values. We should be using moist values with moisture + amrex::Real p_loc = p_0 * std::pow(R_d * rhotheta * (1.0 + R_v/R_d*qv) * ip_0, Gamma); + return p_loc / (R_d * rho * (1.0 + R_v/R_d*qv) ); } /** @@ -27,9 +28,9 @@ amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta) * @params[in] rd0cp ratio of R_d to c_p */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -amrex::Real getThgivenRandT(const amrex::Real rho, const amrex::Real T, const amrex::Real rdOcp) +amrex::Real getThgivenRandT(const amrex::Real rho, const amrex::Real T, const amrex::Real rdOcp, const amrex::Real qv=0.0) { - amrex::Real p_loc = rho * R_d * T; + amrex::Real p_loc = rho * R_d * T * (1.0 + R_v/R_d*qv); return T * std::pow((p_0/p_loc),rdOcp); } @@ -65,9 +66,9 @@ amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv = 0.) * @params[in] rd0cp ratio of R_d to c_p */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -amrex::Real getRhogivenThetaPress (const amrex::Real theta, const amrex::Real p, const amrex::Real rdOcp) +amrex::Real getRhogivenThetaPress (const amrex::Real theta, const amrex::Real p, const amrex::Real rdOcp, const amrex::Real qv=0.0) { - return std::pow(p_0, rdOcp) * std::pow(p, iGamma) / (R_d * theta); + return std::pow(p_0, rdOcp) * std::pow(p, iGamma) / (R_d * theta * (1.0 + R_v/R_d*qv) ); } /** From 4a20716be52d3676381e54a635c71384f1bf7a37 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Tue, 10 Oct 2023 13:29:23 -0700 Subject: [PATCH 02/37] Change all functions in EOS.H to work with moisture --- Source/EOS.H | 26 +++++++++++++++++++------- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/Source/EOS.H b/Source/EOS.H index 189303e7e..886674bba 100644 --- a/Source/EOS.H +++ b/Source/EOS.H @@ -15,8 +15,11 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta, const amrex::Real qv=0.0) { - // rho and rhotheta are dry values. We should be using moist values with moisture + // rho and rhotheta are dry values. We should be using moist value of theta when using moisture + // theta_m = theta * (1 + R_v/R_d*qv) amrex::Real p_loc = p_0 * std::pow(R_d * rhotheta * (1.0 + R_v/R_d*qv) * ip_0, Gamma); + // p = rho_d * R_d * T_v (not T) + // T_v = T * (1 + R_v/R_d*qv) return p_loc / (R_d * rho * (1.0 + R_v/R_d*qv) ); } @@ -30,7 +33,9 @@ amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta, AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getThgivenRandT(const amrex::Real rho, const amrex::Real T, const amrex::Real rdOcp, const amrex::Real qv=0.0) { + // p = rho_d * R_d * T_moist amrex::Real p_loc = rho * R_d * T * (1.0 + R_v/R_d*qv); + // theta_d = T * (p0/p)^(R_d/C_p) return T * std::pow((p_0/p_loc),rdOcp); } @@ -68,6 +73,8 @@ amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv = 0.) AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getRhogivenThetaPress (const amrex::Real theta, const amrex::Real p, const amrex::Real rdOcp, const amrex::Real qv=0.0) { + // We should be using moist value of theta when using moisture + // theta_m = theta * (1 + R_v/R_d*qv) return std::pow(p_0, rdOcp) * std::pow(p, iGamma) / (R_d * theta * (1.0 + R_v/R_d*qv) ); } @@ -79,9 +86,11 @@ amrex::Real getRhogivenThetaPress (const amrex::Real theta, const amrex::Real p, * @params[in] rd0cp ratio of R_d to c_p */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -amrex::Real getdPdRgivenConstantTheta(const amrex::Real rho, const amrex::Real theta) +amrex::Real getdPdRgivenConstantTheta(const amrex::Real rho, const amrex::Real theta, const amrex::Real qv=0.0) { - return Gamma * p_0 * std::pow( (R_d * theta * ip_0), Gamma) * std::pow(rho, Gamma-1.0) ; + // We should be using moist value of theta when using moisture + // theta_m = theta * (1 + R_v/R_d*qv) + return Gamma * p_0 * std::pow( (R_d * theta * (1.0 + R_v/R_d*qv) * ip_0), Gamma) * std::pow(rho, Gamma-1.0) ; } /** @@ -103,10 +112,12 @@ amrex::Real getExnergivenP(const amrex::Real P, const amrex::Real rdOcp) * @params[in] rd0cp ratio of R_d to c_p */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -amrex::Real getExnergivenRTh(const amrex::Real rhotheta, const amrex::Real rdOcp) +amrex::Real getExnergivenRTh(const amrex::Real rhotheta, const amrex::Real rdOcp, const amrex::Real qv=0.0 ) { // Exner function pi in terms of (rho theta) - return std::pow(R_d * rhotheta * ip_0, Gamma * rdOcp); + // We should be using moist value of theta when using moisture + // theta_m = theta * (1 + R_v/R_d*qv) + return std::pow(R_d * rhotheta * (1.0 + R_v/R_d*qv) * ip_0, Gamma * rdOcp); } /** @@ -115,11 +126,12 @@ amrex::Real getExnergivenRTh(const amrex::Real rhotheta, const amrex::Real rdOcp * @params[in] p pressure */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -amrex::Real getRhoThetagivenP(const amrex::Real p) +amrex::Real getRhoThetagivenP(const amrex::Real p, const amrex::Real qv=0.0) { // diagnostic relation for the full pressure // see https://erf.readthedocs.io/en/latest/theory/NavierStokesEquations.html - return std::pow(p*std::pow(p_0, Gamma-1), iGamma) * iR_d; + // For cases with mositure theta = theta_m / (1 + R_v/R_d*qv) + return std::pow(p*std::pow(p_0, Gamma-1), iGamma) * iR_d / (1.0 + R_v/R_d*qv) ; } #endif From 5e6005609d1c19d87255a107a710ac91ab437a0d Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 15 Nov 2023 15:35:29 -0800 Subject: [PATCH 03/37] Adding qp for buoyancy type 1 --- Source/TimeIntegration/ERF_make_buoyancy.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/TimeIntegration/ERF_make_buoyancy.cpp b/Source/TimeIntegration/ERF_make_buoyancy.cpp index 7ad36127c..b3c461a14 100644 --- a/Source/TimeIntegration/ERF_make_buoyancy.cpp +++ b/Source/TimeIntegration/ERF_make_buoyancy.cpp @@ -173,9 +173,9 @@ void make_buoyancy (Vector& S_data, amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Real rhop_hi = cell_data(i,j,k ,Rho_comp) * (1.0 + qv_data(i,j,k ) + qc_data(i,j,k ) - + qi_data(i,j,k )) - r0_arr(i,j,k ); + + qi_data(i,j,k ) + cell_data(i,j,k,RhoQp_comp)/cell_data(i,j,k,Rho_comp)) - r0_arr(i,j,k ); Real rhop_lo = cell_data(i,j,k-1,Rho_comp) * (1.0 + qv_data(i,j,k-1) + qc_data(i,j,k-1) - + qi_data(i,j,k-1)) - r0_arr(i,j,k-1); + + qi_data(i,j,k-1) + + cell_data(i,j,k-1,RhoQp_comp)/cell_data(i,j,k-1,Rho_comp)) - r0_arr(i,j,k-1); buoyancy_fab(i, j, k) = grav_gpu[2] * 0.5 * ( rhop_hi + rhop_lo ); }); } // mfi From 5582efc27d497b6f62dcc42d2627c45c87b15b81 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 15 Nov 2023 15:48:46 -0800 Subject: [PATCH 04/37] Updating docs for buoyancy term --- Docs/sphinx_doc/theory/Buoyancy.rst | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index 9e5f0d819..36cc012f8 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -7,6 +7,28 @@ .. _Buoyancy: +Density of the mixture +======================== + +The total density in a given cell is given by +.. math:: + \rho &=& \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, + +where `m_a` is the mass of dry air, `m_v` is the mass of water vapor, `m_c` is the mass of liquid water, and `m_p` is the mass of precipitate. +From the definitions of the mass mixing ratio (ratio of mass of a component to mass of dry air), we have for any component + +.. math:: + q_i \equiv = \frac{m_i}{m_a}. + +Using this we can write + +.. math:: + \rho = m_a\frac{(1 + q_v + q_c + q_p)}{V} + = \rho_d(1 + q_v + q_c + q_p), + +where `\rho_d \equiv \cfrac{m_a}{V}` is the density of dry air. + + Buoyancy ========= From ee798397b0959701d99f47fbe3f41fc5818ca2b6 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 15 Nov 2023 15:51:04 -0800 Subject: [PATCH 05/37] Updating docs for buoyancy term --- Docs/sphinx_doc/theory/Buoyancy.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index 36cc012f8..2005be0dc 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -11,6 +11,7 @@ Density of the mixture ======================== The total density in a given cell is given by + .. math:: \rho &=& \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, From 38c286126925d1cfbd4697e9a18766735d16e0bb Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 15 Nov 2023 15:55:29 -0800 Subject: [PATCH 06/37] Update docs for buoyancy --- Docs/sphinx_doc/theory/Buoyancy.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index 2005be0dc..1eb960c17 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -15,7 +15,7 @@ The total density in a given cell is given by .. math:: \rho &=& \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, -where `m_a` is the mass of dry air, `m_v` is the mass of water vapor, `m_c` is the mass of liquid water, and `m_p` is the mass of precipitate. +where :math:`m_a` is the mass of dry air, :math:`m_v` is the mass of water vapor, :math:`m_c` is the mass of liquid water, and :math:`m_p` is the mass of precipitate. From the definitions of the mass mixing ratio (ratio of mass of a component to mass of dry air), we have for any component .. math:: @@ -27,7 +27,7 @@ Using this we can write \rho = m_a\frac{(1 + q_v + q_c + q_p)}{V} = \rho_d(1 + q_v + q_c + q_p), -where `\rho_d \equiv \cfrac{m_a}{V}` is the density of dry air. +where :math:`\rho_d \equiv \cfrac{m_a}{V}` is the density of dry air. Buoyancy From fb48acc89f673a07be777772eeb6484fc97e217d Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 15 Nov 2023 16:09:37 -0800 Subject: [PATCH 07/37] Avoiding extra divide --- Source/TimeIntegration/ERF_make_buoyancy.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/TimeIntegration/ERF_make_buoyancy.cpp b/Source/TimeIntegration/ERF_make_buoyancy.cpp index b3c461a14..6e3b7f39e 100644 --- a/Source/TimeIntegration/ERF_make_buoyancy.cpp +++ b/Source/TimeIntegration/ERF_make_buoyancy.cpp @@ -173,9 +173,9 @@ void make_buoyancy (Vector& S_data, amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Real rhop_hi = cell_data(i,j,k ,Rho_comp) * (1.0 + qv_data(i,j,k ) + qc_data(i,j,k ) - + qi_data(i,j,k ) + cell_data(i,j,k,RhoQp_comp)/cell_data(i,j,k,Rho_comp)) - r0_arr(i,j,k ); + + qi_data(i,j,k )) + cell_data(i,j,k,RhoQp_comp) - r0_arr(i,j,k ); Real rhop_lo = cell_data(i,j,k-1,Rho_comp) * (1.0 + qv_data(i,j,k-1) + qc_data(i,j,k-1) - + qi_data(i,j,k-1) + + cell_data(i,j,k-1,RhoQp_comp)/cell_data(i,j,k-1,Rho_comp)) - r0_arr(i,j,k-1); + + qi_data(i,j,k-1)) + cell_data(i,j,k-1,RhoQp_comp) - r0_arr(i,j,k-1); buoyancy_fab(i, j, k) = grav_gpu[2] * 0.5 * ( rhop_hi + rhop_lo ); }); } // mfi From 926ad1ef197ebe9ee702dda3ee586f36b1e4b6b7 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 15 Nov 2023 16:34:23 -0800 Subject: [PATCH 08/37] Update buoyancy docs --- Docs/sphinx_doc/theory/Buoyancy.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index 1eb960c17..204d262e3 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -13,7 +13,7 @@ Density of the mixture The total density in a given cell is given by .. math:: - \rho &=& \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, + \rho = \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, where :math:`m_a` is the mass of dry air, :math:`m_v` is the mass of water vapor, :math:`m_c` is the mass of liquid water, and :math:`m_p` is the mass of precipitate. From the definitions of the mass mixing ratio (ratio of mass of a component to mass of dry air), we have for any component From fe0b6d67c114498e0881550959ff31fa08633972 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 16 Nov 2023 10:26:09 -0800 Subject: [PATCH 09/37] Some corrections to buoyancy docs --- Docs/sphinx_doc/theory/Buoyancy.rst | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index 204d262e3..87df660f0 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -7,19 +7,25 @@ .. _Buoyancy: +Buoyancy +========= + +ERF has three options for how to define the buoyancy force. Even in the absence of moisture these +expressions are not equivalent. + Density of the mixture -======================== +----------------------- -The total density in a given cell is given by +The total density in a cell containing air, water vapor, liquid water and precipitates is given by .. math:: - \rho = \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, + \rho &=& \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, where :math:`m_a` is the mass of dry air, :math:`m_v` is the mass of water vapor, :math:`m_c` is the mass of liquid water, and :math:`m_p` is the mass of precipitate. From the definitions of the mass mixing ratio (ratio of mass of a component to mass of dry air), we have for any component .. math:: - q_i \equiv = \frac{m_i}{m_a}. + q_i \equiv \frac{m_i}{m_a}. Using this we can write @@ -30,12 +36,6 @@ Using this we can write where :math:`\rho_d \equiv \cfrac{m_a}{V}` is the density of dry air. -Buoyancy -========= - -ERF has three options for how to define the buoyancy force. Even in the absence of moisture these -expressions are not equivalent. - Type 1 ------ @@ -47,8 +47,10 @@ One version of the buoyancy force is expressed simply as .. math:: \rho^\prime = \rho_{total} - \rho_0 -where the full density :math:`\rho_{total}` is the sum of dry and moist components and :math:`\rho_0` is the base state density -for dry air only. +where the total density :math:`\rho_{total} = \rho_d(1 + q_v + q_c + q_p)` is the sum of dry and moist components and :math:`\rho_0` is the total density +for the background state. For eg., a usual scenario is that of a background state that contains only air and vapor and no cloud water or precipitates. For such a state, +the total background density :math:`\rho_0 = \rho_{d_0}(1 + q_{v_0})`, where :math:`\rho_{d_0}` and :math:`q_{v_0}` are the background dry density and vapor mixing ratio respectively. +As a check, we observe that :math:`\rho^\prime_0 = 0`, which means that the background state is not buoyant. Type 2 ------ From a89a48f0162e973d50ca4ed2f0a1b330f1404c80 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 16 Nov 2023 10:33:44 -0800 Subject: [PATCH 10/37] Some corrections to Docs in buoyancy --- Docs/sphinx_doc/theory/Buoyancy.rst | 23 ----------------------- 1 file changed, 23 deletions(-) diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index de71f28b9..87df660f0 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -7,29 +7,6 @@ .. _Buoyancy: -Density of the mixture -======================== - -The total density in a given cell is given by - -.. math:: - \rho = \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, - -where :math:`m_a` is the mass of dry air, :math:`m_v` is the mass of water vapor, :math:`m_c` is the mass of liquid water, and :math:`m_p` is the mass of precipitate. -From the definitions of the mass mixing ratio (ratio of mass of a component to mass of dry air), we have for any component - -.. math:: - q_i \equiv = \frac{m_i}{m_a}. - -Using this we can write - -.. math:: - \rho = m_a\frac{(1 + q_v + q_c + q_p)}{V} - = \rho_d(1 + q_v + q_c + q_p), - -where :math:`\rho_d \equiv \cfrac{m_a}{V}` is the density of dry air. - - Buoyancy ========= From ab100b4f3fd1bc30bf048f95b0ba2431f740821b Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 16 Nov 2023 12:20:04 -0800 Subject: [PATCH 11/37] Correcting docs issue --- Docs/sphinx_doc/theory/Buoyancy.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index 87df660f0..6213c379f 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -19,7 +19,7 @@ Density of the mixture The total density in a cell containing air, water vapor, liquid water and precipitates is given by .. math:: - \rho &=& \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, + \rho = \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, where :math:`m_a` is the mass of dry air, :math:`m_v` is the mass of water vapor, :math:`m_c` is the mass of liquid water, and :math:`m_p` is the mass of precipitate. From the definitions of the mass mixing ratio (ratio of mass of a component to mass of dry air), we have for any component From d277d3390729f47103a51baa405b22d3c41874ea Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 16 Nov 2023 12:22:43 -0800 Subject: [PATCH 12/37] Correcting docs --- Docs/sphinx_doc/theory/Buoyancy.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index 87df660f0..6213c379f 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -19,7 +19,7 @@ Density of the mixture The total density in a cell containing air, water vapor, liquid water and precipitates is given by .. math:: - \rho &=& \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, + \rho = \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V}, where :math:`m_a` is the mass of dry air, :math:`m_v` is the mass of water vapor, :math:`m_c` is the mass of liquid water, and :math:`m_p` is the mass of precipitate. From the definitions of the mass mixing ratio (ratio of mass of a component to mass of dry air), we have for any component From 051117a4c3df86f352df0d709c43be1c0c48222d Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 16 Nov 2023 15:07:49 -0800 Subject: [PATCH 13/37] adding Kessler microphysics docs --- Docs/sphinx_doc/theory/Microphysics.rst | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/Docs/sphinx_doc/theory/Microphysics.rst b/Docs/sphinx_doc/theory/Microphysics.rst index 42684f0cf..eadc410e8 100644 --- a/Docs/sphinx_doc/theory/Microphysics.rst +++ b/Docs/sphinx_doc/theory/Microphysics.rst @@ -7,6 +7,21 @@ .. _Microphysics: +Kessler Microphysics model +=========================== +Governing equations for the microphysical quantities for Kessler microphysics from `gabervsek2012dry`_ are + +.. math:: + \frac{\partial q_v}{\partial t} = -C_c + E_c + E_r + \frac{\partial q_c}{\partial t} = C_c - E_c - (A_c + K_c) + \frac{\partial q_p}{\partial t} = \frac{1}{\overline{\rho}}\frac{\partial}{\partial z}(\overline{\rho}Vq_p) + (A_c + K_c) - E_r + \frac{\partial q_t}{\partial t} = \frac{\partial q_v}{\partial t} + \frac{\partial q_c}{\partial t} + = E_r - (A_c + K_c) +where :math:`C_c` is the rate of condensation of water vapor to cloud water, :math:`E_c` is the rate of evaporation of cloud water to water vapor, :math:`A_c` is the autoconversion of cloud water to rain, :math:`K_c` is the accretion of cloud water to rain drops, :math:`E_r` is the evaporation of rain to water vapor and :math:`F_r` is the sedimentation of rain. The parametrization used is given in \cite{klemp1978simulation}, and is given below. Note that in all the equations, $p$ is specified in millibars and $\overline{\rho}$ is specified in g cm$^{-3}$. + +.. _`gabervsek2012dry`: https://journals.ametsoc.org/view/journals/mwre/140/10/mwr-d-11-00144.1.xml + + Single Moment Microphysics Model =================================== The conversion rates among the moist hydrometeors are parameterized assuming that From 16e76d76d1b0e9c591489b5b96a921a2b84950ab Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 16 Nov 2023 15:10:32 -0800 Subject: [PATCH 14/37] adding Kessler microphysics docs --- Docs/sphinx_doc/theory/Microphysics.rst | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/Docs/sphinx_doc/theory/Microphysics.rst b/Docs/sphinx_doc/theory/Microphysics.rst index eadc410e8..d18b21d7f 100644 --- a/Docs/sphinx_doc/theory/Microphysics.rst +++ b/Docs/sphinx_doc/theory/Microphysics.rst @@ -13,10 +13,14 @@ Governing equations for the microphysical quantities for Kessler microphysics fr .. math:: \frac{\partial q_v}{\partial t} = -C_c + E_c + E_r +.. math:: \frac{\partial q_c}{\partial t} = C_c - E_c - (A_c + K_c) +.. math:: \frac{\partial q_p}{\partial t} = \frac{1}{\overline{\rho}}\frac{\partial}{\partial z}(\overline{\rho}Vq_p) + (A_c + K_c) - E_r +.. math:: \frac{\partial q_t}{\partial t} = \frac{\partial q_v}{\partial t} + \frac{\partial q_c}{\partial t} = E_r - (A_c + K_c) + where :math:`C_c` is the rate of condensation of water vapor to cloud water, :math:`E_c` is the rate of evaporation of cloud water to water vapor, :math:`A_c` is the autoconversion of cloud water to rain, :math:`K_c` is the accretion of cloud water to rain drops, :math:`E_r` is the evaporation of rain to water vapor and :math:`F_r` is the sedimentation of rain. The parametrization used is given in \cite{klemp1978simulation}, and is given below. Note that in all the equations, $p$ is specified in millibars and $\overline{\rho}$ is specified in g cm$^{-3}$. .. _`gabervsek2012dry`: https://journals.ametsoc.org/view/journals/mwre/140/10/mwr-d-11-00144.1.xml From f33cdd8060d3bb8e5ed26dc46c554ff8dd6d3baa Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 16 Nov 2023 15:14:35 -0800 Subject: [PATCH 15/37] adding Kessler microphysics docs --- Docs/sphinx_doc/theory/Microphysics.rst | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/Docs/sphinx_doc/theory/Microphysics.rst b/Docs/sphinx_doc/theory/Microphysics.rst index d18b21d7f..f747ace71 100644 --- a/Docs/sphinx_doc/theory/Microphysics.rst +++ b/Docs/sphinx_doc/theory/Microphysics.rst @@ -21,9 +21,13 @@ Governing equations for the microphysical quantities for Kessler microphysics fr \frac{\partial q_t}{\partial t} = \frac{\partial q_v}{\partial t} + \frac{\partial q_c}{\partial t} = E_r - (A_c + K_c) -where :math:`C_c` is the rate of condensation of water vapor to cloud water, :math:`E_c` is the rate of evaporation of cloud water to water vapor, :math:`A_c` is the autoconversion of cloud water to rain, :math:`K_c` is the accretion of cloud water to rain drops, :math:`E_r` is the evaporation of rain to water vapor and :math:`F_r` is the sedimentation of rain. The parametrization used is given in \cite{klemp1978simulation}, and is given below. Note that in all the equations, $p$ is specified in millibars and $\overline{\rho}$ is specified in g cm$^{-3}$. +where :math:`C_c` is the rate of condensation of water vapor to cloud water, :math:`E_c` is the rate of evaporation of cloud water to water vapor, +:math:`A_c` is the autoconversion of cloud water to rain, :math:`K_c` is the accretion of cloud water to rain drops, :math:`E_r` is the evaporation of +rain to water vapor and :math:`F_r` is the sedimentation of rain. The parametrization used is given in `klemp1978simulation`_, and is given +below. Note that in all the equations, :math:`p` is specified in millibars and :math:`\overline{\rho}` is specified in g cm:math:`^{-3}`. .. _`gabervsek2012dry`: https://journals.ametsoc.org/view/journals/mwre/140/10/mwr-d-11-00144.1.xml +.. _`klemp1978simulation`: https://journals.ametsoc.org/view/journals/atsc/35/6/1520-0469_1978_035_1070_tsotdc_2_0_co_2.xml Single Moment Microphysics Model From a0bb74dddf0c5de225692a763a2bb01bcfaf0aec Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 16 Nov 2023 15:16:11 -0800 Subject: [PATCH 16/37] adding Kessler microphysics docs --- Docs/sphinx_doc/theory/Microphysics.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Docs/sphinx_doc/theory/Microphysics.rst b/Docs/sphinx_doc/theory/Microphysics.rst index f747ace71..a644d86d4 100644 --- a/Docs/sphinx_doc/theory/Microphysics.rst +++ b/Docs/sphinx_doc/theory/Microphysics.rst @@ -24,7 +24,7 @@ Governing equations for the microphysical quantities for Kessler microphysics fr where :math:`C_c` is the rate of condensation of water vapor to cloud water, :math:`E_c` is the rate of evaporation of cloud water to water vapor, :math:`A_c` is the autoconversion of cloud water to rain, :math:`K_c` is the accretion of cloud water to rain drops, :math:`E_r` is the evaporation of rain to water vapor and :math:`F_r` is the sedimentation of rain. The parametrization used is given in `klemp1978simulation`_, and is given -below. Note that in all the equations, :math:`p` is specified in millibars and :math:`\overline{\rho}` is specified in g cm:math:`^{-3}`. +below. Note that in all the equations, :math:`p` is specified in millibars and :math:`\overline{\rho}` is specified in g cm :math:`^{-3}`. .. _`gabervsek2012dry`: https://journals.ametsoc.org/view/journals/mwre/140/10/mwr-d-11-00144.1.xml .. _`klemp1978simulation`: https://journals.ametsoc.org/view/journals/atsc/35/6/1520-0469_1978_035_1070_tsotdc_2_0_co_2.xml From 84bd4efb9ef83c36ba5103b3a1613cd19dcc8f62 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 29 Nov 2023 05:00:18 -0800 Subject: [PATCH 17/37] Removing unused vars --- Exec/SquallLine_2D/GNUmakefile | 2 +- Exec/SquallLine_2D/inputs_moisture | 2 +- Exec/SquallLine_2D/prob.H | 11 +-- Exec/SquallLine_2D/prob.cpp | 92 +++++++------------ Source/Initialization/ERF_init1d.cpp | 2 +- .../Microphysics/Kessler/Diagnose_Kessler.cpp | 1 - Source/Microphysics/Kessler/Kessler.cpp | 26 +----- .../Microphysics/Kessler/Update_Kessler.cpp | 2 - Source/prob_common.H | 1 - 9 files changed, 38 insertions(+), 101 deletions(-) diff --git a/Exec/SquallLine_2D/GNUmakefile b/Exec/SquallLine_2D/GNUmakefile index 40f807cd1..96d247391 100644 --- a/Exec/SquallLine_2D/GNUmakefile +++ b/Exec/SquallLine_2D/GNUmakefile @@ -19,7 +19,7 @@ USE_HIP = FALSE USE_SYCL = FALSE # Debugging -DEBUG = FALSE +DEBUG = TRUE TEST = TRUE USE_ASSERTION = TRUE diff --git a/Exec/SquallLine_2D/inputs_moisture b/Exec/SquallLine_2D/inputs_moisture index 1a6d7fb48..05eb31e60 100644 --- a/Exec/SquallLine_2D/inputs_moisture +++ b/Exec/SquallLine_2D/inputs_moisture @@ -1,5 +1,5 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -max_step = 9000 +max_step = 1 stop_time = 90000.0 amrex.fpe_trap_invalid = 1 diff --git a/Exec/SquallLine_2D/prob.H b/Exec/SquallLine_2D/prob.H index 1d5891061..be547c6e9 100644 --- a/Exec/SquallLine_2D/prob.H +++ b/Exec/SquallLine_2D/prob.H @@ -13,16 +13,8 @@ #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; + amrex::Real Theta_0 = 300.0; }; // namespace ProbParm @@ -63,7 +55,6 @@ public: void erf_init_dens_hse_moist (amrex::MultiFab& rho_hse, std::unique_ptr& z_phys_nd, - std::unique_ptr& z_phys_cc, amrex::Geometry const& geom) override; void init_custom_terrain ( diff --git a/Exec/SquallLine_2D/prob.cpp b/Exec/SquallLine_2D/prob.cpp index 8b7f2cc65..94d218fa7 100644 --- a/Exec/SquallLine_2D/prob.cpp +++ b/Exec/SquallLine_2D/prob.cpp @@ -14,23 +14,11 @@ Problem::Problem () { // Parse params amrex::ParmParse pp("prob"); - pp.query("T_0", parms.T_0); - pp.query("U_0", parms.U_0); - pp.query("x_c", parms.x_c); - pp.query("z_c", parms.z_c); - pp.query("x_r", parms.x_r); - pp.query("z_r", parms.z_r); - pp.query("T_pert", parms.T_pert); } - -Real z_tr = 12000.0; -Real height = 1200.0; -Real R_v_by_R_d = R_v/R_d; - Real compute_theta (Real z) { - Real theta_0=300.0, theta_tr=343.0, T_tr=213.0; + Real theta_0=300.0, theta_tr=343.0, T_tr=213.0, z_tr = 12000.0; if(z <= z_tr) { return theta_0 + (theta_tr - theta_0)*std::pow(z/z_tr,1.25); } else { @@ -46,7 +34,7 @@ Real compute_saturation_pressure (const Real T_b) Real p_s = exp(34.494 - 4924.99/(T_b - 273.15 + 237.1))/std::pow(T_b - 273.15 + 105.0,1.57); - Real T = T_b - 273.15; + //Real T = T_b - 273.15; //Real p_s = 0.61121e3*exp((18.678 - T/234.5)*(T/(257.14 + T))); @@ -55,8 +43,11 @@ Real compute_saturation_pressure (const Real T_b) AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE -Real compute_relative_humidity (Real z, Real height, Real z_tr, Real p_b, Real T_b) +Real compute_relative_humidity (Real z, Real p_b, Real T_b) { + Real height = 1200.0; + Real z_tr = 12000.0; + Real p_s = compute_saturation_pressure(T_b); Real q_s = 0.622*p_s/(p_b - p_s); @@ -79,8 +70,9 @@ Real compute_vapor_pressure (Real p_s, Real RH) AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE -Real vapor_mixing_ratio (const Real z, const Real height, const Real p_b, const Real T_b, const Real RH){ - +Real vapor_mixing_ratio (const Real z, const Real p_b, const Real T_b, const Real RH) +{ + Real height = 1200.0; Real p_s = compute_saturation_pressure(T_b); Real p_v = compute_vapor_pressure(p_s, RH); @@ -100,13 +92,13 @@ Real compute_temperature (const Real p_b, const Real theta_b) return theta_b*std::pow(p_b/p_0,R_d/Cp_d); } -Real compute_dewpoint_temperature (const Real z, const Real p_b, const Real T_b, const Real RH) +Real compute_dewpoint_temperature (const Real T_b, const Real RH) { Real T_dp, gamma, T; T = T_b - 273.15; - Real a = 6.11e-3*p_0, b = 18.678, c = 257.14, d = 234.5; + Real b = 18.678, c = 257.14, d = 234.5; gamma = log(RH*exp((b - T/d)*T/(c + T))); T_dp = c*gamma/(b - gamma); @@ -120,11 +112,11 @@ void compute_rho (const Real z, const Real pressure, Real &theta, Real& rho, Rea theta = compute_theta(z); T_b = compute_temperature(pressure, theta); - Real RH = compute_relative_humidity(z, height, z_tr, pressure, T_b); - q_v = vapor_mixing_ratio(z, height, pressure, T_b, RH); + Real RH = compute_relative_humidity(z, pressure, T_b); + q_v = vapor_mixing_ratio(z, pressure, T_b, RH); rho = pressure/(R_d*T_b*(1.0 + (R_v/R_d)*q_v)); rho = rho*(1.0 + q_v); - T_dp = compute_dewpoint_temperature(z, pressure, T_b, RH); + T_dp = compute_dewpoint_temperature(T_b, RH); } double compute_F (const Real p_k, const Real p_k_minus_1, Real &theta_k, Real& rho_k, Real& q_v_k, @@ -168,20 +160,25 @@ init_isentropic_hse_no_terrain(Real *theta, Real* r, Real* p, Real *q_v, const Real& dz, const Real& prob_lo_z, const int& khi) { - Real z, p_b, theta_b, T_b, rho_b, RH, T_dp, rho0, theta0, T0, q_v0; - p_b = p_0; - //T_b = T0; + FILE *file_IC; + file_IC = fopen("input_sounding_probcpp.txt","w"); + Real z, T_b, T_dp; // Compute the quantities at z = 0.5*dz (first cell center) z = prob_lo_z + 0.5*dz; p[0] = p_0; compute_p_k(p[0], p_0, theta[0], r[0], q_v[0], T_dp, T_b, dz, z, 0.0); + fprintf(file_IC, "%0.15g %0.15g %0.15g %0.15g %0.15g %0.15g %0.15g\n", z, T_b-273.15, T_dp, p[0], r[0], theta[0], q_v[0]); + for (int k=1;k<=khi;k++){ z = prob_lo_z + (k+0.5)*dz; p[k] = p[k-1]; compute_p_k(p[k], p[k-1], theta[k], r[k], q_v[k], T_dp, T_b, dz, z, r[k-1]); + fprintf(file_IC, "%0.15g %0.15g %0.15g %0.15g %0.15g %0.15g %0.15g\n", z, T_b-273.15, T_dp, p[k], r[k], theta[k], q_v[k]); } + fclose(file_IC); + r[khi+1] = r[khi]; } @@ -189,7 +186,6 @@ init_isentropic_hse_no_terrain(Real *theta, Real* r, Real* p, Real *q_v, void Problem::erf_init_dens_hse_moist (MultiFab& rho_hse, std::unique_ptr& z_phys_nd, - std::unique_ptr& z_phys_cc, Geometry const& geom) { @@ -197,11 +193,6 @@ Problem::erf_init_dens_hse_moist (MultiFab& rho_hse, const Real dz = geom.CellSize()[2]; const int khi = geom.Domain().bigEnd()[2]; - - const Real T_sfc = 300.; - const Real rho_sfc = p_0 / (R_d*T_sfc); - const Real Thetabar = T_sfc; - // use_terrain = 1 if (z_phys_nd) { @@ -210,7 +201,7 @@ Problem::erf_init_dens_hse_moist (MultiFab& rho_hse, for ( MFIter mfi(rho_hse, TileNoZ()); mfi.isValid(); ++mfi ) { Array4 rho_arr = rho_hse.array(mfi); - Array4 z_cc_arr = z_phys_cc->const_array(mfi); + //Array4 z_cc_arr = z_phys_cc->const_array(mfi); // Create a flat box with same horizontal extent but only one cell in vertical const Box& tbz = mfi.nodaltilebox(2); @@ -219,8 +210,7 @@ Problem::erf_init_dens_hse_moist (MultiFab& rho_hse, b2d.setRange(2,0); ParallelFor(b2d, [=] AMREX_GPU_DEVICE (int i, int j, int) { - Array1D r;; - Array1D p;; + Array1D r; //init_isentropic_hse_terrain(i,j,rho_sfc,Thetabar,&(r(0)),&(p(0)),z_cc_arr,khi); @@ -252,8 +242,6 @@ Problem::erf_init_dens_hse_moist (MultiFab& rho_hse, amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_q_v.begin(), h_q_v.end(), d_q_v.begin()); Real* r = d_r.data(); - Real* q_v = d_q_v.data(); - Real* theta = d_t.data(); #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) @@ -306,19 +294,6 @@ Problem::init_custom_pert ( // This is what we do at k = 0 -- note we assume p = p_0 and T = T_0 at z=0 const amrex::Real& dz = geomdata.CellSize()[2]; const amrex::Real& prob_lo_z = geomdata.ProbLo()[2]; - const amrex::Real& prob_hi_z = geomdata.ProbHi()[2]; - - const amrex::Real rdOcp = sc.rdOcp; - - // const amrex::Real thetatr = 343.0; - // const amrex::Real theta0 = 300.0; - const amrex::Real ztr = 12000.0; - const amrex::Real Ttr = 213.0; - const amrex::Real Ttop = 213.0; - const amrex::Real deltaz = 1000.*0.0; - const amrex::Real zs = 5000.; - const amrex::Real us = 30.; - const amrex::Real uc = 15.; // Call the routine to calculate the 1d background condition @@ -340,17 +315,16 @@ Problem::init_custom_pert ( amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_q_v.begin(), h_q_v.end(), d_q_v.begin()); Real* t = d_t.data(); - Real* r = d_r.data(); - Real* q_v = d_q_v.data(); Real* p = d_p.data(); - Real d_z_tr = 12000.0; - Real d_height = 1200.0; const Real x_c = 0.0, z_c = 1.5e3, x_r = 4.0e3, z_r = 1.5e3, r_c = 1.0, theta_c = 3.0; //const Real x_c = 0.0, z_c = 2.0e3, x_r = 10.0e3, z_r = 1.5e3, r_c = 1.0, theta_c = 3.0; - amrex::ParallelForRNG(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k, const amrex::RandomEngine& engine) noexcept + Real Rd_by_Cp = sc.rdOcp; + Rd_by_Cp = Rd_by_Cp; + + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Geometry (note we must include these here to get the data on device) const auto prob_lo = geomdata.ProbLo(); @@ -375,15 +349,15 @@ Problem::init_custom_pert ( theta_total = t[k] + delta_theta; temperature = compute_temperature(p[k], theta_total); Real T_b = compute_temperature(p[k], t[k]); - RH = compute_relative_humidity(z, d_height, d_z_tr, p[k], T_b); - Real q_v_hot = vapor_mixing_ratio(z, d_height, p[k], T_b, RH); + RH = compute_relative_humidity(z, p[k], T_b); + Real q_v_hot = vapor_mixing_ratio(z, p[k], T_b, RH); rho = p[k]/(R_d*temperature*(1.0 + (R_v/R_d)*q_v_hot)); // Compute background quantities Real temperature_back = compute_temperature(p[k], t[k]); Real T_back = compute_temperature(p[k], t[k]); - Real RH_back = compute_relative_humidity(z, d_height, d_z_tr, p[k], T_back); - Real q_v_back = vapor_mixing_ratio(z, d_height, p[k], T_back, RH_back); + Real RH_back = compute_relative_humidity(z, p[k], T_back); + Real q_v_back = vapor_mixing_ratio(z, p[k], T_back, RH_back); Real rho_back = p[k]/(R_d*temperature_back*(1.0 + (R_v/R_d)*q_v_back)); // This version perturbs rho but not p @@ -470,7 +444,7 @@ Problem::erf_init_rayleigh (amrex::Vector& tau, ubar[k] = 2.0; vbar[k] = 1.0; wbar[k] = 0.0; - //thetabar[k] = parms.Theta_0; + thetabar[k] = parms.Theta_0; } } diff --git a/Source/Initialization/ERF_init1d.cpp b/Source/Initialization/ERF_init1d.cpp index ce47b5f80..2046e3954 100644 --- a/Source/Initialization/ERF_init1d.cpp +++ b/Source/Initialization/ERF_init1d.cpp @@ -132,7 +132,7 @@ ERF::initHSE (int lev) #ifdef ERF_USE_MOISTURE // Initial r_hse may or may not be in HSE -- defined in prob.cpp if(solverChoice.use_moist_background){ - prob->erf_init_dens_hse_moist(r_hse, z_phys_nd[lev], z_phys_cc[lev], geom[lev]); + prob->erf_init_dens_hse_moist(r_hse, z_phys_nd[lev], geom[lev]); }else{ prob->erf_init_dens_hse(r_hse, z_phys_nd[lev], z_phys_cc[lev], geom[lev]); } diff --git a/Source/Microphysics/Kessler/Diagnose_Kessler.cpp b/Source/Microphysics/Kessler/Diagnose_Kessler.cpp index 4b74bcf28..43210f501 100644 --- a/Source/Microphysics/Kessler/Diagnose_Kessler.cpp +++ b/Source/Microphysics/Kessler/Diagnose_Kessler.cpp @@ -25,7 +25,6 @@ void Kessler::Diagnose () { auto qci_array = qci->array(mfi); auto qpl_array = qpl->array(mfi); auto qpi_array = qpi->array(mfi); - auto tabs_array = tabs->array(mfi); const auto& box3d = mfi.tilebox(); diff --git a/Source/Microphysics/Kessler/Kessler.cpp b/Source/Microphysics/Kessler/Kessler.cpp index d12cd566f..529bb03e9 100644 --- a/Source/Microphysics/Kessler/Kessler.cpp +++ b/Source/Microphysics/Kessler/Kessler.cpp @@ -13,27 +13,6 @@ using namespace amrex; */ void Kessler::AdvanceKessler() { - Real powr1 = (3.0 + b_rain) / 4.0; - Real powr2 = (5.0 + b_rain) / 8.0; - Real pows1 = (3.0 + b_snow) / 4.0; - Real pows2 = (5.0 + b_snow) / 8.0; - Real powg1 = (3.0 + b_grau) / 4.0; - Real powg2 = (5.0 + b_grau) / 8.0; - - auto accrrc_t = accrrc.table(); - auto accrsc_t = accrsc.table(); - auto accrsi_t = accrsi.table(); - auto accrgc_t = accrgc.table(); - auto accrgi_t = accrgi.table(); - auto coefice_t = coefice.table(); - auto evapr1_t = evapr1.table(); - auto evapr2_t = evapr2.table(); - auto evaps1_t = evaps1.table(); - auto evaps2_t = evaps2.table(); - auto evapg1_t = evapg1.table(); - auto evapg2_t = evapg2.table(); - auto pres1d_t = pres1d.table(); - auto qt = mic_fab_vars[MicVar_Kess::qt]; auto qp = mic_fab_vars[MicVar_Kess::qp]; auto qn = mic_fab_vars[MicVar_Kess::qn]; @@ -52,7 +31,6 @@ void Kessler::AdvanceKessler() { MultiFab fz; auto ba = tabs->boxArray(); auto dm = tabs->DistributionMap(); - auto ngrow = tabs->nGrowVect(); fz.define(convert(ba, IntVect(0,0,1)), dm, 1, 0); // No ghost cells for ( MFIter mfi(fz, TilingIfNotGPU()); mfi.isValid(); ++mfi ){ @@ -136,7 +114,6 @@ void Kessler::AdvanceKessler() { auto qt_array = mic_fab_vars[MicVar_Kess::qt]->array(mfi); auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi); - auto qcl_array = mic_fab_vars[MicVar_Kess::qcl]->array(mfi); auto theta_array = theta->array(mfi); auto qv_array = qv->array(mfi); auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); @@ -160,8 +137,7 @@ void Kessler::AdvanceKessler() { } //------- Autoconversion/accretion - Real omn, omp, omg, qcc, qii, autor, autos, accrr, qrr, accrcs, accris, - qss, accrcg, accrig, tmp, qgg, dq_clwater_to_rain, dq_rain_to_vapor, dq_clwater_to_vapor, dq_vapor_to_clwater, qsat; + Real qcc, autor, accrr, dq_clwater_to_rain, dq_rain_to_vapor, dq_clwater_to_vapor, dq_vapor_to_clwater, qsat; Real pressure = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k))/100.0; erf_qsatw(tabs_array(i,j,k), pressure, qsat); diff --git a/Source/Microphysics/Kessler/Update_Kessler.cpp b/Source/Microphysics/Kessler/Update_Kessler.cpp index fcef4df9b..a32d44452 100644 --- a/Source/Microphysics/Kessler/Update_Kessler.cpp +++ b/Source/Microphysics/Kessler/Update_Kessler.cpp @@ -33,8 +33,6 @@ void Kessler::Update (amrex::MultiFab& cons, auto theta_arr = mic_fab_vars[MicVar_Kess::theta]->array(mfi); auto qt_arr = mic_fab_vars[MicVar_Kess::qt]->array(mfi); auto qp_arr = mic_fab_vars[MicVar_Kess::qp]->array(mfi); - auto qpl_arr = mic_fab_vars[MicVar_Kess::qpl]->array(mfi); - auto qpi_arr = mic_fab_vars[MicVar_Kess::qpi]->array(mfi); auto qgraup_arr= qgraup_mf.array(mfi); diff --git a/Source/prob_common.H b/Source/prob_common.H index 931f04479..ca365830b 100644 --- a/Source/prob_common.H +++ b/Source/prob_common.H @@ -52,7 +52,6 @@ public: virtual void erf_init_dens_hse_moist (amrex::MultiFab& /*rho_hse*/, std::unique_ptr& /*z_phys_nd*/, - std::unique_ptr& /*z_phys_cc*/, amrex::Geometry const& /*geom*/) { From 32ce1029dbb53da8498fd846c69a226ce9898cac Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 29 Nov 2023 05:40:04 -0800 Subject: [PATCH 18/37] Updating prob for squall line --- Exec/SquallLine_2D/prob.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/Exec/SquallLine_2D/prob.cpp b/Exec/SquallLine_2D/prob.cpp index 94d218fa7..5ed946ff0 100644 --- a/Exec/SquallLine_2D/prob.cpp +++ b/Exec/SquallLine_2D/prob.cpp @@ -92,6 +92,8 @@ Real compute_temperature (const Real p_b, const Real theta_b) return theta_b*std::pow(p_b/p_0,R_d/Cp_d); } +AMREX_FORCE_INLINE +AMREX_GPU_HOST_DEVICE Real compute_dewpoint_temperature (const Real T_b, const Real RH) { @@ -104,7 +106,6 @@ Real compute_dewpoint_temperature (const Real T_b, const Real RH) T_dp = c*gamma/(b - gamma); return T_dp; - } void compute_rho (const Real z, const Real pressure, Real &theta, Real& rho, Real& q_v, Real& T_dp, Real& T_b) @@ -119,8 +120,8 @@ void compute_rho (const Real z, const Real pressure, Real &theta, Real& rho, Rea T_dp = compute_dewpoint_temperature(T_b, RH); } -double compute_F (const Real p_k, const Real p_k_minus_1, Real &theta_k, Real& rho_k, Real& q_v_k, - Real& T_dp, Real& T_b, const Real dz, const Real z, const Real rho_k_minus_1) +double compute_F (const Real& p_k, const Real& p_k_minus_1, Real &theta_k, Real& rho_k, Real& q_v_k, + Real& T_dp, Real& T_b, const Real& dz, const Real& z, const Real& rho_k_minus_1) { Real F; From cc034f5599c81e1f61136326a55842ffd7aef6f9 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 29 Nov 2023 05:41:58 -0800 Subject: [PATCH 19/37] Updating prob for squall line --- Exec/SquallLine_2D/prob.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/Exec/SquallLine_2D/prob.cpp b/Exec/SquallLine_2D/prob.cpp index 5ed946ff0..f6e77845f 100644 --- a/Exec/SquallLine_2D/prob.cpp +++ b/Exec/SquallLine_2D/prob.cpp @@ -45,8 +45,8 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE Real compute_relative_humidity (Real z, Real p_b, Real T_b) { - Real height = 1200.0; - Real z_tr = 12000.0; + Real height = 1200.0; + Real z_tr = 12000.0; Real p_s = compute_saturation_pressure(T_b); @@ -72,7 +72,7 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE Real vapor_mixing_ratio (const Real z, const Real p_b, const Real T_b, const Real RH) { - Real height = 1200.0; + Real height = 1200.0; Real p_s = compute_saturation_pressure(T_b); Real p_v = compute_vapor_pressure(p_s, RH); @@ -161,7 +161,7 @@ init_isentropic_hse_no_terrain(Real *theta, Real* r, Real* p, Real *q_v, const Real& dz, const Real& prob_lo_z, const int& khi) { - FILE *file_IC; + FILE *file_IC; file_IC = fopen("input_sounding_probcpp.txt","w"); Real z, T_b, T_dp; @@ -169,16 +169,16 @@ init_isentropic_hse_no_terrain(Real *theta, Real* r, Real* p, Real *q_v, z = prob_lo_z + 0.5*dz; p[0] = p_0; compute_p_k(p[0], p_0, theta[0], r[0], q_v[0], T_dp, T_b, dz, z, 0.0); - fprintf(file_IC, "%0.15g %0.15g %0.15g %0.15g %0.15g %0.15g %0.15g\n", z, T_b-273.15, T_dp, p[0], r[0], theta[0], q_v[0]); + fprintf(file_IC, "%0.15g %0.15g %0.15g %0.15g %0.15g %0.15g %0.15g\n", z, T_b-273.15, T_dp, p[0], r[0], theta[0], q_v[0]); for (int k=1;k<=khi;k++){ z = prob_lo_z + (k+0.5)*dz; p[k] = p[k-1]; compute_p_k(p[k], p[k-1], theta[k], r[k], q_v[k], T_dp, T_b, dz, z, r[k-1]); - fprintf(file_IC, "%0.15g %0.15g %0.15g %0.15g %0.15g %0.15g %0.15g\n", z, T_b-273.15, T_dp, p[k], r[k], theta[k], q_v[k]); + fprintf(file_IC, "%0.15g %0.15g %0.15g %0.15g %0.15g %0.15g %0.15g\n", z, T_b-273.15, T_dp, p[k], r[k], theta[k], q_v[k]); } - fclose(file_IC); + fclose(file_IC); r[khi+1] = r[khi]; @@ -322,8 +322,8 @@ Problem::init_custom_pert ( const Real x_c = 0.0, z_c = 1.5e3, x_r = 4.0e3, z_r = 1.5e3, r_c = 1.0, theta_c = 3.0; //const Real x_c = 0.0, z_c = 2.0e3, x_r = 10.0e3, z_r = 1.5e3, r_c = 1.0, theta_c = 3.0; - Real Rd_by_Cp = sc.rdOcp; - Rd_by_Cp = Rd_by_Cp; + Real Rd_by_Cp = sc.rdOcp; + Rd_by_Cp = Rd_by_Cp; amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { From 45e2bae0e61139f3ce3de0fa8c0a75e323666e5d Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 29 Nov 2023 06:09:31 -0800 Subject: [PATCH 20/37] Minor change to makefile --- Exec/SquallLine_2D/GNUmakefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Exec/SquallLine_2D/GNUmakefile b/Exec/SquallLine_2D/GNUmakefile index 96d247391..40f807cd1 100644 --- a/Exec/SquallLine_2D/GNUmakefile +++ b/Exec/SquallLine_2D/GNUmakefile @@ -19,7 +19,7 @@ USE_HIP = FALSE USE_SYCL = FALSE # Debugging -DEBUG = TRUE +DEBUG = FALSE TEST = TRUE USE_ASSERTION = TRUE From 83f594b573d9e8394a71bb1ca2c73759f417f334 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 29 Nov 2023 06:10:01 -0800 Subject: [PATCH 21/37] Minor change --- Exec/SquallLine_2D/inputs_moisture | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Exec/SquallLine_2D/inputs_moisture b/Exec/SquallLine_2D/inputs_moisture index 05eb31e60..1a6d7fb48 100644 --- a/Exec/SquallLine_2D/inputs_moisture +++ b/Exec/SquallLine_2D/inputs_moisture @@ -1,5 +1,5 @@ # ------------------ INPUTS TO MAIN PROGRAM ------------------- -max_step = 1 +max_step = 9000 stop_time = 90000.0 amrex.fpe_trap_invalid = 1 From a7ccf8be314ef43a1576b2d812d51666a605fac9 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 29 Nov 2023 07:50:31 -0800 Subject: [PATCH 22/37] Finalizing Kessler --- Exec/SquallLine_2D/inputs_moisture | 17 +++-- Exec/SquallLine_2D/inputs_moisture_outflow | 86 ---------------------- Exec/SquallLine_2D/prob.H | 61 +++++++++++++++ Exec/SquallLine_2D/prob.cpp | 40 ++++++---- 4 files changed, 99 insertions(+), 105 deletions(-) delete mode 100644 Exec/SquallLine_2D/inputs_moisture_outflow diff --git a/Exec/SquallLine_2D/inputs_moisture b/Exec/SquallLine_2D/inputs_moisture index 1a6d7fb48..2e2757fe5 100644 --- a/Exec/SquallLine_2D/inputs_moisture +++ b/Exec/SquallLine_2D/inputs_moisture @@ -45,7 +45,7 @@ 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 +#amr.restart = chk09000 # PLOTFILES erf.plot_file_1 = plt # root name of plotfile @@ -70,7 +70,7 @@ erf.molec_diff_type = "ConstantAlpha" 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.alpha_C = 25.0 erf.moisture_model = "Kessler" erf.use_moist_background = true @@ -79,6 +79,13 @@ 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 +prob.z_tr = 12000.0 +prob.height = 1200.0 +prob.theta_0 = 300.0 +prob.theta_tr = 343.0 +prob.T_tr = 213.0 +prob.x_c = 0.0 +prob.z_c = 1500.0 +prob.x_r = 4000.0 +prob.z_r = 1500.0 +prob.theta_c = 3.0 diff --git a/Exec/SquallLine_2D/inputs_moisture_outflow b/Exec/SquallLine_2D/inputs_moisture_outflow deleted file mode 100644 index 0411708df..000000000 --- a/Exec/SquallLine_2D/inputs_moisture_outflow +++ /dev/null @@ -1,86 +0,0 @@ -# ------------------ 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 index be547c6e9..cfd38c304 100644 --- a/Exec/SquallLine_2D/prob.H +++ b/Exec/SquallLine_2D/prob.H @@ -13,6 +13,18 @@ #include "TileNoZ.H" struct ProbParm : ProbParmDefaults { + amrex::Real z_tr = 12000.0; + amrex::Real height = 1200.0; + amrex::Real theta_0 = 300.0; + amrex::Real theta_tr = 343.0; + amrex::Real T_tr = 213.0; + + amrex::Real x_c = 0.0; + amrex::Real z_c = 1.5e3; + amrex::Real x_r = 4.0e3; + amrex::Real z_r = 1.5e3; + amrex::Real theta_c = 3.0; + amrex::Real C_p = 1004.0; amrex::Real Theta_0 = 300.0; }; // namespace ProbParm @@ -70,6 +82,55 @@ public: amrex::Vector& thetabar, amrex::Geometry const& geom) override; + amrex::Real compute_theta (amrex::Real z); + + amrex::Real compute_relative_humidity (const amrex::Real& z, + const amrex::Real& p_b, + const amrex::Real& T_b); + + amrex::Real vapor_mixing_ratio (const amrex::Real& z, + const amrex::Real& p_b, + const amrex::Real& T_b, + const amrex::Real& RH); + + amrex::Real compute_p_k (amrex::Real& p_k, + const amrex::Real p_k_minus_1, + amrex::Real& theta_k, + amrex::Real& rho_k, + amrex::Real& q_v_k, + amrex::Real& T_dp, + amrex::Real& T_b, + const amrex::Real dz, + const amrex::Real z, + const amrex::Real rho_k_minus_1); + + amrex::Real compute_F (const amrex::Real& p_k, + const amrex::Real& p_k_minus_1, + amrex::Real &theta_k, + amrex::Real& rho_k, + amrex::Real& q_v_k, + amrex::Real& T_dp, + amrex::Real& T_b, + const amrex::Real& dz, + const amrex::Real& z, + const amrex::Real& rho_k_minus_1); + + void compute_rho (const amrex::Real& z, + const amrex::Real& pressure, + amrex::Real &theta, + amrex::Real& rho, + amrex::Real& q_v, + amrex::Real& T_dp, + amrex::Real& T_b); + + void init_isentropic_hse_no_terrain(amrex::Real *theta, + amrex::Real* r, + amrex::Real* p, + amrex::Real *q_v, + const amrex::Real& dz, + const amrex::Real& prob_lo_z, + const int& khi); + protected: std::string name () override { return "Supercell"; } diff --git a/Exec/SquallLine_2D/prob.cpp b/Exec/SquallLine_2D/prob.cpp index f6e77845f..06ef5607a 100644 --- a/Exec/SquallLine_2D/prob.cpp +++ b/Exec/SquallLine_2D/prob.cpp @@ -14,11 +14,21 @@ Problem::Problem () { // Parse params amrex::ParmParse pp("prob"); + pp.query("z_tr", parms.z_tr); + pp.query("height", parms.height); + pp.query("theta_0", parms.theta_0); + pp.query("theta_tr", parms.theta_tr); + pp.query("T_tr", parms.T_tr); + pp.query("x_c", parms.x_c); + pp.query("z_c", parms.z_c); + pp.query("x_r", parms.x_r); + pp.query("z_r", parms.z_r); + pp.query("theta_c", parms.theta_c); } -Real compute_theta (Real z) +Real Problem::compute_theta (const Real z) { - Real theta_0=300.0, theta_tr=343.0, T_tr=213.0, z_tr = 12000.0; + Real theta_0 = parms.theta_0, theta_tr = parms.theta_tr, T_tr = parms.T_tr, z_tr = parms.z_tr; if(z <= z_tr) { return theta_0 + (theta_tr - theta_0)*std::pow(z/z_tr,1.25); } else { @@ -43,10 +53,10 @@ Real compute_saturation_pressure (const Real T_b) AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE -Real compute_relative_humidity (Real z, Real p_b, Real T_b) +Real Problem::compute_relative_humidity (const Real& z, const Real& p_b, const Real& T_b) { - Real height = 1200.0; - Real z_tr = 12000.0; + Real height = parms.height; + Real z_tr = parms.z_tr; Real p_s = compute_saturation_pressure(T_b); @@ -63,14 +73,14 @@ Real compute_relative_humidity (Real z, Real p_b, Real T_b) AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE -Real compute_vapor_pressure (Real p_s, Real RH) +Real compute_vapor_pressure (const Real p_s, const Real RH) { return p_s*RH; } AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE -Real vapor_mixing_ratio (const Real z, const Real p_b, const Real T_b, const Real RH) +Real Problem::vapor_mixing_ratio (const Real& z, const Real& p_b, const Real& T_b, const Real& RH) { Real height = 1200.0; Real p_s = compute_saturation_pressure(T_b); @@ -87,14 +97,14 @@ Real vapor_mixing_ratio (const Real z, const Real p_b, const Real T_b, const Rea AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE -Real compute_temperature (const Real p_b, const Real theta_b) +Real compute_temperature (const Real& p_b, const Real& theta_b) { return theta_b*std::pow(p_b/p_0,R_d/Cp_d); } AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE -Real compute_dewpoint_temperature (const Real T_b, const Real RH) +Real compute_dewpoint_temperature (const Real& T_b, const Real& RH) { Real T_dp, gamma, T; @@ -108,7 +118,7 @@ Real compute_dewpoint_temperature (const Real T_b, const Real RH) return T_dp; } -void compute_rho (const Real z, const Real pressure, Real &theta, Real& rho, Real& q_v, Real& T_dp, Real& T_b) +void Problem::compute_rho (const Real& z, const Real& pressure, Real& theta, Real& rho, Real& q_v, Real& T_dp, Real& T_b) { theta = compute_theta(z); @@ -120,7 +130,7 @@ void compute_rho (const Real z, const Real pressure, Real &theta, Real& rho, Rea T_dp = compute_dewpoint_temperature(T_b, RH); } -double compute_F (const Real& p_k, const Real& p_k_minus_1, Real &theta_k, Real& rho_k, Real& q_v_k, +Real Problem::compute_F (const Real& p_k, const Real& p_k_minus_1, Real &theta_k, Real& rho_k, Real& q_v_k, Real& T_dp, Real& T_b, const Real& dz, const Real& z, const Real& rho_k_minus_1) { @@ -140,7 +150,7 @@ double compute_F (const Real& p_k, const Real& p_k_minus_1, Real &theta_k, Real& return F; } -Real compute_p_k (Real &p_k, const Real p_k_minus_1, Real &theta_k, Real& rho_k, Real& q_v_k, +Real Problem::compute_p_k (Real &p_k, const Real p_k_minus_1, Real &theta_k, Real& rho_k, Real& q_v_k, Real& T_dp, Real& T_b, const Real dz, const Real z, const Real rho_k_minus_1) { @@ -157,10 +167,11 @@ Real compute_p_k (Real &p_k, const Real p_k_minus_1, Real &theta_k, Real& rho_k, } void -init_isentropic_hse_no_terrain(Real *theta, Real* r, Real* p, Real *q_v, +Problem::init_isentropic_hse_no_terrain(Real *theta, Real* r, Real* p, Real *q_v, const Real& dz, const Real& prob_lo_z, const int& khi) { + FILE *file_IC; file_IC = fopen("input_sounding_probcpp.txt","w"); Real z, T_b, T_dp; @@ -190,6 +201,7 @@ Problem::erf_init_dens_hse_moist (MultiFab& rho_hse, Geometry const& geom) { + const Real prob_lo_z = geom.ProbLo()[2]; const Real dz = geom.CellSize()[2]; const int khi = geom.Domain().bigEnd()[2]; @@ -319,7 +331,7 @@ Problem::init_custom_pert ( Real* p = d_p.data(); - const Real x_c = 0.0, z_c = 1.5e3, x_r = 4.0e3, z_r = 1.5e3, r_c = 1.0, theta_c = 3.0; + const Real x_c = parms.x_c, z_c = parms.z_c, x_r = parms.x_r, z_r = parms.z_r, theta_c = parms.theta_c, r_c = 1.0; //const Real x_c = 0.0, z_c = 2.0e3, x_r = 10.0e3, z_r = 1.5e3, r_c = 1.0, theta_c = 3.0; Real Rd_by_Cp = sc.rdOcp; From 43fd3582c681152e827869eb0ac23a3ec0d0f853 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 29 Nov 2023 08:01:30 -0800 Subject: [PATCH 23/37] Commenting out IC file writing --- Exec/SquallLine_2D/prob.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Exec/SquallLine_2D/prob.cpp b/Exec/SquallLine_2D/prob.cpp index 06ef5607a..b24260ebc 100644 --- a/Exec/SquallLine_2D/prob.cpp +++ b/Exec/SquallLine_2D/prob.cpp @@ -172,24 +172,24 @@ Problem::init_isentropic_hse_no_terrain(Real *theta, Real* r, Real* p, Real *q_v const int& khi) { - FILE *file_IC; - file_IC = fopen("input_sounding_probcpp.txt","w"); + //FILE *file_IC; + //file_IC = fopen("input_sounding_probcpp.txt","w"); Real z, T_b, T_dp; // Compute the quantities at z = 0.5*dz (first cell center) z = prob_lo_z + 0.5*dz; p[0] = p_0; compute_p_k(p[0], p_0, theta[0], r[0], q_v[0], T_dp, T_b, dz, z, 0.0); - fprintf(file_IC, "%0.15g %0.15g %0.15g %0.15g %0.15g %0.15g %0.15g\n", z, T_b-273.15, T_dp, p[0], r[0], theta[0], q_v[0]); + //fprintf(file_IC, "%0.15g %0.15g %0.15g %0.15g %0.15g %0.15g %0.15g\n", z, T_b-273.15, T_dp, p[0], r[0], theta[0], q_v[0]); for (int k=1;k<=khi;k++){ z = prob_lo_z + (k+0.5)*dz; p[k] = p[k-1]; compute_p_k(p[k], p[k-1], theta[k], r[k], q_v[k], T_dp, T_b, dz, z, r[k-1]); - fprintf(file_IC, "%0.15g %0.15g %0.15g %0.15g %0.15g %0.15g %0.15g\n", z, T_b-273.15, T_dp, p[k], r[k], theta[k], q_v[k]); + //fprintf(file_IC, "%0.15g %0.15g %0.15g %0.15g %0.15g %0.15g %0.15g\n", z, T_b-273.15, T_dp, p[k], r[k], theta[k], q_v[k]); } - fclose(file_IC); + //fclose(file_IC); r[khi+1] = r[khi]; From 07e01e9565b4f75ec3086359a5f2a350f9585fb9 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 29 Nov 2023 08:37:03 -0800 Subject: [PATCH 24/37] Updating prob.H with host device macros --- Exec/SquallLine_2D/prob.H | 4 ++++ Exec/SquallLine_2D/prob.cpp | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/Exec/SquallLine_2D/prob.H b/Exec/SquallLine_2D/prob.H index cfd38c304..3040122e6 100644 --- a/Exec/SquallLine_2D/prob.H +++ b/Exec/SquallLine_2D/prob.H @@ -84,10 +84,14 @@ public: amrex::Real compute_theta (amrex::Real z); + AMREX_FORCE_INLINE + AMREX_GPU_HOST_DEVICE amrex::Real compute_relative_humidity (const amrex::Real& z, const amrex::Real& p_b, const amrex::Real& T_b); + AMREX_FORCE_INLINE + AMREX_GPU_HOST_DEVICE amrex::Real vapor_mixing_ratio (const amrex::Real& z, const amrex::Real& p_b, const amrex::Real& T_b, diff --git a/Exec/SquallLine_2D/prob.cpp b/Exec/SquallLine_2D/prob.cpp index b24260ebc..2acb16208 100644 --- a/Exec/SquallLine_2D/prob.cpp +++ b/Exec/SquallLine_2D/prob.cpp @@ -82,7 +82,7 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE Real Problem::vapor_mixing_ratio (const Real& z, const Real& p_b, const Real& T_b, const Real& RH) { - Real height = 1200.0; + Real height = parms.height; Real p_s = compute_saturation_pressure(T_b); Real p_v = compute_vapor_pressure(p_s, RH); From f058575fb902f5c21e33a1855fb5330f0cb89465 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 29 Nov 2023 09:12:44 -0800 Subject: [PATCH 25/37] Correcting argument list --- Exec/SquallLine_2D/prob.H | 22 +++++++++++----------- Exec/SquallLine_2D/prob.cpp | 8 ++++---- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/Exec/SquallLine_2D/prob.H b/Exec/SquallLine_2D/prob.H index 3040122e6..0d7a7809d 100644 --- a/Exec/SquallLine_2D/prob.H +++ b/Exec/SquallLine_2D/prob.H @@ -86,23 +86,23 @@ public: AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE - amrex::Real compute_relative_humidity (const amrex::Real& z, - const amrex::Real& p_b, - const amrex::Real& T_b); + amrex::Real compute_relative_humidity (const amrex::Real z, + const amrex::Real p_b, + const amrex::Real T_b); AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE - amrex::Real vapor_mixing_ratio (const amrex::Real& z, - const amrex::Real& p_b, - const amrex::Real& T_b, - const amrex::Real& RH); + amrex::Real vapor_mixing_ratio (const amrex::Real z, + const amrex::Real p_b, + const amrex::Real T_b, + const amrex::Real RH); amrex::Real compute_p_k (amrex::Real& p_k, const amrex::Real p_k_minus_1, amrex::Real& theta_k, amrex::Real& rho_k, amrex::Real& q_v_k, - amrex::Real& T_dp, + amrex::Real& T_dp, amrex::Real& T_b, const amrex::Real dz, const amrex::Real z, @@ -113,7 +113,7 @@ public: amrex::Real &theta_k, amrex::Real& rho_k, amrex::Real& q_v_k, - amrex::Real& T_dp, + amrex::Real& T_dp, amrex::Real& T_b, const amrex::Real& dz, const amrex::Real& z, @@ -131,9 +131,9 @@ public: amrex::Real* r, amrex::Real* p, amrex::Real *q_v, - const amrex::Real& dz, + const amrex::Real& dz, const amrex::Real& prob_lo_z, - const int& khi); + const int& khi); protected: std::string name () override { return "Supercell"; } diff --git a/Exec/SquallLine_2D/prob.cpp b/Exec/SquallLine_2D/prob.cpp index 2acb16208..a8b3a8f54 100644 --- a/Exec/SquallLine_2D/prob.cpp +++ b/Exec/SquallLine_2D/prob.cpp @@ -53,7 +53,7 @@ Real compute_saturation_pressure (const Real T_b) AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE -Real Problem::compute_relative_humidity (const Real& z, const Real& p_b, const Real& T_b) +Real Problem::compute_relative_humidity (const Real z, const Real p_b, const Real T_b) { Real height = parms.height; Real z_tr = parms.z_tr; @@ -80,7 +80,7 @@ Real compute_vapor_pressure (const Real p_s, const Real RH) AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE -Real Problem::vapor_mixing_ratio (const Real& z, const Real& p_b, const Real& T_b, const Real& RH) +Real Problem::vapor_mixing_ratio (const Real z, const Real p_b, const Real T_b, const Real RH) { Real height = parms.height; Real p_s = compute_saturation_pressure(T_b); @@ -97,14 +97,14 @@ Real Problem::vapor_mixing_ratio (const Real& z, const Real& p_b, const Real& T_ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE -Real compute_temperature (const Real& p_b, const Real& theta_b) +Real compute_temperature (const Real p_b, const Real theta_b) { return theta_b*std::pow(p_b/p_0,R_d/Cp_d); } AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE -Real compute_dewpoint_temperature (const Real& T_b, const Real& RH) +Real compute_dewpoint_temperature (const Real T_b, const Real RH) { Real T_dp, gamma, T; From c07e3f8f445435f7a6d51ee035b14f9562e8fbb9 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 29 Nov 2023 09:53:32 -0800 Subject: [PATCH 26/37] Correcting argument list --- Exec/SquallLine_2D/prob.H | 3 +++ Exec/SquallLine_2D/prob.cpp | 20 ++++++++------------ 2 files changed, 11 insertions(+), 12 deletions(-) diff --git a/Exec/SquallLine_2D/prob.H b/Exec/SquallLine_2D/prob.H index 0d7a7809d..cd4556b6b 100644 --- a/Exec/SquallLine_2D/prob.H +++ b/Exec/SquallLine_2D/prob.H @@ -87,12 +87,15 @@ public: AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE amrex::Real compute_relative_humidity (const amrex::Real z, + const amrex::Real height, + const amrex::Real z_tr, const amrex::Real p_b, const amrex::Real T_b); AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE amrex::Real vapor_mixing_ratio (const amrex::Real z, + const amrex::Real height, const amrex::Real p_b, const amrex::Real T_b, const amrex::Real RH); diff --git a/Exec/SquallLine_2D/prob.cpp b/Exec/SquallLine_2D/prob.cpp index a8b3a8f54..23e43a241 100644 --- a/Exec/SquallLine_2D/prob.cpp +++ b/Exec/SquallLine_2D/prob.cpp @@ -53,11 +53,8 @@ Real compute_saturation_pressure (const Real T_b) AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE -Real Problem::compute_relative_humidity (const Real z, const Real p_b, const Real T_b) +Real Problem::compute_relative_humidity (const Real z, const Real height, const Real z_tr, const Real p_b, const Real T_b) { - Real height = parms.height; - Real z_tr = parms.z_tr; - Real p_s = compute_saturation_pressure(T_b); Real q_s = 0.622*p_s/(p_b - p_s); @@ -80,9 +77,8 @@ Real compute_vapor_pressure (const Real p_s, const Real RH) AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE -Real Problem::vapor_mixing_ratio (const Real z, const Real p_b, const Real T_b, const Real RH) +Real Problem::vapor_mixing_ratio (const Real z, const Real height, const Real p_b, const Real T_b, const Real RH) { - Real height = parms.height; Real p_s = compute_saturation_pressure(T_b); Real p_v = compute_vapor_pressure(p_s, RH); @@ -123,8 +119,8 @@ void Problem::compute_rho (const Real& z, const Real& pressure, Real& theta, Rea theta = compute_theta(z); T_b = compute_temperature(pressure, theta); - Real RH = compute_relative_humidity(z, pressure, T_b); - q_v = vapor_mixing_ratio(z, pressure, T_b, RH); + Real RH = compute_relative_humidity(z, parms.height, parms.z_tr, pressure, T_b); + q_v = vapor_mixing_ratio(z, parms.height, pressure, T_b, RH); rho = pressure/(R_d*T_b*(1.0 + (R_v/R_d)*q_v)); rho = rho*(1.0 + q_v); T_dp = compute_dewpoint_temperature(T_b, RH); @@ -362,15 +358,15 @@ Problem::init_custom_pert ( theta_total = t[k] + delta_theta; temperature = compute_temperature(p[k], theta_total); Real T_b = compute_temperature(p[k], t[k]); - RH = compute_relative_humidity(z, p[k], T_b); - Real q_v_hot = vapor_mixing_ratio(z, p[k], T_b, RH); + RH = compute_relative_humidity(z, parms.height, parms.z_tr, p[k], T_b); + Real q_v_hot = vapor_mixing_ratio(z, parms.height, p[k], T_b, RH); rho = p[k]/(R_d*temperature*(1.0 + (R_v/R_d)*q_v_hot)); // Compute background quantities Real temperature_back = compute_temperature(p[k], t[k]); Real T_back = compute_temperature(p[k], t[k]); - Real RH_back = compute_relative_humidity(z, p[k], T_back); - Real q_v_back = vapor_mixing_ratio(z, p[k], T_back, RH_back); + Real RH_back = compute_relative_humidity(z, parms.height, parms.z_tr, p[k], T_back); + Real q_v_back = vapor_mixing_ratio(z, parms.height, p[k], T_back, RH_back); Real rho_back = p[k]/(R_d*temperature_back*(1.0 + (R_v/R_d)*q_v_back)); // This version perturbs rho but not p From 7efb72c9546a1eaedd2dc9084e634d69bab20d4d Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 29 Nov 2023 10:50:32 -0800 Subject: [PATCH 27/37] Correcting CUDA compilation errors --- Exec/SquallLine_2D/prob.cpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/Exec/SquallLine_2D/prob.cpp b/Exec/SquallLine_2D/prob.cpp index 23e43a241..ec9789d26 100644 --- a/Exec/SquallLine_2D/prob.cpp +++ b/Exec/SquallLine_2D/prob.cpp @@ -333,7 +333,9 @@ Problem::init_custom_pert ( Real Rd_by_Cp = sc.rdOcp; Rd_by_Cp = Rd_by_Cp; - amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + auto lparms = parms; + + amrex::ParallelFor(bx, [=, lparms=lparms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Geometry (note we must include these here to get the data on device) const auto prob_lo = geomdata.ProbLo(); @@ -358,15 +360,15 @@ Problem::init_custom_pert ( theta_total = t[k] + delta_theta; temperature = compute_temperature(p[k], theta_total); Real T_b = compute_temperature(p[k], t[k]); - RH = compute_relative_humidity(z, parms.height, parms.z_tr, p[k], T_b); - Real q_v_hot = vapor_mixing_ratio(z, parms.height, p[k], T_b, RH); + RH = compute_relative_humidity(z, lparms.height, lparms.z_tr, p[k], T_b); + Real q_v_hot = vapor_mixing_ratio(z, lparms.height, p[k], T_b, RH); rho = p[k]/(R_d*temperature*(1.0 + (R_v/R_d)*q_v_hot)); // Compute background quantities Real temperature_back = compute_temperature(p[k], t[k]); Real T_back = compute_temperature(p[k], t[k]); - Real RH_back = compute_relative_humidity(z, parms.height, parms.z_tr, p[k], T_back); - Real q_v_back = vapor_mixing_ratio(z, parms.height, p[k], T_back, RH_back); + Real RH_back = compute_relative_humidity(z, lparms.height, lparms.z_tr, p[k], T_back); + Real q_v_back = vapor_mixing_ratio(z, lparms.height, p[k], T_back, RH_back); Real rho_back = p[k]/(R_d*temperature_back*(1.0 + (R_v/R_d)*q_v_back)); // This version perturbs rho but not p From db8b0ba2b4dcb4b25bbbf7108d0b00206e66e7b0 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 29 Nov 2023 11:31:25 -0800 Subject: [PATCH 28/37] Correcting CUDA compilation errors --- Exec/SquallLine_2D/prob.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/Exec/SquallLine_2D/prob.cpp b/Exec/SquallLine_2D/prob.cpp index ec9789d26..85006ff02 100644 --- a/Exec/SquallLine_2D/prob.cpp +++ b/Exec/SquallLine_2D/prob.cpp @@ -333,9 +333,10 @@ Problem::init_custom_pert ( Real Rd_by_Cp = sc.rdOcp; Rd_by_Cp = Rd_by_Cp; - auto lparms = parms; + Real height = parms.height; + Real z_tr = parms.z_tr; - amrex::ParallelFor(bx, [=, lparms=lparms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Geometry (note we must include these here to get the data on device) const auto prob_lo = geomdata.ProbLo(); @@ -360,15 +361,15 @@ Problem::init_custom_pert ( theta_total = t[k] + delta_theta; temperature = compute_temperature(p[k], theta_total); Real T_b = compute_temperature(p[k], t[k]); - RH = compute_relative_humidity(z, lparms.height, lparms.z_tr, p[k], T_b); - Real q_v_hot = vapor_mixing_ratio(z, lparms.height, p[k], T_b, RH); + RH = compute_relative_humidity(z, height, z_tr, p[k], T_b); + Real q_v_hot = vapor_mixing_ratio(z, height, p[k], T_b, RH); rho = p[k]/(R_d*temperature*(1.0 + (R_v/R_d)*q_v_hot)); // Compute background quantities Real temperature_back = compute_temperature(p[k], t[k]); Real T_back = compute_temperature(p[k], t[k]); - Real RH_back = compute_relative_humidity(z, lparms.height, lparms.z_tr, p[k], T_back); - Real q_v_back = vapor_mixing_ratio(z, lparms.height, p[k], T_back, RH_back); + Real RH_back = compute_relative_humidity(z, height, z_tr, p[k], T_back); + Real q_v_back = vapor_mixing_ratio(z, height, p[k], T_back, RH_back); Real rho_back = p[k]/(R_d*temperature_back*(1.0 + (R_v/R_d)*q_v_back)); // This version perturbs rho but not p From b5adfd973073a1f7352033639aae4e9771b92cc5 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 29 Nov 2023 12:08:50 -0800 Subject: [PATCH 29/37] Correcting CUDA compilation errors --- Exec/SquallLine_2D/prob.H | 16 ---------------- Exec/SquallLine_2D/prob.cpp | 4 ++-- 2 files changed, 2 insertions(+), 18 deletions(-) diff --git a/Exec/SquallLine_2D/prob.H b/Exec/SquallLine_2D/prob.H index cd4556b6b..90f8f73f5 100644 --- a/Exec/SquallLine_2D/prob.H +++ b/Exec/SquallLine_2D/prob.H @@ -84,22 +84,6 @@ public: amrex::Real compute_theta (amrex::Real z); - AMREX_FORCE_INLINE - AMREX_GPU_HOST_DEVICE - amrex::Real compute_relative_humidity (const amrex::Real z, - const amrex::Real height, - const amrex::Real z_tr, - const amrex::Real p_b, - const amrex::Real T_b); - - AMREX_FORCE_INLINE - AMREX_GPU_HOST_DEVICE - amrex::Real vapor_mixing_ratio (const amrex::Real z, - const amrex::Real height, - const amrex::Real p_b, - const amrex::Real T_b, - const amrex::Real RH); - amrex::Real compute_p_k (amrex::Real& p_k, const amrex::Real p_k_minus_1, amrex::Real& theta_k, diff --git a/Exec/SquallLine_2D/prob.cpp b/Exec/SquallLine_2D/prob.cpp index 85006ff02..973c924f4 100644 --- a/Exec/SquallLine_2D/prob.cpp +++ b/Exec/SquallLine_2D/prob.cpp @@ -53,7 +53,7 @@ Real compute_saturation_pressure (const Real T_b) AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE -Real Problem::compute_relative_humidity (const Real z, const Real height, const Real z_tr, const Real p_b, const Real T_b) +Real compute_relative_humidity (const Real z, const Real height, const Real z_tr, const Real p_b, const Real T_b) { Real p_s = compute_saturation_pressure(T_b); @@ -77,7 +77,7 @@ Real compute_vapor_pressure (const Real p_s, const Real RH) AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE -Real Problem::vapor_mixing_ratio (const Real z, const Real height, const Real p_b, const Real T_b, const Real RH) +Real vapor_mixing_ratio (const Real z, const Real height, const Real p_b, const Real T_b, const Real RH) { Real p_s = compute_saturation_pressure(T_b); Real p_v = compute_vapor_pressure(p_s, RH); From 17e6253fb637384bfbd7f3e51e32fadeed71d386 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 30 Nov 2023 09:29:18 -0800 Subject: [PATCH 30/37] Some changes to docs --- Docs/sphinx_doc/theory/Buoyancy.rst | 52 ++++++++++++++++++ Docs/sphinx_doc/theory/WetEquations.rst | 70 +++++++++++++++++++++++++ 2 files changed, 122 insertions(+) diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index 6213c379f..4c4238b80 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -109,3 +109,55 @@ which removes the need to compute horizontal averages of these quantities. Thi We note that this version of the buoyancy force matches that given in Marat F. Khairoutdinov and David A. Randall's paper (J. Atm Sciences, 607, 1983) if we neglect :math:`\frac{p^\prime}{\bar{p_0}}`. + +Type 4 +------ +.. math:: + + \begin{equation} + \mathbf{B} = \rho'\mathbf{g} \approx -\rho\Bigg(\frac{T'}{T} + 0.61 q_v' - q_c - q_p - \frac{p'}{p}\Bigg), + \end{equation} + +The derivation follows. The total density is given by :math:`\rho = \rho_d(1 + q_v + q_c + q_p)`, which can be written as + +.. math:: + + \rho = \frac{p (1 + q_v + q_c + q_p)}{R_dT\Bigg(1 + \cfrac{R_v}{R_d}q_v\Bigg)} + +This can be written using binomial expansion as + +.. math:: + + \begin{align*} + \rho &= \frac{p}{R_dT} (1 + q_v + q_c + q_p)\Bigg(1 + \frac{R_v}{R_d}q_v\Bigg)^{-1} \\ + &= \frac{p}{R_dT} (1 + q_v + q_c + q_p)\Bigg(1 - \frac{R_v}{R_d}q_v + O(q_v^2)\Bigg) \\ + &= \frac{p}{R_dT}\Bigg(1 + q_v + q_c + q_p - \frac{R_v}{R_d}q_v + \text{H.O.T. such as } O(q_v^2) + O(q_vq_c)\Bigg) \\ + &\approx \frac{p}{R_dT}\Bigg(1 + q_v + q_c + q_p - \frac{R_v}{R_d}q_v\Bigg) + \end{align*} + +Taking log on both sides, we get + +.. math:: + + \log{\rho} = \log{p} - \log{R_d} - \log{T} + \log(1 - 0.61 q_v + q_c + q_p) + +Taking derivative gives + +.. math:: + + \frac{\rho'}{\rho} = \frac{p'}{p} - \frac{T'}{T} + \frac{(-0.61 q_v' + q_c' + q_p')}{(1 - 0.61 q_v + q_c + q_p)} + +Using :math:`- 0.61 q_v + q_c + q_p \ll 1`, we have + +.. math:: + + \frac{\rho'}{\rho} = \frac{p'}{p} - \frac{T'}{T} + (-0.61 q_v' + q_c' + q_p') + +Since the background values of cloud water and precipitate mass mixing ratios -- :math:`q_c` and :math:`q_p` are zero, we have :math:`q_c' = q_c` and :math:`q_p' = q_p`. Hence, we have + +.. math:: + + \begin{equation} + \rho'\approx -\rho\Bigg(\frac{T'}{T} + 0.61 q_v' - q_c - q_p - \frac{p'}{p}\Bigg), + \end{equation} + diff --git a/Docs/sphinx_doc/theory/WetEquations.rst b/Docs/sphinx_doc/theory/WetEquations.rst index 38288e484..9f231d4af 100644 --- a/Docs/sphinx_doc/theory/WetEquations.rst +++ b/Docs/sphinx_doc/theory/WetEquations.rst @@ -139,3 +139,73 @@ In this set of equations, the subgrid turbulent parameterization effects are inc of water vapor to/from water through condensation/evaporation, which is determined by the microphysics parameterization processes. :math:`\mathbf{B}` is the buoyancy force, which is defined in :ref:`Buoyancy `. +Initial condition with second-order integration of the hydrostatic equation +============================================================================= + +We have the hydrostatic equation given by + +.. math:: + + \frac{\partial p}{\partial z} = -\rho g, + +where :math:`\rho = \rho_d(1 + q_v)`, :math:`\rho_d` is the dry density, and :math:`q_v` is the mass mixing ratio of water vapor. Using an average value of :math:`\rho` for the integration, we get + +.. math:: + + p(k) = p(k-1) - \frac{(\rho(k-1) + \rho(k))}{2} g\Delta z. + +The density at a point is a function of the pressure, potential temperature, and relative humidity. The latter two quantities are computed using user-specified profiles, and hence, for simplicity, we write :math:`\rho(k) = f(p(k))`. Hence + +.. math:: + + p(k) = p(k-1) - \frac{\rho(k-1)}{2}g\Delta z - \frac{f(p(k))}{2}g\Delta z. + +Now, we define + +.. math:: + + F(p(k)) \equiv p(k) - p(k-1) + \frac{\rho(k-1)}{2}g\Delta z + \frac{f(p(k))}{2}g\Delta z = 0. + +This is a non-linear equation in :math:`p(k)`. Consider a Newton-Raphson iteration (where :math:`n` denotes the iteration number) procedure + +.. math:: + + F(p+\delta p) \approx F(p) + \delta p \frac{\partial F}{\partial p} = 0, + +which implies + +.. math:: + + \delta p = -\frac{F}{F'}, + +with the gradient being evaluated as + +.. math:: + + F' = \frac{F(p+\epsilon) - F(p)}{\epsilon}, + +and the iteration update is given by + +.. math:: + + p^{n+1} = p^n + \delta p. + +For the first cell (:math:`k=0`), which is at a height of :math:`z = \frac{\Delta z}{2}` from the base, we have + +.. math:: + + p(k) = p_0 - \rho(k)g\frac{\Delta z}{2}, + +where :math:`p_0 = 1e5 \, \text{N/m}^2` is the pressure at the base. Hence, we define + +.. math:: + + F(p(k)) \equiv p(k) - p_0 + \rho(k)g\frac{\Delta z}{2}, + +and the Newton-Raphson procedure is the same. + + + + + + From 56957bd18dd54149bf21dce22ba35eebb3710874 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Thu, 30 Nov 2023 09:30:58 -0800 Subject: [PATCH 31/37] Updating inputs --- Exec/SquallLine_2D/inputs_moisture | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Exec/SquallLine_2D/inputs_moisture b/Exec/SquallLine_2D/inputs_moisture index 2e2757fe5..7ede7c674 100644 --- a/Exec/SquallLine_2D/inputs_moisture +++ b/Exec/SquallLine_2D/inputs_moisture @@ -68,9 +68,9 @@ erf.les_type = "None" 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.dynamicViscosity = 100.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s erf.alpha_T = 00.0 # [m^2/s] -erf.alpha_C = 25.0 +erf.alpha_C = 100.0 erf.moisture_model = "Kessler" erf.use_moist_background = true From 14073876bc2c1f404de1303a68ed1cd111d9631d Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Sat, 2 Dec 2023 16:38:09 -0800 Subject: [PATCH 32/37] Adding inputs for Gabersek test case: --- Exec/SquallLine_2D/inputs_moisture_Gabersek | 91 +++++++++++++++++++++ 1 file changed, 91 insertions(+) create mode 100644 Exec/SquallLine_2D/inputs_moisture_Gabersek diff --git a/Exec/SquallLine_2D/inputs_moisture_Gabersek b/Exec/SquallLine_2D/inputs_moisture_Gabersek new file mode 100644 index 000000000..ad5b642b4 --- /dev/null +++ b/Exec/SquallLine_2D/inputs_moisture_Gabersek @@ -0,0 +1,91 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 18000 +stop_time = 90000.0 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 2048 1024 2048 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_lo = -60000. 0. 0. +geometry.prob_hi = 60000. 400. 20000. +amr.n_cell = 384 4 128 # 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 = 0.5 # fixed time step [s] -- Straka et al 1993 +erf.fixed_fast_dt = 0.25 # 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 = 1 +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 = 100.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.z_tr = 12000.0 +prob.height = 1200.0 +prob.theta_0 = 300.0 +prob.theta_tr = 343.0 +prob.T_tr = 213.0 +prob.x_c = 0.0 +prob.z_c = 2000.0 +prob.x_r = 10000.0 +prob.z_r = 1500.0 +prob.theta_c = 3.0 From c94b22aaed5f68ff1ef8653b17e65cf0728c4d0a Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Tue, 5 Dec 2023 15:34:20 -0800 Subject: [PATCH 33/37] Updating docs --- Docs/sphinx_doc/index.rst | 1 + Docs/sphinx_doc/theory/Initialization.rst | 155 ++++++++++++++++++++++ Docs/sphinx_doc/theory/Microphysics.rst | 30 +++-- Docs/sphinx_doc/theory/WetEquations.rst | 71 ---------- 4 files changed, 174 insertions(+), 83 deletions(-) create mode 100644 Docs/sphinx_doc/theory/Initialization.rst diff --git a/Docs/sphinx_doc/index.rst b/Docs/sphinx_doc/index.rst index 62e87c7ce..5e57dfa3d 100644 --- a/Docs/sphinx_doc/index.rst +++ b/Docs/sphinx_doc/index.rst @@ -49,6 +49,7 @@ In addition to this documentation, there is API documentation for ERF generated theory/DryEquations.rst theory/WetEquations.rst + theory/Initialization.rst theory/Buoyancy.rst theory/Microphysics.rst theory/DNSvsLES.rst diff --git a/Docs/sphinx_doc/theory/Initialization.rst b/Docs/sphinx_doc/theory/Initialization.rst new file mode 100644 index 000000000..d1fa66279 --- /dev/null +++ b/Docs/sphinx_doc/theory/Initialization.rst @@ -0,0 +1,155 @@ + + .. role:: cpp(code) + :language: c++ + + .. role:: f(code) + :language: fortran + +.. _Initialization: + +Initialization +================== + +To initialize a background (base) state for a simulation with a hydrostatic atmosphere, the hydrostatic equation balancing the pressure gradient +and gravity must be satisfied. This section details the procedure for the initialization of the background state. The procedure is similar for the +cases with and without moisture, the only difference being that the density for the cases with moisture has to be the total density +:math:`\rho = \rho_d(1 + q_t)`, where :math:`\rho_d` is the dry density, and :math:`q_t` is the total mass mixing ratio -- water vapor and liquid water, instead +of the dry density :math:`\rho_d` for cases without moisture. + +Computation of the dry density +------------------------------- +We express the total pressure :math:`p` as: + +.. math:: + p = \rho_d R_d T_v. + +By definition, we have: + +.. math:: + T_v = \theta_m\left(\frac{p}{p_0}\right)^{R_d/C_p}, + +.. math:: + T = \theta_d\left(\frac{p}{p_0}\right)^{R_d/C_p}. + +This gives: + +.. math:: + p = \rho_d R_d \theta_m\left(\frac{p}{p_0}\right)^{R_d/C_p}, + +which, on rearranging, yields: + +.. math:: + p = p_0\left(\frac{\rho_d R_d \theta_m}{p_0}\right)^\gamma. + +To obtain :math:`\theta_m`, consider the density of dry air: + +.. math:: + \rho_d = \frac{p - p_v}{R_d T}. + +Substituting for :math:`\rho_d` from the above equation, we get: + +.. math:: + \frac{p}{T_v} = \frac{p - p_v}{T}, + +which implies: + +.. math:: + \frac{T_v}{T} = \frac{p}{p - p_v} = \frac{1}{1-\left(\cfrac{p_v}{p}\right)}. + +We also have: + +.. math:: + q_v = \frac{\rho_v}{\rho_d} = \frac{p_v}{R_v T}\frac{R_d T}{p-p_v} = \frac{r p_v}{p - p_v}, + +where :math:`p_v` is the partial pressure of water vapor, :math:`r = R_d/R_v \approx 0.622`, and :math:`q_v` is the vapor mass mixing ratio—the ratio of the +mass of vapor to the mass of dry air. Rearranging and using :math:`q_v \ll r`, we get: + +.. math:: + \frac{p_v}{p} = \frac{1}{1 + \left(\cfrac{r}{q_v}\right)} \approx \frac{q_v}{r}, + +which, on substitution in the equation for :math:`\frac{T_v}{T}`, gives: + +.. math:: + \frac{T_v}{T} = \frac{1}{1 - \left(\cfrac{q_v}{r}\right)}. + +As :math:`q_v \ll 1`, a binomial expansion, ignoring higher-order terms, gives: + +.. math:: + T_v = T\left(1 + \frac{R_v}{R_d}q_v\right). + +Hence, the density of dry air is given by: + +.. math:: + \rho_d = \frac{p}{R_d T_v} = \frac{p}{R_d T\left(1 + \cfrac{R_v}{R_d}q_v\right)}. + + +Initialization with a second-order integration of the hydrostatic equation +---------------------------------------------------------------------------- + +We have the hydrostatic equation given by + +.. math:: + + \frac{\partial p}{\partial z} = -\rho g, + +where :math:`\rho = \rho_d(1 + q_t)`, :math:`\rho_d` is the dry density, and :math:`q_t` is the total mass mixing ratio -- water vapor and liquid water. Using an average value of :math:`\rho` for the integration from :math:`z = z(k-1)` to :math:`z(k)`, we get + +.. math:: + + p(k) = p(k-1) - \frac{(\rho(k-1) + \rho(k))}{2} g\Delta z. + +The density at a point is a function of the pressure, potential temperature, and relative humidity. The latter two quantities are computed using user-specified profiles, and hence, for simplicity, we write :math:`\rho(k) = f(p(k))`. Hence + +.. math:: + + p(k) = p(k-1) - \frac{\rho(k-1)}{2}g\Delta z - \frac{f(p(k))}{2}g\Delta z. + +Now, we define + +.. math:: + + F(p(k)) \equiv p(k) - p(k-1) + \frac{\rho(k-1)}{2}g\Delta z + \frac{f(p(k))}{2}g\Delta z = 0. + +This is a non-linear equation in :math:`p(k)`. Consider a Newton-Raphson iteration (where :math:`n` denotes the iteration number) procedure + +.. math:: + + F(p+\delta p) \approx F(p) + \delta p \frac{\partial F}{\partial p} = 0, + +which implies + +.. math:: + + \delta p = -\frac{F}{F'}, + +with the gradient being evaluated as + +.. math:: + + F' = \frac{F(p+\epsilon) - F(p)}{\epsilon}, + +and the iteration update is given by + +.. math:: + + p^{n+1} = p^n + \delta p. + +For the first cell (:math:`k=0`), which is at a height of :math:`z = \frac{\Delta z}{2}` from the base, we have + +.. math:: + + p(0) = p_0 - \rho(0)g\frac{\Delta z}{2}, + +where :math:`p_0 = 1e5 \, \text{N/m}^2` is the pressure at the base. Hence, we define + +.. math:: + + F(p(0)) \equiv p(0) - p_0 + \rho(0)g\frac{\Delta z}{2}, + +and the Newton-Raphson procedure is the same. + + + + + + diff --git a/Docs/sphinx_doc/theory/Microphysics.rst b/Docs/sphinx_doc/theory/Microphysics.rst index c87426065..6ce46670e 100644 --- a/Docs/sphinx_doc/theory/Microphysics.rst +++ b/Docs/sphinx_doc/theory/Microphysics.rst @@ -7,8 +7,12 @@ .. _Microphysics: +Microphysics model +==================== + Kessler Microphysics model -=========================== +--------------------------- +The Kessler microphysics model is a simple version of cloud microphysics which has precipitation only in the form of rain. Hence :math:`q_p = q_r`. Governing equations for the microphysical quantities for Kessler microphysics from `gabervsek2012dry`_ are .. math:: @@ -34,22 +38,24 @@ of the source terms are given below. \newgeometry{left=2cm,right=2cm,top=2cm,bottom=2cm} Rate of condensation of water vapor/evaporation of cloud water --------------------------------------------------------------- +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ From `klemp1978simulation`_, we have the following expressions. .. _`klemp1978simulation`: https://journals.ametsoc.org/view/journals/atsc/35/6/1520-0469_1978_035_1070_tsotdc_2_0_co_2.xml +If the air is not saturated, i.e. :math:`q_v > q_{vs}` + .. math:: C_c = \frac{q_v - q_{vs}}{1 + \cfrac{q_{vs}^*4093L}{C_p(T-36)^2}} -If the air is not saturated, i.e. :math:`q_v < q_s`, then cloud water evaporates to water vapor at the rate +If the air is not saturated, i.e. :math:`q_v < q_{vs}`, then cloud water evaporates to water vapor at the rate .. math:: E_c = \frac{q_{vs} - q_v}{1 + \cfrac{q_{vs}^*4093L}{C_p(T-36)^2}} Rate of autoconversion of cloud water into rain ------------------------------------------------- +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The rate of autoconversion of cloud water into rain is given by @@ -59,7 +65,7 @@ The rate of autoconversion of cloud water into rain is given by where :math:`k_1 = 0.001` s\ :sup:`-1`, :math:`a = 0.001` kg kg\ :sup:`-1`. Rate of accretion of cloud water into rain water drops ------------------------------------------------------- +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The rate of accretion of cloud water into rain water drops is given by @@ -69,7 +75,7 @@ The rate of accretion of cloud water into rain water drops is given by where :math:`k_2= 2.2` s\ :sup:`-1`. The rate of evaporation of rain into water vapor -------------------------------------------------- +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The rate of evaporation of rain into water vapor is given by @@ -79,7 +85,7 @@ The rate of evaporation of rain into water vapor is given by where the ventilation factor :math:`C = 1.6 + 124.9(\overline{\rho}q_r)^{0.2046}`. Terminal fall velocity of rain -------------------------------- +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The terminal fall velocity of rain is given by @@ -93,7 +99,7 @@ The terminal fall velocity of rain is given by Single Moment Microphysics Model -=================================== +---------------------------------- The conversion rates among the moist hydrometeors are parameterized assuming that .. math:: @@ -116,7 +122,7 @@ Marat F. Khairoutdinov and David A. Randall's (J. Atm Sciences, 607, 1983). The implementation of microphysics model in ERF is similar to the that in the SAM code (http://rossby.msrc.sunysb.edu/~marat/SAM.html) Accretion ------------------- +~~~~~~~~~~~~ There are several different type of accretional growth mechanisms that need to be included; these describe the interaction of water vapor and cloud water with rain water. @@ -131,21 +137,21 @@ The accretion of raindrops by accretion of cloud water is Q_{racw} = \frac{\pi E_{RW}n_{0R}\alpha q_{c}\Gamma(3+b)}{4\lambda_{R}^{3+b}}(\frac{\rho_{0}}{\rho})^{1/2} The bergeron Process ------------------------- +~~~~~~~~~~~~~~~~~~~~~~ The cloud water transform to snow by deposition and rimming can be written as .. math:: Q_{sfw} = N_{150}\left(\alpha_{1}m_{150}^{\alpha_{2}}+\pi E_{iw}\rho q_{c}R_{150}^{2}U_{150}\right) Autoconversion ------------------------- +~~~~~~~~~~~~~~~~~~~~~~ The collision and coalescence of cloud water to from raindrops is parameterized as following .. math:: Q_{raut} = \rho\left(q_{c}-q_{c0}\right)^{2}\left[1.2 \times 10^{-4}+{1.569 \times 10^{-12}N_{1}/[D_{0}(q_{c}-q_{c0})]}\right]^{-1} Evaporation ------------------------- +~~~~~~~~~~~~~~~~~~~~~~ The evaporation rate of rain is .. math:: diff --git a/Docs/sphinx_doc/theory/WetEquations.rst b/Docs/sphinx_doc/theory/WetEquations.rst index 9f231d4af..27e75e28f 100644 --- a/Docs/sphinx_doc/theory/WetEquations.rst +++ b/Docs/sphinx_doc/theory/WetEquations.rst @@ -138,74 +138,3 @@ In this set of equations, the subgrid turbulent parameterization effects are inc :math:`\mathbf{F}` stands for the external force, and :math:`Q` and :math:`F_Q` represent the mass and energy transformation of water vapor to/from water through condensation/evaporation, which is determined by the microphysics parameterization processes. :math:`\mathbf{B}` is the buoyancy force, which is defined in :ref:`Buoyancy `. - -Initial condition with second-order integration of the hydrostatic equation -============================================================================= - -We have the hydrostatic equation given by - -.. math:: - - \frac{\partial p}{\partial z} = -\rho g, - -where :math:`\rho = \rho_d(1 + q_v)`, :math:`\rho_d` is the dry density, and :math:`q_v` is the mass mixing ratio of water vapor. Using an average value of :math:`\rho` for the integration, we get - -.. math:: - - p(k) = p(k-1) - \frac{(\rho(k-1) + \rho(k))}{2} g\Delta z. - -The density at a point is a function of the pressure, potential temperature, and relative humidity. The latter two quantities are computed using user-specified profiles, and hence, for simplicity, we write :math:`\rho(k) = f(p(k))`. Hence - -.. math:: - - p(k) = p(k-1) - \frac{\rho(k-1)}{2}g\Delta z - \frac{f(p(k))}{2}g\Delta z. - -Now, we define - -.. math:: - - F(p(k)) \equiv p(k) - p(k-1) + \frac{\rho(k-1)}{2}g\Delta z + \frac{f(p(k))}{2}g\Delta z = 0. - -This is a non-linear equation in :math:`p(k)`. Consider a Newton-Raphson iteration (where :math:`n` denotes the iteration number) procedure - -.. math:: - - F(p+\delta p) \approx F(p) + \delta p \frac{\partial F}{\partial p} = 0, - -which implies - -.. math:: - - \delta p = -\frac{F}{F'}, - -with the gradient being evaluated as - -.. math:: - - F' = \frac{F(p+\epsilon) - F(p)}{\epsilon}, - -and the iteration update is given by - -.. math:: - - p^{n+1} = p^n + \delta p. - -For the first cell (:math:`k=0`), which is at a height of :math:`z = \frac{\Delta z}{2}` from the base, we have - -.. math:: - - p(k) = p_0 - \rho(k)g\frac{\Delta z}{2}, - -where :math:`p_0 = 1e5 \, \text{N/m}^2` is the pressure at the base. Hence, we define - -.. math:: - - F(p(k)) \equiv p(k) - p_0 + \rho(k)g\frac{\Delta z}{2}, - -and the Newton-Raphson procedure is the same. - - - - - - From b5faafc3f6ff346da8b1cde48887813fb96ffa59 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Tue, 5 Dec 2023 15:34:40 -0800 Subject: [PATCH 34/37] Updating docs --- Docs/sphinx_doc/theory/Buoyancy.rst | 12 ++++++------ Docs/sphinx_doc/theory/Initialization.rst | 10 +++++----- Docs/sphinx_doc/theory/Microphysics.rst | 4 ++-- 3 files changed, 13 insertions(+), 13 deletions(-) diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index 4c4238b80..92e987ab5 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -112,13 +112,13 @@ if we neglect :math:`\frac{p^\prime}{\bar{p_0}}`. Type 4 ------ -.. math:: +.. math:: \begin{equation} \mathbf{B} = \rho'\mathbf{g} \approx -\rho\Bigg(\frac{T'}{T} + 0.61 q_v' - q_c - q_p - \frac{p'}{p}\Bigg), \end{equation} -The derivation follows. The total density is given by :math:`\rho = \rho_d(1 + q_v + q_c + q_p)`, which can be written as +The derivation follows. The total density is given by :math:`\rho = \rho_d(1 + q_v + q_c + q_p)`, which can be written as .. math:: @@ -155,9 +155,9 @@ Using :math:`- 0.61 q_v + q_c + q_p \ll 1`, we have Since the background values of cloud water and precipitate mass mixing ratios -- :math:`q_c` and :math:`q_p` are zero, we have :math:`q_c' = q_c` and :math:`q_p' = q_p`. Hence, we have -.. math:: +.. math:: - \begin{equation} - \rho'\approx -\rho\Bigg(\frac{T'}{T} + 0.61 q_v' - q_c - q_p - \frac{p'}{p}\Bigg), - \end{equation} + \begin{equation} + \rho'\approx -\rho\Bigg(\frac{T'}{T} + 0.61 q_v' - q_c - q_p - \frac{p'}{p}\Bigg), + \end{equation} diff --git a/Docs/sphinx_doc/theory/Initialization.rst b/Docs/sphinx_doc/theory/Initialization.rst index d1fa66279..9e732c550 100644 --- a/Docs/sphinx_doc/theory/Initialization.rst +++ b/Docs/sphinx_doc/theory/Initialization.rst @@ -10,10 +10,10 @@ Initialization ================== -To initialize a background (base) state for a simulation with a hydrostatic atmosphere, the hydrostatic equation balancing the pressure gradient -and gravity must be satisfied. This section details the procedure for the initialization of the background state. The procedure is similar for the -cases with and without moisture, the only difference being that the density for the cases with moisture has to be the total density -:math:`\rho = \rho_d(1 + q_t)`, where :math:`\rho_d` is the dry density, and :math:`q_t` is the total mass mixing ratio -- water vapor and liquid water, instead +To initialize a background (base) state for a simulation with a hydrostatic atmosphere, the hydrostatic equation balancing the pressure gradient +and gravity must be satisfied. This section details the procedure for the initialization of the background state. The procedure is similar for the +cases with and without moisture, the only difference being that the density for the cases with moisture has to be the total density +:math:`\rho = \rho_d(1 + q_t)`, where :math:`\rho_d` is the dry density, and :math:`q_t` is the total mass mixing ratio -- water vapor and liquid water, instead of the dry density :math:`\rho_d` for cases without moisture. Computation of the dry density @@ -61,7 +61,7 @@ We also have: .. math:: q_v = \frac{\rho_v}{\rho_d} = \frac{p_v}{R_v T}\frac{R_d T}{p-p_v} = \frac{r p_v}{p - p_v}, -where :math:`p_v` is the partial pressure of water vapor, :math:`r = R_d/R_v \approx 0.622`, and :math:`q_v` is the vapor mass mixing ratio—the ratio of the +where :math:`p_v` is the partial pressure of water vapor, :math:`r = R_d/R_v \approx 0.622`, and :math:`q_v` is the vapor mass mixing ratio—the ratio of the mass of vapor to the mass of dry air. Rearranging and using :math:`q_v \ll r`, we get: .. math:: diff --git a/Docs/sphinx_doc/theory/Microphysics.rst b/Docs/sphinx_doc/theory/Microphysics.rst index 6ce46670e..05607a0b4 100644 --- a/Docs/sphinx_doc/theory/Microphysics.rst +++ b/Docs/sphinx_doc/theory/Microphysics.rst @@ -28,7 +28,7 @@ Governing equations for the microphysical quantities for Kessler microphysics fr where :math:`C_c` is the rate of condensation of water vapor to cloud water, :math:`E_c` is the rate of evaporation of cloud water to water vapor, :math:`A_c` is the autoconversion of cloud water to rain, :math:`K_c` is the accretion of cloud water to rain drops, :math:`E_r` is the evaporation of rain to water vapor and :math:`F_r` is the sedimentation of rain. The parametrization used is given in `klemp1978simulation`_, and is given -below. Note that in all the equations, :math:`p` is specified in millibars and :math:`\overline{\rho}` is specified in g cm :math:`^{-3}`. The parametrization +below. Note that in all the equations, :math:`p` is specified in millibars and :math:`\overline{\rho}` is specified in g cm :math:`^{-3}`. The parametrization of the source terms are given below. .. _`gabervsek2012dry`: https://journals.ametsoc.org/view/journals/mwre/140/10/mwr-d-11-00144.1.xml @@ -40,7 +40,7 @@ of the source terms are given below. Rate of condensation of water vapor/evaporation of cloud water ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -From `klemp1978simulation`_, we have the following expressions. +From `klemp1978simulation`_, we have the following expressions. .. _`klemp1978simulation`: https://journals.ametsoc.org/view/journals/atsc/35/6/1520-0469_1978_035_1070_tsotdc_2_0_co_2.xml From 933cd1a8598e51a723edb0f36985e24a85016585 Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Tue, 5 Dec 2023 15:36:16 -0800 Subject: [PATCH 35/37] Updating docs --- Docs/sphinx_doc/theory/Initialization.rst | 6 ------ 1 file changed, 6 deletions(-) diff --git a/Docs/sphinx_doc/theory/Initialization.rst b/Docs/sphinx_doc/theory/Initialization.rst index 9e732c550..10a7e1a50 100644 --- a/Docs/sphinx_doc/theory/Initialization.rst +++ b/Docs/sphinx_doc/theory/Initialization.rst @@ -147,9 +147,3 @@ where :math:`p_0 = 1e5 \, \text{N/m}^2` is the pressure at the base. Hence, we d F(p(0)) \equiv p(0) - p_0 + \rho(0)g\frac{\Delta z}{2}, and the Newton-Raphson procedure is the same. - - - - - - From 89b3846e1b5643ddf1d2cecdbc9a066d3fd2980e Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 6 Dec 2023 01:27:37 -0800 Subject: [PATCH 36/37] Updating docs --- Docs/sphinx_doc/theory/Buoyancy.rst | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index 92e987ab5..dbf289ac2 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -112,6 +112,11 @@ if we neglect :math:`\frac{p^\prime}{\bar{p_0}}`. Type 4 ------ +This expression for buoyancy is from `khairoutdinov2003cloud`_ and `bryan2002benchmark`_. + +.. _`khairoutdinov2003cloud`: https://journals.ametsoc.org/view/journals/atsc/60/4/1520-0469_2003_060_0607_crmota_2.0.co_2.xml +.. _`bryan2002benchmark`: https://journals.ametsoc.org/view/journals/mwre/130/12/1520-0493_2002_130_2917_absfmn_2.0.co_2.xml + .. math:: \begin{equation} From 59439a708494a34b86ba99a961dd6cb2fc9981ea Mon Sep 17 00:00:00 2001 From: Mahesh Natarajan Date: Wed, 6 Dec 2023 22:22:56 -0800 Subject: [PATCH 37/37] Correcting VisIt plot issue --- Source/Microphysics/Kessler/Kessler.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/Source/Microphysics/Kessler/Kessler.cpp b/Source/Microphysics/Kessler/Kessler.cpp index 529bb03e9..6c32a0ff2 100644 --- a/Source/Microphysics/Kessler/Kessler.cpp +++ b/Source/Microphysics/Kessler/Kessler.cpp @@ -197,8 +197,12 @@ void Kessler::AdvanceKessler() { } + if(std::fabs(fz_array(i,j,k+1)) < 1e-14) fz_array(i,j,k+1) = 0.0; + if(std::fabs(fz_array(i,j,k)) < 1e-14) fz_array(i,j,k) = 0.0; Real dq_sed = 1.0/rho_array(i,j,k)*(fz_array(i,j,k+1) - fz_array(i,j,k))/dz*dtn; - + if(std::fabs(dq_sed) < 1e-14)dq_sed = 0.0; + //dq_sed = 0.0; + qt_array(i,j,k) = qt_array(i,j,k) + dq_rain_to_vapor - dq_clwater_to_rain; qp_array(i,j,k) = qp_array(i,j,k) + dq_sed + dq_clwater_to_rain - dq_rain_to_vapor; qn_array(i,j,k) = qn_array(i,j,k) + dq_vapor_to_clwater - dq_clwater_to_vapor - dq_clwater_to_rain;