ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_VelocityToMomentum.cpp File Reference
#include <AMReX.H>
#include <AMReX_MultiFab.H>
#include <ERF_Utils.H>
#include <AMReX_MultiFabUtil.H>
Include dependency graph for ERF_VelocityToMomentum.cpp:

Functions

void VelocityToMomentum (const MultiFab &xvel_in, const IntVect &xvel_ngrow, const MultiFab &yvel_in, const IntVect &yvel_ngrow, const MultiFab &zvel_in, const IntVect &zvel_ngrow, const MultiFab &density, MultiFab &xmom, MultiFab &ymom, MultiFab &zmom, const Box &domain, const Vector< BCRec > &domain_bcs_type_h, const MultiFab *c_vfrac)
 

Function Documentation

◆ VelocityToMomentum()

void VelocityToMomentum ( const MultiFab &  xvel_in,
const IntVect &  xvel_ngrow,
const MultiFab &  yvel_in,
const IntVect &  yvel_ngrow,
const MultiFab &  zvel_in,
const IntVect &  zvel_ngrow,
const MultiFab &  density,
MultiFab &  xmom,
MultiFab &  ymom,
MultiFab &  zmom,
const Box &  domain,
const Vector< BCRec > &  domain_bcs_type_h,
const MultiFab *  c_vfrac 
)

Convert velocity to momentum by multiplying by density averaged onto faces.

