ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Functions
ERF_EBSlopes.H File Reference
#include <AMReX_EBCellFlag.H>
#include <ERF_IndexDefines.H>
Include dependency graph for ERF_EBSlopes.H:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > erf_calc_slopes_eb ([[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)
 
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::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)
 

Function Documentation

◆ erf_calc_slopes_eb()

AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> erf_calc_slopes_eb ( [[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 
)

Function for computing the slopes in the staggered grids This function now works for only face-centered grids

Parameters
[in]igrid_queryindex of grids where the query point is located (0=CC, 1=U, 2=V, 3=W)
[in]igrid_dataindex of grids where data points are located (0=CC, 1=U, 2=V, 3=W)
[in]bcent_ebcentroid of the EB face
[in]state_ebstate at the EB face (Dirichlet data)
[in]statearray of state variable
[in]ccentarray of cell centroid coordinates
30 {
31  AMREX_ASSERT( igrid_query == igrid_data );
32 
33  // Fitting stencil
34 
35  amrex::RealVect bias{0.,0.,0.};
36  // int ii_lo = -2;
37  // int jj_lo = -2;
38  // int kk_lo = -2;
39  // int ii_hi = 2;
40  // int jj_hi = 2;
41  // int kk_hi = 2;
42  // int dim_a = 125;
43 
44  // const int dim_a = 48; // This is the maximum size of the stencil
45  constexpr int dim_a = 27;
46 
47  int ii_lo = -1;
48  int jj_lo = -1;
49  int kk_lo = -1;
50  int ii_hi = 1;
51  int jj_hi = 1;
52  int kk_hi = 1;
53 
54  amrex::Real A[dim_a][AMREX_SPACEDIM];
55 
56  // Array of distances from the query point
57  int ll=0;
58  for(int kk(kk_lo); kk<=kk_hi; kk++) {
59  for(int jj(jj_lo); jj<=jj_hi; jj++) {
60  for(int ii(ii_lo); ii<=ii_hi; ii++) {
61 
62  if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0)) {
63  A[ll][0] = ( amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0) + bias[0] - bcent_eb[0] ) * dx;
64  A[ll][1] = ( amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1) + bias[1] - bcent_eb[1] ) * dy;
65  A[ll][2] = ( amrex::Real(kk) + ccent(i+ii,j+jj,k+kk,2) + bias[2] - bcent_eb[2] ) * dz;
66  } else {
67  A[ll][0] = amrex::Real(0.0);
68  A[ll][1] = amrex::Real(0.0);
69  A[ll][2] = amrex::Real(0.0);
70  }
71  ll++;
72  }}}
73 
74  //
75  // Calculate the slopes given the matrix A (See amrex_calc_slopes_eb_given_A)
76  //
77  amrex::Real du[dim_a];
78 
79  ll=0;
80  for(int kk(kk_lo); kk<=kk_hi; kk++) {
81  for(int jj(jj_lo); jj<=jj_hi; jj++) {
82  for(int ii(ii_lo); ii<=ii_hi; ii++) {
83 
84  if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0)) {
85  du[ll] = state(i+ii,j+jj,k+kk) - state_eb;
86  } else {
87  du[ll] = amrex::Real(0.0);
88  }
89  ll++;
90  }}}
91 
92  amrex::Real AtA[AMREX_SPACEDIM][AMREX_SPACEDIM];
93  amrex::Real Atb[AMREX_SPACEDIM];
94 
95  for(int jj(0); jj<AMREX_SPACEDIM; ++jj){
96  for(int ii(0); ii<AMREX_SPACEDIM; ++ii){ // NOLINT(modernize-loop-convert)
97  AtA[jj][ii] = amrex::Real(0.0);
98  }
99  Atb[jj] = amrex::Real(0.0);
100  }
101 
102  for(int lc(0); lc < dim_a; ++lc)
103  {
104  AtA[0][0] += A[lc][0]* A[lc][0];
105  AtA[0][1] += A[lc][0]* A[lc][1];
106  AtA[0][2] += A[lc][0]* A[lc][2];
107  AtA[1][1] += A[lc][1]* A[lc][1];
108  AtA[1][2] += A[lc][1]* A[lc][2];
109  AtA[2][2] += A[lc][2]* A[lc][2];
110 
111  Atb[0] += A[lc][0]*du[lc];
112  Atb[1] += A[lc][1]*du[lc];
113  Atb[2] += A[lc][2]*du[lc];
114  }
115 
116  // Fill in symmetric
117  AtA[1][0] = AtA[0][1];
118  AtA[2][0] = AtA[0][2];
119  AtA[2][1] = AtA[1][2];
120 
121  amrex::Real detAtA =
122  AtA[0][0]*(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[2][1]) -
123  AtA[0][1]*(AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
124  AtA[0][2]*(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
125 
126  amrex::Real detAtA_x =
127  Atb[0] *(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[1][2]) -
128  AtA[0][1]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) +
129  AtA[0][2]*(Atb[1] * AtA[2][1] - AtA[1][1]*Atb[2] );
130 
131  // Slope at centroid of (i,j,k)
132  amrex::Real xslope = detAtA_x / detAtA;
133 
134  amrex::Real detAtA_y =
135  AtA[0][0]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) -
136  Atb[0] * (AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
137  AtA[0][2]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]);
138 
139  // Slope at centroid of (i,j,k)
140  amrex::Real yslope = detAtA_y / detAtA;
141 
142  amrex::Real detAtA_z =
143  AtA[0][0]*(AtA[1][1]*Atb[2] - Atb[1] *AtA[1][2]) -
144  AtA[0][1]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]) +
145  Atb[0] *(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
146 
147  // Slope at centroid of (i,j,k)
148  amrex::Real zslope = detAtA_z / detAtA;
149 
150  return {xslope,yslope,zslope};
151 }

