diff --git a/Source/BoundaryConditions/MOSTStress.H b/Source/BoundaryConditions/MOSTStress.H index 8cc6f27cc..2e4a45649 100644 --- a/Source/BoundaryConditions/MOSTStress.H +++ b/Source/BoundaryConditions/MOSTStress.H @@ -1174,7 +1174,7 @@ struct custom_flux const amrex::Array4& /*t_surf_arr*/, const amrex::Array4& dest_arr) const { - amrex::Real rho, eta, qv; + amrex::Real rho, eta; int ic, jc; ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i; @@ -1183,7 +1183,6 @@ struct custom_flux jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc; rho = cons_arr(ic,jc,zlo,Rho_comp); - qv = cons_arr(ic,jc,zlo,RhoQ1_comp) / rho; amrex::Real qstar = q_star_arr(ic,jc,zlo); amrex::Real moflux = (std::abs(qstar) > eps) ? -rho * qstar : 0.0; @@ -1220,12 +1219,12 @@ struct custom_flux const amrex::Array4& /*vely_arr*/, const amrex::Array4& /*umm_arr*/, const amrex::Array4& /*tm_arr*/, - const amrex::Array4& u_star_arr, + const amrex::Array4& /*u_star_arr*/, const amrex::Array4& t_star_arr, const amrex::Array4& /*t_surf_arr*/, const amrex::Array4& dest_arr) const { - amrex::Real rho, eta, theta; + amrex::Real rho, eta; int ic, jc; ic = i < lbound(cons_arr).x ? lbound(cons_arr).x : i; @@ -1234,9 +1233,7 @@ struct custom_flux jc = jc > ubound(cons_arr).y ? ubound(cons_arr).y : jc; rho = cons_arr(ic,jc,zlo,Rho_comp); - theta = cons_arr(ic,jc,zlo,RhoTheta_comp) / rho; - amrex::Real ustar = u_star_arr(ic,jc,zlo); amrex::Real tstar = t_star_arr(ic,jc,zlo); amrex::Real moflux = (std::abs(tstar) > eps) ? -rho * tstar : 0.0; amrex::Real deltaz = dz * (zlo - k); diff --git a/Source/Microphysics/SAM/Cloud_SAM.cpp b/Source/Microphysics/SAM/Cloud_SAM.cpp index 7783cb5bf..84d51549f 100644 --- a/Source/Microphysics/SAM/Cloud_SAM.cpp +++ b/Source/Microphysics/SAM/Cloud_SAM.cpp @@ -18,8 +18,6 @@ void SAM::Cloud (const SolverChoice& solverChoice) { Real fac_fus = m_fac_fus; Real rdOcp = m_rdOcp; - Real tol = 1.0e-4; - SolverChoice sc = solverChoice; int SAM_moisture_type = 1; @@ -47,16 +45,12 @@ void SAM::Cloud (const SolverChoice& solverChoice) { ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // Saturation moisture fractions - Real omn, domn; - Real qsat, dqsat; - Real qsatw, dqsatw; - Real qsati, dqsati; + Real omn; + Real qsat; + Real qsatw; + Real qsati; // Newton iteration vars - int niter; - Real fff, dfff, dtabs; - Real lstar, dlstar; - Real lstarw, lstari; Real delta_qv, delta_qc, delta_qi; // NOTE: Conversion before iterations is necessary to @@ -118,122 +112,40 @@ void SAM::Cloud (const SolverChoice& solverChoice) { pres_array(i,j,k) *= 0.01; } - // Initial guess for temperature & pressure - Real tabs = tabs_array(i,j,k); - Real pres = pres_array(i,j,k); - // Saturation moisture fractions - erf_qsatw(tabs, pres, qsatw); - erf_qsati(tabs, pres, qsati); + erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsatw); + erf_qsati(tabs_array(i,j,k), pres_array(i,j,k), qsati); qsat = omn * qsatw + (1.0-omn) * qsati; // We have enough total moisture to relax to equilibrium if (qt_array(i,j,k) > qsat) { - niter = 0; - dtabs = 1; - //================================================== - // Newton iteration to qv=qsat (cloud phase only) - //================================================== - do { - // Latent heats and their derivatives wrt to T - lstarw = fac_cond; - lstari = fac_fus; - domn = 0.0; - - // Saturation moisture fractions - erf_qsatw(tabs, pres, qsatw); - erf_qsati(tabs, pres, qsati); - erf_dtqsatw(tabs, pres, dqsatw); - erf_dtqsati(tabs, pres, dqsati); - - if (SAM_moisture_type == 1) { - // Cloud ice not permitted (condensation & fusion) - if(tabs >= tbgmax) { - omn = 1.0; - } - // Cloud water not permitted (sublimation & fusion) - else if(tabs <= tbgmin) { - omn = 0.0; - lstarw = fac_sub; - } - // Mixed cloud phase (condensation & fusion) - else { - omn = an*tabs-bn; - domn = an; - } - } else if (SAM_moisture_type == 2) { - omn = 1.0; - domn = 0.0; - } - - // Linear combination of each component - qsat = omn * qsatw + (1.0-omn) * qsati; - dqsat = omn * dqsatw + (1.0-omn) * dqsati - + domn * qsatw - domn * qsati; - lstar = omn * lstarw + (1.0-omn) * lstari; - dlstar = domn * lstarw - domn * lstari; - - // Function for root finding: - // 0 = -T_new + T_old + L_eff/C_p * (qv - qsat) - fff = -tabs + tabs_array(i,j,k) + lstar*(qv_array(i,j,k) - qsat); - - // Derivative of function (T_new iterated on) - dfff = -1.0 + dlstar*(qv_array(i,j,k) - qsat) - lstar*dqsat; - - // Update the temperature - dtabs = -fff/dfff; - tabs = tabs+dtabs; - - // For now at least we perform this iteration at constant pressure - // For the moist bubble case, the results are indistinguisable - // between running with this used vs commented out - // pres = rho_array(i,j,k) * R_d * tabs - // * (1.0 + R_v/R_d * qsat) * 0.01; - - // Update iteration - niter = niter+1; - } while(std::abs(dtabs) > tol && niter < 20); - - // Update qsat from last iteration (dq = dq/dt * dt) - qsat += dqsat*dtabs; - - // Changes in each component - delta_qv = qv_array(i,j,k) - qsat; - delta_qc = std::max(-qcl_array(i,j,k), delta_qv * omn); - delta_qi = std::max(-qci_array(i,j,k), delta_qv * (1.0-omn)); - - // Partition the change in non-precipitating q - qv_array(i,j,k) = qsat; - qcl_array(i,j,k) += delta_qc; - qci_array(i,j,k) += delta_qi; - qn_array(i,j,k) = qcl_array(i,j,k) + qci_array(i,j,k); - qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k); // Update temperature - tabs_array(i,j,k) = tabs; + tabs_array(i,j,k) = NewtonIterSat(i, j, k , SAM_moisture_type , + fac_cond , fac_fus , fac_sub , + an , bn , + tabs_array, pres_array, + qv_array , qcl_array , qci_array, + qn_array , qt_array); // Update theta from temperature (it is essential to do this BEFORE the pressure is updated) // This would be inconsistent with updating the pressure as part of the iteration above. // Empirically based on the moist bubble rise case, getting the correct theta here // depends on using the old (unchanged by the phase changes) pressure. - theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); // Update pressure to be consistent with updated theta_array pres_array(i,j,k) = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k)); - // This was used in the earlier implmentation when we updated theta using this new pressure - // pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) - // * (1.0 + R_v/R_d * qv_array(i,j,k)); - // Pressure unit conversion pres_array(i,j,k) *= 0.01; // - // We cannot blindly relax to qsat, but we can convert qc/qi -> qv - // The concept here is that if we assume the temperature change due to this conversion - // doesn't change qsat, we can safely put all the moisture into qv without reaching qv > qsat - // because we are in the "else" part of the test on (qt_array(i,j,k) > qsat) + // We cannot blindly relax to qsat, but we can convert qc/qi -> qv. + // The concept here is that if we put all the moisture into qv and modify + // the temperature, we can then check if qv > qsat occurs (for final T/P/qv). + // If the reduction in T/qsat and increase in qv does trigger the + // aforementioned condition, we can do Newton iteration to drive qv = qsat. // } else { // Changes in each component @@ -251,15 +163,45 @@ void SAM::Cloud (const SolverChoice& solverChoice) { // Update temperature (endothermic since we evap/sublime) tabs_array(i,j,k) -= fac_cond * delta_qc + fac_sub * delta_qi; - // Update pressure - pres_array(i,j,k) = rho_array(i,j,k) * R_d * tabs_array(i,j,k) - * (1.0 + R_v/R_d * qv_array(i,j,k)); + // Update theta from temperature (it is essential to do this BEFORE the pressure is updated) + // This would be inconsistent with updating the pressure as part of the iteration above. + // Empirically based on the moist bubble rise case, getting the correct theta here + // depends on using the old (unchanged by the phase changes) pressure. + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); - // Update theta from temperature - theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp); + // Update pressure to be consistent with updated theta_array + pres_array(i,j,k) = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k)); // Pressure unit conversion pres_array(i,j,k) *= 0.01; + + // Verify assumption that qv > qsat does not occur + erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsatw); + erf_qsati(tabs_array(i,j,k), pres_array(i,j,k), qsati); + qsat = omn * qsatw + (1.0-omn) * qsati; + if (qt_array(i,j,k) > qsat) { + + // Update temperature + tabs_array(i,j,k) = NewtonIterSat(i, j, k , SAM_moisture_type , + fac_cond , fac_fus , fac_sub , + an , bn , + tabs_array, pres_array, + qv_array , qcl_array , qci_array, + qn_array , qt_array); + + // Update theta from temperature (it is essential to do this BEFORE the pressure is updated) + // This would be inconsistent with updating the pressure as part of the iteration above. + // Empirically based on the moist bubble rise case, getting the correct theta here + // depends on using the old (unchanged by the phase changes) pressure. + theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp); + + // Update pressure to be consistent with updated theta_array + pres_array(i,j,k) = getPgivenRTh(rho_array(i,j,k)*theta_array(i,j,k),qv_array(i,j,k)); + + // Pressure unit conversion + pres_array(i,j,k) *= 0.01; + + } } }); } // mfi diff --git a/Source/Microphysics/SAM/SAM.H b/Source/Microphysics/SAM/SAM.H index 3b4d9f280..8a88866a6 100644 --- a/Source/Microphysics/SAM/SAM.H +++ b/Source/Microphysics/SAM/SAM.H @@ -149,6 +149,122 @@ public: int Qstate_Size () override { return SAM::m_qstate_size; } + AMREX_GPU_HOST_DEVICE + AMREX_FORCE_INLINE + static amrex::Real + NewtonIterSat (int& i, int& j, int& k, + const int& SAM_moisture_type, + const amrex::Real& fac_cond, + const amrex::Real& fac_fus, + const amrex::Real& fac_sub, + const amrex::Real& an, + const amrex::Real& bn, + const amrex::Array4& tabs_array, + const amrex::Array4& pres_array, + const amrex::Array4& qv_array, + const amrex::Array4& qc_array, + const amrex::Array4& qi_array, + const amrex::Array4& qn_array, + const amrex::Array4& qt_array) + { + // Solution tolerance + amrex::Real tol = 1.0e-4; + + // Saturation moisture fractions + amrex::Real omn, domn; + amrex::Real qsat, dqsat; + amrex::Real qsatw, dqsatw; + amrex::Real qsati, dqsati; + + // Newton iteration vars + int niter; + amrex::Real fff, dfff, dtabs; + amrex::Real lstar, dlstar; + amrex::Real lstarw, lstari; + amrex::Real delta_qv, delta_qc, delta_qi; + + // Initial guess for temperature & pressure + amrex::Real tabs = tabs_array(i,j,k); + amrex::Real pres = pres_array(i,j,k); + + niter = 0; + dtabs = 1; + //================================================== + // Newton iteration to qv=qsat (cloud phase only) + //================================================== + do { + // Latent heats and their derivatives wrt to T + lstarw = fac_cond; + lstari = fac_fus; + domn = 0.0; + + // Saturation moisture fractions + erf_qsatw(tabs, pres, qsatw); + erf_qsati(tabs, pres, qsati); + erf_dtqsatw(tabs, pres, dqsatw); + erf_dtqsati(tabs, pres, dqsati); + + if (SAM_moisture_type == 1) { + // Cloud ice not permitted (condensation & fusion) + if(tabs >= tbgmax) { + omn = 1.0; + } + // Cloud water not permitted (sublimation & fusion) + else if(tabs <= tbgmin) { + omn = 0.0; + lstarw = fac_sub; + } + // Mixed cloud phase (condensation & fusion) + else { + omn = an*tabs-bn; + domn = an; + } + } else if (SAM_moisture_type == 2) { + omn = 1.0; + domn = 0.0; + } + + // Linear combination of each component + qsat = omn * qsatw + (1.0-omn) * qsati; + dqsat = omn * dqsatw + (1.0-omn) * dqsati + + domn * qsatw - domn * qsati; + lstar = omn * lstarw + (1.0-omn) * lstari; + dlstar = domn * lstarw - domn * lstari; + + // Function for root finding: + // 0 = -T_new + T_old + L_eff/C_p * (qv - qsat) + fff = -tabs + tabs_array(i,j,k) + lstar*(qv_array(i,j,k) - qsat); + + // Derivative of function (T_new iterated on) + dfff = -1.0 + dlstar*(qv_array(i,j,k) - qsat) - lstar*dqsat; + + // Update the temperature + dtabs = -fff/dfff; + tabs += dtabs; + + // Update iteration + niter = niter+1; + } while(std::abs(dtabs) > tol && niter < 20); + + // Update qsat from last iteration (dq = dq/dt * dt) + qsat += dqsat*dtabs; + + // Changes in each component + delta_qv = qv_array(i,j,k) - qsat; + delta_qc = std::max(-qc_array(i,j,k), delta_qv * omn); + delta_qi = std::max(-qi_array(i,j,k), delta_qv * (1.0-omn)); + + // Partition the change in non-precipitating q + qv_array(i,j,k) = qsat; + qc_array(i,j,k) += delta_qc; + qi_array(i,j,k) += delta_qi; + qn_array(i,j,k) = qc_array(i,j,k) + qi_array(i,j,k); + qt_array(i,j,k) = qv_array(i,j,k) + qn_array(i,j,k); + + // Return to temperature + return tabs; + } + private: // Number of qmoist variables (qt, qv, qcl, qci, qp, qpr, qps, qpg) int m_qmoist_size = 11; diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index 6b33f0e9f..059923b56 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -202,12 +202,6 @@ void erf_slow_rhs_pre (int level, int finest_level, // Define updates and fluxes in the current RK stage // ***************************************************************************** - // Open bc will be imposed upon all vars (we only access cons here for simplicity) - const bool xlo_open = (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::open); - const bool xhi_open = (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::open); - const bool ylo_open = (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::open); - const bool yhi_open = (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::open); - #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif