REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_BoundaryConditions_xvel.cpp
Go to the documentation of this file.
1#include "AMReX_PhysBCFunct.H"
2#include <REMORA_PhysBCFunct.H>
3
4using namespace amrex;
5
6/**
7 * @param[inout] dest_arr data on which to apply BCs
8 * @param[in ] bx box to update on
9 * @param[in ] domain domain box
10 * @param[in ] dxInv pm or pn
11 * @param[in ] msku land-sea mask on u-points
12 * @param[in ] calc_arr_in data to use in the RHS of calculations
13 * @param[in ] time current time
14 * @param[in ] bccomp index into both domain_bcs_type_bcr and bc_extdir_vals for icomp=0
15 */
16void REMORAPhysBCFunct::impose_xvel_bcs (const Array4<Real>& dest_arr, const Box& bx, const Box& domain,
17 const GpuArray<Real,AMREX_SPACEDIM> /*dxInv*/, const Array4<const Real>& msku,
18 const Array4<const Real>& calc_arr_in,
19 Real /*time*/, int bccomp)
20{
21 BL_PROFILE_VAR("impose_xvel_bcs()",impose_xvel_bcs);
22 const auto& dom_lo = amrex::lbound(domain);
23 const auto& dom_hi = amrex::ubound(domain);
24
25 // Based on BCRec for the domain, we need to make BCRec for this Box
26 // bccomp is used as starting index for m_domain_bcs_type
27 // 0 is used as starting index for bcrs
28 int ncomp = 1;
29 Vector<BCRec> bcrs(ncomp);
30 amrex::setBC(bx, domain, bccomp, 0, ncomp, m_domain_bcs_type, bcrs);
31
32 // xlo: ori = 0
33 // ylo: ori = 1
34 // zlo: ori = 2
35 // xhi: ori = 3
36 // yhi: ori = 4
37 // zhi: ori = 5
38
39 amrex::Gpu::DeviceVector<BCRec> bcrs_d(ncomp);
40#ifdef AMREX_USE_GPU
41 Gpu::htod_memcpy_async(bcrs_d.data(), bcrs.data(), sizeof(BCRec)*ncomp);
42#else
43 std::memcpy(bcrs_d.data(), bcrs.data(), sizeof(BCRec)*ncomp);
44#endif
45 const amrex::BCRec* bc_ptr = bcrs_d.data();
46 const auto* bc_extdir_vals_ptr = m_bc_extdir_vals_d.data();
47
48 GeometryData const& geomdata = m_geom.data();
49 bool is_periodic_in_x = geomdata.isPeriodic(0);
50 bool is_periodic_in_y = geomdata.isPeriodic(1);
51 const Real eps= Real(1.0e-20);
52 // If calc_arr isn't passed in, we actually just want to use dest_arr values instead.
53 // That matches ROMS.
54 Array4<const Real> calc_arr = calc_arr_in.ok() ? calc_arr_in : Array4<const Real>(dest_arr);
55
56 Box dest_arr_box = growHi(convert(Box(dest_arr),IntVect(1,0,0)),0,-1);
57 // First do all ext_dir bcs
58 if (!is_periodic_in_x or bccomp == BCVars::foextrap_bc(m_ncons))
59 {
60 Box bx_xlo(bx); bx_xlo.setBig (0,dom_lo.x-1);
61 Box bx_xhi(bx); bx_xhi.setSmall(0,dom_hi.x+2);
62 Box bx_xlo_face(bx); bx_xlo_face.setSmall(0,dom_lo.x ); bx_xlo_face.setBig(0,dom_lo.x );
63 Box bx_xhi_face(bx); bx_xhi_face.setSmall(0,dom_hi.x+1); bx_xhi_face.setBig(0,dom_hi.x+1);
64 ParallelFor(
65 // We only set the values on the domain faces themselves if EXT_DIR or actual outflow
66 grow(bx_xlo_face,IntVect(0,-1,0)) & dest_arr_box, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
67 if (bc_ptr[n].lo(0) == REMORABCType::ext_dir) {
68 dest_arr(i,j,k) = bc_extdir_vals_ptr[bccomp+n][0]*msku(i,j,0);
69 } else if (bc_ptr[n].lo(0) == REMORABCType::foextrap) {
70 dest_arr(i,j,k) = dest_arr(dom_lo.x+1,j,k)*msku(i,j,0);
71 } else if (bc_ptr[n].lo(0) == REMORABCType::orlanski_rad) {
72 Real grad_lo_ip1 = calc_arr(dom_lo.x+1,j ,k) - calc_arr(dom_lo.x+1,j-1,k);
73 Real grad_lo_ijp1 = calc_arr(dom_lo.x+1,j+1,k) - calc_arr(dom_lo.x+1,j ,k);
74 Real dUdt = calc_arr(dom_lo.x+1,j,k) - dest_arr(dom_lo.x+1,j,k);
75 Real dUdx = dest_arr(dom_lo.x+1,j,k) - dest_arr(dom_lo.x+2,j,k);
76 if (dUdt * dUdx < zero) dUdt = zero;
77 Real dUde = (dUdt * (grad_lo_ip1 + grad_lo_ijp1) > zero) ? grad_lo_ip1 : grad_lo_ijp1;
78 Real cff = std::max(dUdx*dUdx+dUde*dUde,eps);
79 Real Cx = dUdt*dUdx;
80 dest_arr(i,j,k) = (cff * calc_arr(i,j,k) + Cx * dest_arr(dom_lo.x+1,j,k)) * msku(i,j,0) / (cff+Cx);
81 }
82 });
83 ParallelFor(
84 grow(bx_xlo,IntVect(0,-1,0)) & dest_arr_box, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
85 int inner = (bc_ptr[n].lo(0) == REMORABCType::foextrap) ? 1 : 0;
86 int iflip = dom_lo.x - i;
87 if (bc_ptr[n].lo(0) == REMORABCType::ext_dir) {
88 dest_arr(i,j,k) = bc_extdir_vals_ptr[bccomp+n][0]*msku(i,j,0);
89 } else if (bc_ptr[n].lo(0) == REMORABCType::foextrap || bc_ptr[n].lo(0) == REMORABCType::clamped ||
90 bc_ptr[n].lo(0) == REMORABCType::orlanski_rad) {
91 dest_arr(i,j,k) = dest_arr(dom_lo.x+inner,j,k)*msku(i,j,0);
92 } else if (bc_ptr[n].lo(0) == REMORABCType::reflect_even) {
93 dest_arr(i,j,k) = dest_arr(iflip,j,k)*msku(i,j,0);
94 } else if (bc_ptr[n].lo(0) == REMORABCType::reflect_odd) {
95 dest_arr(i,j,k) = -dest_arr(iflip,j,k)*msku(i,j,0);
96 }
97 });
98 // We only set the values on the domain faces themselves if EXT_DIR or actual outflow
99 ParallelFor(
100 grow(bx_xhi_face,IntVect(0,-1,0)) & dest_arr_box, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
101 if (bc_ptr[n].hi(0) == REMORABCType::ext_dir) {
102 dest_arr(i,j,k) = bc_extdir_vals_ptr[bccomp+n][3]*msku(i,j,0);
103 } else if (bc_ptr[n].hi(0) == REMORABCType::foextrap) {
104 dest_arr(i,j,k) = dest_arr(dom_hi.x,j,k)*msku(i,j,0);
105 } else if (bc_ptr[n].hi(0) == REMORABCType::orlanski_rad) {
106 Real grad_hi = calc_arr(dom_hi.x ,j ,k) - calc_arr(dom_hi.x ,j-1,k);
107 Real grad_hi_jp1 = calc_arr(dom_hi.x ,j+1,k) - calc_arr(dom_hi.x ,j ,k);
108 Real dUdt = calc_arr(dom_hi.x,j,k) - dest_arr(dom_hi.x ,j,k);
109 Real dUdx = dest_arr(dom_hi.x,j,k) - dest_arr(dom_hi.x-1,j,k);
110 if (dUdt * dUdx < zero) dUdt = zero;
111 Real dUde = (dUdt * (grad_hi + grad_hi_jp1) > zero) ? grad_hi : grad_hi_jp1;
112 Real cff = std::max(dUdx*dUdx+dUde*dUde,eps);
113 Real Cx = dUdt * dUdx;
114 dest_arr(i,j,k) = (cff * calc_arr(dom_hi.x+1,j,k) + Cx * dest_arr(dom_hi.x,j,k)) * msku(i,j,0) / (cff + Cx);
115 }
116 });
117 ParallelFor(
118 grow(bx_xhi,IntVect(0,-1,0)) & dest_arr_box, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
119 int iflip = 2*(dom_hi.x + 1) - i;
120 int inner = (bc_ptr[n].hi(0) == REMORABCType::foextrap) ? 1 : 0;
121 if (bc_ptr[n].hi(0) == REMORABCType::ext_dir) {
122 dest_arr(i,j,k) = bc_extdir_vals_ptr[bccomp+n][3]*msku(i,j,0);
123 } else if (bc_ptr[n].hi(0) == REMORABCType::foextrap || bc_ptr[n].hi(0) == REMORABCType::clamped ||
124 bc_ptr[n].hi(0) == REMORABCType::orlanski_rad) {
125 dest_arr(i,j,k) = dest_arr(dom_hi.x+1-inner,j,k)*msku(i,j,0);
126 } else if (bc_ptr[n].hi(0) == REMORABCType::reflect_even) {
127 dest_arr(i,j,k) = dest_arr(iflip,j,k)*msku(i,j,0);
128 } else if (bc_ptr[n].hi(0) == REMORABCType::reflect_odd) {
129 dest_arr(i,j,k) = -dest_arr(iflip,j,k)*msku(i,j,0);
130 }
131 });
132 } // not periodic in x
133
134 if (!is_periodic_in_y or bccomp == BCVars::foextrap_bc(m_ncons))
135 {
136 // Populate ghost cells on lo-y and hi-y domain boundaries
137 Box bx_ylo(bx); bx_ylo.setBig (1,dom_lo.y-1);
138 Box bx_yhi(bx); bx_yhi.setSmall(1,dom_hi.y+1);
139 ParallelFor(
140 grow(bx_ylo,IntVect(-1,0,0)) & dest_arr_box, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
141 int jflip = dom_lo.y - 1 - j;
142 if (bc_ptr[n].lo(1) == REMORABCType::ext_dir) {
143 dest_arr(i,j,k) = bc_extdir_vals_ptr[bccomp+n][1]*msku(i,j,0);
144 } else if (bc_ptr[n].lo(1) == REMORABCType::foextrap || bc_ptr[n].lo(1) == REMORABCType::clamped) {
145 dest_arr(i,j,k) = dest_arr(i,dom_lo.y,k)*msku(i,j,0);
146 } else if (bc_ptr[n].lo(1) == REMORABCType::orlanski_rad) {
147 Real grad_lo = calc_arr(i+1,dom_lo.y ,k) - calc_arr(i ,dom_lo.y ,k);
148 Real grad_lo_im1 = calc_arr(i ,dom_lo.y ,k) - calc_arr(i-1,dom_lo.y ,k);
149 Real dUdt = calc_arr(i,dom_lo.y,k) - dest_arr(i,dom_lo.y ,k);
150 Real dUde = dest_arr(i,dom_lo.y,k) - dest_arr(i,dom_lo.y+1,k);
151 if (dUdt * dUde < zero) dUdt = zero;
152 Real dUdx = (dUdt * (grad_lo_im1 + grad_lo) > zero) ? grad_lo_im1 : grad_lo;
153 Real cff = std::max(dUdx * dUdx + dUde * dUde, eps);
154 Real Ce = dUdt * dUde;
155 dest_arr(i,j,k) = (cff * calc_arr(i,dom_lo.y-1,k) + Ce * dest_arr(i,dom_lo.y,k)) * msku(i,j,0) / (cff + Ce);
156 } else if (bc_ptr[n].lo(1) == REMORABCType::reflect_even) {
157 dest_arr(i,j,k) = dest_arr(i,jflip,k)*msku(i,j,0);
158 } else if (bc_ptr[n].lo(1) == REMORABCType::reflect_odd) {
159 dest_arr(i,j,k) = -dest_arr(i,jflip,k)*msku(i,j,0);
160 }
161 },
162 grow(bx_yhi,IntVect(-1,0,0)) & dest_arr_box, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
163 int jflip = 2*dom_hi.y + 1 - j;
164 if (bc_ptr[n].hi(1) == REMORABCType::ext_dir) {
165 dest_arr(i,j,k) = bc_extdir_vals_ptr[bccomp+n][4]*msku(i,j,0);
166 } else if (bc_ptr[n].hi(1) == REMORABCType::foextrap || bc_ptr[n].hi(1) == REMORABCType::clamped) {
167 dest_arr(i,j,k) = dest_arr(i,dom_hi.y,k)*msku(i,j,0);
168 } else if (bc_ptr[n].hi(1) == REMORABCType::orlanski_rad) {
169 Real grad_hi = calc_arr(i+1,dom_hi.y ,k) - calc_arr(i ,dom_hi.y ,k);
170 Real grad_hi_im1 = calc_arr(i ,dom_hi.y ,k) - calc_arr(i-1,dom_hi.y ,k);
171 Real dUdt = calc_arr(i,dom_hi.y,k) - dest_arr(i,dom_hi.y ,k);
172 Real dUde = dest_arr(i,dom_hi.y,k) - dest_arr(i,dom_hi.y-1,k);
173 if (dUdt * dUde < zero) dUdt = zero;
174 Real dUdx = (dUdt * (grad_hi_im1 + grad_hi) > zero) ? grad_hi_im1 : grad_hi;
175 Real cff = std::max(dUdx*dUdx+dUde*dUde,eps);
176 Real Ce = dUdt * dUde;
177 dest_arr(i,j,k) = (cff * calc_arr(i,dom_hi.y+1,k) + Ce * dest_arr(i,dom_hi.y,k)) * msku(i,j,0) / (cff+Ce);
178 } else if (bc_ptr[n].hi(1) == REMORABCType::reflect_even) {
179 dest_arr(i,j,k) = dest_arr(i,jflip,k)*msku(i,j,0);
180 } else if (bc_ptr[n].hi(1) == REMORABCType::reflect_odd) {
181 dest_arr(i,j,k) = -dest_arr(i,jflip,k)*msku(i,j,0);
182 }
183 }
184 );
185 } // not periodic in y
186
187 {
188 // Populate ghost cells on lo-z and hi-z domain boundaries
189 Box bx_zlo(bx); bx_zlo.setBig (2,dom_lo.z-1);
190 Box bx_zhi(bx); bx_zhi.setSmall(2,dom_hi.z+1);
191 ParallelFor(
192 bx_zlo & dest_arr_box, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
193 int kflip = dom_lo.z - 1 - k;
194 if (bc_ptr[n].lo(2) == REMORABCType::ext_dir) {
195 dest_arr(i,j,k) = bc_extdir_vals_ptr[bccomp+n][2]*msku(i,j,0);
196 } else if (bc_ptr[n].lo(2) == REMORABCType::foextrap) {
197 dest_arr(i,j,k) = dest_arr(i,j,dom_lo.z)*msku(i,j,0);
198 } else if (bc_ptr[n].lo(2) == REMORABCType::reflect_even) {
199 dest_arr(i,j,k) = dest_arr(i,j,kflip)*msku(i,j,0);
200 } else if (bc_ptr[n].lo(2) == REMORABCType::reflect_odd) {
201 dest_arr(i,j,k) = -dest_arr(i,j,kflip)*msku(i,j,0);
202 }
203 },
204 bx_zhi & dest_arr_box, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) {
205 int kflip = 2*dom_hi.z + 1 - k;
206 if (bc_ptr[n].hi(2) == REMORABCType::ext_dir) {
207 dest_arr(i,j,k) = bc_extdir_vals_ptr[bccomp+n][5]*msku(i,j,0);
208 } else if (bc_ptr[n].hi(2) == REMORABCType::foextrap) {
209 dest_arr(i,j,k) = dest_arr(i,j,dom_hi.z)*msku(i,j,0);
210 } else if (bc_ptr[n].hi(2) == REMORABCType::reflect_even) {
211 dest_arr(i,j,k) = dest_arr(i,j,kflip)*msku(i,j,0);
212 } else if (bc_ptr[n].hi(2) == REMORABCType::reflect_odd) {
213 dest_arr(i,j,k) = -dest_arr(i,j,kflip)*msku(i,j,0);
214 }
215 }
216 );
217 } // z
218
219 if ((!is_periodic_in_x or bccomp == BCVars::foextrap_bc(m_ncons)) and
220 (!is_periodic_in_y or bccomp == BCVars::foextrap_bc(m_ncons))) {
221 Box xlo(bx); xlo.setBig (0,dom_lo.x );
222 Box xhi(bx); xhi.setSmall(0,dom_hi.x+1);
223 Box ylo(bx); ylo.setBig (1,dom_lo.y-1);
224 Box yhi(bx); yhi.setSmall(1,dom_hi.y+1);
225 Box xlo_ylo = xlo & ylo;
226 Box xlo_yhi = xlo & yhi;
227 Box xhi_ylo = xhi & ylo;
228 Box xhi_yhi = xhi & yhi;
229
230 const bool clamp_west = m_domain_bcs_type[bccomp].lo(0) == REMORABCType::clamped;
231 const bool clamp_east = m_domain_bcs_type[bccomp].hi(0) == REMORABCType::clamped;
232 const bool clamp_south = m_domain_bcs_type[bccomp].lo(1) == REMORABCType::clamped;
233 const bool clamp_north = m_domain_bcs_type[bccomp].hi(1) == REMORABCType::clamped;
234
235 if (!clamp_west && !clamp_south) {
236 ParallelFor(xlo_ylo & dest_arr_box, [=] AMREX_GPU_DEVICE (int i, int j, int k)
237 {
238 dest_arr(i,j,k) = Real(0.5) * (dest_arr(i,dom_lo.y,k) + dest_arr(dom_lo.x+1,j,k));
239 });
240 }
241 if (!clamp_west && !clamp_north) {
242 ParallelFor(xlo_yhi & dest_arr_box, [=] AMREX_GPU_DEVICE (int i, int j, int k)
243 {
244 dest_arr(i,j,k) = Real(0.5) * (dest_arr(i,dom_hi.y,k) + dest_arr(dom_lo.x+1,j,k));
245 });
246 }
247 if (!clamp_east && !clamp_south) {
248 ParallelFor(xhi_ylo & dest_arr_box, [=] AMREX_GPU_DEVICE (int i, int j, int k)
249 {
250 dest_arr(i,j,k) = Real(0.5) * (dest_arr(i,dom_lo.y,k) + dest_arr(dom_hi.x,j,k));
251 });
252 }
253 if (!clamp_east && !clamp_north) {
254 ParallelFor(xhi_yhi & dest_arr_box, [=] AMREX_GPU_DEVICE (int i, int j, int k)
255 {
256 dest_arr(i,j,k) = Real(0.5) * (dest_arr(i,dom_hi.y,k) + dest_arr(dom_hi.x,j,k));
257 });
258 }
259 }
260
261
262 Gpu::streamSynchronize();
263}
zero
constexpr amrex::Real zero
Definition REMORA_Constants.H:6
setVal
mf_h setVal(geomdata.ProbHi(2))
geomdata
auto geomdata
Definition REMORA_InitAnalyticGridScale_BoundaryLayer.H:2
REMORA_PhysBCFunct.H
REMORAPhysBCFunct::impose_xvel_bcs
void impose_xvel_bcs(const amrex::Array4< amrex::Real > &dest_arr, const amrex::Box &bx, const amrex::Box &domain, const amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > dxInv, const amrex::Array4< const amrex::Real > &msku, const amrex::Array4< const amrex::Real > &calc_arr, amrex::Real time, int bccomp)
apply x-velocity type boundary conditions
Definition REMORA_BoundaryConditions_xvel.cpp:16
REMORAPhysBCFunct::m_domain_bcs_type
amrex::Vector< amrex::BCRec > m_domain_bcs_type
Definition REMORA_PhysBCFunct.H:79
REMORAPhysBCFunct::m_geom
amrex::Geometry m_geom
Definition REMORA_PhysBCFunct.H:78
REMORAPhysBCFunct::m_bc_extdir_vals_d
amrex::Gpu::DeviceVector< amrex::GpuArray< amrex::Real, AMREX_SPACEDIM *2 > > m_bc_extdir_vals_d
Definition REMORA_PhysBCFunct.H:82
BCVars::foextrap_bc
int foextrap_bc(int ncons) noexcept
Definition REMORA_IndexDefines.H:49