Referenced by DiffusionSrcForMom_EB().

Here is the caller graph for this function:

◆ erf_calc_slopes_eb_staggered()

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::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 
)

Function for computing the slopes in the staggered grids This function now works for only face-centered grids

Parameters
[in]igrid_queryindex of grids where the query point is located (0=CC, 1=U, 2=V, 3=W)
[in]igrid_dataindex of grids where data points are located (0=CC, 1=U, 2=V, 3=W)
[in]bcent_ebcentroid of the EB face
[in]state_ebstate at the EB face (Dirichlet data)
[in]statearray of state variable
[in]ccentarray of cell centroid coordinates
177 {
178  // Fitting stencil
179 
180  amrex::RealVect bias{0.,0.,0.};
181 
182  constexpr int dim_a = 48;
183 
184  int ii_lo = -1;
185  int jj_lo = -1;
186  int kk_lo = -1;
187  int ii_hi = 1;
188  int jj_hi = 1;
189  int kk_hi = 1;
190 
191  // Determine which directions are involved
192  bool x_involved = (igrid_query == Vars::xvel || igrid_data == Vars::xvel);
193  bool y_involved = (igrid_query == Vars::yvel || igrid_data == Vars::yvel);
194  bool z_involved = (igrid_query == Vars::zvel || igrid_data == Vars::zvel);
195 
196  if (x_involved) ii_hi = 2;
197  if (y_involved) jj_hi = 2;
198  if (z_involved) kk_hi = 2;
199 
200  // Set bias in the index space between the two staggered grids
201  bias[igrid_query-1] = 0.5;
202  bias[igrid_data-1] = -0.5;
203 
204  amrex::Real A[dim_a][AMREX_SPACEDIM];
205 
206  // Array of distances from the query point
207  int ll=0;
208  for(int kk(kk_lo); kk<=kk_hi; kk++) {
209  for(int jj(jj_lo); jj<=jj_hi; jj++) {
210  for(int ii(ii_lo); ii<=ii_hi; ii++) {
211 
212  if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0)) {
213  A[ll][0] = ( amrex::Real(ii) + ccent(i+ii,j+jj,k+kk,0) + bias[0] - bcent_eb[0] ) * dx;
214  A[ll][1] = ( amrex::Real(jj) + ccent(i+ii,j+jj,k+kk,1) + bias[1] - bcent_eb[1] ) * dy;
215  A[ll][2] = ( amrex::Real(kk) + ccent(i+ii,j+jj,k+kk,2) + bias[2] - bcent_eb[2] ) * dz;
216  } else {
217  A[ll][0] = amrex::Real(0.0);
218  A[ll][1] = amrex::Real(0.0);
219  A[ll][2] = amrex::Real(0.0);
220  }
221  ll++;
222  }}}
223 
224  //
225  // Calculate the slopes given the matrix A (See amrex_calc_slopes_eb_given_A)
226  //
227  amrex::Real du[dim_a];
228 
229  ll=0;
230  for(int kk(kk_lo); kk<=kk_hi; kk++) {
231  for(int jj(jj_lo); jj<=jj_hi; jj++) {
232  for(int ii(ii_lo); ii<=ii_hi; ii++) {
233 
234  if (!flag(i+ii,j+jj,k+kk).isCovered() && !(ii==0 && jj==0 && kk==0)) {
235  du[ll] = state(i+ii,j+jj,k+kk) - state_eb;
236  } else {
237  du[ll] = amrex::Real(0.0);
238  }
239  ll++;
240  }}}
241 
242  amrex::Real AtA[AMREX_SPACEDIM][AMREX_SPACEDIM];
243  amrex::Real Atb[AMREX_SPACEDIM];
244 
245  for(int jj(0); jj<AMREX_SPACEDIM; ++jj){
246  for(int ii(0); ii<AMREX_SPACEDIM; ++ii){ // NOLINT(modernize-loop-convert)
247  AtA[jj][ii] = amrex::Real(0.0);
248  }
249  Atb[jj] = amrex::Real(0.0);
250  }
251 
252  for(int lc(0); lc < dim_a; ++lc)
253  {
254  AtA[0][0] += A[lc][0]* A[lc][0];
255  AtA[0][1] += A[lc][0]* A[lc][1];
256  AtA[0][2] += A[lc][0]* A[lc][2];
257  AtA[1][1] += A[lc][1]* A[lc][1];
258  AtA[1][2] += A[lc][1]* A[lc][2];
259  AtA[2][2] += A[lc][2]* A[lc][2];
260 
261  Atb[0] += A[lc][0]*du[lc];
262  Atb[1] += A[lc][1]*du[lc];
263  Atb[2] += A[lc][2]*du[lc];
264  }
265 
266  // Fill in symmetric
267  AtA[1][0] = AtA[0][1];
268  AtA[2][0] = AtA[0][2];
269  AtA[2][1] = AtA[1][2];
270 
271  amrex::Real detAtA =
272  AtA[0][0]*(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[2][1]) -
273  AtA[0][1]*(AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
274  AtA[0][2]*(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
275 
276  amrex::Real detAtA_x =
277  Atb[0] *(AtA[1][1]*AtA[2][2] - AtA[1][2]*AtA[1][2]) -
278  AtA[0][1]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) +
279  AtA[0][2]*(Atb[1] * AtA[2][1] - AtA[1][1]*Atb[2] );
280 
281  // Slope at centroid of (i,j,k)
282  amrex::Real xslope = detAtA_x / detAtA;
283 
284  amrex::Real detAtA_y =
285  AtA[0][0]*(Atb[1] * AtA[2][2] - AtA[1][2]*Atb[2] ) -
286  Atb[0] * (AtA[1][0]*AtA[2][2] - AtA[1][2]*AtA[2][0]) +
287  AtA[0][2]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]);
288 
289  // Slope at centroid of (i,j,k)
290  amrex::Real yslope = detAtA_y / detAtA;
291 
292  amrex::Real detAtA_z =
293  AtA[0][0]*(AtA[1][1]*Atb[2] - Atb[1] *AtA[1][2]) -
294  AtA[0][1]*(AtA[1][0]*Atb[2] - Atb[1] *AtA[2][0]) +
295  Atb[0] *(AtA[1][0]*AtA[2][1] - AtA[1][1]*AtA[2][0]);
296 
297  // Slope at centroid of (i,j,k)
298  amrex::Real zslope = detAtA_z / detAtA;
299 
300  return {xslope,yslope,zslope};
301 }
Vars::xvel
@ xvel
Definition: ERF_IndexDefines.H:141
Vars::zvel
@ zvel
Definition: ERF_IndexDefines.H:143
Vars::yvel
@ yvel
Definition: ERF_IndexDefines.H:142

Referenced by DiffusionSrcForMom_EB().

Here is the caller graph for this function: