Block-Structured AMR Software Framework
AMReX_EB_LeastSquares_3D_K.H
Go to the documentation of this file.
1 #ifndef AMREX_MLEB_LEASTSQUARES_K_H_
2 #define AMREX_MLEB_LEASTSQUARES_K_H_
3 #include <AMReX_Config.H>
4 #include <AMReX_EBCellFlag.H>
5 
6 namespace amrex {
7 
8 // This is the 3D version (i.e. for 10x10 (A^T A) matrix)
9 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
10 void decomp_chol_np10(Array2D<Real,0,9,0,9>& aa)
11 {
12  int neq = 10;
13 
14  Real p[10];
15  Real sum1;
16  int ising;
17 
18  for (int ii = 0; ii < neq; ii++)
19  {
20  ising = 0;
21 
22  for (int jj = ii; jj < neq; jj++)
23  {
24  sum1 = aa(ii,jj);
25 
26  for (int kk = ii-1; kk >= 0; kk--)
27  {
28  sum1 = sum1 - aa(ii,kk)*aa(jj,kk);
29  }
30 
31  if (ii == jj)
32  {
33  if (sum1 <= Real(0.0))
34  {
35  p[ii] = Real(0.0);
36  ising = 1;
37  } else {
38  p[ii] = std::sqrt(sum1);
39  }
40  } else {
41  if (ising == 0) {
42  aa(jj,ii) = sum1 / p[ii];
43  } else {
44  aa(jj,ii) = Real(0.0);
45  }
46  }
47  }
48  }
49 
50  for (int ii = 0; ii < neq; ii++)
51  {
52  for (int jj = ii+1; jj < neq; jj++)
53  {
54  aa(ii,jj) = Real(0.0); // Zero upper triangle
55  aa(ii,jj) = aa(jj,ii); // Zero upper triangle
56  }
57 
58  aa(ii,ii) = p[ii];
59  }
60 }
61 
62 // This is the 3D version (i.e. for 10x10 (A^T A) matrix)
63 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
64 void cholsol_np10(Array2D<Real,0,35,0,9>& Amatrix, Array1D<Real,0,9>& b)
65 {
66  int neq = 10;
67 
68  Array2D<Real,0,9,0,9> AtA;
69 
70  for (int irow = 0; irow < neq; irow++) {
71  for (int icol = 0; icol < neq; icol++) {
72  AtA(irow,icol) = Real(0.0);
73  }
74  }
75 
76  for (int irow = 0; irow < 36; irow++)
77  {
78  AtA(0,0) += Amatrix(irow,0)*Amatrix(irow,0); // e^T e
79  AtA(0,1) += Amatrix(irow,0)*Amatrix(irow,1); // e^T x
80  AtA(0,2) += Amatrix(irow,0)*Amatrix(irow,2); // e^T y
81  AtA(0,3) += Amatrix(irow,0)*Amatrix(irow,3); // e^T z
82  AtA(0,4) += Amatrix(irow,0)*Amatrix(irow,4); // e^T x^2
83  AtA(0,5) += Amatrix(irow,0)*Amatrix(irow,5); // e^T x*y
84  AtA(0,6) += Amatrix(irow,0)*Amatrix(irow,6); // e^T y^2
85  AtA(0,7) += Amatrix(irow,0)*Amatrix(irow,7); // e^T x*z
86  AtA(0,8) += Amatrix(irow,0)*Amatrix(irow,8); // e^T y*z
87  AtA(0,9) += Amatrix(irow,0)*Amatrix(irow,9); // e^T z^2
88 
89  AtA(1,1) += Amatrix(irow,1)*Amatrix(irow,1); // x^T x
90  AtA(1,2) += Amatrix(irow,1)*Amatrix(irow,2); // x^T y
91  AtA(1,3) += Amatrix(irow,1)*Amatrix(irow,3); // x^T y
92  AtA(1,4) += Amatrix(irow,1)*Amatrix(irow,4); // x^T (x^2)
93  AtA(1,5) += Amatrix(irow,1)*Amatrix(irow,5); // x^T (xy)
94  AtA(1,6) += Amatrix(irow,1)*Amatrix(irow,6); // x^T (y^2)
95  AtA(1,7) += Amatrix(irow,1)*Amatrix(irow,7); // x^T x*z
96  AtA(1,8) += Amatrix(irow,1)*Amatrix(irow,8); // x^T y*z
97  AtA(1,9) += Amatrix(irow,1)*Amatrix(irow,9); // x^T z^2
98 
99  AtA(2,2) += Amatrix(irow,2)*Amatrix(irow,2); // y^T y
100  AtA(2,3) += Amatrix(irow,2)*Amatrix(irow,3); // y^T z
101  AtA(2,4) += Amatrix(irow,2)*Amatrix(irow,4); // y^T (x^2)
102  AtA(2,5) += Amatrix(irow,2)*Amatrix(irow,5); // y^T (xy)
103  AtA(2,6) += Amatrix(irow,2)*Amatrix(irow,6); // y^T (y^2)
104  AtA(2,7) += Amatrix(irow,2)*Amatrix(irow,7); // y^T x*z
105  AtA(2,8) += Amatrix(irow,2)*Amatrix(irow,8); // y^T y*z
106  AtA(2,9) += Amatrix(irow,2)*Amatrix(irow,9); // y^T z^2
107 
108  AtA(3,3) += Amatrix(irow,3)*Amatrix(irow,3); // z^T z
109  AtA(3,4) += Amatrix(irow,3)*Amatrix(irow,4); // z^T (x^2)
110  AtA(3,5) += Amatrix(irow,3)*Amatrix(irow,5); // z^T (xy)
111  AtA(3,6) += Amatrix(irow,3)*Amatrix(irow,6); // z^T (y^2)
112  AtA(3,7) += Amatrix(irow,3)*Amatrix(irow,7); // z^T x*z
113  AtA(3,8) += Amatrix(irow,3)*Amatrix(irow,8); // z^T y*z
114  AtA(3,9) += Amatrix(irow,3)*Amatrix(irow,9); // z^T z^2
115 
116  AtA(4,4) += Amatrix(irow,4)*Amatrix(irow,4); // (x^2)^T (x^2)
117  AtA(4,5) += Amatrix(irow,4)*Amatrix(irow,5); // (x^2)^T (xy)
118  AtA(4,6) += Amatrix(irow,4)*Amatrix(irow,6); // (x^2)^T (y^2)
119  AtA(4,7) += Amatrix(irow,4)*Amatrix(irow,7); // (x^2)^T x*z
120  AtA(4,8) += Amatrix(irow,4)*Amatrix(irow,8); // (x^2)^T y*z
121  AtA(4,9) += Amatrix(irow,4)*Amatrix(irow,9); // (x^2)^T z^2
122 
123  AtA(5,5) += Amatrix(irow,5)*Amatrix(irow,5); // (xy)^T (xy)
124  AtA(5,6) += Amatrix(irow,5)*Amatrix(irow,6); // (xy)^T (y^2)
125  AtA(5,7) += Amatrix(irow,5)*Amatrix(irow,7); // (xy)^T x*z
126  AtA(5,8) += Amatrix(irow,5)*Amatrix(irow,8); // (xy)^T y*z
127  AtA(5,9) += Amatrix(irow,5)*Amatrix(irow,9); // (xy)^T z^2
128 
129  AtA(6,6) += Amatrix(irow,6)*Amatrix(irow,6); // (y^2)^T (y^2)
130  AtA(6,7) += Amatrix(irow,6)*Amatrix(irow,7); // (y^2)^T x*z
131  AtA(6,8) += Amatrix(irow,6)*Amatrix(irow,8); // (y^2)^T y*z
132  AtA(6,9) += Amatrix(irow,6)*Amatrix(irow,9); // (y^2)^T z^2
133 
134  AtA(7,7) += Amatrix(irow,7)*Amatrix(irow,7); // (xz)^T x*z
135  AtA(7,8) += Amatrix(irow,7)*Amatrix(irow,8); // (xz)^T y*z
136  AtA(7,9) += Amatrix(irow,7)*Amatrix(irow,9); // (xz)^T z^2
137 
138  AtA(8,8) += Amatrix(irow,8)*Amatrix(irow,8); // (yz)^T y*z
139  AtA(8,9) += Amatrix(irow,8)*Amatrix(irow,9); // (yz)^T z^2
140 
141  AtA(9,9) += Amatrix(irow,9)*Amatrix(irow,9); // (z^2)^T z^2
142  }
143 
144  for (int irow = 0; irow < neq-1; irow++) {
145  for (int icol = irow+1; icol < neq; icol++) {
146  AtA(icol,irow) = AtA(irow,icol);
147  }
148  }
149 
150  decomp_chol_np10(AtA);
151 
152  if (AtA(0,0) > 0.) {
153  b(0) = b(0) / AtA(0,0);
154  } else {
155  b(0) = 0.;
156  }
157 
158  for (int ii = 1; ii < neq; ii++)
159  {
160  if (AtA(ii,ii) > 0.)
161  {
162  for (int jj = 0; jj < ii; jj++) {
163  b(ii) = b(ii) - AtA(ii,jj)*b(jj);
164  }
165 
166  b(ii) = b(ii) / AtA(ii,ii);
167  }
168  else
169  {
170  b(ii) = Real(0.0);
171  }
172  }
173 
174  if (AtA(neq-1,neq-1) > 0.) {
175  b(neq-1) = b(neq-1) / AtA(neq-1,neq-1);
176  } else {
177  b(neq-1) = Real(0.0);
178  }
179 
180  for (int ii = neq-2; ii >= 0; ii--)
181  {
182  if (AtA(ii,ii) > 0.)
183  {
184  for (int jj = ii+1; jj < neq; jj++) {
185  b(ii) = b(ii) - AtA(ii,jj)*b(jj);
186  }
187 
188  b(ii) = b(ii) / AtA(ii,ii);
189  }
190  else
191  {
192  b(ii) = Real(0.0);
193  }
194  }
195 }
196 
197 // This is the 3D version (i.e. for 10x10 (A^T A) matrix) but for A = 54 x 10
198 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
199 void cholsol_for_eb(Array2D<Real,0,53,0,9>& Amatrix, Array1D<Real,0,9>& b)
200 {
201  int neq = 10;
202 
203  Array2D<Real,0,9,0,9> AtA;
204 
205  for (int irow = 0; irow < neq; irow++) {
206  for (int icol = 0; icol < neq; icol++) {
207  AtA(irow,icol) = Real(0.0);
208  }
209  }
210 
211  for (int irow = 0; irow < 54; irow++)
212  {
213  AtA(0,0) += Amatrix(irow,0)*Amatrix(irow,0); // e^T e
214  AtA(0,1) += Amatrix(irow,0)*Amatrix(irow,1); // e^T x
215  AtA(0,2) += Amatrix(irow,0)*Amatrix(irow,2); // e^T y
216  AtA(0,3) += Amatrix(irow,0)*Amatrix(irow,3); // e^T z
217  AtA(0,4) += Amatrix(irow,0)*Amatrix(irow,4); // e^T x^2
218  AtA(0,5) += Amatrix(irow,0)*Amatrix(irow,5); // e^T x*y
219  AtA(0,6) += Amatrix(irow,0)*Amatrix(irow,6); // e^T y^2
220  AtA(0,7) += Amatrix(irow,0)*Amatrix(irow,7); // e^T x*z
221  AtA(0,8) += Amatrix(irow,0)*Amatrix(irow,8); // e^T y*z
222  AtA(0,9) += Amatrix(irow,0)*Amatrix(irow,9); // e^T z^2
223 
224  AtA(1,1) += Amatrix(irow,1)*Amatrix(irow,1); // x^T x
225  AtA(1,2) += Amatrix(irow,1)*Amatrix(irow,2); // x^T y
226  AtA(1,3) += Amatrix(irow,1)*Amatrix(irow,3); // x^T y
227  AtA(1,4) += Amatrix(irow,1)*Amatrix(irow,4); // x^T (x^2)
228  AtA(1,5) += Amatrix(irow,1)*Amatrix(irow,5); // x^T (xy)
229  AtA(1,6) += Amatrix(irow,1)*Amatrix(irow,6); // x^T (y^2)
230  AtA(1,7) += Amatrix(irow,1)*Amatrix(irow,7); // x^T x*z
231  AtA(1,8) += Amatrix(irow,1)*Amatrix(irow,8); // x^T y*z
232  AtA(1,9) += Amatrix(irow,1)*Amatrix(irow,9); // x^T z^2
233 
234  AtA(2,2) += Amatrix(irow,2)*Amatrix(irow,2); // y^T y
235  AtA(2,3) += Amatrix(irow,2)*Amatrix(irow,3); // y^T z
236  AtA(2,4) += Amatrix(irow,2)*Amatrix(irow,4); // y^T (x^2)
237  AtA(2,5) += Amatrix(irow,2)*Amatrix(irow,5); // y^T (xy)
238  AtA(2,6) += Amatrix(irow,2)*Amatrix(irow,6); // y^T (y^2)
239  AtA(2,7) += Amatrix(irow,2)*Amatrix(irow,7); // y^T x*z
240  AtA(2,8) += Amatrix(irow,2)*Amatrix(irow,8); // y^T y*z
241  AtA(2,9) += Amatrix(irow,2)*Amatrix(irow,9); // y^T z^2
242 
243  AtA(3,3) += Amatrix(irow,3)*Amatrix(irow,3); // z^T z
244  AtA(3,4) += Amatrix(irow,3)*Amatrix(irow,4); // z^T (x^2)
245  AtA(3,5) += Amatrix(irow,3)*Amatrix(irow,5); // z^T (xy)
246  AtA(3,6) += Amatrix(irow,3)*Amatrix(irow,6); // z^T (y^2)
247  AtA(3,7) += Amatrix(irow,3)*Amatrix(irow,7); // z^T x*z
248  AtA(3,8) += Amatrix(irow,3)*Amatrix(irow,8); // z^T y*z
249  AtA(3,9) += Amatrix(irow,3)*Amatrix(irow,9); // z^T z^2
250 
251  AtA(4,4) += Amatrix(irow,4)*Amatrix(irow,4); // (x^2)^T (x^2)
252  AtA(4,5) += Amatrix(irow,4)*Amatrix(irow,5); // (x^2)^T (xy)
253  AtA(4,6) += Amatrix(irow,4)*Amatrix(irow,6); // (x^2)^T (y^2)
254  AtA(4,7) += Amatrix(irow,4)*Amatrix(irow,7); // (x^2)^T x*z
255  AtA(4,8) += Amatrix(irow,4)*Amatrix(irow,8); // (x^2)^T y*z
256  AtA(4,9) += Amatrix(irow,4)*Amatrix(irow,9); // (x^2)^T z^2
257 
258  AtA(5,5) += Amatrix(irow,5)*Amatrix(irow,5); // (xy)^T (xy)
259  AtA(5,6) += Amatrix(irow,5)*Amatrix(irow,6); // (xy)^T (y^2)
260  AtA(5,7) += Amatrix(irow,5)*Amatrix(irow,7); // (xy)^T x*z
261  AtA(5,8) += Amatrix(irow,5)*Amatrix(irow,8); // (xy)^T y*z
262  AtA(5,9) += Amatrix(irow,5)*Amatrix(irow,9); // (xy)^T z^2
263 
264  AtA(6,6) += Amatrix(irow,6)*Amatrix(irow,6); // (y^2)^T (y^2)
265  AtA(6,7) += Amatrix(irow,6)*Amatrix(irow,7); // (y^2)^T x*z
266  AtA(6,8) += Amatrix(irow,6)*Amatrix(irow,8); // (y^2)^T y*z
267  AtA(6,9) += Amatrix(irow,6)*Amatrix(irow,9); // (y^2)^T z^2
268 
269  AtA(7,7) += Amatrix(irow,7)*Amatrix(irow,7); // (xz)^T x*z
270  AtA(7,8) += Amatrix(irow,7)*Amatrix(irow,8); // (xz)^T y*z
271  AtA(7,9) += Amatrix(irow,7)*Amatrix(irow,9); // (xz)^T z^2
272 
273  AtA(8,8) += Amatrix(irow,8)*Amatrix(irow,8); // (yz)^T y*z
274  AtA(8,9) += Amatrix(irow,8)*Amatrix(irow,9); // (yz)^T z^2
275 
276  AtA(9,9) += Amatrix(irow,9)*Amatrix(irow,9); // (z^2)^T z^2
277  }
278 
279  for (int irow = 0; irow < neq-1; irow++) {
280  for (int icol = irow+1; icol < neq; icol++) {
281  AtA(icol,irow) = AtA(irow,icol);
282  }
283  }
284 
285  decomp_chol_np10(AtA);
286 
287  if (AtA(0,0) > 0.) {
288  b(0) = b(0) / AtA(0,0);
289  } else {
290  b(0) = 0.;
291  }
292 
293  for (int ii = 1; ii < neq; ii++)
294  {
295  if (AtA(ii,ii) > 0.)
296  {
297  for (int jj = 0; jj < ii; jj++) {
298  b(ii) = b(ii) - AtA(ii,jj)*b(jj);
299  }
300 
301  b(ii) = b(ii) / AtA(ii,ii);
302  }
303  else
304  {
305  b(ii) = Real(0.0);
306  }
307  }
308 
309  if (AtA(neq-1,neq-1) > 0.) {
310  b(neq-1) = b(neq-1) / AtA(neq-1,neq-1);
311  } else {
312  b(neq-1) = Real(0.0);
313  }
314 
315  for (int ii = neq-2; ii >= 0; ii--)
316  {
317  if (AtA(ii,ii) > 0.)
318  {
319  for (int jj = ii+1; jj < neq; jj++) {
320  b(ii) = b(ii) - AtA(ii,jj)*b(jj);
321  }
322 
323  b(ii) = b(ii) / AtA(ii,ii);
324  }
325  else
326  {
327  b(ii) = Real(0.0);
328  }
329  }
330 }
331 
332 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
333 Real grad_x_of_phi_on_centroids(int i,int j,int k,int n,
334  Array4<Real const> const& phi,
335  Array4<Real const> const& phieb,
336  Array4<EBCellFlag const> const& flag,
337  Array4<Real const> const& ccent,
338  Array4<Real const> const& bcent,
339  Real& yloc_on_xface, Real& zloc_on_xface,
340  bool is_eb_dirichlet, bool is_eb_inhomog)
341 {
342  AMREX_ALWAYS_ASSERT(is_eb_dirichlet);
343  Array2D<Real,0,35,0,9> Amatrix;
344  Array1D<Real,0, 9 > rhs;
345 
346  // Order of column -- first 9 are cell centroids, second 9 are EB centroids
347 
348  for (int irow = 0; irow < 36; irow++) {
349  for (int icol = 0; icol < 10; icol++) {
350  Amatrix(irow,icol) = Real(0.0);
351  }
352  }
353 
354  // Columns 0-3: [e x y z]
355  for (int ii = i-1; ii <= i; ii++) { // Normal to face
356  for (int kk = k-1; kk <= k+1; kk++) { // Tangential to face
357  for (int jj = j-1; jj <= j+1; jj++) { // Tangential to face
358  if (!flag(ii,jj,kk).isCovered())
359  {
360  int a_ind = (jj-(j-1)) + 3*(kk-(k-1)) + 9*(ii-(i-1));
361 
362  Real x_off = static_cast<Real>(ii-i) + Real(0.5);
363  Real y_off = static_cast<Real>(jj-j);
364  Real z_off = static_cast<Real>(kk-k);
365 
366  Amatrix(a_ind,0) = Real(1.0);
367  Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0);
368  Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1) - yloc_on_xface;
369  Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2) - zloc_on_xface;
370 
371  if (!flag(ii,jj,kk).isRegular())
372  {
373  Amatrix(a_ind+18,0) = Real(1.0);
374  Amatrix(a_ind+18,1) = x_off + bcent(ii,jj,kk,0);
375  Amatrix(a_ind+18,2) = y_off + bcent(ii,jj,kk,1) - yloc_on_xface;
376  Amatrix(a_ind+18,3) = z_off + bcent(ii,jj,kk,2) - zloc_on_xface;
377  }
378  }
379  }
380  }
381  }
382 
383  // Columns 4-9 : [x*x x*y y*y x*z y*z z*z]
384  for (int irow = 0; irow < 36; irow++)
385  {
386  Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1);
387  Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2);
388  Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2);
389  Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3);
390  Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3);
391  Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3);
392  }
393 
394  // Make the RHS = A^T v
395  for (int irow = 0; irow < 10; irow++)
396  {
397  rhs(irow) = 0.;
398 
399  for (int ii = i-1; ii <= i; ii++) { // Normal to face
400  for (int kk = k-1; kk <= k+1; kk++) { // Tangential to face
401  for (int jj = j-1; jj <= j+1; jj++) { // Tangential to face
402  if (!flag(ii,jj,kk).isCovered())
403  {
404  int a_ind = (jj-(j-1)) + 3*(kk-(k-1)) + 9*(ii-(i-1));
405 
406  rhs(irow) += Amatrix(a_ind,irow)* phi(ii,jj,kk,n);
407 
408  if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog) {
409  rhs(irow) += Amatrix(a_ind+18,irow)*phieb(ii,jj,kk,n);
410  }
411  }
412  }
413  }
414  }
415  }
416 
417  cholsol_np10(Amatrix, rhs);
418 
419  return rhs(1);
420 }
421 
422 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
423 Real grad_y_of_phi_on_centroids(int i,int j,int k,int n,
424  Array4<Real const> const& phi,
425  Array4<Real const> const& phieb,
426  Array4<EBCellFlag const> const& flag,
427  Array4<Real const> const& ccent,
428  Array4<Real const> const& bcent,
429  Real& xloc_on_yface, Real& zloc_on_yface,
430  bool is_eb_dirichlet, bool is_eb_inhomog)
431 {
432  AMREX_ALWAYS_ASSERT(is_eb_dirichlet);
433  Array2D<Real,0,35,0,9> Amatrix;
434  Array1D<Real,0, 9 > rhs;
435 
436  // Order of column -- first 9 are cell centroids, second 9 are EB centroids
437 
438  for (int irow = 0; irow < 36; irow++) {
439  for (int icol = 0; icol < 10; icol++) {
440  Amatrix(irow,icol) = Real(0.0);
441  }
442  }
443 
444  // Columns 0-2: [e x y]
445  for (int jj = j-1; jj <= j; jj++) { // Normal to face
446  for (int kk = k-1; kk <= k+1; kk++) { // Tangential to face
447  for (int ii = i-1; ii <= i+1; ii++) { // Tangential to face
448  if (!flag(ii,jj,kk).isCovered())
449  {
450  int a_ind = (ii-(i-1)) + 3*(kk-(k-1)) + 9*(jj-(j-1));
451 
452  Real x_off = static_cast<Real>(ii-i);
453  Real y_off = static_cast<Real>(jj-j) + Real(0.5);
454  Real z_off = static_cast<Real>(kk-k);
455 
456  Amatrix(a_ind,0) = Real(1.0);
457  Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0) - xloc_on_yface;
458  Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1);
459  Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2) - zloc_on_yface;
460 
461  if (!flag(ii,jj,kk).isRegular())
462  {
463  Amatrix(a_ind+18,0) = Real(1.0);
464  Amatrix(a_ind+18,1) = x_off + bcent(ii,jj,kk,0) - xloc_on_yface;
465  Amatrix(a_ind+18,2) = y_off + bcent(ii,jj,kk,1);
466  Amatrix(a_ind+18,3) = z_off + bcent(ii,jj,kk,2) - zloc_on_yface;
467  }
468  }
469  }
470  }
471  }
472 
473  // Columns 4-9 : [x*x x*y y*y x*z y*z z*z]
474  for (int irow = 0; irow < 36; irow++)
475  {
476  Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1);
477  Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2);
478  Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2);
479  Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3);
480  Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3);
481  Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3);
482  }
483 
484  // Make the RHS = A^T v
485  for (int irow = 0; irow < 10; irow++)
486  {
487  rhs(irow) = 0.;
488 
489  for (int jj = j-1; jj <= j; jj++) { // Normal to face
490  for (int kk = k-1; kk <= k+1; kk++) { // Tangential to face
491  for (int ii = i-1; ii <= i+1; ii++) { // Tangential to face
492  if (!flag(ii,jj,kk).isCovered())
493  {
494  int a_ind = (ii-(i-1)) + 3*(kk-(k-1)) + 9*(jj-(j-1));
495 
496  rhs(irow) += Amatrix(a_ind,irow)* phi(ii,jj,kk,n);
497 
498  if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog) {
499  rhs(irow) += Amatrix(a_ind+18,irow)*phieb(ii,jj,kk,n);
500  }
501  }
502  }
503  }
504  }
505  }
506 
507  cholsol_np10(Amatrix, rhs);
508 
509  return rhs(2);
510 }
511 
512 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
513 Real grad_z_of_phi_on_centroids(int i,int j,int k,int n,
514  Array4<Real const> const& phi,
515  Array4<Real const> const& phieb,
516  Array4<EBCellFlag const> const& flag,
517  Array4<Real const> const& ccent,
518  Array4<Real const> const& bcent,
519  Real& xloc_on_zface, Real& yloc_on_zface,
520  bool is_eb_dirichlet, bool is_eb_inhomog)
521 {
522  AMREX_ALWAYS_ASSERT(is_eb_dirichlet);
523  Array2D<Real,0,35,0,9> Amatrix;
524  Array1D<Real,0, 9 > rhs;
525 
526  // Order of column -- first 9 are cell centroids, second 9 are EB centroids
527 
528  for (int irow = 0; irow < 36; irow++) {
529  for (int icol = 0; icol < 10; icol++) {
530  Amatrix(irow,icol) = Real(0.0);
531  }
532  }
533 
534  // Columns 0-3: [e x y z]
535  for (int kk = k-1; kk <= k; kk++) // Normal to face
536  {
537  for (int jj = j-1; jj <= j+1; jj++) { // Tangential to face
538  for (int ii = i-1; ii <= i+1; ii++) // Tangential to face
539  {
540  if (!flag(ii,jj,kk).isCovered())
541  {
542  int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-1));
543 
544  Real x_off = static_cast<Real>(ii-i);
545  Real y_off = static_cast<Real>(jj-j);
546  Real z_off = static_cast<Real>(kk-k) + Real(0.5);
547 
548  Amatrix(a_ind,0) = Real(1.0);
549  Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0) - xloc_on_zface;
550  Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1) - yloc_on_zface;
551  Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2);
552 
553  if (!flag(ii,jj,kk).isRegular())
554  {
555  Amatrix(a_ind+18,0) = Real(1.0);
556  Amatrix(a_ind+18,1) = x_off + bcent(ii,jj,kk,0) - xloc_on_zface;
557  Amatrix(a_ind+18,2) = y_off + bcent(ii,jj,kk,1) - yloc_on_zface;
558  Amatrix(a_ind+18,3) = z_off + bcent(ii,jj,kk,2);
559  }
560  }
561  }
562  }
563  }
564 
565  // Columns 4-9 : [x*x x*y y*y x*z y*z z*z]
566  for (int irow = 0; irow < 36; irow++)
567  {
568  Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1);
569  Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2);
570  Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2);
571  Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3);
572  Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3);
573  Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3);
574  }
575 
576  // Make the RHS = A^T v
577  for (int irow = 0; irow < 10; irow++)
578  {
579  rhs(irow) = 0.;
580 
581  for (int kk = k-1; kk <= k; kk++) { // Normal to face
582  for (int jj = j-1; jj <= j+1; jj++) { // Tangential to face
583  for (int ii = i-1; ii <= i+1; ii++) { // Tangential to face
584  if (!flag(ii,jj,kk).isCovered())
585  {
586  int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-1));
587 
588  rhs(irow) += Amatrix(a_ind,irow)* phi(ii,jj,kk,n);
589 
590  if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog) {
591  rhs(irow) += Amatrix(a_ind+18,irow)*phieb(ii,jj,kk,n);
592  }
593  }
594  }
595  }
596  }
597  }
598 
599  cholsol_np10(Amatrix, rhs);
600 
601  return rhs(3);
602 }
603 
604 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
605 Real grad_eb_of_phi_on_centroids(int i,int j,int k,int n,
606  Array4<Real const> const& phi,
607  Array4<Real const> const& phieb,
608  Array4<EBCellFlag const> const& flag,
609  Array4<Real const> const& ccent,
610  Array4<Real const> const& bcent,
611  Real& nrmx, Real& nrmy, Real& nrmz,
612  bool is_eb_inhomog)
613 {
614  Array2D<Real,0,53,0,9> Amatrix;
615  Array1D<Real,0, 9 > rhs;
616 
617  // Order of column -- first 27 are cell centroids, second 27 are EB centroids
618 
619  for (int irow = 0; irow < 54; irow++) {
620  for (int icol = 0; icol < 10; icol++) {
621  Amatrix(irow,icol) = Real(0.0);
622  }
623  }
624 
625  // Columns 0-3: [e x y z]
626  for (int kk = k-1; kk <= k+1; kk++) {
627  for (int jj = j-1; jj <= j+1; jj++) {
628  for (int ii = i-1; ii <= i+1; ii++)
629  {
630  if (!flag(ii,jj,kk).isCovered())
631  {
632  int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-1));
633 
634  Real x_off = static_cast<Real>(ii-i) - bcent(i,j,k,0);
635  Real y_off = static_cast<Real>(jj-j) - bcent(i,j,k,1);
636  Real z_off = static_cast<Real>(kk-k) - bcent(i,j,k,2);
637 
638  Amatrix(a_ind,0) = Real(1.0);
639  Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0);
640  Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1);
641  Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2);
642 
643  if (!flag(ii,jj,kk).isRegular())
644  {
645  Amatrix(a_ind+27,0) = Real(1.0);
646  Amatrix(a_ind+27,1) = x_off + bcent(ii,jj,kk,0);
647  Amatrix(a_ind+27,2) = y_off + bcent(ii,jj,kk,1);
648  Amatrix(a_ind+27,3) = z_off + bcent(ii,jj,kk,2);
649  }
650  }
651  }
652  }
653  }
654 
655  // Columns 4-9 : [x*x x*y y*y x*z y*z z*z]
656  for (int irow = 0; irow < 54; irow++)
657  {
658  Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1);
659  Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2);
660  Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2);
661  Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3);
662  Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3);
663  Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3);
664  }
665 
666  // Make the RHS = A^T v
667  for (int irow = 0; irow < 10; irow++)
668  {
669  rhs(irow) = 0.;
670 
671  for (int kk = k-1; kk <= k+1; kk++) {
672  for (int jj = j-1; jj <= j+1; jj++) {
673  for (int ii = i-1; ii <= i+1; ii++) {
674  if (!flag(ii,jj,kk).isCovered())
675  {
676  int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-1));
677 
678  rhs(irow) += Amatrix(a_ind,irow)* phi(ii,jj,kk,n);
679 
680  if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog) {
681  rhs(irow) += Amatrix(a_ind+27,irow)*phieb(ii,jj,kk,n);
682  }
683  }
684  }
685  }
686  }
687  }
688 
689  cholsol_for_eb(Amatrix, rhs);
690 
691  Real dphidn = rhs(1)*nrmx + rhs(2)*nrmy + rhs(3)*nrmz;
692  return dphidn;
693 }
694 
695 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
696 Real grad_x_of_phi_on_centroids_extdir(int i,int j,int k,int n,
697  Array4<Real const> const& phi,
698  Array4<Real const> const& phieb,
699  Array4<EBCellFlag const> const& flag,
700  Array4<Real const> const& ccent,
701  Array4<Real const> const& bcent,
702  Array4<Real const> const& vfrac,
703  Real& yloc_on_xface, Real& zloc_on_xface,
704  bool is_eb_dirichlet, bool is_eb_inhomog,
705  const bool on_x_face, const int domlo_x, const int domhi_x,
706  const bool on_y_face, const int domlo_y, const int domhi_y,
707  const bool on_z_face, const int domlo_z, const int domhi_z)
708 {
709  AMREX_ALWAYS_ASSERT(is_eb_dirichlet);
710 
711  Array2D<Real,0,53,0,9> Amatrix;
712  Array1D<Real,0, 9 > rhs;
713 
714  // Order of column -- first 9 are cell centroids, second 9 are EB centroids
715 
716  for (int irow = 0; irow < 54; irow++) {
717  for (int icol = 0; icol < 10; icol++) {
718  Amatrix(irow,icol) = Real(0.0);
719  }
720  }
721 
722  const int im = (i > domhi_x) ? 2 : 1;
723  const int ip = 2 - im;
724 
725  // Columns 0-3: [e x y z]
726  for (int ii = i-im; ii <= i+ip; ii++) { // Normal to face
727  for (int kk = k-1; kk <= k+1; kk++) { // Tangential to face
728  for (int jj = j-1; jj <= j+1; jj++) // Tangential to face
729  {
730  // Don't include corner cells. Could make this even more strict
731  // by removing the on_?_face restrictions.
732  if (((on_x_face && ii < domlo_x) && (on_y_face && jj < domlo_y)) ||
733  ((on_x_face && ii < domlo_x) && (on_y_face && jj > domhi_y)) ||
734  ((on_x_face && ii > domhi_x) && (on_y_face && jj < domlo_y)) ||
735  ((on_x_face && ii > domhi_x) && (on_y_face && jj > domhi_y)) ||
736  ((on_y_face && jj < domlo_y) && (on_z_face && kk < domlo_z)) ||
737  ((on_y_face && jj < domlo_y) && (on_z_face && kk > domhi_z)) ||
738  ((on_y_face && jj > domhi_y) && (on_z_face && kk < domlo_z)) ||
739  ((on_y_face && jj > domhi_y) && (on_z_face && kk > domhi_z)) ||
740  ((on_x_face && ii < domlo_x) && (on_z_face && kk < domlo_z)) ||
741  ((on_x_face && ii < domlo_x) && (on_z_face && kk > domhi_z)) ||
742  ((on_x_face && ii > domhi_x) && (on_z_face && kk < domlo_z)) ||
743  ((on_x_face && ii > domhi_x) && (on_z_face && kk > domhi_z))) {
744  continue;
745  }
746 
747  if ( !phi.contains(ii,jj,kk) ) {
748  continue;
749  }
750 
751  if (!flag(ii,jj,kk).isCovered())
752  {
753 
754  int a_ind = (jj-(j-1)) + 3*(kk-(k-1)) + 9*(ii-(i-im));
755 
756  Real x_off = static_cast<Real>(ii-i) + Real(0.5);
757  Real y_off = static_cast<Real>(jj-j);
758  Real z_off = static_cast<Real>(kk-k);
759 
760  if(on_x_face){
761  if (ii < domlo_x && (vfrac(ii+1,jj,kk) != Real(1.0) || vfrac(ii+2,jj,kk) != Real(1.0)) ) {
762  continue;
763  }
764  if (ii > domhi_x && (vfrac(ii-1,jj,kk) != Real(1.0) || vfrac(ii-2,jj,kk) != Real(1.0))) {
765  continue;
766  }
767  }
768 
769  if(on_y_face){
770  if (jj < domlo_y && (vfrac(ii,jj+1,kk) != Real(1.0) || vfrac(ii,jj+2,kk) != Real(1.0)) ) {
771  continue;
772  }
773  if (jj > domhi_y && (vfrac(ii,jj-1,kk) != Real(1.0) || vfrac(ii,jj-2,kk) != Real(1.0)) ) {
774  continue;
775  }
776  }
777 
778  if(on_z_face){
779  if (kk < domlo_z && (vfrac(ii,jj,kk+1) != Real(1.0) || vfrac(ii,jj,kk+2) != Real(1.0)) ) {
780  continue;
781  }
782  if (kk > domhi_z && (vfrac(ii,jj,kk-1) != Real(1.0) || vfrac(ii,jj,kk-2) != Real(1.0)) ) {
783  continue;
784  }
785  }
786 
787  Amatrix(a_ind,0) = Real(1.0);
788  Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0);
789  Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1) - yloc_on_xface;
790  Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2) - zloc_on_xface;
791 
792  // Add in information about the location of the EB. Exclude
793  // EBs that are outside the domain.
794  if (flag(ii,jj,kk).isSingleValued() &&
795  domlo_x <= ii && ii <= domhi_x &&
796  domlo_y <= jj && jj <= domhi_y &&
797  domlo_z <= kk && kk <= domhi_z)
798  {
799  Amatrix(a_ind+27,0) = Real(1.0);
800  Amatrix(a_ind+27,1) = x_off + bcent(ii,jj,kk,0);
801  Amatrix(a_ind+27,2) = y_off + bcent(ii,jj,kk,1) - yloc_on_xface;
802  Amatrix(a_ind+27,3) = z_off + bcent(ii,jj,kk,2) - zloc_on_xface;
803  }
804  }
805  }
806  }
807  }
808 
809  // Columns 4-9 : [x*x x*y y*y x*z y*z z*z]
810  for (int irow = 0; irow < 54; irow++)
811  {
812  Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1);
813  Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2);
814  Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2);
815  Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3);
816  Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3);
817  Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3);
818  }
819 
820  // Make the RHS = A^T v
821  for (int irow = 0; irow < 10; irow++)
822  {
823  rhs(irow) = 0.;
824 
825  for (int ii = i-im; ii <= i+ip; ii++) { // Normal to face
826  for (int kk = k-1; kk <= k+1; kk++) { // Tangential to face
827  for (int jj = j-1; jj <= j+1; jj++) { // Tangential to face
828  if ( !phi.contains(ii,jj,kk) ) {
829  continue;
830  }
831 
832  if (!flag(ii,jj,kk).isCovered())
833  {
834  int a_ind = (jj-(j-1)) + 3*(kk-(k-1)) + 9*(ii-(i-im));
835  Real phi_val = Amatrix(a_ind,0)*phi(ii,jj,kk,n);
836 
837  rhs(irow) += Amatrix(a_ind,irow)* phi_val;
838 
839  if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog && Amatrix(a_ind+27,irow) != Real(0.0)) {
840  rhs(irow) += Amatrix(a_ind+27,irow)*phieb(ii,jj,kk,n);
841  }
842  }
843  }
844  }
845  }
846  }
847 
848  cholsol_for_eb(Amatrix, rhs);
849 
850  return rhs(1);
851 }
852 
853 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
854 Real grad_y_of_phi_on_centroids_extdir(int i,int j,int k,int n,
855  Array4<Real const> const& phi,
856  Array4<Real const> const& phieb,
857  Array4<EBCellFlag const> const& flag,
858  Array4<Real const> const& ccent,
859  Array4<Real const> const& bcent,
860  Array4<Real const> const& vfrac,
861  Real& xloc_on_yface, Real& zloc_on_yface,
862  bool is_eb_dirichlet, bool is_eb_inhomog,
863  const bool on_x_face, const int domlo_x, const int domhi_x,
864  const bool on_y_face, const int domlo_y, const int domhi_y,
865  const bool on_z_face, const int domlo_z, const int domhi_z)
866 {
867  AMREX_ALWAYS_ASSERT(is_eb_dirichlet);
868  Array2D<Real,0,53,0,9> Amatrix;
869  Array1D<Real,0, 9 > rhs;
870 
871  // Order of column -- first 9 are cell centroids, second 9 are EB centroids
872 
873  for (int irow = 0; irow < 54; irow++) {
874  for (int icol = 0; icol < 10; icol++) {
875  Amatrix(irow,icol) = Real(0.0);
876  }
877  }
878 
879  const int jm = (j > domhi_y) ? 2 : 1;
880  const int jp = 2 - jm;
881 
882  // Columns 0-2: [e x y]
883  for (int jj = j-jm; jj <= j+jp; jj++) { // Normal to face
884  for (int kk = k-1; kk <= k+1; kk++) { // Tangential to face
885  for (int ii = i-1; ii <= i+1; ii++) // Tangential to face
886  {
887  // Don't include corner cells. Could make this even more strict
888  // by removing the on_?_face restrictions.
889  if (((on_x_face && ii < domlo_x) && (on_y_face && jj < domlo_y)) ||
890  ((on_x_face && ii < domlo_x) && (on_y_face && jj > domhi_y)) ||
891  ((on_x_face && ii > domhi_x) && (on_y_face && jj < domlo_y)) ||
892  ((on_x_face && ii > domhi_x) && (on_y_face && jj > domhi_y)) ||
893  ((on_y_face && jj < domlo_y) && (on_z_face && kk < domlo_z)) ||
894  ((on_y_face && jj < domlo_y) && (on_z_face && kk > domhi_z)) ||
895  ((on_y_face && jj > domhi_y) && (on_z_face && kk < domlo_z)) ||
896  ((on_y_face && jj > domhi_y) && (on_z_face && kk > domhi_z)) ||
897  ((on_x_face && ii < domlo_x) && (on_z_face && kk < domlo_z)) ||
898  ((on_x_face && ii < domlo_x) && (on_z_face && kk > domhi_z)) ||
899  ((on_x_face && ii > domhi_x) && (on_z_face && kk < domlo_z)) ||
900  ((on_x_face && ii > domhi_x) && (on_z_face && kk > domhi_z))) {
901  continue;
902  }
903 
904  if ( !phi.contains(ii,jj,kk) ) {
905  continue;
906  }
907 
908  if (!flag(ii,jj,kk).isCovered())
909  {
910  int a_ind = (ii-(i-1)) + 3*(kk-(k-1)) + 9*(jj-(j-jm));
911 
912  Real x_off = static_cast<Real>(ii-i);
913  Real y_off = static_cast<Real>(jj-j) + Real(0.5);
914  Real z_off = static_cast<Real>(kk-k);
915 
916  if(on_x_face){
917  if (ii < domlo_x && (vfrac(ii+1,jj,kk) != Real(1.0) || vfrac(ii+2,jj,kk) != Real(1.0)) ) {
918  continue;
919  }
920  if (ii > domhi_x && (vfrac(ii-1,jj,kk) != Real(1.0) || vfrac(ii-2,jj,kk) != Real(1.0))) {
921  continue;
922  }
923  }
924 
925  if(on_y_face){
926  if (jj < domlo_y && (vfrac(ii,jj+1,kk) != Real(1.0) || vfrac(ii,jj+2,kk) != Real(1.0)) ) {
927  continue;
928  }
929  if (jj > domhi_y && (vfrac(ii,jj-1,kk) != Real(1.0) || vfrac(ii,jj-2,kk) != Real(1.0)) ) {
930  continue;
931  }
932  }
933 
934  if(on_z_face){
935  if (kk < domlo_z && (vfrac(ii,jj,kk+1) != Real(1.0) || vfrac(ii,jj,kk+2) != Real(1.0)) ) {
936  continue;
937  }
938  if (kk > domhi_z && (vfrac(ii,jj,kk-1) != Real(1.0) || vfrac(ii,jj,kk-2) != Real(1.0)) ) {
939  continue;
940  }
941  }
942 
943  Amatrix(a_ind,0) = Real(1.0);
944  Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0) - xloc_on_yface;
945  Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1);
946  Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2) - zloc_on_yface;
947 
948  // Add in information about the location of the EB. Exclude
949  // EBs that are outside the domain.
950  if (flag(ii,jj,kk).isSingleValued() &&
951  domlo_x <= ii && ii <= domhi_x &&
952  domlo_y <= jj && jj <= domhi_y &&
953  domlo_z <= kk && kk <= domhi_z)
954  {
955  Amatrix(a_ind+27,0) = Real(1.0);
956  Amatrix(a_ind+27,1) = x_off + bcent(ii,jj,kk,0) - xloc_on_yface;
957  Amatrix(a_ind+27,2) = y_off + bcent(ii,jj,kk,1);
958  Amatrix(a_ind+27,3) = z_off + bcent(ii,jj,kk,2) - zloc_on_yface;
959  }
960  }
961  }
962  }
963  }
964 
965  // Columns 4-9 : [x*x x*y y*y x*z y*z z*z]
966  for (int irow = 0; irow < 54; irow++)
967  {
968  Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1);
969  Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2);
970  Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2);
971  Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3);
972  Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3);
973  Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3);
974  }
975 
976  // Make the RHS = A^T v
977  for (int irow = 0; irow < 10; irow++)
978  {
979  rhs(irow) = 0.;
980 
981  for (int jj = j-jm; jj <= j+jp; jj++) { // Normal to face
982  for (int kk = k-1; kk <= k+1; kk++) { // Tangential to face
983  for (int ii = i-1; ii <= i+1; ii++) {// Tangential to face
984  if ( !phi.contains(ii,jj,kk) ) {
985  continue;
986  }
987 
988  if (!flag(ii,jj,kk).isCovered())
989  {
990  int a_ind = (ii-(i-1)) + 3*(kk-(k-1)) + 9*(jj-(j-jm));
991  Real phi_val = Amatrix(a_ind,0) * phi(ii,jj,kk,n);
992 
993  rhs(irow) += Amatrix(a_ind,irow)* phi_val;
994 
995  if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog && Amatrix(a_ind+27,irow) != Real(0.0)) {
996  rhs(irow) += Amatrix(a_ind+27,irow)*phieb(ii,jj,kk,n);
997  }
998  }
999  }
1000  }
1001  }
1002  }
1003 
1004  cholsol_for_eb(Amatrix, rhs);
1005 
1006  return rhs(2);
1007 }
1008 
1009 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1010 Real grad_z_of_phi_on_centroids_extdir(int i,int j,int k,int n,
1011  Array4<Real const> const& phi,
1012  Array4<Real const> const& phieb,
1013  Array4<EBCellFlag const> const& flag,
1014  Array4<Real const> const& ccent,
1015  Array4<Real const> const& bcent,
1016  Array4<Real const> const& vfrac,
1017  Real& xloc_on_zface, Real& yloc_on_zface,
1018  bool is_eb_dirichlet, bool is_eb_inhomog,
1019  const bool on_x_face, const int domlo_x, const int domhi_x,
1020  const bool on_y_face, const int domlo_y, const int domhi_y,
1021  const bool on_z_face, const int domlo_z, const int domhi_z)
1022 {
1023  AMREX_ALWAYS_ASSERT(is_eb_dirichlet);
1024  Array2D<Real,0,53,0,9> Amatrix;
1025  Array1D<Real,0, 9 > rhs;
1026 
1027  // Order of column -- first 9 are cell centroids, second 9 are EB centroids
1028 
1029  for (int irow = 0; irow < 54; irow++) {
1030  for (int icol = 0; icol < 10; icol++) {
1031  Amatrix(irow,icol) = Real(0.0);
1032  }
1033  }
1034 
1035  const int km = (k > domhi_z) ? 2 : 1;
1036  const int kp = 2 - km;
1037 
1038  // Columns 0-3: [e x y z]
1039  for (int kk = k-km; kk <= k+kp; kk++) // Normal to face
1040  {
1041  for (int jj = j-1; jj <= j+1; jj++) { // Tangential to face
1042  for (int ii = i-1; ii <= i+1; ii++) // Tangential to face
1043  {
1044  // Don't include corner cells. Could make this even more strict
1045  // by removing the on_?_face restrictions.
1046  if (((on_x_face && ii < domlo_x) && (on_y_face && jj < domlo_y)) ||
1047  ((on_x_face && ii < domlo_x) && (on_y_face && jj > domhi_y)) ||
1048  ((on_x_face && ii > domhi_x) && (on_y_face && jj < domlo_y)) ||
1049  ((on_x_face && ii > domhi_x) && (on_y_face && jj > domhi_y)) ||
1050  ((on_y_face && jj < domlo_y) && (on_z_face && kk < domlo_z)) ||
1051  ((on_y_face && jj < domlo_y) && (on_z_face && kk > domhi_z)) ||
1052  ((on_y_face && jj > domhi_y) && (on_z_face && kk < domlo_z)) ||
1053  ((on_y_face && jj > domhi_y) && (on_z_face && kk > domhi_z)) ||
1054  ((on_x_face && ii < domlo_x) && (on_z_face && kk < domlo_z)) ||
1055  ((on_x_face && ii < domlo_x) && (on_z_face && kk > domhi_z)) ||
1056  ((on_x_face && ii > domhi_x) && (on_z_face && kk < domlo_z)) ||
1057  ((on_x_face && ii > domhi_x) && (on_z_face && kk > domhi_z))) {
1058  continue;
1059  }
1060 
1061  if (!phi.contains(ii,jj,kk)) {
1062  continue;
1063  }
1064 
1065  if (!flag(ii,jj,kk).isCovered())
1066  {
1067  int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-km));
1068 
1069  Real x_off = static_cast<Real>(ii-i);
1070  Real y_off = static_cast<Real>(jj-j);
1071  Real z_off = static_cast<Real>(kk-k) + Real(0.5);
1072 
1073  if(on_x_face){
1074  if (ii < domlo_x && (vfrac(ii+1,jj,kk) != Real(1.0) || vfrac(ii+2,jj,kk) != Real(1.0)) ) {
1075  continue;
1076  }
1077  if (ii > domhi_x && (vfrac(ii-1,jj,kk) != Real(1.0) || vfrac(ii-2,jj,kk) != Real(1.0))) {
1078  continue;
1079  }
1080  }
1081 
1082  if(on_y_face){
1083  if (jj < domlo_y && (vfrac(ii,jj+1,kk) != Real(1.0) || vfrac(ii,jj+2,kk) != Real(1.0)) ) {
1084  continue;
1085  }
1086  if (jj > domhi_y && (vfrac(ii,jj-1,kk) != Real(1.0) || vfrac(ii,jj-2,kk) != Real(1.0)) ) {
1087  continue;
1088  }
1089  }
1090 
1091  if(on_z_face){
1092  if (kk < domlo_z && (vfrac(ii,jj,kk+1) != Real(1.0) || vfrac(ii,jj,kk+2) != Real(1.0)) ) {
1093  continue;
1094  }
1095  if (kk > domhi_z && (vfrac(ii,jj,kk-1) != Real(1.0) || vfrac(ii,jj,kk-2) != Real(1.0)) ) {
1096  continue;
1097  }
1098  }
1099 
1100  Amatrix(a_ind,0) = Real(1.0);
1101  Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0) - xloc_on_zface;
1102  Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1) - yloc_on_zface ;
1103  Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2);
1104 
1105  // Add in information about the location of the EB. Exclude
1106  // EBs that are outside the domain.
1107  if (flag(ii,jj,kk).isSingleValued() &&
1108  domlo_x <= ii && ii <= domhi_x &&
1109  domlo_y <= jj && jj <= domhi_y &&
1110  domlo_z <= kk && kk <= domhi_z)
1111  {
1112  Amatrix(a_ind+27,0) = Real(1.0);
1113  Amatrix(a_ind+27,1) = x_off + bcent(ii,jj,kk,0) - xloc_on_zface;
1114  Amatrix(a_ind+27,2) = y_off + bcent(ii,jj,kk,1) - yloc_on_zface;
1115  Amatrix(a_ind+27,3) = z_off + bcent(ii,jj,kk,2);
1116  }
1117  }
1118  }
1119  }
1120  }
1121 
1122  // Columns 4-9 : [x*x x*y y*y x*z y*z z*z]
1123  for (int irow = 0; irow < 54; irow++)
1124  {
1125  Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1);
1126  Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2);
1127  Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2);
1128  Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3);
1129  Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3);
1130  Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3);
1131  }
1132 
1133  // Make the RHS = A^T v
1134  for (int irow = 0; irow < 10; irow++)
1135  {
1136  rhs(irow) = 0.;
1137 
1138  for (int kk = k-km; kk <= k+kp; kk++) { // Normal to face
1139  for (int jj = j-1; jj <= j+1; jj++) { // Tangential to face
1140  for (int ii = i-1; ii <= i+1; ii++) {// Tangential to face
1141  if ( !phi.contains(ii,jj,kk) ) {
1142  continue;
1143  }
1144 
1145  if (!flag(ii,jj,kk).isCovered())
1146  {
1147  int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-km));
1148  Real phi_val = Amatrix(a_ind,0) * phi(ii,jj,kk,n);
1149 
1150  rhs(irow) += Amatrix(a_ind,irow)* phi_val;
1151 
1152  if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog && Amatrix(a_ind+27,irow) != Real(0.0)) {
1153  rhs(irow) += Amatrix(a_ind+27,irow)*phieb(ii,jj,kk,n);
1154  }
1155  }
1156  }
1157  }
1158  }
1159  }
1160 
1161  cholsol_for_eb(Amatrix, rhs);
1162 
1163  return rhs(3);
1164 }
1165 
1166 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1167 Real grad_eb_of_phi_on_centroids_extdir(int i,int j,int k,int n,
1168  Array4<Real const> const& phi,
1169  Array4<Real const> const& phieb,
1170  Array4<EBCellFlag const> const& flag,
1171  Array4<Real const> const& ccent,
1172  Array4<Real const> const& bcent,
1173  Array4<Real const> const& vfrac,
1174  Real& nrmx, Real& nrmy, Real& nrmz,
1175  bool is_eb_inhomog,
1176  const bool on_x_face, const int domlo_x, const int domhi_x,
1177  const bool on_y_face, const int domlo_y, const int domhi_y,
1178  const bool on_z_face, const int domlo_z, const int domhi_z)
1179 {
1180  Array2D<Real,0,53,0,9> Amatrix;
1181  Array1D<Real,0, 9 > rhs;
1182 
1183  // Order of column -- first 27 are cell centroids, second 27 are EB centroids
1184 
1185  for (int irow = 0; irow < 54; irow++) {
1186  for (int icol = 0; icol < 10; icol++) {
1187  Amatrix(irow,icol) = Real(0.0);
1188  }
1189  }
1190 
1191  // Columns 0-3: [e x y z]
1192  for (int kk = k-1; kk <= k+1; kk++) {
1193  for (int jj = j-1; jj <= j+1; jj++) {
1194  for (int ii = i-1; ii <= i+1; ii++)
1195  {
1196  // This is likely overkill for EB grads.
1197  // Don't include corner cells. Could make this even more strict
1198  // by removing the on_?_face restrictions.
1199  if (((on_x_face && ii < domlo_x) && (on_y_face && jj < domlo_y)) ||
1200  ((on_x_face && ii < domlo_x) && (on_y_face && jj > domhi_y)) ||
1201  ((on_x_face && ii > domhi_x) && (on_y_face && jj < domlo_y)) ||
1202  ((on_x_face && ii > domhi_x) && (on_y_face && jj > domhi_y)) ||
1203  ((on_y_face && jj < domlo_y) && (on_z_face && kk < domlo_z)) ||
1204  ((on_y_face && jj < domlo_y) && (on_z_face && kk > domhi_z)) ||
1205  ((on_y_face && jj > domhi_y) && (on_z_face && kk < domlo_z)) ||
1206  ((on_y_face && jj > domhi_y) && (on_z_face && kk > domhi_z)) ||
1207  ((on_x_face && ii < domlo_x) && (on_z_face && kk < domlo_z)) ||
1208  ((on_x_face && ii < domlo_x) && (on_z_face && kk > domhi_z)) ||
1209  ((on_x_face && ii > domhi_x) && (on_z_face && kk < domlo_z)) ||
1210  ((on_x_face && ii > domhi_x) && (on_z_face && kk > domhi_z))) {
1211  continue;
1212  }
1213 
1214  if ( !phi.contains(ii,jj,kk) ) {
1215  continue;
1216  }
1217 
1218  if (!flag(ii,jj,kk).isCovered())
1219  {
1220  int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-1));
1221 
1222  Real x_off = static_cast<Real>(ii-i) - bcent(i,j,k,0);
1223  Real y_off = static_cast<Real>(jj-j) - bcent(i,j,k,1);
1224  Real z_off = static_cast<Real>(kk-k) - bcent(i,j,k,2);
1225 
1226  if(on_x_face){
1227  if (ii < domlo_x && (vfrac(ii+1,jj,kk) != Real(1.0) || vfrac(ii+2,jj,kk) != Real(1.0)) ) {
1228  continue;
1229  }
1230  if (ii > domhi_x && (vfrac(ii-1,jj,kk) != Real(1.0) || vfrac(ii-2,jj,kk) != Real(1.0))) {
1231  continue;
1232  }
1233  }
1234 
1235  if(on_y_face){
1236  if (jj < domlo_y && (vfrac(ii,jj+1,kk) != Real(1.0) || vfrac(ii,jj+2,kk) != Real(1.0)) ) {
1237  continue;
1238  }
1239  if (jj > domhi_y && (vfrac(ii,jj-1,kk) != Real(1.0) || vfrac(ii,jj-2,kk) != Real(1.0)) ) {
1240  continue;
1241  }
1242  }
1243 
1244  if(on_z_face){
1245  if (kk < domlo_z && (vfrac(ii,jj,kk+1) != Real(1.0) || vfrac(ii,jj,kk+2) != Real(1.0)) ) {
1246  continue;
1247  }
1248  if (kk > domhi_z && (vfrac(ii,jj,kk-1) != Real(1.0) || vfrac(ii,jj,kk-2) != Real(1.0)) ) {
1249  continue;
1250  }
1251  }
1252 
1253  Amatrix(a_ind,0) = Real(1.0);
1254  Amatrix(a_ind,1) = x_off + ccent(ii,jj,kk,0);
1255  Amatrix(a_ind,2) = y_off + ccent(ii,jj,kk,1);
1256  Amatrix(a_ind,3) = z_off + ccent(ii,jj,kk,2);
1257 
1258  if (flag(ii,jj,kk).isSingleValued() &&
1259  domlo_x <= ii && ii <= domhi_x &&
1260  domlo_y <= jj && jj <= domhi_y &&
1261  domlo_z <= kk && kk <= domhi_z)
1262  {
1263  Amatrix(a_ind+27,0) = Real(1.0);
1264  Amatrix(a_ind+27,1) = x_off + bcent(ii,jj,kk,0);
1265  Amatrix(a_ind+27,2) = y_off + bcent(ii,jj,kk,1);
1266  Amatrix(a_ind+27,3) = z_off + bcent(ii,jj,kk,2);
1267  }
1268  }
1269  }
1270  }
1271  }
1272 
1273  // Columns 4-9 : [x*x x*y y*y x*z y*z z*z]
1274  for (int irow = 0; irow < 54; irow++)
1275  {
1276  Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,1);
1277  Amatrix(irow,5) = Amatrix(irow,1) * Amatrix(irow,2);
1278  Amatrix(irow,6) = Amatrix(irow,2) * Amatrix(irow,2);
1279  Amatrix(irow,7) = Amatrix(irow,1) * Amatrix(irow,3);
1280  Amatrix(irow,8) = Amatrix(irow,2) * Amatrix(irow,3);
1281  Amatrix(irow,9) = Amatrix(irow,3) * Amatrix(irow,3);
1282  }
1283 
1284  // Make the RHS = A^T v
1285  for (int irow = 0; irow < 10; irow++)
1286  {
1287  rhs(irow) = 0.;
1288 
1289  for (int kk = k-1; kk <= k+1; kk++) {
1290  for (int jj = j-1; jj <= j+1; jj++) {
1291  for (int ii = i-1; ii <= i+1; ii++) {
1292  if ( !phi.contains(ii,jj,kk) ) {
1293  continue;
1294  }
1295 
1296  if (!flag(ii,jj,kk).isCovered())
1297  {
1298  int a_ind = (ii-(i-1)) + 3*(jj-(j-1)) + 9*(kk-(k-1));
1299 
1300  Real phi_val = Amatrix(a_ind,0) * phi(ii,jj,kk,n);
1301 
1302  rhs(irow) += Amatrix(a_ind,irow)* phi_val;
1303 
1304  if (flag(ii,jj,kk).isSingleValued() && is_eb_inhomog && Amatrix(a_ind+27,irow) != Real(0.0)) {
1305  rhs(irow) += Amatrix(a_ind+27,irow)*phieb(ii,jj,kk,n);
1306  }
1307  }
1308  }
1309  }
1310  }
1311  }
1312 
1313  cholsol_for_eb(Amatrix, rhs);
1314 
1315  Real dphidn = rhs(1)*nrmx + rhs(2)*nrmy + rhs(3)*nrmz;
1316  return dphidn;
1317 }
1318 
1319 }
1320 
1321 #endif
AMREX_ALWAYS_ASSERT
#define AMREX_ALWAYS_ASSERT(EX)
Definition: AMReX_BLassert.H:50
AMReX_EBCellFlag.H
AMREX_FORCE_INLINE
#define AMREX_FORCE_INLINE
Definition: AMReX_Extension.H:119
AMREX_GPU_HOST_DEVICE
#define AMREX_GPU_HOST_DEVICE
Definition: AMReX_GpuQualifiers.H:20
amrex
Definition: AMReX_Amr.cpp:49
amrex::grad_y_of_phi_on_centroids
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_y_of_phi_on_centroids(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Real &xloc_on_yface, bool is_eb_dirichlet, bool is_eb_inhomog)
Definition: AMReX_EB_LeastSquares_2D_K.H:331
amrex::cholsol_for_eb
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void cholsol_for_eb(Array2D< Real, 0, 17, 0, 5 > &Amatrix, Array1D< Real, 0, 5 > &b)
Definition: AMReX_EB_LeastSquares_2D_K.H:155
amrex::grad_eb_of_phi_on_centroids
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_eb_of_phi_on_centroids(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Real &nrmx, Real &nrmy, bool is_eb_inhomog)
Definition: AMReX_EB_LeastSquares_2D_K.H:412
amrex::grad_x_of_phi_on_centroids_extdir
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_x_of_phi_on_centroids_extdir(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Array4< Real const > const &vfrac, Real &yloc_on_xface, bool is_eb_dirichlet, bool is_eb_inhomog, const bool on_x_face, const int domlo_x, const int domhi_x, const bool on_y_face, const int domlo_y, const int domhi_y)
Definition: AMReX_EB_LeastSquares_2D_K.H:495
amrex::grad_z_of_phi_on_centroids
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_z_of_phi_on_centroids(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Real &xloc_on_zface, Real &yloc_on_zface, bool is_eb_dirichlet, bool is_eb_inhomog)
Definition: AMReX_EB_LeastSquares_3D_K.H:513
amrex::grad_eb_of_phi_on_centroids_extdir
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_eb_of_phi_on_centroids_extdir(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Array4< Real const > const &vfrac, Real &nrmx, Real &nrmy, bool is_eb_inhomog, const bool on_x_face, const int domlo_x, const int domhi_x, const bool on_y_face, const int domlo_y, const int domhi_y)
Definition: AMReX_EB_LeastSquares_2D_K.H:765
amrex::grad_z_of_phi_on_centroids_extdir
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_z_of_phi_on_centroids_extdir(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Array4< Real const > const &vfrac, Real &xloc_on_zface, Real &yloc_on_zface, bool is_eb_dirichlet, bool is_eb_inhomog, const bool on_x_face, const int domlo_x, const int domhi_x, const bool on_y_face, const int domlo_y, const int domhi_y, const bool on_z_face, const int domlo_z, const int domhi_z)
Definition: AMReX_EB_LeastSquares_3D_K.H:1010
amrex::grad_y_of_phi_on_centroids_extdir
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_y_of_phi_on_centroids_extdir(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Array4< Real const > const &vfrac, Real &xloc_on_yface, bool is_eb_dirichlet, bool is_eb_inhomog, const bool on_x_face, const int domlo_x, const int domhi_x, const bool on_y_face, const int domlo_y, const int domhi_y)
Definition: AMReX_EB_LeastSquares_2D_K.H:632
amrex::cholsol_np10
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void cholsol_np10(Array2D< Real, 0, 35, 0, 9 > &Amatrix, Array1D< Real, 0, 9 > &b)
Definition: AMReX_EB_LeastSquares_3D_K.H:64
amrex::grad_x_of_phi_on_centroids
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real grad_x_of_phi_on_centroids(int i, int j, int k, int n, Array4< Real const > const &phi, Array4< Real const > const &phieb, Array4< EBCellFlag const > const &flag, Array4< Real const > const &ccent, Array4< Real const > const &bcent, Real &yloc_on_xface, bool is_eb_dirichlet, bool is_eb_inhomog)
Definition: AMReX_EB_LeastSquares_2D_K.H:246
amrex::decomp_chol_np10
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void decomp_chol_np10(Array2D< Real, 0, 9, 0, 9 > &aa)
Definition: AMReX_EB_LeastSquares_3D_K.H:10
amrex::sqrt
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > sqrt(const GpuComplex< T > &a_z) noexcept
Return the square root of a complex number.
Definition: AMReX_GpuComplex.H:373
amrex::Array1D
Definition: AMReX_Array.H:164
amrex::Array2D
Definition: AMReX_Array.H:285
amrex::Array4
Definition: AMReX_Array4.H:61
amrex::Array4::contains
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool contains(int i, int j, int k) const noexcept
Definition: AMReX_Array4.H:251