From 7b812deee4221ed9c7a692ace5302989c21ce9cf Mon Sep 17 00:00:00 2001 From: Eliot Quon Date: Wed, 31 Jul 2024 22:45:45 -0600 Subject: [PATCH] Generalize geostrophic forcing (#1718) * Add coriolis_3d flag (default: true) * Remove hard-coded coriolis factor when using custom geo wind profiles --- Source/DataStructs/DataStruct.H | 9 +++++++-- Source/SourceTerms/ERF_make_mom_sources.cpp | 12 +++++------- 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H index e2ee13a6e..4f9f68a63 100644 --- a/Source/DataStructs/DataStruct.H +++ b/Source/DataStructs/DataStruct.H @@ -430,10 +430,14 @@ struct SolverChoice { amrex::Real latitude = 90.0; pp.query("latitude", latitude); + pp.query("coriolis_3d", coriolis_3d); + // Convert to radians latitude *= (PI/180.); sinphi = std::sin(latitude); - cosphi = std::cos(latitude); + if (coriolis_3d) { + cosphi = std::cos(latitude); + } amrex::Print() << "Coriolis frequency, f = " << coriolis_factor * sinphi << " 1/s" << std::endl; @@ -478,6 +482,7 @@ struct SolverChoice { // Specify what additional physics/forcing modules we use bool use_gravity = false; bool use_coriolis = false; + bool coriolis_3d = true; bool rayleigh_damp_U = false; bool rayleigh_damp_V = false; @@ -503,7 +508,7 @@ struct SolverChoice { // Coriolis forcing amrex::Real coriolis_factor = 0.0; - amrex::Real cosphi = 0.0 ; + amrex::Real cosphi = 0.0; amrex::Real sinphi = 0.0; // User-specified forcings in problem definition diff --git a/Source/SourceTerms/ERF_make_mom_sources.cpp b/Source/SourceTerms/ERF_make_mom_sources.cpp index b51d33f69..9c6e574d9 100644 --- a/Source/SourceTerms/ERF_make_mom_sources.cpp +++ b/Source/SourceTerms/ERF_make_mom_sources.cpp @@ -93,7 +93,7 @@ void make_mom_sources (int /*level*/, // Flag for Geostrophic forcing // ***************************************************************************** auto abl_geo_forcing = solverChoice.abl_geo_forcing; - auto geostrophic_wind = solverChoice.custom_geostrophic_profile; + auto geo_wind_profile = solverChoice.custom_geostrophic_profile; // ***************************************************************************** // Data for Rayleigh damping @@ -286,20 +286,18 @@ void make_mom_sources (int /*level*/, // ***************************************************************************** // Add height-dependent GEOSTROPHIC forcing // ***************************************************************************** - if (geostrophic_wind) { + if (geo_wind_profile) { ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) ); - Real rho_v_loc = 0.25 * (rho_v(i,j+1,k) + rho_v(i,j,k) + rho_v(i-1,j+1,k) + rho_v(i-1,j,k)); - xmom_src_arr(i, j, k) += -0.376E-4 * (rho_on_u_face * dptr_v_geos[k] - rho_v_loc); + xmom_src_arr(i, j, k) -= coriolis_factor * rho_on_u_face * dptr_v_geos[k] * sinphi; }); ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) { Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) ); - Real rho_u_loc = 0.25 * (rho_u(i+1,j,k) + rho_u(i,j,k) + rho_u(i+1,j-1,k) + rho_u(i,j-1,k)); - ymom_src_arr(i, j, k) += 0.376E-4 * (rho_on_v_face * dptr_u_geos[k] - rho_u_loc); + ymom_src_arr(i, j, k) += coriolis_factor * rho_on_v_face * dptr_u_geos[k] * sinphi; }); - } // geostrophic_wind + } // geo_wind_profile // ***************************************************************************** // Add custom SUBSIDENCE terms