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

Functions

void MomentumToVelocity (MultiFab &xvel, MultiFab &yvel, MultiFab &zvel, const MultiFab &density, const MultiFab &xmom_in, const MultiFab &ymom_in, const MultiFab &zmom_in, const Box &domain, const Vector< BCRec > &domain_bcs_type_h, const MultiFab *c_vfrac)
 

Function Documentation

◆ MomentumToVelocity()

void MomentumToVelocity ( MultiFab &  xvel,
MultiFab &  yvel,
MultiFab &  zvel,
const MultiFab &  density,
const MultiFab &  xmom_in,
const MultiFab &  ymom_in,
const MultiFab &  zmom_in,
const Box &  domain,
const Vector< BCRec > &  domain_bcs_type_h,
const MultiFab *  c_vfrac 
)

Convert momentum to velocity by dividing by density averaged onto faces

Parameters
[out]xvelx-component of velocity
[out]yvely-component of velocity
[out]zvelz-component of velocity
[in]densitydensity at cell centers
[in]xmom_inx-component of momentum
[in]ymom_iny-component of momentum
[in]zmom_inz-component of momentum
[in]domainDomain at this level
[in]domain_bcs_type_hhost vector for domain boundary conditions
31  {
32  BL_PROFILE_VAR("MomentumToVelocity()",MomentumToVelocity);
33 
34  const BCRec* bc_ptr_h = domain_bcs_type_h.data();
35 
36 #ifdef _OPENMP
37 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
38 #endif
39  for ( MFIter mfi(density,TilingIfNotGPU()); mfi.isValid(); ++mfi)
40  {
41  // We need velocity in the interior ghost cells (init == real)
42  Box bx = mfi.tilebox();
43 
44  const Box& tbx = surroundingNodes(bx,0);
45  const Box& tby = surroundingNodes(bx,1);
46  const Box& tbz = surroundingNodes(bx,2);
47 
48  // Conserved variables on cell centers -- we use this for density
49  const Array4<const Real>& dens_arr = density.array(mfi);
50 
51  // Momentum on faces
52  Array4<Real const> const& momx = xmom_in.const_array(mfi);
53  Array4<Real const> const& momy = ymom_in.const_array(mfi);
54  Array4<Real const> const& momz = zmom_in.const_array(mfi);
55 
56  // Velocity on faces
57  const Array4<Real>& velx = xvel.array(mfi);
58  const Array4<Real>& vely = yvel.array(mfi);
59  const Array4<Real>& velz = zvel.array(mfi);
60 
61  if (c_vfrac==nullptr) {
62  ParallelFor(tbx, tby, tbz,
63  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
64  velx(i,j,k) = momx(i,j,k) * two / (dens_arr(i,j,k,Rho_comp) + dens_arr(i-1,j,k,Rho_comp));
65  },
66  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
67  vely(i,j,k) = momy(i,j,k) * two / (dens_arr(i,j,k,Rho_comp) + dens_arr(i,j-1,k,Rho_comp));
68  },
69  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
70  velz(i,j,k) = momz(i,j,k) * two / (dens_arr(i,j,k,Rho_comp) + dens_arr(i,j,k-1,Rho_comp));
71  });
72  } else {
73  // EB
74  const Array4<const Real>& c_vfrac_arr = c_vfrac->const_array(mfi);
75 
76  ParallelFor(tbx, tby, tbz,
77  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
78  Real vfrac_i = c_vfrac_arr(i,j,k);
79  Real vfrac_im1 = c_vfrac_arr(i-1,j,k);
80  Real vfrac_sum = vfrac_i + vfrac_im1;
81  if (vfrac_sum > zero) {
82  Real rho = (vfrac_i > zero ? vfrac_i * dens_arr(i,j,k,Rho_comp) : zero)
83  + (vfrac_im1 > zero ? vfrac_im1 * dens_arr(i-1,j,k,Rho_comp) : zero);
84  rho /= vfrac_sum;
85  velx(i,j,k) = momx(i,j,k) / rho;
86  }
87  },
88  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
89  Real vfrac_j = c_vfrac_arr(i,j,k);
90  Real vfrac_jm1 = c_vfrac_arr(i,j-1,k);
91  Real vfrac_sum = vfrac_j + vfrac_jm1;
92  if (vfrac_sum > zero) {
93  Real rho = (vfrac_j > zero ? vfrac_j * dens_arr(i,j,k,Rho_comp) : zero)
94  + (vfrac_jm1 > zero ? vfrac_jm1 * dens_arr(i,j-1,k,Rho_comp) : zero);
95  rho /= vfrac_sum;
96  vely(i,j,k) = momy(i,j,k) / rho;
97  }
98  },
99  [=] AMREX_GPU_DEVICE (int i, int j, int k) {
100  Real vfrac_k = c_vfrac_arr(i,j,k);
101  Real vfrac_km1 = c_vfrac_arr(i,j,k-1);
102  Real vfrac_sum = vfrac_k + vfrac_km1;
103  if (vfrac_sum > zero) {
104  Real rho = (vfrac_k > zero ? vfrac_k * dens_arr(i,j,k,Rho_comp) : zero)
105  + (vfrac_km1 > zero ? vfrac_km1 * dens_arr(i,j,k-1,Rho_comp) : zero);
106  rho /= vfrac_sum;
107  velz(i,j,k) = momz(i,j,k) / rho;
108  }
109  });
110  }
111 
112  if (bx.smallEnd(0) == domain.smallEnd(0)) {
113  if (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::ext_dir)
114  {
115  ParallelFor(makeSlab(tbx,0,domain.smallEnd(0)), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
116  velx(i,j,k) = momx(i,j,k) / dens_arr(i-1,j,k,Rho_comp);
117  });
118  }
119  else if (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::ext_dir_upwind)
120  {
121  ParallelFor(makeSlab(tbx,0,domain.smallEnd(0)), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
122  if (momx(i,j,k) >= zero) {
123  velx(i,j,k) = momx(i,j,k) / dens_arr(i-1,j,k,Rho_comp);
124  }
125  });
126  }
127  }
128 
129  if (bx.bigEnd(0) == domain.bigEnd(0)) {
130  if (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::ext_dir)
131  {
132  ParallelFor(makeSlab(tbx,0,domain.bigEnd(0)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
133  velx(i,j,k) = momx(i,j,k) / dens_arr(i,j,k,Rho_comp);
134  });
135  }
136  else if (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::ext_dir_upwind)
137  {
138  ParallelFor(makeSlab(tbx,0,domain.bigEnd(0)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
139  if (momx(i,j,k) <= zero) {
140  velx(i,j,k) = momx(i,j,k) / dens_arr(i,j,k,Rho_comp);
141  }
142  });
143  }
144  }
145 
146  if (bx.smallEnd(1) == domain.smallEnd(1)) {
147  if (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::ext_dir)
148  {
149  ParallelFor(makeSlab(tby,1,domain.smallEnd(1)), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
150  vely(i,j,k) = momy(i,j,k) / dens_arr(i,j-1,k,Rho_comp);
151  });
152  }
153  else if (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::ext_dir_upwind)
154  {
155  ParallelFor(makeSlab(tby,1,domain.smallEnd(1)), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
156  if (momy(i,j,k) >= zero) {
157  vely(i,j,k) = momy(i,j,k) / dens_arr(i,j-1,k,Rho_comp);
158  }
159  });
160  }
161  }
162 
163  if (bx.bigEnd(1) == domain.bigEnd(1)) {
164  if (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::ext_dir)
165  {
166  ParallelFor(makeSlab(tby,1,domain.bigEnd(1)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
167  vely(i,j,k) = momy(i,j,k) / dens_arr(i,j,k,Rho_comp);
168  });
169  }
170  else if (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::ext_dir_upwind)
171  {
172  ParallelFor(makeSlab(tby,1,domain.bigEnd(1)+1), [=] AMREX_GPU_DEVICE (int i, int j, int k) {
173  if (momy(i,j,k) <= zero) {
174  vely(i,j,k) = momy(i,j,k) / dens_arr(i,j,k,Rho_comp);
175  }
176  });
177  }
178  }
179  } // end MFIter
180 }
OpenSide::hi
@ hi
OpenSide::lo
@ lo
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:39
MomentumToVelocity
void MomentumToVelocity(MultiFab &xvel, MultiFab &yvel, MultiFab &zvel, const MultiFab &density, const MultiFab &xmom_in, const MultiFab &ymom_in, const MultiFab &zmom_in, const Box &domain, const Vector< BCRec > &domain_bcs_type_h, const MultiFab *c_vfrac)
Definition: ERF_MomentumToVelocity.cpp:25
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);})
two
constexpr amrex::Real two
Definition: ERF_NumericalConstants.H:31
zero
constexpr amrex::Real zero
Definition: ERF_NumericalConstants.H:29
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
BCVars::cons_bc
@ cons_bc
Definition: ERF_IndexDefines.H:89
ERFBCType::ext_dir
@ ext_dir
Definition: ERF_IndexDefines.H:297
ERFBCType::ext_dir_upwind
@ ext_dir_upwind
Definition: ERF_IndexDefines.H:305
MicVar_Kess::rho
@ rho
Definition: ERF_Kessler.H:25
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:215
Vars::zvel
@ zvel
Definition: ERF_IndexDefines.H:217
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:216

Referenced by ERF::AverageDownTo(), ERF::FillCoarsePatch(), ERF::FillIntermediatePatch(), ERF::FillPatchFineLevel(), and ERF::project_initial_velocity().

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