Block-Structured AMR Software Framework
AMReX_HypreMLABecLap_K.H
Go to the documentation of this file.
1 #ifndef AMREX_HYPRE_ML_ABECLAP_K_H_
2 #define AMREX_HYPRE_ML_ABECLAP_K_H_
3 #include <AMReX_Config.H>
4 
5 #include <AMReX_Array4.H>
6 #include <AMReX_LO_BCTYPES.H>
7 #include <AMReX_LOUtil_K.H>
8 #include <AMReX_REAL.H>
9 
10 #include <HYPRE_utilities.h>
11 
12 namespace amrex {
13 
14 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
15 void hypmlabeclap_mat (GpuArray<Real,2*AMREX_SPACEDIM+1>& sten, int i, int j, int k,
16  Dim3 const& boxlo, Dim3 const& boxhi,
17  Real sa, Array4<Real const> const& a,
18  Real sb, GpuArray<Real,AMREX_SPACEDIM> const& dx,
19  GpuArray<Array4<Real const>, AMREX_SPACEDIM> const& b,
20  GpuArray<int,AMREX_SPACEDIM*2> const& bctype,
21  GpuArray<Real,AMREX_SPACEDIM*2> const& bcl,
22  GpuArray<Array4<int const>, AMREX_SPACEDIM*2> const& bcmsk,
23  GpuArray<Array4<Real const>, AMREX_SPACEDIM*2> const& bcval,
24  GpuArray<Array4<Real>, AMREX_SPACEDIM*2> const& bcrhs,
25  int level, IntVect const& fixed_pt)
26 {
27  Real bxm = b[0] ? b[0](i ,j ,k ) : Real(1.0);
28  Real bxp = b[0] ? b[0](i+1,j ,k ) : Real(1.0);
29  Real bym = b[1] ? b[1](i ,j ,k ) : Real(1.0);
30  Real byp = b[1] ? b[1](i ,j+1,k ) : Real(1.0);
31 #if (AMREX_SPACEDIM > 2)
32  Real bzm = b[2] ? b[2](i ,j ,k ) : Real(1.0);
33  Real bzp = b[2] ? b[2](i ,j ,k+1) : Real(1.0);
34 #endif
35  Real ac = a ? a(i,j,k) : Real(0.0);
36 
37  sten[1] = -(sb / (dx[0]*dx[0])) * bxm;
38  sten[2] = -(sb / (dx[0]*dx[0])) * bxp;
39  sten[3] = -(sb / (dx[1]*dx[1])) * bym;
40  sten[4] = -(sb / (dx[1]*dx[1])) * byp;
41 #if (AMREX_SPACEDIM == 2)
42  sten[0] = -(sten[1] + sten[2] + sten[3] + sten[4]) + sa*ac;
43 #else
44  sten[5] = -(sb / (dx[2]*dx[2])) * bzm;
45  sten[6] = -(sb / (dx[2]*dx[2])) * bzp;
46  sten[0] = -(sten[1] + sten[2] + sten[3] + sten[4] + sten[5] + sten[6]) + sa*ac;
47 #endif
48 
49  // xlo
50  if (i == boxlo.x) {
51  int const cdir = Orientation(Direction::x, Orientation::low);
52  int const bcmk = bcmsk[cdir](i-1,j,k);
53  if (bcmk > 0) {
54  int const bct = bctype[cdir];
55  Real cc[3];
56  if (bct == AMREX_LO_DIRICHLET) {
57  Real xx[3] = {-bcl[cdir], dx[0]*Real(0.5), dx[0]*Real(1.5)};
58  poly_interp_coeff<3>(dx[0]*Real(-0.5), xx, cc);
59  } else { // Neumann
60  cc[0] = Real(0.0);
61  cc[1] = Real(1.0);
62  cc[2] = Real(0.0);
63  }
64  Real fac = (sb / (dx[0]*dx[0])) * bxm;
65  if (bct == AMREX_LO_DIRICHLET && (level == 0 || bcmk == 2)) {
66  // bcmk == 2 means outside the domain.
67  // We need to modify RHS at external Dirichlet boundaries.
68  bcrhs[cdir](i,j,k) = fac * cc[0] * bcval[cdir](i-1,j,k);
69  } else {
70  bcrhs[cdir](i,j,k) = Real(0.0);
71  }
72  sten[0] -= fac * cc[1];
73  sten[1] = Real(0.0);
74  sten[2] -= fac * cc[2];
75  }
76  }
77 
78  // xhi
79  if (i == boxhi.x) {
80  int const cdir = Orientation(Direction::x, Orientation::high);
81  int const bcmk = bcmsk[cdir](i+1,j,k);
82  if (bcmk > 0) {
83  int const bct = bctype[cdir];
84  Real cc[3];
85  if (bct == AMREX_LO_DIRICHLET) {
86  Real xx[3] = {-bcl[cdir], dx[0]*Real(0.5), dx[0]*Real(1.5)};
87  poly_interp_coeff<3>(dx[0]*Real(-0.5), xx, cc);
88  } else { // Neumann
89  cc[0] = Real(0.0);
90  cc[1] = Real(1.0);
91  cc[2] = Real(0.0);
92  }
93  Real fac = (sb / (dx[0]*dx[0])) * bxp;
94  if (bct == AMREX_LO_DIRICHLET && (level == 0 || bcmk == 2)) {
95  // bcmk == 2 means outside the domain.
96  // We need to modify RHS at external Dirichlet boundaries.
97  bcrhs[cdir](i,j,k) = fac * cc[0] * bcval[cdir](i+1,j,k);
98  } else {
99  bcrhs[cdir](i,j,k) = Real(0.0);
100  }
101  sten[0] -= fac * cc[1];
102  sten[1] -= fac * cc[2];
103  sten[2] = Real(0.0);
104  }
105  }
106 
107  // ylo
108  if (j == boxlo.y) {
109  int const cdir = Orientation(Direction::y, Orientation::low);
110  int const bcmk = bcmsk[cdir](i,j-1,k);
111  if (bcmk > 0) {
112  int const bct = bctype[cdir];
113  Real cc[3];
114  if (bct == AMREX_LO_DIRICHLET) {
115  Real xx[3] = {-bcl[cdir], dx[1]*Real(0.5), dx[1]*Real(1.5)};
116  poly_interp_coeff<3>(dx[1]*Real(-0.5), xx, cc);
117  } else { // Neumann
118  cc[0] = Real(0.0);
119  cc[1] = Real(1.0);
120  cc[2] = Real(0.0);
121  }
122  Real fac = (sb / (dx[1]*dx[1])) * bym;
123  if (bct == AMREX_LO_DIRICHLET && (level == 0 || bcmk == 2)) {
124  // bcmk == 2 means outside the domain.
125  // We need to modify RHS at external Dirichlet boundaries.
126  bcrhs[cdir](i,j,k) = fac * cc[0] * bcval[cdir](i,j-1,k);
127  } else {
128  bcrhs[cdir](i,j,k) = Real(0.0);
129  }
130  sten[0] -= fac * cc[1];
131  sten[3] = Real(0.0);
132  sten[4] -= fac * cc[2];
133  }
134  }
135 
136  // yhi
137  if (j == boxhi.y) {
138  int const cdir = Orientation(Direction::y, Orientation::high);
139  int const bcmk = bcmsk[cdir](i,j+1,k);
140  if (bcmk > 0) {
141  int const bct = bctype[cdir];
142  Real cc[3];
143  if (bct == AMREX_LO_DIRICHLET) {
144  Real xx[3] = {-bcl[cdir], dx[1]*Real(0.5), dx[1]*Real(1.5)};
145  poly_interp_coeff<3>(dx[1]*Real(-0.5), xx, cc);
146  } else { // Neumann
147  cc[0] = Real(0.0);
148  cc[1] = Real(1.0);
149  cc[2] = Real(0.0);
150  }
151  Real fac = (sb / (dx[1]*dx[1])) * byp;
152  if (bct == AMREX_LO_DIRICHLET && (level == 0 || bcmk == 2)) {
153  // bcmk == 2 means outside the domain.
154  // We need to modify RHS at external Dirichlet boundaries.
155  bcrhs[cdir](i,j,k) = fac * cc[0] * bcval[cdir](i,j+1,k);
156  } else {
157  bcrhs[cdir](i,j,k) = Real(0.0);
158  }
159  sten[0] -= fac * cc[1];
160  sten[3] -= fac * cc[2];
161  sten[4] = Real(0.0);
162  }
163  }
164 
165 #if (AMREX_SPACEDIM > 2)
166 
167  // zlo
168  if (k == boxlo.z) {
169  int const cdir = Orientation(Direction::z, Orientation::low);
170  int const bcmk = bcmsk[cdir](i,j,k-1);
171  if (bcmk > 0) {
172  int const bct = bctype[cdir];
173  Real cc[3];
174  if (bct == AMREX_LO_DIRICHLET) {
175  Real xx[3] = {-bcl[cdir], dx[2]*Real(0.5), dx[2]*Real(1.5)};
176  poly_interp_coeff<3>(dx[2]*Real(-0.5), xx, cc);
177  } else { // Neumann
178  cc[0] = Real(0.0);
179  cc[1] = Real(1.0);
180  cc[2] = Real(0.0);
181  }
182  Real fac = (sb / (dx[2]*dx[2])) * bzm;
183  if (bct == AMREX_LO_DIRICHLET && (level == 0 || bcmk == 2)) {
184  // bcmk == 2 means outside the domain.
185  // We need to modify RHS at external Dirichlet boundaries.
186  bcrhs[cdir](i,j,k) = fac * cc[0] * bcval[cdir](i,j,k-1);
187  } else {
188  bcrhs[cdir](i,j,k) = Real(0.0);
189  }
190  sten[0] -= fac * cc[1];
191  sten[5] = Real(0.0);
192  sten[6] -= fac * cc[2];
193  }
194  }
195 
196  // zhi
197  if (k == boxhi.z) {
198  int const cdir = Orientation(Direction::z, Orientation::high);
199  int const bcmk = bcmsk[cdir](i,j,k+1);
200  if (bcmk > 0) {
201  int const bct = bctype[cdir];
202  Real cc[3];
203  if (bct == AMREX_LO_DIRICHLET) {
204  Real xx[3] = {-bcl[cdir], dx[2]*Real(0.5), dx[2]*Real(1.5)};
205  poly_interp_coeff<3>(dx[2]*Real(-0.5), xx, cc);
206  } else { // Neumann
207  cc[0] = Real(0.0);
208  cc[1] = Real(1.0);
209  cc[2] = Real(0.0);
210  }
211  Real fac = (sb / (dx[2]*dx[2])) * bzp;
212  if (bct == AMREX_LO_DIRICHLET && (level == 0 || bcmk == 2)) {
213  // bcmk == 2 means outside the domain.
214  // We need to modify RHS at external Dirichlet boundaries.
215  bcrhs[cdir](i,j,k) = fac * cc[0] * bcval[cdir](i,j,k+1);
216  } else {
217  bcrhs[cdir](i,j,k) = Real(0.0);
218  }
219  sten[0] -= fac * cc[1];
220  sten[5] -= fac * cc[2];
221  sten[6] = Real(0.0);
222  }
223  }
224 
225 #endif
226 
227  if (fixed_pt == IntVect(AMREX_D_DECL(i,j,k))) {
228  for (int n = 1; n < 2*AMREX_SPACEDIM+1; ++n) {
229  sten[n] = Real(0.0);
230  }
231  }
232 }
233 
234 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
235 void hypmlabeclap_rhs (int i, int j, int k, Dim3 const& boxlo, Dim3 const& boxhi,
236  Array4<Real> const& rhs1,
237  Array4<Real const> const& rhs0,
238  GpuArray<Array4<int const>, AMREX_SPACEDIM*2> const& bcmsk,
239  GpuArray<Array4<Real const>, AMREX_SPACEDIM*2> const& bcrhs)
240 {
241  rhs1(i,j,k) = rhs0(i,j,k);
242 
243  // xlo
244  if (i == boxlo.x) {
245  int cdir = Orientation(Direction::x, Orientation::low);
246  if (bcmsk[cdir](i-1,j,k) > 0) {
247  rhs1(i,j,k) += bcrhs[cdir](i,j,k);
248  }
249  }
250 
251  // xhi
252  if (i == boxhi.x) {
253  int cdir = Orientation(Direction::x, Orientation::high);
254  if (bcmsk[cdir](i+1,j,k) > 0) {
255  rhs1(i,j,k) += bcrhs[cdir](i,j,k);
256  }
257  }
258 
259  // ylo
260  if (j == boxlo.y) {
261  int cdir = Orientation(Direction::y, Orientation::low);
262  if (bcmsk[cdir](i,j-1,k) > 0) {
263  rhs1(i,j,k) += bcrhs[cdir](i,j,k);
264  }
265  }
266 
267  // yhi
268  if (j == boxhi.y) {
269  int cdir = Orientation(Direction::y, Orientation::high);
270  if (bcmsk[cdir](i,j+1,k) > 0) {
271  rhs1(i,j,k) += bcrhs[cdir](i,j,k);
272  }
273  }
274 
275 #if (AMREX_SPACEDIM > 2)
276 
277  // zlo
278  if (k == boxlo.z) {
279  int cdir = Orientation(Direction::z, Orientation::low);
280  if (bcmsk[cdir](i,j,k-1) > 0) {
281  rhs1(i,j,k) += bcrhs[cdir](i,j,k);
282  }
283  }
284 
285  // zhi
286  if (k == boxhi.z) {
287  int cdir = Orientation(Direction::z, Orientation::high);
288  if (bcmsk[cdir](i,j,k+1) > 0) {
289  rhs1(i,j,k) += bcrhs[cdir](i,j,k);
290  }
291  }
292 
293 #endif
294 }
295 
296 }
297 
298 #if (AMREX_SPACEDIM == 2)
299 #include <AMReX_HypreMLABecLap_2D_K.H>
300 #else
301 #include <AMReX_HypreMLABecLap_3D_K.H>
302 #endif
303 
304 #endif
AMReX_Array4.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_HypreMLABecLap_2D_K.H
AMReX_HypreMLABecLap_3D_K.H
AMReX_LOUtil_K.H
AMReX_LO_BCTYPES.H
AMREX_LO_DIRICHLET
#define AMREX_LO_DIRICHLET
Definition: AMReX_LO_BCTYPES.H:5
AMReX_REAL.H
AMREX_D_DECL
#define AMREX_D_DECL(a, b, c)
Definition: AMReX_SPACE.H:104
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
Definition: AMReX_Amr.cpp:49
amrex::hypmlabeclap_rhs
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void hypmlabeclap_rhs(int i, int j, int k, Dim3 const &boxlo, Dim3 const &boxhi, Array4< Real > const &rhs1, Array4< Real const > const &rhs0, GpuArray< Array4< int const >, AMREX_SPACEDIM *2 > const &bcmsk, GpuArray< Array4< Real const >, AMREX_SPACEDIM *2 > const &bcrhs)
Definition: AMReX_HypreMLABecLap_K.H:235
amrex::IntVect
IntVectND< AMREX_SPACEDIM > IntVect
Definition: AMReX_BaseFwd.H:30
amrex::hypmlabeclap_mat
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void hypmlabeclap_mat(GpuArray< Real, 2 *AMREX_SPACEDIM+1 > &sten, int i, int j, int k, Dim3 const &boxlo, Dim3 const &boxhi, Real sa, Array4< Real const > const &a, Real sb, GpuArray< Real, AMREX_SPACEDIM > const &dx, GpuArray< Array4< Real const >, AMREX_SPACEDIM > const &b, GpuArray< int, AMREX_SPACEDIM *2 > const &bctype, GpuArray< Real, AMREX_SPACEDIM *2 > const &bcl, GpuArray< Array4< int const >, AMREX_SPACEDIM *2 > const &bcmsk, GpuArray< Array4< Real const >, AMREX_SPACEDIM *2 > const &bcval, GpuArray< Array4< Real >, AMREX_SPACEDIM *2 > const &bcrhs, int level, IntVect const &fixed_pt)
Definition: AMReX_HypreMLABecLap_K.H:15
amrex::Array4
Definition: AMReX_Array4.H:61
amrex::Dim3
Definition: AMReX_Dim3.H:12
amrex::Dim3::x
int x
Definition: AMReX_Dim3.H:12
amrex::Dim3::z
int z
Definition: AMReX_Dim3.H:12
amrex::Dim3::y
int y
Definition: AMReX_Dim3.H:12
amrex::GpuArray
Definition: AMReX_Array.H:33