Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
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
6namespace amrex {
7
8// This is the 2D version (i.e. for 6x6 (A^T A) matrix)
9AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
10void 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)
63AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
64void 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
154AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
155void 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
245AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
246Real 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
330AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
331Real 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
411AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
412Real 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
494AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
495Real 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
631AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
632Real 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
764AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
765Real 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::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