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]levlevel
[in]S_datacurrent solution
[in]S_primprimitive variables (i.e. conserved variables divided by density)
[in]qttotal water in the state
[out]buoyancybuoyancy term computed here
[in]geomContainer for geometric information
[in]solverChoiceContainer for solver parameters
[in]base_statebase state
[in]n_qstateNumber of moist variables used by the current model
[in]ebfactContainer of EB information
[in]anelasticAre we solving the anelastic equations (1 if yes, 0 if no)
43 {
44  BL_PROFILE("make_buoyancy()");
45 
46  const Array<Real,AMREX_SPACEDIM> grav{zero, zero, -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  MultiFab r0 (base_state, make_alias, BaseState::r0_comp , 1);
53  MultiFab p0 (base_state, make_alias, BaseState::p0_comp , 1);
54  MultiFab th0(base_state, make_alias, BaseState::th0_comp, 1);
55  MultiFab qv0(base_state, make_alias, BaseState::qv0_comp, 1);
56 
57 #ifdef _OPENMP
58 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
59 #endif
60  for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi)
61  {
62  Box tbz = mfi.tilebox();
63 
64  // We don't compute a source term for z-momentum on the bottom or top boundary
65  if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1);
66  if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1);
67 
68  const Array4<const Real> & cell_data = S_data[IntVars::cons].array(mfi);
69  const Array4<const Real> & cell_prim = S_prim.array(mfi);
70  const Array4<const Real> & qt_arr = qt.array(mfi);
71  const Array4< Real> & buoyancy_fab = buoyancy.array(mfi);
72 
73  // Base state density and pressure
74  const Array4<const Real>& r0_arr = r0.const_array(mfi);
75  const Array4<const Real>& p0_arr = p0.const_array(mfi);
76  const Array4<const Real>& th0_arr = th0.const_array(mfi);
77  const Array4<const Real>& qv0_arr = qv0.const_array(mfi);
78 
79  if (solverChoice.terrain_type != TerrainType::EB) {
80 
81  if ( anelastic && (solverChoice.moisture_type == MoistureType::None) )
82  {
83  // ******************************************************************************************
84  // Dry anelastic
85  // ******************************************************************************************
86  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
87  {
88  //
89  // Return -rho0 g (thetaprime / theta0)
90  //
91  buoyancy_fab(i, j, k) = buoyancy_dry_anelastic(i,j,k,grav_gpu[2],
92  r0_arr,th0_arr,cell_data);
93  });
94  }
95  else if ( anelastic && (solverChoice.moisture_type != MoistureType::None) )
96  {
97  // ******************************************************************************************
98  // Moist anelastic
99  // ******************************************************************************************
100  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
101  {
102  //
103  // Return -rho0 g (thetaprime / theta0)
104  //
105  //buoyancy_fab(i, j, k) = buoyancy_moist_anelastic(i,j,k,grav_gpu[2],RvoRd_d,
106  // r0_arr,th0_arr,qv0_arr,cell_data,qt_arr);
107 
108  // NOTE: Using the type 4, which we formally derived.
109  // The above has errors and needs rederiving.
110  buoyancy_fab(i, j, k) = buoyancy_moist_Thpert(i,j,k,n_qstate,grav_gpu[2],
111  r0_arr,th0_arr,qv0_arr,cell_prim,qt_arr);
112  });
113  }
114  else if ( !anelastic && (solverChoice.moisture_type == MoistureType::None) )
115  {
116  // ******************************************************************************************
117  // Dry compressible
118  // ******************************************************************************************
119  if (solverChoice.buoyancy_type[lev] == 1) {
120 
121  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
122  {
123  //
124  // Return -rho0 g (thetaprime / theta0)
125  //
126  buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,grav_gpu[2],
127  r0_arr,qv0_arr,cell_data,qt_arr);
128  });
129  }
130  else if (solverChoice.buoyancy_type[lev] == 2 || solverChoice.buoyancy_type[lev] == 3)
131  {
132  ParallelFor(tbz, [=,rdOcp_d=solverChoice.rdOcp] AMREX_GPU_DEVICE (int i, int j, int k)
133  {
134  //
135  // Return -rho0 g (Tprime / T0)
136  //
137  buoyancy_fab(i, j, k) = buoyancy_dry_Tpert(i,j,k,grav_gpu[2],rdOcp_d,
138  r0_arr,p0_arr,th0_arr,cell_data);
139  });
140  }
141  else if (solverChoice.buoyancy_type[lev] == 4)
142  {
143  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
144  {
145  //
146  // Return -rho0 g (Theta_prime / Theta_0)
147  //
148  buoyancy_fab(i, j, k) = buoyancy_dry_Thpert(i,j,k,grav_gpu[2],
149  r0_arr,th0_arr,cell_data);
150  });
151  } // buoyancy_type for dry compressible
152  }
153  else // if ( !anelastic && (solverChoice.moisture_type != MoistureType::None) )
154  {
155  // ******************************************************************************************
156  // Moist compressible
157  // ******************************************************************************************
158 
159  if ( (solverChoice.moisture_type == MoistureType::Kessler_NoRain) ||
160  (solverChoice.moisture_type == MoistureType::SAM) ||
161  (solverChoice.moisture_type == MoistureType::SAM_NoPrecip_NoIce) )
162  {
163  AMREX_ALWAYS_ASSERT(solverChoice.buoyancy_type[lev] == 1);
164  }
165 
166  if (solverChoice.buoyancy_type[lev] == 1)
167  {
168  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
169  {
170  buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,grav_gpu[2],
171  r0_arr,qv0_arr,cell_data,qt_arr);
172  });
173  }
174  else if (solverChoice.buoyancy_type[lev] == 2 || solverChoice.buoyancy_type[lev] == 3)
175  {
176 
177  ParallelFor(tbz, [=,RdoCp_d=RdoCp] AMREX_GPU_DEVICE (int i, int j, int k)
178  {
179  buoyancy_fab(i, j, k) = buoyancy_moist_Tpert(i,j,k,n_qstate,grav_gpu[2],RdoCp_d,
180  r0_arr,th0_arr,qv0_arr,p0_arr,
181  cell_prim,cell_data,qt_arr);
182  });
183  }
184  else if (solverChoice.buoyancy_type[lev] == 4)
185  {
186  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
187  {
188  buoyancy_fab(i, j, k) = buoyancy_moist_Thpert(i,j,k,n_qstate,grav_gpu[2],
189  r0_arr,th0_arr,qv0_arr,cell_prim,qt_arr);
190  });
191  }
192  } // moist compressible
193 
194  } else {
195 
196  if ( anelastic && (solverChoice.moisture_type == MoistureType::None) ) {
197 
198  if (grav_gpu[2]==0) {
199  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
200  {
201  buoyancy_fab(i, j, k) = zero;
202  });
203  } else {
204 
205  Array4<const EBCellFlag> cellflg = (ebfact.get_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
206 
207  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
208  {
209  buoyancy_fab(i, j, k) = buoyancy_dry_anelastic_eb(i,j,k,grav_gpu[2],
210  r0_arr,th0_arr,cell_data,cellflg);
211  });
212  }
213  }
214  else
215  {
216  if (grav_gpu[2]==0) {
217  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
218  {
219  buoyancy_fab(i, j, k) = zero;
220  });
221  } else {
222  // Currently, only dry compressible is supported
223  AMREX_ASSERT( !anelastic && (solverChoice.moisture_type == MoistureType::None) && solverChoice.buoyancy_type[lev] == 1 );
224 
225  Array4<const EBCellFlag> cellflg = (ebfact.get_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
226 
227  ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
228  {
229  buoyancy_fab(i, j, k) = buoyancy_rhopert_eb(i,j,k,grav_gpu[2],
230  r0_arr,qv0_arr,cell_data,qt_arr,cellflg);
231  });
232  }
233  }
234  } // TerrainType::EB
235  } // mfi
236 }
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_dry_Thpert(int &i, int &j, int &k, const amrex::Real &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:203
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_dry_Tpert(int &i, int &j, int &k, const amrex::Real &grav_gpu, const amrex::Real &, 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:177
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_rhopert(int &i, int &j, int &k, const amrex::Real &grav_gpu, const amrex::Array4< const amrex::Real > &r0_arr, const amrex::Array4< const amrex::Real > &qv0_arr, const amrex::Array4< const amrex::Real > &cell_data, const amrex::Array4< const amrex::Real > &qt_arr)
Definition: ERF_BuoyancyUtils.H:116
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_moist_Tpert(int &i, int &j, int &k, const int &n_qstate, const amrex::Real &grav_gpu, const amrex::Real &, 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:221
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_moist_Thpert(int &i, int &j, int &k, const int &n_qstate, const amrex::Real &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:258
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_dry_anelastic_eb(int &i, int &j, int &k, const amrex::Real &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, amrex::Array4< amrex::EBCellFlag const > const &flag)
Definition: ERF_BuoyancyUtils.H:30
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_rhopert_eb(int &i, int &j, int &k, const amrex::Real &grav_gpu, const amrex::Array4< const amrex::Real > &r0_arr, const amrex::Array4< const amrex::Real > &qv0_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:133
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real buoyancy_dry_anelastic(int &i, int &j, int &k, const amrex::Real &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
constexpr amrex::Real zero
Definition: ERF_Constants.H:8
constexpr amrex::Real RdoCp
Definition: ERF_Constants.H:54
const int khi
Definition: ERF_InitCustomPert_Bubble.H:21
AMREX_ALWAYS_ASSERT(bx.length()[2]==khi+1)
ParallelFor(grown_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) { qrcuten_arr(i, j, k)=Real(0);qscuten_arr(i, j, k)=Real(0);qicuten_arr(i, j, k)=Real(0);})
const std::unique_ptr< amrex::EBFArrayBoxFactory > & get_const_factory() const noexcept
Definition: ERF_EB.H:46
@ qv0_comp
Definition: ERF_IndexDefines.H:77
@ p0_comp
Definition: ERF_IndexDefines.H:74
@ th0_comp
Definition: ERF_IndexDefines.H:76
@ r0_comp
Definition: ERF_IndexDefines.H:73
@ cons
Definition: ERF_IndexDefines.H:193
@ qt
Definition: ERF_Kessler.H:29
real(c_double), parameter p0
Definition: ERF_module_model_constants.F90:40
real(kind=kind_phys), parameter, private r0
Definition: ERF_module_mp_wsm6.F90:21
amrex::Real rdOcp
Definition: ERF_DataStruct.H:1329
amrex::Vector< int > buoyancy_type
Definition: ERF_DataStruct.H:1260
amrex::Real gravity
Definition: ERF_DataStruct.H:1327
MoistureType moisture_type
Definition: ERF_DataStruct.H:1424
static TerrainType terrain_type
Definition: ERF_DataStruct.H:1230

Referenced by ERF::Write3DPlotFile().

Here is the call graph for this function:
Here is the caller graph for this function: