Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
AMReX_HypreMLABecLap_1D_K.H
Go to the documentation of this file.
1#ifndef AMREX_HYPRE_ML_ABECLAP_2D_K_H_
2#define AMREX_HYPRE_ML_ABECLAP_2D_K_H_
3
4#include <AMReX_Array.H>
5#include <AMReX_Orientation.H>
6
7namespace amrex {
8
9AMREX_GPU_DEVICE AMREX_FORCE_INLINE
10void hypmlabeclap_f2c_set_values (IntVect const& cell, Real* values,
11 GpuArray<Real,AMREX_SPACEDIM> const& dx, Real sb,
12 GpuArray<Array4<Real const>,AMREX_SPACEDIM> const& b,
13 GpuArray<Array4<int const>,AMREX_SPACEDIM*2> const& bmask,
14 IntVect const& refratio, int not_covered)
15{
16 Array2D<Real,-1,1,-1,1> tmp;
17 for (auto& x : tmp) { x = Real(0.0); }
18
19 Array2D<int,-1,1,-1,1> used;
20 for (auto& x : used) { x = 0; }
21
22 for (OrientationIter ori; ori; ++ori) {
23 auto const face = ori();
24 int const idir = face.coordDir();
25 int const idir1 = 1-idir; // the transverse direction
26 IntVect offset(0);
27 offset[idir] = face.isLow() ? -1 : +1;
28 IntVect const cell_out = cell + offset;
29 auto const& msk = bmask[face];
30 if (msk.contains(cell_out) && msk(cell_out) == not_covered) {
31 // There is a coarse cell on the other side of the face. There
32 // are three cases for the coarse cells involved. (1)
33 // Interpolation using 3 coarse cells. (2) Upward biased
34 // interpolation using 2 coarse cells. (3) Doward biased
35 // interpolation using 2 coarse cells. (Here up and down means
36 // the y-dirction, if we assume idir is the x-direction.)
37 IntVect offset_t(0);
38 offset_t[idir1] = refratio[idir1];
39 Real bcoeff = b[idir] ? b[idir](face.isLow() ? cell : cell_out) : Real(1.0); // b is on face
40 Real poly_coef[3];
41 {
42 Real xx[3] = {Real(-0.5)*Real(refratio[idir]), Real(0.5), Real(1.5)};
43 poly_interp_coeff<3>(Real(-0.5), xx, poly_coef);
44 }
45 Real fac = -(sb / (dx[idir]*dx[idir])) * bcoeff * poly_coef[0];
46 int const rr1 = refratio[idir1];
47 int const i1 = cell[idir1];
48 int const i1c = amrex::coarsen(i1, rr1);
49 Real xInt = Real(-0.5) + (i1-i1c*rr1+Real(0.5))/Real(rr1);
50 Real xc[] = {Real(-1.0), Real(0.0), Real(1.0)};
51 Real c[] = {Real(0.0), Real(0.0), Real(0.0)};
52 int cc[] = {0, 0, 0};
53 if (msk(cell_out-offset_t) == not_covered &&
54 msk(cell_out+offset_t) == not_covered)
55 {
56 poly_interp_coeff<3>(xInt, xc, c);
57 cc[0] = cc[1] = cc[2] = 1;
58 } else if (msk(cell_out+offset_t) == not_covered) {
59 poly_interp_coeff<2>(xInt, &(xc[1]), &(c[1]));
60 cc[1] = cc[2] = 1;
61 } else {
62 poly_interp_coeff<2>(xInt, xc, c);
63 cc[0] = cc[1] = 1;
64 }
65 if (face == Orientation(0, Orientation::low)) {
66 for (int m = 0; m < 3; ++m) {
67 tmp(-1,m-1) += c[m] * fac;
68 used(-1,m-1) += cc[m];
69 }
70 } else if (face == Orientation(0, Orientation::high)) {
71 for (int m = 0; m < 3; ++m) {
72 tmp(1,m-1) += c[m] * fac;
73 used(1,m-1) += cc[m];
74 }
75 } else if (face == Orientation(1, Orientation::low)) {
76 for (int m = 0; m < 3; ++m) {
77 tmp(m-1,-1) += c[m] * fac;
78 used(m-1,-1) += cc[m];
79 }
80 } else if (face == Orientation(1, Orientation::high)) {
81 for (int m = 0; m < 3; ++m) {
82 tmp(m-1,1) += c[m] * fac;
83 used(m-1,1) += cc[m];
84 }
85 }
86 }
87 }
88
89 auto const* ptmp = tmp.begin();
90 auto const* pused = used.begin();
91 for (int m = 0; m < 9; ++m) {
92 if (pused[m]) {
93 (*values) += ptmp[m];
94 ++values;
95 }
96 }
97}
98
99AMREX_GPU_DEVICE AMREX_FORCE_INLINE
100void hypmlabeclap_c2f (int i, int j, int k,
101 Array4<GpuArray<Real,2*AMREX_SPACEDIM+1>> const& stencil,
102 GpuArray<HYPRE_Int,AMREX_SPACEDIM>* civ, HYPRE_Int* nentries,
103 int* entry_offset, Real* entry_values,
104 Array4<int const> const& offset_from,
105 Array4<int const> const& nentries_to,
106 Array4<int const> const& offset_to,
107 GpuArray<Real,AMREX_SPACEDIM> const& dx, Real sb,
108 Array4<int const> const& offset_bx,
109 Array4<int const> const& offset_by,
110 Real const* bx, Real const* by,
111 Array4<int const> const& fine_mask,
112 IntVect const& rr, Array4<int const> const& crse_mask)
113{
114 if (fine_mask(i,j,k)) {
115 // Let's set off-diagonal elements to zero
116 for (int m = 1; m < 2*AMREX_SPACEDIM+1; ++m) {
117 stencil(i,j,k)[m] = Real(0.0);
118 }
119 } else if (nentries_to(i,j,k) > 0) {
120 int const fromoff = offset_from(i,j,k);
121 civ[fromoff][0] = i;
122 civ[fromoff][1] = j;
123 nentries[fromoff] = nentries_to(i,j,k);
124 int foff = offset_to(i,j,k);
125 entry_offset[fromoff] = foff;
126
127 // We must iterate the faces in the lexicographical order of fine
128 // neighbor cells, because that's the order when non-stencil entries
129 // were added to Hypre's graph. Also note that a coarse cell will
130 // not have entries to fine cells at both ends of a direction. Thus
131 // we do not have to worry about the order between fine cells at the
132 // small and big ends of the same direction.
133
134 if (fine_mask(i,j-1,k)) {
135 stencil(i,j,k)[0] += stencil(i,j,k)[3];
136 stencil(i,j,k)[3] = Real(0.0);
137 // Reflux: sb/h^2*by(i,j,k)*(phi(i,j,k)-phi(i,j-1,k)) is replaced by
138 // sb/h*sum_{fine faces}(dphi/dy*by)/n_fine_faces
139 Real dyf = dx[1] / Real(rr[1]);
140 Real dycinv = Real(1.0) / dx[1];
141 Real dyfinv = Real(1.0) / dyf;
142 Real cc[3];
143 Real yy[3] = {dx[1]*Real(-0.5), dyf*Real(0.5), dyf*Real(1.5)};
144 poly_interp_coeff<3>(dyf*Real(-0.5), yy, cc);
145 for (int irx = 0; irx < rr[0]; ++irx) {
146 Real bym = by ? by[offset_by(i,j,k)+irx] : Real(1.0);
147 Real fac = sb*dycinv*dyfinv*bym/Real(rr[0]);
148 // int ii = i*rr[0] + irx
149 // int jj = j*rr[1]
150 // dphi/dy = (phi_interp - phi_fine(jj-1)) / dy_fine
151 // = (phi_coarse*cc[0] + phi_fine(jj-1)*(cc[1]-1)
152 // + phi_fine(jj-2)* cc[2]) / dy_fine
153 // So the entry for fine cell (jj-1) is (cc[1]-1)*fac
154 // fine cell (jj-2) is cc[2] *fac
155 entry_values[foff+irx ] += fac* cc[2];
156 entry_values[foff+irx+rr[0]] += fac*(cc[1]-Real(1.0));
157
158 // The coarse cell's stencils need updates too.
159 Real xInt = Real(-0.5) + (irx+Real(0.5))/Real(rr[0]);
160 Real xc[3] = {Real(-1.0), Real(0.0), Real(1.0)};
161 Real ct[3] = {Real(0.0), Real(0.0), Real(0.0)};
162 if ( fine_mask(i-1,j,k) ||
163 !crse_mask(i-1,j,k))
164 {
165 poly_interp_coeff<2>(xInt, &(xc[1]), &(ct[1]));
166 } else if ( fine_mask(i+1,j,k) ||
167 !crse_mask(i+1,j,k))
168 {
169 poly_interp_coeff<2>(xInt, xc, ct);
170 } else {
171 poly_interp_coeff<3>(xInt, xc, ct);
172 }
173 // phi_coarse = ct[0]*phi(i-1) + ct[1]*phi(i) + ct[2]*phi(i+1)
174 stencil(i,j,k)[0] += (fac*cc[0])*ct[1];
175 stencil(i,j,k)[1] += (fac*cc[0])*ct[0];
176 stencil(i,j,k)[2] += (fac*cc[0])*ct[2];
177 }
178 foff += 2*rr[0];
179 }
180
181 if (fine_mask(i-1,j,k)) {
182 stencil(i,j,k)[0] += stencil(i,j,k)[1];
183 stencil(i,j,k)[1] = Real(0.0);
184 // Reflux: sb/h^2*bx(i,j,k)*(phi(i,j,k)-phi(i-1,j,k)) is replaced by
185 // sb/h*sum_{fine faces}(dphi/dx*bx).
186 Real dxf = dx[0] / Real(rr[0]);
187 Real dxcinv = Real(1.0) / dx[0];
188 Real dxfinv = Real(1.0) / dxf;
189 Real cc[3];
190 Real xx[3] = {dx[0]*Real(-0.5), dxf*Real(0.5), dxf*Real(1.5)};
191 poly_interp_coeff<3>(dxf*Real(-0.5), xx, cc);
192 for (int iry = 0; iry < rr[1]; ++iry) {
193 Real bxm = bx ? bx[offset_bx(i,j,k)+iry] : Real(1.0);
194 Real fac = sb*dxcinv*dxfinv*bxm/Real(rr[1]);
195 // int ii = i*rr[0]
196 // int jj = j*rr[1] + iry
197 // dphi/dx = (phi_interp - phi_fine(ii-1)) / dx_fine
198 // = (phi_coarse*cc[0] + phi_fine(ii-1)*(cc[1]-1)
199 // + phi_fine(ii-2)* cc[2]) / dx_fine
200 // So the entry for fine cell(ii-1) is (cc[1]-1)*fac
201 // fine cell(ii-2) is cc[2] *fac
202 entry_values[foff++] = fac* cc[2];
203 entry_values[foff++] = fac*(cc[1] - Real(1.0));
204
205 // The coarse cell's stencils need updates too.
206 Real yInt = Real(-0.5) + (iry+Real(0.5))/Real(rr[1]);
207 Real yc[3] = {Real(-1.0), Real(0.0), Real(1.0)};
208 Real ct[3] = {Real(0.0), Real(0.0), Real(0.0)};
209 if ( fine_mask(i,j-1,k) ||
210 !crse_mask(i,j-1,k))
211 {
212 poly_interp_coeff<2>(yInt, &(yc[1]), &(ct[1]));
213 } else if ( fine_mask(i,j+1,k) ||
214 !crse_mask(i,j+1,k))
215 {
216 poly_interp_coeff<2>(yInt, yc, ct);
217 } else {
218 poly_interp_coeff<3>(yInt, yc, ct);
219 }
220 // phi_coarse = ct[0]*phi(j-1) + ct[1]*phi(j) + ct[2]*phi(j+1)
221 stencil(i,j,k)[0] += (fac*cc[0])*ct[1];
222 stencil(i,j,k)[3] += (fac*cc[0])*ct[0];
223 stencil(i,j,k)[4] += (fac*cc[0])*ct[2];
224 }
225 }
226
227 if (fine_mask(i+1,j,k)) {
228 stencil(i,j,k)[0] += stencil(i,j,k)[2];
229 stencil(i,j,k)[2] = Real(0.0);
230 // Reflux: sb/h^2*bx(i+1,j,k)*(phi(i,j,k)-phi(i+1,j,k)) is replaced by
231 // sb/h*sum_{fine faces}(-dphi/dx*bx).
232 Real dxf = dx[0] / Real(rr[0]);
233 Real dxcinv = Real(1.0) / dx[0];
234 Real dxfinv = Real(1.0) / dxf;
235 Real cc[3];
236 Real xx[3] = {dx[0]*Real(-0.5), dxf*Real(0.5), dxf*Real(1.5)};
237 poly_interp_coeff<3>(dxf*Real(-0.5), xx, cc);
238 for (int iry = 0; iry < rr[1]; ++iry) {
239 Real bxp = bx ? bx[offset_bx(i+1,j,k)+iry] : Real(1.0);
240 Real fac = sb*dxcinv*dxfinv*bxp/Real(rr[1]);
241 // int ii = i*rr[0] + (rr[0]-1)
242 // int jj = j*rr[1] + iry
243 // -dphi/dx = (phi_interp - phi_fine(ii+1)) / dx_fine
244 // = (phi_coarse*cc[0] + phi_fine(ii+1)*(cc[1]-1)
245 // + phi_fine(ii+2)* cc[2]) / dx_fine
246 // So the entry for fine cell(ii+1) is (cc[1]-1)*fac
247 // fine cell(ii+2) is cc[2] *fac
248 entry_values[foff++] = fac*(cc[1] - Real(1.0));
249 entry_values[foff++] = fac* cc[2];
250
251 // The coarse cell's stencils need updates too.
252 Real yInt = Real(-0.5) + (iry+Real(0.5))/Real(rr[1]);
253 Real yc[3] = {Real(-1.0), Real(0.0), Real(1.0)};
254 Real ct[3] = {Real(0.0), Real(0.0), Real(0.0)};
255 if ( fine_mask(i,j-1,k) ||
256 !crse_mask(i,j-1,k))
257 {
258 poly_interp_coeff<2>(yInt, &(yc[1]), &(ct[1]));
259 } else if ( fine_mask(i,j+1,k) ||
260 !crse_mask(i,j+1,k))
261 {
262 poly_interp_coeff<2>(yInt, yc, ct);
263 } else {
264 poly_interp_coeff<3>(yInt, yc, ct);
265 }
266 // phi_coarse = ct[0]*phi(j-1) + ct[1]*phi(j) + ct[2]*phi(j+1)
267 stencil(i,j,k)[0] += (fac*cc[0])*ct[1];
268 stencil(i,j,k)[3] += (fac*cc[0])*ct[0];
269 stencil(i,j,k)[4] += (fac*cc[0])*ct[2];
270 }
271 }
272
273 if (fine_mask(i,j+1,k)) {
274 stencil(i,j,k)[0] += stencil(i,j,k)[4];
275 stencil(i,j,k)[4] = Real(0.0);
276 // Reflux: sb/h^2*by(i,j+1,k)*(phi(i,j,k)-phi(i,j+1,k)) is replaced by
277 // sb/h*sum_{fine faces}(-dphi/dy*by)
278 Real dyf = dx[1] / Real(rr[1]);
279 Real dycinv = Real(1.0) / dx[1];
280 Real dyfinv = Real(1.0) / dyf;
281 Real cc[3];
282 Real yy[3] = {dx[1]*Real(-0.5), dyf*Real(0.5), dyf*Real(1.5)};
283 poly_interp_coeff<3>(dyf*Real(-0.5), yy, cc);
284 for (int irx = 0; irx < rr[0]; ++irx) {
285 Real byp = by ? by[offset_by(i,j+1,k)+irx] : Real(1.0);
286 Real fac = sb*dycinv*dyfinv*byp/Real(rr[0]);
287 // int ii = i*rr[0] + irx
288 // int jj = j*rr[1] + (rr[1]-1)
289 // -dphi/dy = (phi_interp - phi_fine(jj+1)) / dy_fine
290 // = (phi_coarse*cc[0] + phi_fine(jj+1)*(cc[1]-1)
291 // + phi_fine(jj+2)* cc[2]) / dy_fine
292 // So the entry for fine cell (jj+1) is (cc[1]-1)*fac
293 // fine cell (jj+2) is cc[2] *fac
294 entry_values[foff+irx ] += fac*(cc[1]-Real(1.0));
295 entry_values[foff+irx+rr[0]] += fac* cc[2];
296
297 // The coarse cell's stencils need updates too.
298 Real xInt = Real(-0.5) + (irx+Real(0.5))/Real(rr[0]);
299 Real xc[3] = {Real(-1.0), Real(0.0), Real(1.0)};
300 Real ct[3] = {Real(0.0), Real(0.0), Real(0.0)};
301 if ( fine_mask(i-1,j,k) ||
302 !crse_mask(i-1,j,k))
303 {
304 poly_interp_coeff<2>(xInt, &(xc[1]), &(ct[1]));
305 } else if ( fine_mask(i+1,j,k) ||
306 !crse_mask(i+1,j,k))
307 {
308 poly_interp_coeff<2>(xInt, xc, ct);
309 } else {
310 poly_interp_coeff<3>(xInt, xc, ct);
311 }
312 // phi_coarse = ct[0]*phi(i-1) + ct[1]*phi(i) + ct[2]*phi(i+1)
313 stencil(i,j,k)[0] += (fac*cc[0])*ct[1];
314 stencil(i,j,k)[1] += (fac*cc[0])*ct[0];
315 stencil(i,j,k)[2] += (fac*cc[0])*ct[2];
316 }
317 // not needed: foff += 2*rr[0];
318 }
319 }
320}
321
322}
323
324#endif
AMReX_Array.H
AMREX_FORCE_INLINE
#define AMREX_FORCE_INLINE
Definition AMReX_Extension.H:119
AMREX_GPU_DEVICE
#define AMREX_GPU_DEVICE
Definition AMReX_GpuQualifiers.H:18
idir
int idir
Definition AMReX_HypreMLABecLap.cpp:1093
offset
Array4< int const > offset
Definition AMReX_HypreMLABecLap.cpp:1089
AMReX_Orientation.H
amrex::OrientationIter
An Iterator over the Orientation of Faces of a Box.
Definition AMReX_Orientation.H:135
amrex::Orientation
Encapsulation of the Orientation of the Faces of a Box.
Definition AMReX_Orientation.H:29
amrex::Orientation::low
@ low
Definition AMReX_Orientation.H:34
amrex::Orientation::high
@ high
Definition AMReX_Orientation.H:34
amrex::Real
amrex_real Real
Floating Point Type for Fields.
Definition AMReX_REAL.H:79
amrex::coarsen
__host__ __device__ BoxND< dim > coarsen(const BoxND< dim > &b, int ref_ratio) noexcept
Coarsen BoxND by given (positive) coarsening ratio.
Definition AMReX_Box.H:1409
amrex
Definition AMReX_Amr.cpp:49
amrex::hypmlabeclap_c2f
__device__ void hypmlabeclap_c2f(int i, int j, int k, Array4< GpuArray< Real, 2 *3+1 > > const &stencil, GpuArray< HYPRE_Int, 3 > *civ, HYPRE_Int *nentries, int *entry_offset, Real *entry_values, Array4< int const > const &offset_from, Array4< int const > const &nentries_to, Array4< int const > const &offset_to, GpuArray< Real, 3 > const &dx, Real sb, Array4< int const > const &offset_bx, Array4< int const > const &offset_by, Real const *bx, Real const *by, Array4< int const > const &fine_mask, IntVect const &rr, Array4< int const > const &crse_mask)
Definition AMReX_HypreMLABecLap_1D_K.H:100
amrex::hypmlabeclap_f2c_set_values
__device__ void hypmlabeclap_f2c_set_values(IntVect const &cell, Real *values, GpuArray< Real, 3 > const &dx, Real sb, GpuArray< Array4< Real const >, 3 > const &b, GpuArray< Array4< int const >, 3 *2 > const &bmask, IntVect const &refratio, int not_covered)
Definition AMReX_HypreMLABecLap_1D_K.H:10
amrex::Array2D
Definition AMReX_Array.H:347
amrex::Array2D::begin
__host__ __device__ const T * begin() const noexcept
Definition AMReX_Array.H:406
amrex::Array4
Definition AMReX_Array4.H:61
amrex::GpuArray
Fixed-size array that can be used on GPU.
Definition AMReX_Array.H:40