Parameters
[in]xvel_inx-component of velocity
[in]xvel_ngrowhow many cells to grow the tilebox for the x-momentum
[in]yvel_iny-component of velocity
[in]yvel_ngrowhow many cells to grow the tilebox for the y-momentum
[in]zvel_inz-component of velocity
[in]zvel_ngrowhow many cells to grow the tilebox for the z-momentum
[in]densitydensity at cell centers
[out]xmomx-component of momentum
[out]ymomy-component of momentum
[out]zmomz-component of momentum
[in]domainDomain at this level
[in]domain_bcs_type_hhost vector for domain boundary conditions
[in]l_use_ndiffflag describing whether we will later add explicit numerical diffusion
[in]u_vfracu-face volume fraction (optional, only needed for EB)
[in]v_vfracv-face volume fraction (optional, only needed for EB)
[in]w_vfracw-face volume fraction (optional, only needed for EB)
42  {
43  BL_PROFILE_VAR("VelocityToMomentum()",VelocityToMomentum);
44 
45  const BCRec* bc_ptr_h = domain_bcs_type_h.data();
46 
47  // Loop over boxes = valid boxes grown by ngrow
48 #ifdef _OPENMP
49 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
50 #endif
51  for ( MFIter mfi(density,TilingIfNotGPU()); mfi.isValid(); ++mfi)
52  {
53  // We need momentum in the interior ghost cells (init == real)
54  const Box& bx = mfi.tilebox();
55  Box tbx, tby, tbz;
56 
57  tbx = mfi.tilebox(IntVect(1,0,0),xvel_ngrow);
58  tby = mfi.tilebox(IntVect(0,1,0),yvel_ngrow);
59  tbz = mfi.tilebox(IntVect(0,0,1),zvel_ngrow);
60 
61  // Don't actually try to fill w above or below the domain
62  if (tbz.smallEnd(2) < domain.smallEnd(2)) tbz.setSmall(2,domain.smallEnd(2));
63  if (tbz.bigEnd(2) > domain.bigEnd(2)+1) tbz.setBig(2,domain.bigEnd(2)+1);
64 
65  // Conserved/state variables on cell centers -- we use this for density
66  const Array4<const Real>& dens_arr = density.array(mfi);
67 
68  // Momentum on faces, to be computed
69  Array4<Real> const& momx = xmom.array(mfi);
70  Array4<Real> const& momy = ymom.array(mfi);
71  Array4<Real> const& momz = zmom.array(mfi);
72 
73  // Velocity on faces, used in computation
74  const Array4<Real const>& velx = xvel_in.const_array(mfi);
75  const Array4<Real const>& vely = yvel_in.const_array(mfi);
76  const Array4<Real const>& velz = zvel_in.const_array(mfi);
77 
78  // ********************************************************************************************
79 
80  if (c_vfrac==nullptr) {
81  ParallelFor(tbx, tby, tbz,
82  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
83  momx(i,j,k) = velx(i,j,k) * myhalf * (dens_arr(i,j,k,Rho_comp) + dens_arr(i-1,j,k,Rho_comp));
84  },
85  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
86  momy(i,j,k) = vely(i,j,k) * myhalf * (dens_arr(i,j,k,Rho_comp) + dens_arr(i,j-1,k,Rho_comp));
87  },
88  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
89  momz(i,j,k) = velz(i,j,k) * myhalf * (dens_arr(i,j,k,Rho_comp) + dens_arr(i,j,k-1,Rho_comp));
90  });
91  } else {
92  // EB
93  const Array4<const Real>& c_vfrac_arr = c_vfrac->const_array(mfi);
94 
95  ParallelFor(tbx, tby, tbz,
96  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
97  Real vfrac_i = c_vfrac_arr(i,j,k);
98  Real vfrac_im1 = c_vfrac_arr(i-1,j,k);
99  Real vfrac_sum = vfrac_i + vfrac_im1;
100  if (vfrac_sum > zero) {
101  Real rho = (vfrac_i > zero ? vfrac_i * dens_arr(i,j,k,Rho_comp) : zero)
102  + (vfrac_im1 > zero ? vfrac_im1 * dens_arr(i-1,j,k,Rho_comp) : zero);
103  rho /= vfrac_sum;
104  momx(i,j,k) = velx(i,j,k) * rho;
105  } else {
106  momx(i,j,k) = zero;
107  }
108  },
109  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
110  Real vfrac_j = c_vfrac_arr(i,j,k);
111  Real vfrac_jm1 = c_vfrac_arr(i,j-1,k);
112  Real vfrac_sum = vfrac_j + vfrac_jm1;
113  if (vfrac_sum > zero) {
114  Real rho = (vfrac_j > zero ? vfrac_j * dens_arr(i,j,k,Rho_comp) : zero)
115  + (vfrac_jm1 > zero ? vfrac_jm1 * dens_arr(i,j-1,k,Rho_comp) : zero);
116  rho /= vfrac_sum;
117  momy(i,j,k) = vely(i,j,k) * rho;
118  } else {
119  momy(i,j,k) = zero;
120  }
121  },
122  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
123  Real vfrac_k = c_vfrac_arr(i,j,k);
124  Real vfrac_km1 = c_vfrac_arr(i,j,k-1);
125  Real vfrac_sum = vfrac_k + vfrac_km1;
126  if (vfrac_sum > zero) {
127  Real rho = (vfrac_k > zero ? vfrac_k * dens_arr(i,j,k,Rho_comp) : zero)
128  + (vfrac_km1 > zero ? vfrac_km1 * dens_arr(i,j,k-1,Rho_comp) : zero);
129  rho /= vfrac_sum;
130  momz(i,j,k) = velz(i,j,k) * rho;
131  } else {
132  momz(i,j,k) = zero;
133  }
134  });
135  }
136  // ********************************************************************************************
137 
138  if (bx.smallEnd(0) == domain.smallEnd(0)) {
139  if (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::ext_dir)
140  {
141  ParallelFor(makeSlab(tbx,0,domain.smallEnd(0)), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
142  momx(i,j,k) = velx(i,j,k) * dens_arr(i-1,j,k,Rho_comp);
143  });
144  }
145  else if (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::ext_dir_upwind)
146  {
147  ParallelFor(makeSlab(tbx,0,domain.smallEnd(0)), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
148  if (velx(i,j,k) >= zero) {
149  momx(i,j,k) = velx(i,j,k) * dens_arr(i-1,j,k,Rho_comp);
150  }
151  });
152  }
153  }
154 
155  if (bx.bigEnd(0) == domain.bigEnd(0)) {
156  if (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::ext_dir)
157  {
158  ParallelFor(makeSlab(tbx,0,domain.bigEnd(0)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
159  momx(i,j,k) = velx(i,j,k) * dens_arr(i,j,k,Rho_comp);
160  });
161  }
162  else if (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::ext_dir_upwind)
163  {
164  ParallelFor(makeSlab(tbx,0,domain.bigEnd(0)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
165  if (velx(i,j,k) <= zero) {
166  momx(i,j,k) = velx(i,j,k) * dens_arr(i,j,k,Rho_comp);
167  }
168  });
169  }
170  }
171 
172  if (bx.smallEnd(1) == domain.smallEnd(1)) {
173  if (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::ext_dir)
174  {
175  ParallelFor(makeSlab(tby,1,domain.smallEnd(1)), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
176  momy(i,j,k) = vely(i,j,k) * dens_arr(i,j-1,k,Rho_comp);
177  });
178  }
179  else if (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::ext_dir_upwind)
180  {
181  ParallelFor(makeSlab(tby,1,domain.smallEnd(1)), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
182  if (vely(i,j,k) >= zero) {
183  momy(i,j,k) = vely(i,j,k) * dens_arr(i,j-1,k,Rho_comp);
184  }
185  });
186  }
187  }
188 
189  if (bx.bigEnd(1) == domain.bigEnd(1)) {
190  if (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::ext_dir)
191  {
192  ParallelFor(makeSlab(tby,1,domain.bigEnd(1)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
193  momy(i,j,k) = vely(i,j,k) * dens_arr(i,j,k,Rho_comp);
194  });
195  }
196  else if (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::ext_dir_upwind)
197  {
198  ParallelFor(makeSlab(tby,1,domain.bigEnd(1)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
199  if (vely(i,j,k) <= zero) {
200  momy(i,j,k) = vely(i,j,k) * dens_arr(i,j,k,Rho_comp);
201  }
202  });
203  }
204  }
205 
206  // Extrapolate momenta at top of domain (in case they need to be used in OmegaFromW)
207  if (bx.bigEnd(2) == domain.bigEnd(2))
208  {
209  if (bc_ptr_h[BCVars::xvel_bc].hi(2) != ERFBCType::ext_dir) {
210  ParallelFor(makeSlab(tbx,2,domain.bigEnd(2)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
211  momx(i,j,k) = momx(i,j,k-1);
212  });
213  }
214 
215  if (bc_ptr_h[BCVars::yvel_bc].hi(2) != ERFBCType::ext_dir) {
216  ParallelFor(makeSlab(tby,2,domain.bigEnd(2)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
217  momy(i,j,k) = momy(i,j,k-1);
218  });
219  }
220  }
221  } // end MFIter
222 }
OpenSide::hi
@ hi
OpenSide::lo
@ lo
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:39
ParallelFor
ParallelFor(fab_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);})
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:29
myhalf
constexpr amrex::Real myhalf
Definition: ERF_NumericalConstants.H:34
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
VelocityToMomentum
void VelocityToMomentum(const MultiFab &xvel_in, const IntVect &xvel_ngrow, const MultiFab &yvel_in, const IntVect &yvel_ngrow, const MultiFab &zvel_in, const IntVect &zvel_ngrow, const MultiFab &density, MultiFab &xmom, MultiFab &ymom, MultiFab &zmom, const Box &domain, const Vector< BCRec > &domain_bcs_type_h, const MultiFab *c_vfrac)
Definition: ERF_VelocityToMomentum.cpp:31
BCVars::yvel_bc
@ yvel_bc
Definition: ERF_IndexDefines.H:106
BCVars::cons_bc
@ cons_bc
Definition: ERF_IndexDefines.H:89
BCVars::xvel_bc
@ xvel_bc
Definition: ERF_IndexDefines.H:105
ERFBCType::ext_dir
@ ext_dir
Definition: ERF_IndexDefines.H:297
ERFBCType::ext_dir_upwind
@ ext_dir_upwind
Definition: ERF_IndexDefines.H:305
IntVars::ymom
@ ymom
Definition: ERF_IndexDefines.H:234
IntVars::zmom
@ zmom
Definition: ERF_IndexDefines.H:235
IntVars::xmom
@ xmom
Definition: ERF_IndexDefines.H:233
MicVar_Kess::rho
@ rho
Definition: ERF_Kessler.H:25
Here is the call graph for this function: