ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERF_MakeBuoyancy.cpp File Reference
#include <AMReX_MultiFab.H>
#include <AMReX_ArrayLim.H>
#include <AMReX_GpuContainers.H>
#include <ERF_Constants.H>
#include <ERF_EOS.H>
#include <ERF_IndexDefines.H>
#include <ERF_SrcHeaders.H>
#include <ERF_BuoyancyUtils.H>
#include <ERF_EB.H>
Include dependency graph for ERF_MakeBuoyancy.cpp:

Functions

void make_buoyancy (int lev, const Vector< MultiFab > &S_data, const MultiFab &S_prim, const MultiFab &qt, MultiFab &buoyancy, const Geometry geom, const SolverChoice &solverChoice, const MultiFab &base_state, const int n_qstate, const eb_ &ebfact, const int anelastic)
 

Function Documentation

◆ make_buoyancy()

void make_buoyancy ( int  lev,
const Vector< MultiFab > &  S_data,
const MultiFab &  S_prim,
const MultiFab &  qt,
MultiFab &  buoyancy,
const Geometry  geom,
const SolverChoice solverChoice,
const MultiFab &  base_state,
const int  n_qstate,
const eb_ ebfact,
const int  anelastic 
)

Function for computing the buoyancy term to be used in the evolution equation for the z-component of momentum in the slow integrator. There are three options for how buoyancy is computed (two are the same in the absence of moisture).

