ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Functions
ERF_DiffusionSrcForMom_EB.cpp File Reference
#include <AMReX.H>
#include <AMReX_EB_Slopes_K.H>
#include <ERF_EB.H>
#include <ERF_Diffusion.H>
#include <ERF_IndexDefines.H>
#include <ERF_EBSlopes.H>
#include <ERF_DiffStruct.H>
#include <ERF_EBStruct.H>
Include dependency graph for ERF_DiffusionSrcForMom_EB.cpp:

Functions

void DiffusionSrcForMom_EB (const MFIter &mfi, [[maybe_unused]] const Box &domain, const Box &bxx, const Box &bxy, const Box &bxz, const Array4< Real > &rho_u_rhs, const Array4< Real > &rho_v_rhs, const Array4< Real > &rho_w_rhs, const Array4< const Real > &u_arr, const Array4< const Real > &v_arr, const Array4< const Real > &w_arr, const Array4< const Real > &tau11, const Array4< const Real > &tau22, const Array4< const Real > &tau33, const Array4< const Real > &tau12, const Array4< const Real > &tau13, const Array4< const Real > &tau23, const Real *dx_arr, const GpuArray< Real, AMREX_SPACEDIM > &dxInv, const Array4< const Real > &mf_mx, const Array4< const Real > &mf_ux, const Array4< const Real > &mf_vx, const Array4< const Real > &mf_my, const Array4< const Real > &mf_uy, const Array4< const Real > &mf_vy, const SolverChoice &solverChoice, const eb_ &ebfact, [[maybe_unused]] const BCRec *d_bcrec_ptr)
 

Function Documentation

◆ DiffusionSrcForMom_EB()

void DiffusionSrcForMom_EB ( const MFIter &  mfi,
[[maybe_unused] ] const Box &  domain,
const Box &  bxx,
const Box &  bxy,
const Box &  bxz,
const Array4< Real > &  rho_u_rhs,
const Array4< Real > &  rho_v_rhs,
const Array4< Real > &  rho_w_rhs,
const Array4< const Real > &  u_arr,
const Array4< const Real > &  v_arr,
const Array4< const Real > &  w_arr,
const Array4< const Real > &  tau11,
const Array4< const Real > &  tau22,
const Array4< const Real > &  tau33,
const Array4< const Real > &  tau12,
const Array4< const Real > &  tau13,
const Array4< const Real > &  tau23,
const Real *  dx_arr,
const GpuArray< Real, AMREX_SPACEDIM > &  dxInv,
const Array4< const Real > &  mf_mx,
const Array4< const Real > &  mf_ux,
const Array4< const Real > &  mf_vx,
const Array4< const Real > &  mf_my,
const Array4< const Real > &  mf_uy,
const Array4< const Real > &  mf_vy,
const SolverChoice &  solverChoice,
const eb_ &  ebfact,
[[maybe_unused] ] const BCRec *  d_bcrec_ptr 
)

Function for computing the momentum RHS for diffusion operator without terrain.

Parameters
[in]mfiMultiFab Iterator
[in]domaincomputational domain
[in]bxxnodal x box for x-mom
[in]bxynodal y box for y-mom
[in]bxznodal z box for z-mom
[out]rho_u_rhsRHS for x-mom
[out]rho_v_rhsRHS for y-mom
[out]rho_w_rhsRHS for z-mom
[in]tau1111 stress
[in]tau2222 stress
[in]tau3333 stress
[in]tau1212 stress
[in]tau1313 stress
[in]tau2323 stress
[in]dxInvinverse cell size array
[in]mf_mmap factor at cell center
[in]ebfactEB factories for cell- and face-centered variables
60 {
61  BL_PROFILE_VAR("DiffusionSrcForMom_EB()",DiffusionSrcForMom_EB);
62 
63  DiffChoice dc = solverChoice.diffChoice;
64  const bool l_use_constAlpha = ( dc.molec_diff_type == MolecDiffType::ConstantAlpha );
65  Real mu_eff = (l_use_constAlpha) ? 2.0 * dc.dynamic_viscosity / dc.rho0_trans
66  : 2.0 * dc.dynamic_viscosity;
67 
68  auto dxinv = dxInv[0], dyinv = dxInv[1], dzinv = dxInv[2];
69  Real dx = dx_arr[0], dy = dx_arr[1], dz = dx_arr[2];
70  Real vol = dx * dy * dz;
71 
72  EBChoice ebChoice = solverChoice.ebChoice;
73  const bool l_no_slip = (ebChoice.eb_boundary_type == EBBoundaryType::NoSlipWall);
74 
75  const bool l_constraint_x = solverChoice.diffChoice.eb_diff_constraint_x;
76  const bool l_constraint_y = solverChoice.diffChoice.eb_diff_constraint_y;
77  const bool l_constraint_z = solverChoice.diffChoice.eb_diff_constraint_z;
78 
79  // int n = 0;
80 
81  // EB u-factory
82  Array4<const EBCellFlag> u_cellflg = (ebfact.get_u_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
83  Array4<const Real > u_volfrac = (ebfact.get_u_const_factory())->getVolFrac().const_array(mfi);
84  Array4<const Real > u_volcent = (ebfact.get_u_const_factory())->getCentroid().const_array(mfi);
85  Array4<const Real > u_afrac_x = (ebfact.get_u_const_factory())->getAreaFrac()[0]->const_array(mfi);
86  Array4<const Real > u_afrac_y = (ebfact.get_u_const_factory())->getAreaFrac()[1]->const_array(mfi);
87  Array4<const Real > u_afrac_z = (ebfact.get_u_const_factory())->getAreaFrac()[2]->const_array(mfi);
88  Array4<const Real > u_bcent = (ebfact.get_u_const_factory())->getBndryCent().const_array(mfi);
89  Array4<const Real > u_bnorm = (ebfact.get_u_const_factory())->getBndryNorm().const_array(mfi);
90 
91  // EB v-factory
92  Array4<const EBCellFlag> v_cellflg = (ebfact.get_v_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
93  Array4<const Real > v_volfrac = (ebfact.get_v_const_factory())->getVolFrac().const_array(mfi);
94  Array4<const Real > v_volcent = (ebfact.get_v_const_factory())->getCentroid().const_array(mfi);
95  Array4<const Real > v_afrac_x = (ebfact.get_v_const_factory())->getAreaFrac()[0]->const_array(mfi);
96  Array4<const Real > v_afrac_y = (ebfact.get_v_const_factory())->getAreaFrac()[1]->const_array(mfi);
97  Array4<const Real > v_afrac_z = (ebfact.get_v_const_factory())->getAreaFrac()[2]->const_array(mfi);
98  Array4<const Real > v_bcent = (ebfact.get_v_const_factory())->getBndryCent().const_array(mfi);
99  Array4<const Real > v_bnorm = (ebfact.get_v_const_factory())->getBndryNorm().const_array(mfi);
100 
101  // EB w-factory
102  Array4<const EBCellFlag> w_cellflg = (ebfact.get_w_const_factory())->getMultiEBCellFlagFab()[mfi].const_array();
103  Array4<const Real > w_volfrac = (ebfact.get_w_const_factory())->getVolFrac().const_array(mfi);
104  Array4<const Real > w_volcent = (ebfact.get_w_const_factory())->getCentroid().const_array(mfi);
105  Array4<const Real > w_afrac_x = (ebfact.get_w_const_factory())->getAreaFrac()[0]->const_array(mfi);
106  Array4<const Real > w_afrac_y = (ebfact.get_w_const_factory())->getAreaFrac()[1]->const_array(mfi);
107  Array4<const Real > w_afrac_z = (ebfact.get_w_const_factory())->getAreaFrac()[2]->const_array(mfi);
108  Array4<const Real > w_bcent = (ebfact.get_w_const_factory())->getBndryCent().const_array(mfi);
109  Array4<const Real > w_bnorm = (ebfact.get_w_const_factory())->getBndryNorm().const_array(mfi);
110 
111  ParallelFor(bxx, bxy, bxz,
112  [=] AMREX_GPU_DEVICE (int i, int j, int k)
113  {
114  if (u_volfrac(i,j,k)>0.) {
115 
116  // Inv Jacobian
117  Real mfsq = mf_ux(i,j,0) * mf_uy(i,j,0);
118 
119  Real diffContrib = ( (tau11(i , j , k ) * u_afrac_x(i+1,j ,k )
120  - tau11(i-1, j , k ) * u_afrac_x(i ,j ,k ) ) * dxinv * mfsq
121  + (tau12(i , j+1, k ) * u_afrac_y(i ,j+1,k )
122  - tau12(i , j , k ) * u_afrac_y(i ,j ,k ) ) * dyinv * mfsq
123  + (tau13(i , j , k+1) * u_afrac_z(i ,j ,k+1)
124  - tau13(i , j , k ) * u_afrac_z(i ,j ,k )) * dzinv );
125  diffContrib /= u_volfrac(i,j,k);
126 
127  rho_u_rhs(i,j,k) -= diffContrib;
128 
129  if (!l_constraint_x && l_no_slip && u_cellflg(i,j,k).isSingleValued()) {
130 
131  Real axm = u_afrac_x(i ,j ,k );
132  Real axp = u_afrac_x(i+1,j ,k );
133  Real aym = u_afrac_y(i ,j ,k );
134  Real ayp = u_afrac_y(i ,j+1,k );
135  Real azm = u_afrac_z(i ,j ,k );
136  Real azp = u_afrac_z(i ,j ,k+1);
137 
138  Real adx = (axm-axp) * dy * dz;
139  Real ady = (aym-ayp) * dx * dz;
140  Real adz = (azm-azp) * dx * dy;
141 
142  Real barea = std::sqrt(adx*adx + ady*ady + adz*adz);
143 
144  const RealVect bcent_eb {u_bcent(i,j,k,0), u_bcent(i,j,k,1), u_bcent(i,j,k,2)};
145  const Real Dirichlet_u {0.};
146  const Real Dirichlet_v {0.};
147  const Real Dirichlet_w {0.};
148 
149  GpuArray<Real,AMREX_SPACEDIM> slopes_u;
150  GpuArray<Real,AMREX_SPACEDIM> slopes_v;
151  GpuArray<Real,AMREX_SPACEDIM> slopes_w;
152 
153  slopes_u = erf_calc_slopes_eb_Dirichlet ( dx, dy, dz, i, j, k, bcent_eb, Dirichlet_u, u_arr, u_volcent, u_cellflg);
154  slopes_v = erf_calc_slopes_eb_Dirichlet_staggered( Vars::xvel, Vars::yvel, dx, dy, dz, i, j, k, bcent_eb, Dirichlet_v, v_arr, v_volcent, v_cellflg);
155  slopes_w = erf_calc_slopes_eb_Dirichlet_staggered( Vars::xvel, Vars::zvel, dx, dy, dz, i, j, k, bcent_eb, Dirichlet_w, w_arr, w_volcent, w_cellflg);
156 
157  Real dudx = slopes_u[0];
158  Real dudy = slopes_u[1];
159  Real dudz = slopes_u[2];
160  Real dvdx = slopes_v[0];
161  Real dvdy = slopes_v[1];
162  Real dwdx = slopes_w[0];
163  Real dwdz = slopes_w[2];
164 
165  Real tau11_eb = ( dudx - ( dudx + dvdy + dwdz ) / 3. );
166  Real tau12_eb = 0.5 * (dudy + dvdx);
167  Real tau13_eb = 0.5 * (dudz + dwdx);
168 
169  Real dudn = -(u_bnorm(i,j,k,0) * tau11_eb + u_bnorm(i,j,k,1) * tau12_eb + u_bnorm(i,j,k,2) * tau13_eb);
170 
171  rho_u_rhs(i,j,k) -= mu_eff * barea * dudn / (vol * u_volfrac(i,j,k));
172  }
173  }
174 
175  },
176  [=] AMREX_GPU_DEVICE (int i, int j, int k)
177  {
178  if (v_volfrac(i,j,k)>0.) {
179 
180  // Inv Jacobian
181  Real mfsq = mf_vx(i,j,0) * mf_vy(i,j,0);
182 
183  Real diffContrib = ( (tau12(i+1, j , k ) * v_afrac_x(i+1,j ,k )
184  - tau12(i , j , k ) * v_afrac_x(i ,j ,k ) ) * dxinv * mfsq
185  + (tau22(i , j , k ) * v_afrac_y(i ,j+1,k )
186  - tau22(i , j-1, k ) * v_afrac_y(i ,j ,k ) ) * dyinv * mfsq
187  + (tau23(i , j , k+1) * v_afrac_z(i ,j ,k+1)
188  - tau23(i , j , k ) * v_afrac_z(i ,j ,k ) ) * dzinv );
189  diffContrib /= v_volfrac(i,j,k);
190 
191  rho_v_rhs(i,j,k) -= diffContrib;
192 
193  if (!l_constraint_y && l_no_slip && v_cellflg(i,j,k).isSingleValued()) {
194 
195  Real axm = v_afrac_x(i ,j ,k );
196  Real axp = v_afrac_x(i+1,j ,k );
197  Real aym = v_afrac_y(i ,j ,k );
198  Real ayp = v_afrac_y(i ,j+1,k );
199  Real azm = v_afrac_z(i ,j ,k );
200  Real azp = v_afrac_z(i ,j ,k+1);
201 
202  Real adx = (axm-axp) * dy * dz;
203  Real ady = (aym-ayp) * dx * dz;
204  Real adz = (azm-azp) * dx * dy;
205 
206  Real barea = std::sqrt(adx*adx + ady*ady + adz*adz);
207 
208  const RealVect bcent_eb {v_bcent(i,j,k,0), v_bcent(i,j,k,1), v_bcent(i,j,k,2)};
209  const Real Dirichlet_u {0.};
210  const Real Dirichlet_v {0.};
211  const Real Dirichlet_w {0.};
212 
213  GpuArray<Real,AMREX_SPACEDIM> slopes_u;
214  GpuArray<Real,AMREX_SPACEDIM> slopes_v;
215  GpuArray<Real,AMREX_SPACEDIM> slopes_w;
216 
217  slopes_u = erf_calc_slopes_eb_Dirichlet_staggered( Vars::yvel, Vars::xvel, dx, dy, dz, i, j, k, bcent_eb, Dirichlet_u, u_arr, u_volcent, u_cellflg);
218  slopes_v = erf_calc_slopes_eb_Dirichlet ( dx, dy, dz, i, j, k, bcent_eb, Dirichlet_v, v_arr, v_volcent, v_cellflg);
219  slopes_w = erf_calc_slopes_eb_Dirichlet_staggered( Vars::yvel, Vars::zvel, dx, dy, dz, i, j, k, bcent_eb, Dirichlet_w, w_arr, w_volcent, w_cellflg);
220 
221  Real dudx = slopes_u[0];
222  Real dudy = slopes_u[1];
223  Real dvdx = slopes_v[0];
224  Real dvdy = slopes_v[1];
225  Real dvdz = slopes_v[2];
226  Real dwdy = slopes_w[1];
227  Real dwdz = slopes_w[2];
228 
229  Real tau22_eb = ( dvdy - ( dudx + dvdy + dwdz ) / 3. );
230  Real tau12_eb = 0.5 * (dudy + dvdx);
231  Real tau23_eb = 0.5 * (dvdz + dwdy);
232 
233  Real dvdn = -(v_bnorm(i,j,k,0) * tau12_eb + v_bnorm(i,j,k,1) * tau22_eb + v_bnorm(i,j,k,2) * tau23_eb);
234 
235  rho_v_rhs(i,j,k) -= mu_eff * barea * dvdn / (vol * v_volfrac(i,j,k));
236  }
237  }
238  },
239  [=] AMREX_GPU_DEVICE (int i, int j, int k)
240  {
241  if (w_volfrac(i,j,k)>0.) {
242 
243  // Inv Jacobian
244  Real mfsq = mf_mx(i,j,0) * mf_my(i,j,0);
245 
246  Real diffContrib = ( (tau13(i+1, j , k ) * w_afrac_x(i+1,j ,k )
247  - tau13(i , j , k ) * w_afrac_x(i ,j ,k ) ) * dxinv * mfsq
248  + (tau23(i , j+1, k ) * w_afrac_y(i ,j+1,k )
249  - tau23(i , j , k ) * w_afrac_y(i ,j ,k ) ) * dyinv * mfsq
250  + (tau33(i , j , k ) * w_afrac_z(i ,j ,k+1)
251  - tau33(i , j , k-1) * w_afrac_z(i ,j ,k ) ) * dzinv );
252  diffContrib /= w_volfrac(i,j,k);
253  rho_w_rhs(i,j,k) -= diffContrib;
254 
255  if (!l_constraint_z && l_no_slip && w_cellflg(i,j,k).isSingleValued()) {
256 
257  Real axm = w_afrac_x(i ,j ,k );
258  Real axp = w_afrac_x(i+1,j ,k );
259  Real aym = w_afrac_y(i ,j ,k );
260  Real ayp = w_afrac_y(i ,j+1,k );
261  Real azm = w_afrac_z(i ,j ,k );
262  Real azp = w_afrac_z(i ,j ,k+1);
263 
264  Real adx = (axm-axp) * dy * dz;
265  Real ady = (aym-ayp) * dx * dz;
266  Real adz = (azm-azp) * dx * dy;
267 
268  Real barea = std::sqrt(adx*adx + ady*ady + adz*adz);
269 
270  const RealVect bcent_eb {w_bcent(i,j,k,0), w_bcent(i,j,k,1), w_bcent(i,j,k,2)};
271  const Real Dirichlet_u {0.};
272  const Real Dirichlet_v {0.};
273  const Real Dirichlet_w {0.};
274 
275  GpuArray<Real,AMREX_SPACEDIM> slopes_u;
276  GpuArray<Real,AMREX_SPACEDIM> slopes_v;
277  GpuArray<Real,AMREX_SPACEDIM> slopes_w;
278 
279  slopes_u = erf_calc_slopes_eb_Dirichlet_staggered( Vars::zvel, Vars::xvel, dx, dy, dz, i, j, k, bcent_eb, Dirichlet_u, u_arr, u_volcent, u_cellflg);
280  slopes_v = erf_calc_slopes_eb_Dirichlet_staggered( Vars::zvel, Vars::yvel, dx, dy, dz, i, j, k, bcent_eb, Dirichlet_v, v_arr, v_volcent, v_cellflg);
281  slopes_w = erf_calc_slopes_eb_Dirichlet ( dx, dy, dz, i, j, k, bcent_eb, Dirichlet_w, w_arr, w_volcent, w_cellflg);
282 
283  Real dudx = slopes_u[0];
284  Real dudz = slopes_u[2];
285  Real dvdy = slopes_v[1];
286  Real dvdz = slopes_v[2];
287  Real dwdx = slopes_w[0];
288  Real dwdy = slopes_w[1];
289  Real dwdz = slopes_w[2];
290 
291  Real tau33_eb = ( dwdz - ( dudx + dvdy + dwdz ) / 3. );
292  Real tau13_eb = 0.5 * (dudz + dwdx);
293  Real tau23_eb = 0.5 * (dvdz + dwdy);
294 
295  Real dwdn = -(w_bnorm(i,j,k,0) * tau13_eb + w_bnorm(i,j,k,1) * tau23_eb + w_bnorm(i,j,k,2) * tau33_eb);
296 
297  rho_w_rhs(i,j,k) -= mu_eff * barea * dwdn / (vol * w_volfrac(i,j,k));
298  }
299  }
300  });
301 
302 }
tau12
@ tau12
Definition: ERF_DataStruct.H:32
tau23
@ tau23
Definition: ERF_DataStruct.H:32
tau33
@ tau33
Definition: ERF_DataStruct.H:32
tau22
@ tau22
Definition: ERF_DataStruct.H:32
tau11
@ tau11
Definition: ERF_DataStruct.H:32
tau13
@ tau13
Definition: ERF_DataStruct.H:32
MolecDiffType::ConstantAlpha
@ ConstantAlpha
DiffusionSrcForMom_EB
void DiffusionSrcForMom_EB(const MFIter &mfi, [[maybe_unused]] const Box &domain, const Box &bxx, const Box &bxy, const Box &bxz, const Array4< Real > &rho_u_rhs, const Array4< Real > &rho_v_rhs, const Array4< Real > &rho_w_rhs, const Array4< const Real > &u_arr, const Array4< const Real > &v_arr, const Array4< const Real > &w_arr, const Array4< const Real > &tau11, const Array4< const Real > &tau22, const Array4< const Real > &tau33, const Array4< const Real > &tau12, const Array4< const Real > &tau13, const Array4< const Real > &tau23, const Real *dx_arr, const GpuArray< Real, AMREX_SPACEDIM > &dxInv, const Array4< const Real > &mf_mx, const Array4< const Real > &mf_ux, const Array4< const Real > &mf_vx, const Array4< const Real > &mf_my, const Array4< const Real > &mf_uy, const Array4< const Real > &mf_vy, const SolverChoice &solverChoice, const eb_ &ebfact, [[maybe_unused]] const BCRec *d_bcrec_ptr)
Definition: ERF_DiffusionSrcForMom_EB.cpp:34
erf_calc_slopes_eb_Dirichlet
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > erf_calc_slopes_eb_Dirichlet(amrex::Real dx, amrex::Real dy, amrex::Real dz, int i, int j, int k, amrex::RealVect const &bcent_eb, amrex::Real const state_eb, amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::Real const > const &ccent, amrex::Array4< amrex::EBCellFlag const > const &flag)
Definition: ERF_EBSlopes.H:20
erf_calc_slopes_eb_Dirichlet_staggered
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > erf_calc_slopes_eb_Dirichlet_staggered(int igrid_query, int igrid_data, amrex::Real dx, amrex::Real dy, amrex::Real dz, int i, int j, int k, amrex::RealVect const &bcent_eb, amrex::Real const state_eb, amrex::Array4< amrex::Real const > const &state, amrex::Array4< amrex::Real const > const &ccent, amrex::Array4< amrex::EBCellFlag const > const &flag)
Definition: ERF_EBSlopes.H:146
EBBoundaryType::NoSlipWall
@ NoSlipWall
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
eb_::get_w_const_factory
eb_aux_ const * get_w_const_factory() const noexcept
Definition: ERF_EB.H:52
eb_::get_v_const_factory
eb_aux_ const * get_v_const_factory() const noexcept
Definition: ERF_EB.H:51
eb_::get_u_const_factory
eb_aux_ const * get_u_const_factory() const noexcept
Definition: ERF_EB.H:50
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:141
Vars::zvel
@ zvel
Definition: ERF_IndexDefines.H:143
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:142
DiffChoice
Definition: ERF_DiffStruct.H:19
DiffChoice::rho0_trans
amrex::Real rho0_trans
Definition: ERF_DiffStruct.H:91
DiffChoice::eb_diff_constraint_z
bool eb_diff_constraint_z
Definition: ERF_DiffStruct.H:99
DiffChoice::eb_diff_constraint_x
bool eb_diff_constraint_x
Definition: ERF_DiffStruct.H:97
DiffChoice::eb_diff_constraint_y
bool eb_diff_constraint_y
Definition: ERF_DiffStruct.H:98
DiffChoice::molec_diff_type
MolecDiffType molec_diff_type
Definition: ERF_DiffStruct.H:84
DiffChoice::dynamic_viscosity
amrex::Real dynamic_viscosity
Definition: ERF_DiffStruct.H:96
EBChoice
Definition: ERF_EBStruct.H:19
SolverChoice::diffChoice
DiffChoice diffChoice
Definition: ERF_DataStruct.H:991
SolverChoice::ebChoice
EBChoice ebChoice
Definition: ERF_DataStruct.H:995
Here is the call graph for this function: