Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
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
6namespace amrex {
7
8// This is the 3D version (i.e. for 10x10 (A^T A) matrix)
9AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
10void 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)
63AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
64void 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
198AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
199void 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
332AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
333Real 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
422AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
423Real 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
512AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
513Real 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
604AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
605Real 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
695AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
696Real 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
853AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
854Real 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
1009AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1010Real 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
1166AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
1167Real 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::Array1D
Definition AMReX_Array.H:161
amrex::Array2D
Definition AMReX_Array.H:282
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