Block-Structured AMR Software Framework
 
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
Loading...
Searching...
No Matches
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
12namespace amrex {
13
14AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
15void 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
234AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
235void 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::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:34