Parameters
[in]S_datacurrent solution
[in]S_primprimitive variables (i.e. conserved variables divided by density)
[out]buoyancybuoyancy term computed here
[in]qmoistmoisture variables (in order: qv, qc, qi, ...)
[in]qv_dlateral average of cloud vapor
[in]qc_dlateral average of cloud vapor
[in]qd_dlateral average of cloud vapor
[in]geomContainer for geometric information
[in]solverChoiceContainer for solver parameters
[in]r0Reference (hydrostatically stratified) density
[in]n_qstateNumber of moist variables used by the current model
43 {
44  BL_PROFILE("make_buoyancy()");
45 
46  const Array<Real,AMREX_SPACEDIM> grav{0.0, 0.0, -solverChoice.gravity};
47  const GpuArray<Real,AMREX_SPACEDIM> grav_gpu{grav[0], grav[1], grav[2]};
48 
49  const int klo = geom.Domain().smallEnd()[2];
50  const int khi = geom.Domain().bigEnd()[2] + 1;
51 
52  Real rd_over_cp = solverChoice.rdOcp;
53  //Real rv_over_rd = R_v/R_d;
54 
55  MultiFab r0 (base_state, make_alias, BaseState::r0_comp , 1);
56  MultiFab p0 (base_state, make_alias, BaseState::p0_comp , 1);
57  MultiFab th0(base_state, make_alias, BaseState::th0_comp, 1);
58  MultiFab qv0(base_state, make_alias, BaseState::qv0_comp, 1);
59 
60 #ifdef _OPENMP
61 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
62 #endif
63  for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi)
64  {
65  Box tbz = mfi.tilebox();
66 
67  // We don't compute a source term for z-momentum on the bottom or top boundary
68  if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1);
69  if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1);
70 
71  const Array4<const Real> & cell_data = S_data[IntVars::cons].array(mfi);
72  const Array4<const Real> & cell_prim = S_prim.array(mfi);
73  const Array4<const Real> & qt_arr = qt.array(mfi);
74  const Array4< Real> & buoyancy_fab = buoyancy.array(mfi);
75 
76  // Base state density and pressure
77  const Array4<const Real>& r0_arr = r0.const_array(mfi);
78  const Array4<const Real>& p0_arr = p0.const_array(mfi);
79  const Array4<const Real>& th0_arr = th0.const_array(mfi);
80  const Array4<const Real>& qv0_arr = qv0.const_array(mfi);
81 
82  if (solverChoice.terrain_type != TerrainType::EB) {
83 
84  if ( anelastic && (solverChoice.moisture_type == MoistureType::None) )
85  {
86  // ******************************************************************************************
87  // Dry anelastic
88  // ******************************************************************************************
89  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
90  {
91  //
92  // Return -rho0 g (thetaprime / theta0)
93  //
94  buoyancy_fab(i, j, k) = buoyancy_dry_anelastic(i,j,k,grav_gpu[2],
95  r0_arr,th0_arr,cell_data);
96  });
97  }
98  else if ( anelastic && (solverChoice.moisture_type != MoistureType::None) )
99  {
100  // ******************************************************************************************
101  // Moist anelastic
102  // ******************************************************************************************
103  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
104  {
105  //
106  // Return -rho0 g (thetaprime / theta0)
107  //
108  //buoyancy_fab(i, j, k) = buoyancy_moist_anelastic(i,j,k,grav_gpu[2],rv_over_rd,
109  // r0_arr,th0_arr,qv0_arr,cell_data,qt_arr);
110 
111  // NOTE: Using the type 4, which we formally derived.
112  // The above has errors and needs rederiving.
113  buoyancy_fab(i, j, k) = buoyancy_moist_Thpert(i,j,k,n_qstate,grav_gpu[2],
114  r0_arr,th0_arr,qv0_arr,cell_prim,qt_arr);
115  });
116  }
117  else if ( !anelastic && (solverChoice.moisture_type == MoistureType::None) )
118  {
119  // ******************************************************************************************
120  // Dry compressible
121  // ******************************************************************************************
122  if (solverChoice.buoyancy_type[lev] == 1) {
123 
124  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
125  {
126  //
127  // Return -rho0 g (thetaprime / theta0)
128  //
129  buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,grav_gpu[2],
130  r0_arr,cell_data,qt_arr);
131  });
132  }
133  else if (solverChoice.buoyancy_type[lev] == 2 || solverChoice.buoyancy_type[lev] == 3)
134  {
135  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
136  {
137  //
138  // Return -rho0 g (Tprime / T0)
139  //
140  buoyancy_fab(i, j, k) = buoyancy_dry_Tpert(i,j,k,grav_gpu[2],rd_over_cp,
141  r0_arr,p0_arr,th0_arr,cell_data);
142  });
143  }
144  else if (solverChoice.buoyancy_type[lev] == 4)
145  {
146  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
147  {
148  //
149  // Return -rho0 g (Theta_prime / Theta_0)
150  //
151  buoyancy_fab(i, j, k) = buoyancy_dry_Thpert(i,j,k,grav_gpu[2],
152  r0_arr,th0_arr,cell_data);
153  });
154  } // buoyancy_type for dry compressible
155  }
156  else // if ( !anelastic && (solverChoice.moisture_type != MoistureType::None) )
157  {
158  // ******************************************************************************************
159  // Moist compressible
160  // ******************************************************************************************
161 
162  if ( (solverChoice.moisture_type == MoistureType::Kessler_NoRain) ||
163  (solverChoice.moisture_type == MoistureType::SAM) ||
164  (solverChoice.moisture_type == MoistureType::SAM_NoPrecip_NoIce) )
165  {
166  AMREX_ALWAYS_ASSERT(solverChoice.buoyancy_type[lev] == 1);
167  }
168 
169  if (solverChoice.buoyancy_type[lev] == 1)
170  {
171  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
172  {
173  buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,grav_gpu[2],
174  r0_arr,cell_data,qt_arr);
175  });
176  }
177  else if (solverChoice.buoyancy_type[lev] == 2 || solverChoice.buoyancy_type[lev] == 3)
178  {
179 
180  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
181  {
182  buoyancy_fab(i, j, k) = buoyancy_moist_Tpert(i,j,k,n_qstate,grav_gpu[2],rd_over_cp,
183  r0_arr,th0_arr,qv0_arr,p0_arr,
184  cell_prim,cell_data,qt_arr);
185  });
186  }
187  else if (solverChoice.buoyancy_type[lev] == 4)
188  {
189  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
190  {
191  buoyancy_fab(i, j, k) = buoyancy_moist_Thpert(i,j,k,n_qstate,grav_gpu[2],
192  r0_arr,th0_arr,qv0_arr,cell_prim,qt_arr);
193  });
194  }
195  } // moist compressible
196 
197  } else {
198 
199  // Currently, only dry compressible is supported
200  AMREX_ASSERT( !anelastic && (solverChoice.moisture_type == MoistureType::None) && solverChoice.buoyancy_type[lev] == 1 );
201 
202  Array4<const EBCellFlag> cellflg = (ebfact.get_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
203 
204  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
205  {
206  buoyancy_fab(i, j, k) = buoyancy_rhopert_eb(i,j,k,grav_gpu[2],r0_arr,cell_data,qt_arr,cellflg);
207  });
208  } // TerrainType::EB
209  } // mfi
210 }
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_moist_Tpert(int &i, int &j, int &k, const int &n_qstate, amrex::Real const &grav_gpu, amrex::Real const &rd_over_cp, const amrex::Array4< const amrex::Real > &r0_arr, const amrex::Array4< const amrex::Real > &th0_arr, const amrex::Array4< const amrex::Real > &qv0_arr, const amrex::Array4< const amrex::Real > &p0_arr, const amrex::Array4< const amrex::Real > &cell_prim, const amrex::Array4< const amrex::Real > &cell_data, const amrex::Array4< const amrex::Real > &qt_arr)
Definition: ERF_BuoyancyUtils.H:165
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_moist_Thpert(int &i, int &j, int &k, const int &n_qstate, amrex::Real const &grav_gpu, const amrex::Array4< const amrex::Real > &r0_arr, const amrex::Array4< const amrex::Real > &th0_arr, const amrex::Array4< const amrex::Real > &qv0_arr, const amrex::Array4< const amrex::Real > &cell_prim, const amrex::Array4< const amrex::Real > &qt_arr)
Definition: ERF_BuoyancyUtils.H:202
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_dry_anelastic(int &i, int &j, int &k, amrex::Real const &grav_gpu, const amrex::Array4< const amrex::Real > &r0_arr, const amrex::Array4< const amrex::Real > &th0_arr, const amrex::Array4< const amrex::Real > &cell_data)
Definition: ERF_BuoyancyUtils.H:10
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_dry_Thpert(int &i, int &j, int &k, amrex::Real const &grav_gpu, const amrex::Array4< const amrex::Real > &r0_arr, const amrex::Array4< const amrex::Real > &th0_arr, const amrex::Array4< const amrex::Real > &cell_prim)
Definition: ERF_BuoyancyUtils.H:147
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_rhopert_eb(int &i, int &j, int &k, amrex::Real const &grav_gpu, const amrex::Array4< const amrex::Real > &r0_arr, const amrex::Array4< const amrex::Real > &cell_data, const amrex::Array4< const amrex::Real > &qt_arr, amrex::Array4< amrex::EBCellFlag const > const &flag)
Definition: ERF_BuoyancyUtils.H:86
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_dry_Tpert(int &i, int &j, int &k, amrex::Real const &grav_gpu, amrex::Real const &rd_over_cp, const amrex::Array4< const amrex::Real > &r0_arr, const amrex::Array4< const amrex::Real > &p0_arr, const amrex::Array4< const amrex::Real > &th0_arr, const amrex::Array4< const amrex::Real > &cell_data)
Definition: ERF_BuoyancyUtils.H:121
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_rhopert(int &i, int &j, int &k, amrex::Real const &grav_gpu, const amrex::Array4< const amrex::Real > &r0_arr, const amrex::Array4< const amrex::Real > &cell_data, const amrex::Array4< const amrex::Real > &qt_arr)
Definition: ERF_BuoyancyUtils.H:72
amrex::Real Real
Definition: ERF_ShocInterface.H:19
const std::unique_ptr< amrex::EBFArrayBoxFactory > & get_const_factory() const noexcept
Definition: ERF_EB.H:46
@ qv0_comp
Definition: ERF_IndexDefines.H:67
@ p0_comp
Definition: ERF_IndexDefines.H:64
@ th0_comp
Definition: ERF_IndexDefines.H:66
@ r0_comp
Definition: ERF_IndexDefines.H:63
@ cons
Definition: ERF_IndexDefines.H:158
@ qt
Definition: ERF_Kessler.H:27
real(c_double), parameter p0
Definition: ERF_module_model_constants.F90:40
amrex::Real rdOcp
Definition: ERF_DataStruct.H:1058
amrex::Vector< int > buoyancy_type
Definition: ERF_DataStruct.H:1003
amrex::Real gravity
Definition: ERF_DataStruct.H:1056
MoistureType moisture_type
Definition: ERF_DataStruct.H:1101
static TerrainType terrain_type
Definition: ERF_DataStruct.H:971
Here is the call graph for this function: