1 #ifndef ERF_EBSLOPES_H_
2 #define ERF_EBSLOPES_H_
4 #include <AMReX_EBCellFlag.H>
18 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
19 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM>
21 [[maybe_unused]]
int igrid_data,
24 amrex::RealVect
const& bcent_eb,
26 amrex::Array4<amrex::Real const>
const& state,
27 amrex::Array4<amrex::Real const>
const& ccent,
28 amrex::Array4<amrex::EBCellFlag const>
const& flag)
30 AMREX_ASSERT( igrid_query == igrid_data );
34 amrex::RealVect bias{0.,0.,0.};
36 constexpr
int dim_a = 27;
49 for(
int kk(kk_lo); kk<=kk_hi; kk++) {
50 for(
int jj(jj_lo); jj<=jj_hi; jj++) {
51 for(
int ii(ii_lo); ii<=ii_hi; ii++) {
53 if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0)) {
54 A[ll][0] = (
amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0) + bias[0] - bcent_eb[0] ) * dx;
55 A[ll][1] = (
amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1) + bias[1] - bcent_eb[1] ) * dy;
56 A[ll][2] = (
amrex::Real(kk) + ccent(i+ii,j+jj,k+kk,2) + bias[2] - bcent_eb[2] ) * dz;
71 for(
int kk(kk_lo); kk<=kk_hi; kk++) {
72 for(
int jj(jj_lo); jj<=jj_hi; jj++) {
73 for(
int ii(ii_lo); ii<=ii_hi; ii++) {
75 if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0)) {
76 du[ll] = state(i+ii,j+jj,k+kk) - state_eb;
86 for(
int jj(0); jj<AMREX_SPACEDIM; ++jj){
87 for(
int ii(0); ii<AMREX_SPACEDIM; ++ii){
93 for(
int lc(0); lc < dim_a; ++lc)
95 AtA[0][0] += A[lc][0]* A[lc][0];
96 AtA[0][1] += A[lc][0]* A[lc][1];
97 AtA[0][2] += A[lc][0]* A[lc][2];
98 AtA[1][1] += A[lc][1]* A[lc][1];
99 AtA[1][2] += A[lc][1]* A[lc][2];
100 AtA[2][2] += A[lc][2]* A[lc][2];
102 Atb[0] += A[lc][0]*du[lc];
103 Atb[1] += A[lc][1]*du[lc];
104 Atb[2] += A[lc][2]*du[lc];
108 AtA[1][0] = AtA[0][1];
109 AtA[2][0] = AtA[0][2];
110 AtA[2][1] = AtA[1][2];
113 AtA[0][0]*(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[2][1]) -
114 AtA[0][1]*(AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
115 AtA[0][2]*(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
118 Atb[0] *(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[1][2]) -
119 AtA[0][1]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) +
120 AtA[0][2]*(Atb[1] * AtA[2][1] - AtA[1][1]*Atb[2] );
126 AtA[0][0]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) -
127 Atb[0] * (AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
128 AtA[0][2]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]);
134 AtA[0][0]*(AtA[1][1]*Atb[2] - Atb[1] *AtA[1][2]) -
135 AtA[0][1]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]) +
136 Atb[0] *(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
141 return {xslope,yslope,zslope};
155 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
156 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM>
161 amrex::RealVect
const& bcent_eb,
163 amrex::Array4<amrex::Real const>
const& state,
164 amrex::Array4<amrex::Real const>
const& ccent,
165 amrex::Array4<amrex::EBCellFlag const>
const& flag)
169 constexpr
int dim_a = 48;
195 amrex::RealVect bias{0.,0.,0.};
196 bias[igrid_query-1] = 0.5;
197 bias[igrid_data-1] = -0.5;
203 for(
int kk(kk_lo); kk<=kk_hi; kk++) {
204 for(
int jj(jj_lo); jj<=jj_hi; jj++) {
205 for(
int ii(ii_lo); ii<=ii_hi; ii++) {
207 if (!flag(i+ii,j+jj,k+kk).isCovered()) {
208 A[ll][0] = (
amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0) + bias[0] - bcent_eb[0] ) * dx;
209 A[ll][1] = (
amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1) + bias[1] - bcent_eb[1] ) * dy;
210 A[ll][2] = (
amrex::Real(kk) + ccent(i+ii,j+jj,k+kk,2) + bias[2] - bcent_eb[2] ) * dz;
225 for(
int kk(kk_lo); kk<=kk_hi; kk++) {
226 for(
int jj(jj_lo); jj<=jj_hi; jj++) {
227 for(
int ii(ii_lo); ii<=ii_hi; ii++) {
229 if (!flag(i+ii,j+jj,k+kk).isCovered()) {
230 du[ll] = state(i+ii,j+jj,k+kk) - state_eb;
240 for(
int jj(0); jj<AMREX_SPACEDIM; ++jj){
241 for(
int ii(0); ii<AMREX_SPACEDIM; ++ii){
247 for(
int lc(0); lc < dim_a; ++lc)
249 AtA[0][0] += A[lc][0]* A[lc][0];
250 AtA[0][1] += A[lc][0]* A[lc][1];
251 AtA[0][2] += A[lc][0]* A[lc][2];
252 AtA[1][1] += A[lc][1]* A[lc][1];
253 AtA[1][2] += A[lc][1]* A[lc][2];
254 AtA[2][2] += A[lc][2]* A[lc][2];
256 Atb[0] += A[lc][0]*du[lc];
257 Atb[1] += A[lc][1]*du[lc];
258 Atb[2] += A[lc][2]*du[lc];
262 AtA[1][0] = AtA[0][1];
263 AtA[2][0] = AtA[0][2];
264 AtA[2][1] = AtA[1][2];
267 AtA[0][0]*(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[2][1]) -
268 AtA[0][1]*(AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
269 AtA[0][2]*(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
272 Atb[0] *(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[1][2]) -
273 AtA[0][1]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) +
274 AtA[0][2]*(Atb[1] * AtA[2][1] - AtA[1][1]*Atb[2] );
280 AtA[0][0]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) -
281 Atb[0] * (AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
282 AtA[0][2]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]);
288 AtA[0][0]*(AtA[1][1]*Atb[2] - Atb[1] *AtA[1][2]) -
289 AtA[0][1]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]) +
290 Atb[0] *(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
295 return {xslope,yslope,zslope};
298 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
299 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM>
304 amrex::Array4<amrex::Real const>
const& state,
305 amrex::Array4<amrex::Real const>
const& ccent,
306 amrex::Array4<amrex::EBCellFlag const>
const& flag)
312 constexpr
int dim_a = 36;
338 amrex::RealVect bias{0.,0.,0.};
339 bias[igrid_query-1] = 0.5;
345 for(
int kk(kk_lo); kk<=kk_hi; kk++) {
346 for(
int jj(jj_lo); jj<=jj_hi; jj++) {
347 for(
int ii(ii_lo); ii<=ii_hi; ii++) {
350 if (!flag(i+ii,j+jj,k+kk).isCovered()) {
351 A[ll][0] = (
amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0) + bias[0] - ccent(i,j,k,0) ) * dx;
352 A[ll][1] = (
amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1) + bias[1] - ccent(i,j,k,1) ) * dy;
353 A[ll][2] = (
amrex::Real(kk) + ccent(i+ii,j+jj,k+kk,2) + bias[2] - ccent(i,j,k,2) ) * dz;
368 for(
int kk(kk_lo); kk<=kk_hi; kk++) {
369 for(
int jj(jj_lo); jj<=jj_hi; jj++) {
370 for(
int ii(ii_lo); ii<=ii_hi; ii++) {
373 if (!flag(i+ii,j+jj,k+kk).isCovered()) {
374 du[ll] = state(i+ii,j+jj,k+kk) - state(i,j,k);
384 for(
int jj(0); jj<AMREX_SPACEDIM; ++jj){
385 for(
int ii(0); ii<AMREX_SPACEDIM; ++ii){
391 for(
int lc(0); lc < dim_a; ++lc)
393 AtA[0][0] += A[lc][0]* A[lc][0];
394 AtA[0][1] += A[lc][0]* A[lc][1];
395 AtA[0][2] += A[lc][0]* A[lc][2];
396 AtA[1][1] += A[lc][1]* A[lc][1];
397 AtA[1][2] += A[lc][1]* A[lc][2];
398 AtA[2][2] += A[lc][2]* A[lc][2];
400 Atb[0] += A[lc][0]*du[lc];
401 Atb[1] += A[lc][1]*du[lc];
402 Atb[2] += A[lc][2]*du[lc];
406 AtA[1][0] = AtA[0][1];
407 AtA[2][0] = AtA[0][2];
408 AtA[2][1] = AtA[1][2];
411 AtA[0][0]*(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[2][1]) -
412 AtA[0][1]*(AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
413 AtA[0][2]*(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
416 Atb[0] *(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[1][2]) -
417 AtA[0][1]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) +
418 AtA[0][2]*(Atb[1] * AtA[2][1] - AtA[1][1]*Atb[2] );
424 AtA[0][0]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) -
425 Atb[0] * (AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
426 AtA[0][2]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]);
432 AtA[0][0]*(AtA[1][1]*Atb[2] - Atb[1] *AtA[1][2]) -
433 AtA[0][1]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]) +
434 Atb[0] *(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
439 return {xslope,yslope,zslope};
442 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
443 amrex::GpuArray<amrex::Real,AMREX_SPACEDIM>
445 [[maybe_unused]]
int igrid_data,
448 const amrex::Array4<const amrex::Real>& vel_arr,
449 amrex::Array4<amrex::Real const>
const& state,
450 amrex::Array4<amrex::Real const>
const& ccent,
451 amrex::Array4<amrex::EBCellFlag const>
const& flag)
464 if (igrid_query ==
Vars::xvel && vel_arr(i,j,k)>0.0) {
468 if (igrid_query ==
Vars::yvel && vel_arr(i,j,k)>0.0) {
472 if (igrid_query ==
Vars::zvel && vel_arr(i,j,k)>0.0) {
478 amrex::RealVect bias{0.,0.,0.};
479 bias[igrid_query-1] = 0.5;
481 constexpr
int dim_a = 27;
487 for(
int kk(kk_lo); kk<=kk_hi; kk++) {
488 for(
int jj(jj_lo); jj<=jj_hi; jj++) {
489 for(
int ii(ii_lo); ii<=ii_hi; ii++) {
491 if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0)) {
492 A[ll][0] = (
amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0) + bias[0] - ccent(i,j,k,0) ) * dx;
493 A[ll][1] = (
amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1) + bias[1] - ccent(i,j,k,1) ) * dy;
494 A[ll][2] = (
amrex::Real(kk) + ccent(i+ii,j+jj,k+kk,2) + bias[2] - ccent(i,j,k,2) ) * dz;
509 for(
int kk(kk_lo); kk<=kk_hi; kk++) {
510 for(
int jj(jj_lo); jj<=jj_hi; jj++) {
511 for(
int ii(ii_lo); ii<=ii_hi; ii++) {
513 if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0)) {
514 du[ll] = state(i+ii,j+jj,k+kk) - state(i,j,k);
524 for(
int jj(0); jj<AMREX_SPACEDIM; ++jj){
525 for(
int ii(0); ii<AMREX_SPACEDIM; ++ii){
531 for(
int lc(0); lc < dim_a; ++lc)
533 AtA[0][0] += A[lc][0]* A[lc][0];
534 AtA[0][1] += A[lc][0]* A[lc][1];
535 AtA[0][2] += A[lc][0]* A[lc][2];
536 AtA[1][1] += A[lc][1]* A[lc][1];
537 AtA[1][2] += A[lc][1]* A[lc][2];
538 AtA[2][2] += A[lc][2]* A[lc][2];
540 Atb[0] += A[lc][0]*du[lc];
541 Atb[1] += A[lc][1]*du[lc];
542 Atb[2] += A[lc][2]*du[lc];
546 AtA[1][0] = AtA[0][1];
547 AtA[2][0] = AtA[0][2];
548 AtA[2][1] = AtA[1][2];
551 AtA[0][0]*(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[2][1]) -
552 AtA[0][1]*(AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
553 AtA[0][2]*(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
556 Atb[0] *(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[1][2]) -
557 AtA[0][1]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) +
558 AtA[0][2]*(Atb[1] * AtA[2][1] - AtA[1][1]*Atb[2] );
564 AtA[0][0]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) -
565 Atb[0] * (AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
566 AtA[0][2]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]);
572 AtA[0][0]*(AtA[1][1]*Atb[2] - Atb[1] *AtA[1][2]) -
573 AtA[0][1]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]) +
574 Atb[0] *(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
579 return {xslope,yslope,zslope};
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:157
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > erf_calc_slopes_eb_Dirichlet([[maybe_unused]] int igrid_query, [[maybe_unused]] 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:20
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > erf_calc_slopes_eb_staggered_upwind(int igrid_query, [[maybe_unused]] int igrid_data, amrex::Real dx, amrex::Real dy, amrex::Real dz, int i, int j, int k, const amrex::Array4< const amrex::Real > &vel_arr, 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:444
AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > erf_calc_slopes_eb_staggered(int igrid_query, int igrid_data, amrex::Real dx, amrex::Real dy, amrex::Real dz, int i, int j, int k, 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:300
amrex::Real Real
Definition: ERF_ShocInterface.H:16
@ xvel
Definition: ERF_IndexDefines.H:141
@ cons
Definition: ERF_IndexDefines.H:140
@ zvel
Definition: ERF_IndexDefines.H:143
@ yvel
Definition: ERF_IndexDefines.H:142