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) * 0.5 * (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) * 0.5 * (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) * 0.5 * (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  if (c_vfrac_arr(i,j,k) > 0.0 || c_vfrac_arr(i-1,j,k) > 0.0) {
98  Real rho = (c_vfrac_arr(i,j,k) * dens_arr(i,j,k,Rho_comp)
99  + c_vfrac_arr(i-1,j,k) * dens_arr(i-1,j,k,Rho_comp))
100  / (c_vfrac_arr(i,j,k) + c_vfrac_arr(i-1,j,k));
101  momx(i,j,k) = velx(i,j,k) * rho;
102  } else {
103  momx(i,j,k) = velx(i,j,k) * 0.5 * (dens_arr(i,j,k,Rho_comp) + dens_arr(i-1,j,k,Rho_comp));
104  }
105  },
106  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
107  if (c_vfrac_arr(i,j,k) > 0.0 || c_vfrac_arr(i,j-1,k) > 0.0) {
108  Real rho = (c_vfrac_arr(i,j,k) * dens_arr(i,j,k,Rho_comp)
109  + c_vfrac_arr(i,j-1,k) * dens_arr(i,j-1,k,Rho_comp))
110  / (c_vfrac_arr(i,j,k) + c_vfrac_arr(i,j-1,k));
111  momy(i,j,k) = vely(i,j,k) * rho;
112  } else {
113  momy(i,j,k) = vely(i,j,k) * 0.5 * (dens_arr(i,j,k,Rho_comp) + dens_arr(i,j-1,k,Rho_comp));
114  }
115  },
116  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
117  if (c_vfrac_arr(i,j,k) > 0.0 || c_vfrac_arr(i,j,k-1) > 0.0) {
118  Real rho = (c_vfrac_arr(i,j,k) * dens_arr(i,j,k,Rho_comp)
119  + c_vfrac_arr(i,j,k-1) * dens_arr(i,j,k-1,Rho_comp))
120  / (c_vfrac_arr(i,j,k) + c_vfrac_arr(i,j,k-1));
121  momz(i,j,k) = velz(i,j,k) * rho;
122  } else {
123  momz(i,j,k) = velz(i,j,k) * 0.5 * (dens_arr(i,j,k,Rho_comp) + dens_arr(i,j,k-1,Rho_comp));
124  }
125  });
126  }
127  // ********************************************************************************************
128 
129  if (bx.smallEnd(0) == domain.smallEnd(0)) {
130  if (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::ext_dir)
131  {
132  ParallelFor(makeSlab(tbx,0,domain.smallEnd(0)), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
133  momx(i,j,k) = velx(i,j,k) * dens_arr(i-1,j,k,Rho_comp);
134  });
135  }
136  else if (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::ext_dir_upwind)
137  {
138  ParallelFor(makeSlab(tbx,0,domain.smallEnd(0)), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
139  if (velx(i,j,k) >= 0.) {
140  momx(i,j,k) = velx(i,j,k) * dens_arr(i-1,j,k,Rho_comp);
141  }
142  });
143  }
144  }
145 
146  if (bx.bigEnd(0) == domain.bigEnd(0)) {
147  if (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::ext_dir)
148  {
149  ParallelFor(makeSlab(tbx,0,domain.bigEnd(0)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
150  momx(i,j,k) = velx(i,j,k) * dens_arr(i,j,k,Rho_comp);
151  });
152  }
153  else if (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::ext_dir_upwind)
154  {
155  ParallelFor(makeSlab(tbx,0,domain.smallEnd(0)), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
156  if (velx(i,j,k) <= 0.) {
157  momx(i,j,k) = velx(i,j,k) * dens_arr(i,j,k,Rho_comp);
158  }
159  });
160  }
161  }
162 
163  if (bx.smallEnd(1) == domain.smallEnd(1)) {
164  if (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::ext_dir)
165  {
166  ParallelFor(makeSlab(tby,1,domain.smallEnd(1)), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
167  momy(i,j,k) = vely(i,j,k) / dens_arr(i,j-1,k,Rho_comp);
168  });
169  }
170  else if (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::ext_dir_upwind)
171  {
172  ParallelFor(makeSlab(tby,1,domain.smallEnd(1)), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
173  if (vely(i,j,k) >= 0.) {
174  momy(i,j,k) = vely(i,j,k) / dens_arr(i,j-1,k,Rho_comp);
175  }
176  });
177  }
178  }
179 
180  if (bx.bigEnd(1) == domain.bigEnd(1)) {
181  if (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::ext_dir)
182  {
183  ParallelFor(makeSlab(tby,1,domain.bigEnd(1)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
184  momy(i,j,k) = vely(i,j,k) / dens_arr(i,j,k,Rho_comp);
185  });
186  }
187  else if (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::ext_dir_upwind)
188  {
189  ParallelFor(makeSlab(tby,1,domain.smallEnd(1)), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
190  if (vely(i,j,k) <= 0.) {
191  momy(i,j,k) = vely(i,j,k) / dens_arr(i,j,k,Rho_comp);
192  }
193  });
194  }
195  }
196 
197  // Extrapolate momenta at top of domain (in case they need to be used in OmegaFromW)
198  if (bx.bigEnd(2) == domain.bigEnd(2))
199  {
200  ParallelFor(makeSlab(tbx,2,domain.bigEnd(2)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
201  momx(i,j,k) = momx(i,j,k-1);
202  });
203  ParallelFor(makeSlab(tby,2,domain.bigEnd(2)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
204  momy(i,j,k) = momy(i,j,k-1);
205  });
206  }
207  } // end MFIter
208 }
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
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::cons_bc
@ cons_bc
Definition: ERF_IndexDefines.H:76
ERFBCType::ext_dir
@ ext_dir
Definition: ERF_IndexDefines.H:209
ERFBCType::ext_dir_upwind
@ ext_dir_upwind
Definition: ERF_IndexDefines.H:217
IntVars::ymom
@ ymom
Definition: ERF_IndexDefines.H:160
IntVars::zmom
@ zmom
Definition: ERF_IndexDefines.H:161
IntVars::xmom
@ xmom
Definition: ERF_IndexDefines.H:159
MicVar_Kess::rho
@ rho
Definition: ERF_Kessler.H:22