Block-Structured AMR Software Framework
AMReX_EB_LeastSquares_2D_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 2D version (i.e. for 6x6 (A^T A) matrix)
9 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
10 void decomp_chol_np6(Array2D<Real,0,5,0,5>& aa)
11 {
12  int neq = 6;
13 
14  Real p[6];
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 <= 0.)
34  {
35  p[ii] = 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) = 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) = 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 2D version (i.e. for 6x6 (A^T A) matrix)
63 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
64 void cholsol_np6(Array2D<Real,0,11,0,5>& Amatrix, Array1D<Real,0,5>& b)
65 {
66  int neq = 6;
67 
68  Array2D<Real,0,5,0,5> AtA;
69 
70  for (int irow = 0; irow < neq; irow++) {
71  for (int icol = 0; icol < neq; icol++) {
72  AtA(irow,icol) = 0.0;
73  }
74  }
75 
76  for (int irow = 0; irow < 12; 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 x^2
82  AtA(0,4) += Amatrix(irow,0)*Amatrix(irow,4); // e^T x*y
83  AtA(0,5) += Amatrix(irow,0)*Amatrix(irow,5); // e^T y^2
84  AtA(1,1) += Amatrix(irow,1)*Amatrix(irow,1); // x^T x
85  AtA(1,2) += Amatrix(irow,1)*Amatrix(irow,2); // x^T y
86  AtA(1,3) += Amatrix(irow,1)*Amatrix(irow,3); // x^T (x^2)
87  AtA(1,4) += Amatrix(irow,1)*Amatrix(irow,4); // x^T (xy)
88  AtA(1,5) += Amatrix(irow,1)*Amatrix(irow,5); // x^T (y^2)
89  AtA(2,2) += Amatrix(irow,2)*Amatrix(irow,2); // y^T y
90  AtA(2,3) += Amatrix(irow,2)*Amatrix(irow,3); // y^T (x^2)
91  AtA(2,4) += Amatrix(irow,2)*Amatrix(irow,4); // y^T (xy)
92  AtA(2,5) += Amatrix(irow,2)*Amatrix(irow,5); // y^T (y^2)
93  AtA(3,3) += Amatrix(irow,3)*Amatrix(irow,3); // (x^2)^T (x^2)
94  AtA(3,4) += Amatrix(irow,3)*Amatrix(irow,4); // (x^2)^T (x*y)
95  AtA(3,5) += Amatrix(irow,3)*Amatrix(irow,5); // (x^2)^T (y^2)
96  AtA(4,4) += Amatrix(irow,4)*Amatrix(irow,4); // (x*y)^T (x*y)
97  AtA(4,5) += Amatrix(irow,4)*Amatrix(irow,5); // (x*y)^T (y^2)
98  AtA(5,5) += Amatrix(irow,5)*Amatrix(irow,5); // (y^2)^T (y^2)
99  }
100 
101  for (int irow = 0; irow < neq-1; irow++) {
102  for (int icol = irow+1; icol < neq; icol++) {
103  AtA(icol,irow) = AtA(irow,icol);
104  }
105  }
106 
107  decomp_chol_np6(AtA);
108 
109  if (AtA(0,0) > 0.) {
110  b(0) = b(0) / AtA(0,0);
111  } else {
112  b(0) = 0.;
113  }
114 
115  for (int ii = 1; ii < neq; ii++)
116  {
117  if (AtA(ii,ii) > 0.)
118  {
119  for (int jj = 0; jj < ii; jj++) {
120  b(ii) = b(ii) - AtA(ii,jj)*b(jj);
121  }
122 
123  b(ii) = b(ii) / AtA(ii,ii);
124  }
125  else
126  {
127  b(ii) = 0.0;
128  }
129  }
130 
131  if (AtA(neq-1,neq-1) > 0.) {
132  b(neq-1) = b(neq-1) / AtA(neq-1,neq-1);
133  } else {
134  b(neq-1) = 0.0;
135  }
136 
137  for (int ii = neq-2; ii >= 0; ii--)
138  {
139  if (AtA(ii,ii) > 0.)
140  {
141  for (int jj = ii+1; jj < neq; jj++) {
142  b(ii) = b(ii) - AtA(ii,jj)*b(jj);
143  }
144 
145  b(ii) = b(ii) / AtA(ii,ii);
146  }
147  else
148  {
149  b(ii) = 0.0;
150  }
151  }
152 }
153 
154 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
155 void cholsol_for_eb(Array2D<Real,0,17,0,5>& Amatrix, Array1D<Real,0,5>& b)
156 {
157  int neq = 6;
158 
159  Array2D<Real,0,5,0,5> AtA;
160 
161  for (int irow = 0; irow < neq; irow++) {
162  for (int icol = 0; icol < neq; icol++) {
163  AtA(irow,icol) = 0.0;
164  }
165  }
166 
167  for (int irow = 0; irow < 18; irow++)
168  {
169  AtA(0,0) += Amatrix(irow,0)*Amatrix(irow,0); // e^T e
170  AtA(0,1) += Amatrix(irow,0)*Amatrix(irow,1); // e^T x
171  AtA(0,2) += Amatrix(irow,0)*Amatrix(irow,2); // e^T y
172  AtA(0,3) += Amatrix(irow,0)*Amatrix(irow,3); // e^T x^2
173  AtA(0,4) += Amatrix(irow,0)*Amatrix(irow,4); // e^T x*y
174  AtA(0,5) += Amatrix(irow,0)*Amatrix(irow,5); // e^T y^2
175  AtA(1,1) += Amatrix(irow,1)*Amatrix(irow,1); // x^T x
176  AtA(1,2) += Amatrix(irow,1)*Amatrix(irow,2); // x^T y
177  AtA(1,3) += Amatrix(irow,1)*Amatrix(irow,3); // x^T (x^2)
178  AtA(1,4) += Amatrix(irow,1)*Amatrix(irow,4); // x^T (xy)
179  AtA(1,5) += Amatrix(irow,1)*Amatrix(irow,5); // x^T (y^2)
180  AtA(2,2) += Amatrix(irow,2)*Amatrix(irow,2); // y^T y
181  AtA(2,3) += Amatrix(irow,2)*Amatrix(irow,3); // y^T (x^2)
182  AtA(2,4) += Amatrix(irow,2)*Amatrix(irow,4); // y^T (xy)
183  AtA(2,5) += Amatrix(irow,2)*Amatrix(irow,5); // y^T (y^2)
184  AtA(3,3) += Amatrix(irow,3)*Amatrix(irow,3); // (x^2)^T (x^2)
185  AtA(3,4) += Amatrix(irow,3)*Amatrix(irow,4); // (x^2)^T (x*y)
186  AtA(3,5) += Amatrix(irow,3)*Amatrix(irow,5); // (x^2)^T (y^2)
187  AtA(4,4) += Amatrix(irow,4)*Amatrix(irow,4); // (x*y)^T (x*y)
188  AtA(4,5) += Amatrix(irow,4)*Amatrix(irow,5); // (x*y)^T (y^2)
189  AtA(5,5) += Amatrix(irow,5)*Amatrix(irow,5); // (y^2)^T (y^2)
190  }
191 
192  for (int irow = 0; irow < neq-1; irow++) {
193  for (int icol = irow+1; icol < neq; icol++) {
194  AtA(icol,irow) = AtA(irow,icol);
195  }
196  }
197 
198  decomp_chol_np6(AtA);
199 
200  if (AtA(0,0) > 0.) {
201  b(0) = b(0) / AtA(0,0);
202  } else {
203  b(0) = 0.;
204  }
205 
206  for (int ii = 1; ii < neq; ii++)
207  {
208  if (AtA(ii,ii) > 0.)
209  {
210  for (int jj = 0; jj < ii; jj++) {
211  b(ii) = b(ii) - AtA(ii,jj)*b(jj);
212  }
213 
214  b(ii) = b(ii) / AtA(ii,ii);
215  }
216  else
217  {
218  b(ii) = 0.0;
219  }
220  }
221 
222  if (AtA(neq-1,neq-1) > 0.) {
223  b(neq-1) = b(neq-1) / AtA(neq-1,neq-1);
224  } else {
225  b(neq-1) = 0.0;
226  }
227 
228  for (int ii = neq-2; ii >= 0; ii--)
229  {
230  if (AtA(ii,ii) > 0.)
231  {
232  for (int jj = ii+1; jj < neq; jj++) {
233  b(ii) = b(ii) - AtA(ii,jj)*b(jj);
234  }
235 
236  b(ii) = b(ii) / AtA(ii,ii);
237  }
238  else
239  {
240  b(ii) = 0.0;
241  }
242  }
243 }
244 
245 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
246 Real grad_x_of_phi_on_centroids(int i,int j,int k,int n,
247  Array4<Real const> const& phi,
248  Array4<Real const> const& phieb,
249  Array4<EBCellFlag const> const& flag,
250  Array4<Real const> const& ccent,
251  Array4<Real const> const& bcent,
252  Real& yloc_on_xface,
253  bool is_eb_dirichlet, bool is_eb_inhomog)
254 {
255  AMREX_ALWAYS_ASSERT(is_eb_dirichlet);
256  Array2D<Real,0,11,0,5> Amatrix;
257  Array1D<Real,0, 5 > rhs;
258 
259  // Order of column -- first six are cell centroids: (i,j-1) (i,j) (i,j+1) (i-1,j-1) (i-1,j) (i-1,j+1)
260  // Order of column -- second six are EB centroids: (i,j-1) (i,j) (i,j+1) (i-1,j-1) (i-1,j) (i-1,j+1)
261 
262  for (int irow = 0; irow < 12; irow++) {
263  for (int icol = 0; icol < 6; icol++) {
264  Amatrix(irow,icol) = 0.0;
265  }
266  }
267 
268  // Columns: [e x y x*x x*y y*y]
269  for (int ii = i-1; ii <= i; ii++) { // Normal to face
270  for (int jj = j-1; jj <= j+1; jj++) // Tangential to face
271  {
272  if (!flag(ii,jj,k).isCovered())
273  {
274  int a_ind = (jj-(j-1)) + 3*(ii-(i-1));
275 
276  // This shifts the x-offset by a normalized half-dx so that
277  // x=0 is the x-coordinate of the FACE centroid.
278  Real x_off = static_cast<Real>(ii-i) + 0.5;
279  Real y_off = static_cast<Real>(jj-j);
280 
281  Amatrix(a_ind,0) = 1.0;
282  Amatrix(a_ind,1) = x_off + ccent(ii,jj,k,0);
283  Amatrix(a_ind,2) = y_off + ccent(ii,jj,k,1) - yloc_on_xface;
284  Amatrix(a_ind,3) = Amatrix(a_ind,1) * Amatrix(a_ind,1);
285  Amatrix(a_ind,4) = Amatrix(a_ind,1) * Amatrix(a_ind,2);
286  Amatrix(a_ind,5) = Amatrix(a_ind,2) * Amatrix(a_ind,2);
287 
288  if (!flag(ii,jj,k).isRegular())
289  {
290  int a_ind_eb = a_ind + 6;
291 
292  Amatrix(a_ind_eb,0) = 1.0;
293  Amatrix(a_ind_eb,1) = x_off + bcent(ii,jj,k,0);
294  Amatrix(a_ind_eb,2) = y_off + bcent(ii,jj,k,1) - yloc_on_xface;
295  Amatrix(a_ind_eb,3) = Amatrix(a_ind_eb,1) * Amatrix(a_ind_eb,1);
296  Amatrix(a_ind_eb,4) = Amatrix(a_ind_eb,1) * Amatrix(a_ind_eb,2);
297  Amatrix(a_ind_eb,5) = Amatrix(a_ind_eb,2) * Amatrix(a_ind_eb,2);
298  }
299  }
300  }
301  }
302 
303  // Make the RHS = A^T v
304  for (int irow = 0; irow < 6; irow++)
305  {
306  rhs(irow) = 0.; // Only non-zero when inhomogeneous Dirichlet
307 
308  for (int ii = i-1; ii <= i; ii++) { // Normal to face
309  for (int jj = j-1; jj <= j+1; jj++) // Tangential to face
310  {
311  if (!flag(ii,jj,k).isCovered())
312  {
313  int a_ind = (jj-(j-1)) + 3*(ii-(i-1));
314  rhs(irow) += Amatrix(a_ind ,irow)* phi(ii,jj,k,n);
315 
316  if (flag(ii,jj,k).isSingleValued() &&
317  is_eb_dirichlet && is_eb_inhomog) {
318  rhs(irow) += Amatrix(a_ind+6,irow)*phieb(ii,jj,k,n);
319  }
320  }
321  }
322  }
323  }
324 
325  cholsol_np6(Amatrix, rhs);
326 
327  return rhs(1);
328 }
329 
330 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
331 Real grad_y_of_phi_on_centroids(int i,int j,int k,int n,
332  Array4<Real const> const& phi,
333  Array4<Real const> const& phieb,
334  Array4<EBCellFlag const> const& flag,
335  Array4<Real const> const& ccent,
336  Array4<Real const> const& bcent,
337  Real& xloc_on_yface,
338  bool is_eb_dirichlet, bool is_eb_inhomog)
339 {
340  AMREX_ALWAYS_ASSERT(is_eb_dirichlet);
341  Array2D<Real,0,11,0,5> Amatrix;
342  Array1D<Real,0, 5 > rhs;
343 
344  // Order of column -- first six are cell centroids: (i-1,j) (i,j) (i+1,j) (i-1,j-1) (i,j-1) (i+1,j-1)
345  // Order of column -- second six are EB centroids: (i-1,j) (i,j) (i+1,j) (i-1,j-1) (i,j-1) (i+1,j-1)
346 
347  for (int irow = 0; irow < 12; irow++) {
348  for (int icol = 0; icol < 6; icol++) {
349  Amatrix(irow,icol) = 0.0;
350  }
351  }
352 
353  // Columns: [e x y x*x x*y y*y]
354  for (int jj = j-1; jj <= j; jj++) { // Normal to face
355  for (int ii = i-1; ii <= i+1; ii++) // Tangential to face
356  {
357  if (!flag(ii,jj,k).isCovered())
358  {
359  int a_ind = (ii-(i-1)) + 3*(jj-(j-1));
360 
361  Real x_off = static_cast<Real>(ii-i);
362  Real y_off = static_cast<Real>(jj-j) + 0.5;
363 
364  Amatrix(a_ind,0) = 1.0;
365  Amatrix(a_ind,1) = x_off + ccent(ii,jj,k,0) - xloc_on_yface;
366  Amatrix(a_ind,2) = y_off + ccent(ii,jj,k,1);
367  Amatrix(a_ind,3) = Amatrix(a_ind,1) * Amatrix(a_ind,1);
368  Amatrix(a_ind,4) = Amatrix(a_ind,1) * Amatrix(a_ind,2);
369  Amatrix(a_ind,5) = Amatrix(a_ind,2) * Amatrix(a_ind,2);
370 
371  if (!flag(ii,jj,k).isRegular())
372  {
373  int a_ind_eb = a_ind + 6;
374  Amatrix(a_ind_eb,0) = 1.0;
375  Amatrix(a_ind_eb,1) = x_off + bcent(ii,jj,k,0) - xloc_on_yface;
376  Amatrix(a_ind_eb,2) = y_off + bcent(ii,jj,k,1);
377  Amatrix(a_ind_eb,3) = Amatrix(a_ind_eb,1) * Amatrix(a_ind_eb,1);
378  Amatrix(a_ind_eb,4) = Amatrix(a_ind_eb,1) * Amatrix(a_ind_eb,2);
379  Amatrix(a_ind_eb,5) = Amatrix(a_ind_eb,2) * Amatrix(a_ind_eb,2);
380  }
381  }
382  }
383  }
384 
385  // Make the RHS = A^T v
386  for (int irow = 0; irow < 6; irow++)
387  {
388  rhs(irow) = 0.; // Only non-zero when inhomogeneous Dirichlet
389 
390  for (int jj = j-1; jj <= j; jj++) { // Normal to face
391  for (int ii = i-1; ii <= i+1; ii++) { // Tangential to face
392  if (!flag(ii,jj,k).isCovered())
393  {
394  int a_ind = (ii-(i-1)) + 3*(jj-(j-1));
395  rhs(irow) += Amatrix(a_ind ,irow)* phi(ii,jj,k,n);
396 
397  if (flag(ii,jj,k).isSingleValued() &&
398  is_eb_dirichlet && is_eb_inhomog) {
399  rhs(irow) += Amatrix(a_ind+6,irow)*phieb(ii,jj,k,n);
400  }
401  }
402  }
403  }
404  }
405 
406  cholsol_np6(Amatrix, rhs);
407 
408  return rhs(2);
409 }
410 
411 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
412 Real grad_eb_of_phi_on_centroids(int i,int j,int k,int n,
413  Array4<Real const> const& phi,
414  Array4<Real const> const& phieb,
415  Array4<EBCellFlag const> const& flag,
416  Array4<Real const> const& ccent,
417  Array4<Real const> const& bcent,
418  Real& nrmx, Real& nrmy,
419  bool is_eb_inhomog)
420 {
421  Array2D<Real,0,17,0,5> Amatrix;
422  Array1D<Real,0, 5 > rhs;
423 
424  // Order of column -- first 9 are cell centroids, next 9 are EB centroids
425 
426  for (int irow = 0; irow < 18; irow++) {
427  for (int icol = 0; icol < 6; icol++) {
428  Amatrix(irow,icol) = 0.0;
429  }
430  }
431 
432  // Column 0-2: [e x y]
433  for (int ii = i-1; ii <= i+1; ii++) {
434  for (int jj = j-1; jj <= j+1; jj++)
435  {
436  if (!flag(ii,jj,k).isCovered())
437  {
438  int a_ind = (jj-(j-1)) + 3*(ii-(i-1));
439 
440  Real x_off = static_cast<Real>(ii-i);
441  Real y_off = static_cast<Real>(jj-j);
442 
443  if (flag(i,j,k).isConnected((ii-i),(jj-j),0))
444  {
445  Amatrix(a_ind,0) = 1.0;
446  Amatrix(a_ind,1) = x_off + ccent(ii,jj,k,0) - bcent(i,j,k,0);
447  Amatrix(a_ind,2) = y_off + ccent(ii,jj,k,1) - bcent(i,j,k,1);
448  }
449 
450  if (flag(i,j,k).isConnected((ii-i),(jj-j),0) && !flag(ii,jj,k).isRegular())
451  {
452  Amatrix(a_ind+9,0) = 1.0;
453  Amatrix(a_ind+9,1) = x_off + bcent(ii,jj,k,0) - bcent(i,j,k,0);
454  Amatrix(a_ind+9,2) = y_off + bcent(ii,jj,k,1) - bcent(i,j,k,1);
455  }
456  }
457  }
458  }
459 
460  // Columns 3 : [x*x x*y y*y]
461 
462  for (int irow = 0; irow < 18; irow++)
463  {
464  Amatrix(irow,3) = Amatrix(irow,1) * Amatrix(irow,1);
465  Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,2);
466  Amatrix(irow,5) = Amatrix(irow,2) * Amatrix(irow,2);
467  }
468 
469  // Make the RHS = A^T v
470  for (int irow = 0; irow < 6; irow++)
471  {
472  rhs(irow) = 0.;
473 
474  for (int ii = i-1; ii <= i+1; ii++) {
475  for (int jj = j-1; jj <= j+1; jj++) {
476  if (!flag(ii,jj,k).isCovered())
477  {
478  int a_ind = (jj-(j-1)) + 3*(ii-(i-1));
479  rhs(irow) += Amatrix(a_ind,irow) * phi(ii,jj,k,n);
480  if (flag(ii,jj,k).isSingleValued() && is_eb_inhomog) {
481  rhs(irow) += Amatrix(a_ind+9,irow)*phieb(ii,jj,k,n);
482  }
483  }
484  }
485  }
486  }
487 
488  cholsol_for_eb(Amatrix, rhs);
489 
490  Real dphidn = rhs(1)*nrmx + rhs(2)*nrmy;
491  return dphidn;
492 }
493 
494 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
495 Real grad_x_of_phi_on_centroids_extdir(int i,int j,int k,int n,
496  Array4<Real const> const& phi,
497  Array4<Real const> const& phieb,
498  Array4<EBCellFlag const> const& flag,
499  Array4<Real const> const& ccent,
500  Array4<Real const> const& bcent,
501  Array4<Real const> const& vfrac,
502  Real& yloc_on_xface,
503  bool is_eb_dirichlet, bool is_eb_inhomog,
504  const bool on_x_face, const int domlo_x, const int domhi_x,
505  const bool on_y_face, const int domlo_y, const int domhi_y)
506 {
507  AMREX_ALWAYS_ASSERT(is_eb_dirichlet);
508  Array2D<Real,0,17,0,5> Amatrix;
509  Array1D<Real,0, 5 > rhs;
510 
511  // Order of column -- first six are cell centroids: (i,j-1) (i,j) (i,j+1) (i-1,j-1) (i-1,j) (i-1,j+1)
512  // Order of column -- second six are EB centroids: (i,j-1) (i,j) (i,j+1) (i-1,j-1) (i-1,j) (i-1,j+1)
513 
514  for (int irow = 0; irow < 18; irow++) {
515  for (int icol = 0; icol < 6; icol++) {
516  Amatrix(irow,icol) = 0.0;
517  }
518  }
519 
520  const int im = (i > domhi_x) ? 2 : 1;
521  const int ip = 2 - im;
522 
523  // Columns: [e x y x*x x*y y*y]
524  for (int ii = i-im; ii <= i+ip; ii++) { // Normal to face
525  for (int jj = j-1; jj <= j+1; jj++) // Tangential to face
526  {
527 
528  // Don't include corner cells. Could make this even more strict
529  // by removing the on_?_face restrictions.
530  if (((on_x_face && ii < domlo_x) && (on_y_face && jj < domlo_y)) ||
531  ((on_x_face && ii > domhi_x) && (on_y_face && jj > domhi_y)) ||
532  ((on_x_face && ii < domlo_x) && (on_y_face && jj > domhi_y)) ||
533  ((on_x_face && ii > domhi_x) && (on_y_face && jj < domlo_y))) {
534  continue;
535  }
536 
537  if ( !phi.contains(ii,jj,k) ) {
538  continue;
539  }
540 
541  if (!flag(ii,jj,k).isCovered())
542  {
543 
544  int a_ind = (jj-(j-1)) + 3*(ii-(i-im));
545  AMREX_ASSERT(0<= a_ind && a_ind <= 8);
546 
547  // This shifts the x-offset by a normalized half-dx so that
548  // x=0 is the x-coordinate of the FACE centroid.
549  Real x_off = static_cast<Real>(ii-i) + 0.5;
550  Real y_off = static_cast<Real>(jj-j);
551 
552  if(on_x_face){
553  if (ii < domlo_x && (vfrac(ii+1,jj,k) != 1.0 || vfrac(ii+2,jj,k) != 1.0) ) {
554  continue;
555  }
556  if (ii > domhi_x && (vfrac(ii-1,jj,k) != 1.0 || vfrac(ii-2,jj,k) != 1.0)) {
557  continue;
558  }
559  }
560 
561  if(on_y_face){
562  if (jj < domlo_y && (vfrac(ii,jj+1,k) != 1.0 || vfrac(ii,jj+2,k) != 1.0) ) {
563  continue;
564  }
565  if (jj > domhi_y && (vfrac(ii,jj-1,k) != 1.0 || vfrac(ii,jj-2,k) != 1.0) ) {
566  continue;
567  }
568  }
569 
570  Amatrix(a_ind,0) = 1.0;
571  Amatrix(a_ind,1) = x_off + ccent(ii,jj,k,0);
572  Amatrix(a_ind,2) = y_off + ccent(ii,jj,k,1) - yloc_on_xface;
573  Amatrix(a_ind,3) = Amatrix(a_ind,1) * Amatrix(a_ind,1);
574  Amatrix(a_ind,4) = Amatrix(a_ind,1) * Amatrix(a_ind,2);
575  Amatrix(a_ind,5) = Amatrix(a_ind,2) * Amatrix(a_ind,2);
576 
577  // Add in information about the location of the EB. Exclude
578  // EBs that are outside the domain.
579  if (flag(ii,jj,k).isSingleValued() &&
580  domlo_x <= ii && ii <= domhi_x &&
581  domlo_y <= jj && jj <= domhi_y)
582  {
583  int a_ind_eb = a_ind + 9;
584  AMREX_ASSERT(9<= a_ind_eb && a_ind_eb <= 17);
585 
586  Amatrix(a_ind_eb,0) = 1.0;
587  Amatrix(a_ind_eb,1) = x_off + bcent(ii,jj,k,0);
588  Amatrix(a_ind_eb,2) = y_off + bcent(ii,jj,k,1) - yloc_on_xface;
589  Amatrix(a_ind_eb,3) = Amatrix(a_ind_eb,1) * Amatrix(a_ind_eb,1);
590  Amatrix(a_ind_eb,4) = Amatrix(a_ind_eb,1) * Amatrix(a_ind_eb,2);
591  Amatrix(a_ind_eb,5) = Amatrix(a_ind_eb,2) * Amatrix(a_ind_eb,2);
592  }
593  }
594  }
595  }
596 
597  for (int irow = 0; irow < 6; irow++)
598  {
599  rhs(irow) = 0.; // Only non-zero when inhomogeneous Dirichlet
600 
601  for (int ii = i-im; ii <= i+ip; ii++) { // Normal to face
602  for (int jj = j-1; jj <= j+1; jj++) // Tangential to face
603  {
604  if ( !phi.contains(ii,jj,k) ) {
605  continue;
606  }
607 
608  if (!flag(ii,jj,k).isCovered())
609  {
610  int a_ind = (jj-(j-1)) + 3*(ii-(i-im));
611  AMREX_ASSERT(0<=a_ind && a_ind <9);
612 
613  Real phi_val = Amatrix(a_ind,0) * phi(ii,jj,k,n);
614 
615  rhs(irow) += Amatrix(a_ind ,irow)*phi_val;
616 
617  if (flag(ii,jj,k).isSingleValued() &&
618  is_eb_dirichlet && is_eb_inhomog && Amatrix(a_ind+9,irow) != 0.0) {
619  rhs(irow) += Amatrix(a_ind+9,irow)*phieb(ii,jj,k,n);
620  }
621  }
622  }
623  }
624  }
625 
626  cholsol_for_eb(Amatrix, rhs);
627 
628  return rhs(1);
629 }
630 
631 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
632 Real grad_y_of_phi_on_centroids_extdir(int i,int j,int k,int n,
633  Array4<Real const> const& phi,
634  Array4<Real const> const& phieb,
635  Array4<EBCellFlag const> const& flag,
636  Array4<Real const> const& ccent,
637  Array4<Real const> const& bcent,
638  Array4<Real const> const& vfrac,
639  Real& xloc_on_yface,
640  bool is_eb_dirichlet, bool is_eb_inhomog,
641  const bool on_x_face, const int domlo_x, const int domhi_x,
642  const bool on_y_face, const int domlo_y, const int domhi_y)
643 {
644  AMREX_ALWAYS_ASSERT(is_eb_dirichlet);
645  Array2D<Real,0,17,0,5> Amatrix;
646  Array1D<Real,0, 5 > rhs;
647 
648  // Order of column -- first six are cell centroids: (i-1,j) (i,j) (i+1,j) (i-1,j-1) (i,j-1) (i+1,j-1)
649  // Order of column -- second six are EB centroids: (i-1,j) (i,j) (i+1,j) (i-1,j-1) (i,j-1) (i+1,j-1)
650 
651  for (int irow = 0; irow < 18; irow++) {
652  for (int icol = 0; icol < 6; icol++) {
653  Amatrix(irow,icol) = 0.0;
654  }
655  }
656 
657  const int jm = (j > domhi_y) ? 2 : 1;
658  const int jp = 2 - jm;
659 
660  // Columns: [e x y x*x x*y y*y]
661  for (int jj = j-jm; jj <= j+jp; jj++) { // Normal to face
662  for (int ii = i-1; ii <= i+1; ii++) // Tangential to face
663  {
664 
665  // Don't include corner cells.
666  if (((on_x_face && ii < domlo_x) && (on_y_face && jj < domlo_y)) ||
667  ((on_x_face && ii > domhi_x) && (on_y_face && jj > domhi_y)) ||
668  ((on_x_face && ii < domlo_x) && (on_y_face && jj > domhi_y)) ||
669  ((on_x_face && ii > domhi_x) && (on_y_face && jj < domlo_y))) {
670  continue;
671  }
672 
673  if ( !phi.contains(ii,jj,k) ) {
674  continue;
675  }
676 
677  if (!flag(ii,jj,k).isCovered())
678  {
679  int a_ind = (ii-(i-1)) + 3*(jj-(j-jm));
680  AMREX_ASSERT(0<= a_ind && a_ind <= 8);
681 
682  Real x_off = static_cast<Real>(ii-i);
683  Real y_off = static_cast<Real>(jj-j) + 0.5;
684 
685  if(on_x_face){
686  if (ii < domlo_x && (vfrac(ii+1,jj,k) != 1.0 || vfrac(ii+2,jj,k) != 1.0) ) {
687  continue;
688  }
689  if (ii > domhi_x && (vfrac(ii-1,jj,k) != 1.0 || vfrac(ii-2,jj,k) != 1.0)) {
690  continue;
691  }
692  }
693 
694  if(on_y_face){
695  if (jj < domlo_y && (vfrac(ii,jj+1,k) != 1.0 || vfrac(ii,jj+2,k) != 1.0) ) {
696  continue;
697  }
698  if (jj > domhi_y && (vfrac(ii,jj-1,k) != 1.0 || vfrac(ii,jj-2,k) != 1.0) ) {
699  continue;
700  }
701  }
702 
703  Amatrix(a_ind,0) = 1.0;
704  Amatrix(a_ind,1) = x_off + ccent(ii,jj,k,0) - xloc_on_yface;
705  Amatrix(a_ind,2) = y_off + ccent(ii,jj,k,1);
706  Amatrix(a_ind,3) = Amatrix(a_ind,1) * Amatrix(a_ind,1);
707  Amatrix(a_ind,4) = Amatrix(a_ind,1) * Amatrix(a_ind,2);
708  Amatrix(a_ind,5) = Amatrix(a_ind,2) * Amatrix(a_ind,2);
709 
710  if (flag(ii,jj,k).isSingleValued() &&
711  domlo_x <= ii && ii <= domhi_x &&
712  domlo_y <= jj && jj <= domhi_y)
713  {
714 
715  int a_ind_eb = a_ind + 9;
716  AMREX_ASSERT(9<= a_ind_eb && a_ind_eb <= 17);
717 
718  Amatrix(a_ind_eb,0) = 1.0;
719  Amatrix(a_ind_eb,1) = x_off + bcent(ii,jj,k,0) - xloc_on_yface;
720  Amatrix(a_ind_eb,2) = y_off + bcent(ii,jj,k,1);
721  Amatrix(a_ind_eb,3) = Amatrix(a_ind_eb,1) * Amatrix(a_ind_eb,1);
722  Amatrix(a_ind_eb,4) = Amatrix(a_ind_eb,1) * Amatrix(a_ind_eb,2);
723  Amatrix(a_ind_eb,5) = Amatrix(a_ind_eb,2) * Amatrix(a_ind_eb,2);
724  }
725  }
726  }
727  }
728 
729 
730  // Make the RHS = A^T v
731  for (int irow = 0; irow < 6; irow++)
732  {
733  rhs(irow) = 0.; // Only non-zero when inhomogeneous Dirichlet
734 
735  for (int jj = j-jm; jj <= j+jp; jj++) { // Normal to face
736  for (int ii = i-1; ii <= i+1; ii++) {// Tangential to face
737  if ( !phi.contains(ii,jj,k) ) {
738  continue;
739  }
740  if (!flag(ii,jj,k).isCovered())
741  {
742  int a_ind = (ii-(i-1)) + 3*(jj-(j-jm));
743  AMREX_ASSERT(0<=irow && irow < 9);
744 
745  Real phi_val = Amatrix(a_ind,0)*phi(ii,jj,k,n);
746 
747  rhs(irow) += Amatrix(a_ind,irow)*phi_val;
748 
749  if (flag(ii,jj,k).isSingleValued() &&
750  is_eb_dirichlet && is_eb_inhomog && Amatrix(a_ind+9,irow) != Real(0.0)){
751  rhs(irow) += Amatrix(a_ind+9,irow)*phieb(ii,jj,k,n);
752  }
753  }
754  }
755  }
756 
757  }
758 
759  cholsol_for_eb(Amatrix, rhs);
760 
761  return rhs(2);
762 }
763 
764 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
765 Real grad_eb_of_phi_on_centroids_extdir(int i,int j,int k,int n,
766  Array4<Real const> const& phi,
767  Array4<Real const> const& phieb,
768  Array4<EBCellFlag const> const& flag,
769  Array4<Real const> const& ccent,
770  Array4<Real const> const& bcent,
771  Array4<Real const> const& vfrac,
772  Real& nrmx, Real& nrmy,
773  bool is_eb_inhomog,
774  const bool on_x_face, const int domlo_x, const int domhi_x,
775  const bool on_y_face, const int domlo_y, const int domhi_y)
776 {
777  Array2D<Real,0,17,0,5> Amatrix;
778  Array1D<Real,0, 5 > rhs;
779 
780  // Order of column -- first 9 are cell centroids, next 9 are EB centroids
781 
782  for (int irow = 0; irow < 18; irow++) {
783  for (int icol = 0; icol < 6; icol++) {
784  Amatrix(irow,icol) = 0.0;
785  }
786  }
787 
788  // Column 0-2: [e x y]
789  for (int ii = i-1; ii <= i+1; ii++) {
790  for (int jj = j-1; jj <= j+1; jj++)
791  {
792 
793  // This is likely overkill for EB grads.
794  // Don't include corner cells. Could make this even more strict
795  // by removing the on_?_face restrictions.
796  if (((on_x_face && ii < domlo_x) && (on_y_face && jj < domlo_y)) ||
797  ((on_x_face && ii > domhi_x) && (on_y_face && jj > domhi_y))){
798  continue;
799  }
800 
801 
802  if (!flag(ii,jj,k).isCovered())
803  {
804  int a_ind = (jj-(j-1)) + 3*(ii-(i-1));
805 
806  Real x_off = static_cast<Real>(ii-i);
807  Real y_off = static_cast<Real>(jj-j);
808 
809  if(on_x_face){
810  if (ii < domlo_x && (vfrac(ii+1,jj,k) != 1.0 || vfrac(ii+2,jj,k) != 1.0) ) {
811  continue;
812  }
813  if (ii > domhi_x && (vfrac(ii-1,jj,k) != 1.0 || vfrac(ii-2,jj,k) != 1.0)) {
814  continue;
815  }
816  }
817 
818  if(on_y_face){
819  if (jj < domlo_y && (vfrac(ii,jj+1,k) != 1.0 || vfrac(ii,jj+2,k) != 1.0) ) {
820  continue;
821  }
822  if (jj > domhi_y && (vfrac(ii,jj-1,k) != 1.0 || vfrac(ii,jj-2,k) != 1.0) ) {
823  continue;
824  }
825  }
826 
827 
828  if (flag(i,j,k).isConnected((ii-i),(jj-j),0))
829  {
830  Amatrix(a_ind,0) = 1.0;
831  Amatrix(a_ind,1) = x_off + ccent(ii,jj,k,0) - bcent(i,j,k,0);
832  Amatrix(a_ind,2) = y_off + ccent(ii,jj,k,1) - bcent(i,j,k,1);
833  }
834 
835  // Include information from the EB unless it is outside the domain
836  if (flag(i,j,k).isConnected((ii-i),(jj-j),0) && flag(ii,jj,k).isSingleValued() &&
837  domlo_x <= ii && ii <= domhi_x && domlo_y <= jj && jj <= domhi_y)
838  {
839  Amatrix(a_ind+9,0) = 1.0;
840  Amatrix(a_ind+9,1) = x_off + bcent(ii,jj,k,0) - bcent(i,j,k,0);
841  Amatrix(a_ind+9,2) = y_off + bcent(ii,jj,k,1) - bcent(i,j,k,1);
842  }
843  }
844  }
845  }
846 
847  // Columns 3 : [x*x x*y y*y]
848 
849  for (int irow = 0; irow < 18; irow++)
850  {
851  Amatrix(irow,3) = Amatrix(irow,1) * Amatrix(irow,1);
852  Amatrix(irow,4) = Amatrix(irow,1) * Amatrix(irow,2);
853  Amatrix(irow,5) = Amatrix(irow,2) * Amatrix(irow,2);
854  }
855 
856  // Make the RHS = A^T v
857  for (int irow = 0; irow < 6; irow++)
858  {
859  rhs(irow) = 0.;
860 
861  for (int ii = i-1; ii <= i+1; ii++) {
862  for (int jj = j-1; jj <= j+1; jj++) {
863  if ( !phi.contains(ii,jj,k) ) {
864  continue;
865  }
866  if (!flag(ii,jj,k).isCovered())
867  {
868  int a_ind = (jj-(j-1)) + 3*(ii-(i-1));
869  AMREX_ASSERT(0<=a_ind && a_ind <9);
870 
871  // A(a_ind,0) is 1 or 0. Use it to mask phi.
872  Real phi_val = Amatrix(a_ind,0)*phi(ii,jj,k,n);
873 
874  rhs(irow) += Amatrix(a_ind,irow) * phi_val;
875 
876  if (flag(ii,jj,k).isSingleValued() && is_eb_inhomog && Amatrix(a_ind+9,irow) != 0.0) {
877  rhs(irow) += Amatrix(a_ind+9,irow)*phieb(ii,jj,k,n);
878  }
879  }
880  }
881  }
882  }
883 
884  cholsol_for_eb(Amatrix, rhs);
885 
886  Real dphidn = rhs(1)*nrmx + rhs(2)*nrmy;
887  return dphidn;
888 }
889 
890 }
891 #endif
AMREX_ASSERT
#define AMREX_ASSERT(EX)
Definition: AMReX_BLassert.H:38
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::decomp_chol_np6
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void decomp_chol_np6(Array2D< Real, 0, 5, 0, 5 > &aa)
Definition: AMReX_EB_LeastSquares_2D_K.H:10
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_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_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_np6
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void cholsol_np6(Array2D< Real, 0, 11, 0, 5 > &Amatrix, Array1D< Real, 0, 5 > &b)
Definition: AMReX_EB_LeastSquares_2D_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::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