REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_BoundaryConditions_netcdf.cpp
Go to the documentation of this file.
1#include "REMORA.H"
2
3using namespace amrex;
4
5#ifdef REMORA_USE_NETCDF
6/*
7 * @param[inout] mf_to_fill data on which to apply BCs
8 * @param[in ] mf_mask land-sea mask
9 * @param[in ] time current time
10 * @param[in ] bccomp index into both domain_bcs_type_bcr and bc_extdir_vals for icomp=0
11 * @param[in ] bdy_var_type which netcdf boundary data to fill from
12 * @param[in ] icomp_to_fill component to update
13 * @param[in ] icomp_calc component to reference from on RHS
14 * @param[in ] mf_calc data for RHS of calculation
15 * @param[in ] dt_calc time step for the calculation
16 */
17
18void
19REMORA::fill_from_bdyfiles (MultiFab& mf_to_fill, const MultiFab& mf_mask, const Real time, const int bccomp, const int bdy_var_type, const int icomp_to_fill, const int icomp_calc, const MultiFab& mf_calc, const Real dt_calc)
20{
21 int lev = 0;
22
23 // Which variable are we filling
24 int ivar = bdy_var_type;
25
26 //
27 // Note that "domain" is mapped onto the type of box the data is in
28 //
29 Box domain = geom[lev].Domain();
30
31 const auto& mf_index_type = mf_to_fill.boxArray().ixType();
32 domain.convert(mf_index_type);
33
34 const auto& dom_lo = amrex::lbound(domain);
35 const auto& dom_hi = amrex::ubound(domain);
36
37 int ncomp;
38
39 // If we are doing the scalars then do salt as well as temp
40 if (ivar == BdyVars::t) {
41 ncomp = 2;
42 } else {
43 ncomp = 1;
44 }
45
46 // This must be true for the logic below to work
47 AMREX_ALWAYS_ASSERT(Temp_comp == 0);
48 AMREX_ALWAYS_ASSERT(Salt_comp == 1);
49
50 const Real eps= 1.0e-20_rt;
51 const bool null_mf_calc = (!mf_calc.ok());
52
53 for (int icomp = 0; icomp < ncomp; icomp++) // This is to do both temp and salt if doing scalars
54 {
55 boundary_series[ivar+icomp]->update_interpolated_to_time(time);
56 const auto& bdatxlo = boundary_series[ivar+icomp]->xlo_dat_interp.const_array();
57 const auto& bdatxhi = boundary_series[ivar+icomp]->xhi_dat_interp.const_array();
58 const auto& bdatylo = boundary_series[ivar+icomp]->ylo_dat_interp.const_array();
59 const auto& bdatyhi = boundary_series[ivar+icomp]->yhi_dat_interp.const_array();
60
61 const auto& bx_bdatxlo = boundary_series[ivar+icomp]->xlo_dat_interp.box();
62 const auto& bx_bdatxhi = boundary_series[ivar+icomp]->xhi_dat_interp.box();
63 const auto& bx_bdatylo = boundary_series[ivar+icomp]->ylo_dat_interp.box();
64 const auto& bx_bdatyhi = boundary_series[ivar+icomp]->yhi_dat_interp.box();
65
66 if (domain_bcs_type[bccomp+icomp].lo(0) == REMORABCType::flather ||
67 domain_bcs_type[bccomp+icomp].hi(0) == REMORABCType::flather ||
68 domain_bcs_type[bccomp+icomp].lo(1) == REMORABCType::flather ||
69 domain_bcs_type[bccomp+icomp].hi(1) == REMORABCType::flather) {
70 boundary_series[BdyVars::zeta]->update_interpolated_to_time(time);
71 }
72 const auto& bdatxlo_zeta = domain_bcs_type[bccomp+icomp].lo(0) == REMORABCType::flather ?
73 boundary_series[BdyVars::zeta]->xlo_dat_interp.const_array() : Array4<Real>();
74 const auto& bdatxhi_zeta = domain_bcs_type[bccomp+icomp].hi(0) == REMORABCType::flather ?
75 boundary_series[BdyVars::zeta]->xhi_dat_interp.const_array() : Array4<Real>();
76 const auto& bdatylo_zeta = domain_bcs_type[bccomp+icomp].lo(1) == REMORABCType::flather ?
77 boundary_series[BdyVars::zeta]->ylo_dat_interp.const_array() : Array4<Real>();
78 const auto& bdatyhi_zeta = domain_bcs_type[bccomp+icomp].hi(1) == REMORABCType::flather ?
79 boundary_series[BdyVars::zeta]->yhi_dat_interp.const_array() : Array4<Real>();
80
81 const bool apply_west = (domain_bcs_type[bccomp+icomp].lo(0) == REMORABCType::clamped) ||
82 (domain_bcs_type[bccomp+icomp].lo(0) == REMORABCType::flather) ||
83 (domain_bcs_type[bccomp+icomp].lo(0) == REMORABCType::chapman) ||
84 (domain_bcs_type[bccomp+icomp].lo(0) == REMORABCType::orlanski_rad_nudge);
85 const bool apply_east = (domain_bcs_type[bccomp+icomp].hi(0) == REMORABCType::clamped) ||
86 (domain_bcs_type[bccomp+icomp].hi(0) == REMORABCType::flather) ||
87 (domain_bcs_type[bccomp+icomp].hi(0) == REMORABCType::chapman) ||
88 (domain_bcs_type[bccomp+icomp].hi(0) == REMORABCType::orlanski_rad_nudge);
89 const bool apply_north = (domain_bcs_type[bccomp+icomp].lo(1) == REMORABCType::clamped) ||
90 (domain_bcs_type[bccomp+icomp].lo(1) == REMORABCType::flather) ||
91 (domain_bcs_type[bccomp+icomp].lo(1) == REMORABCType::chapman) ||
92 (domain_bcs_type[bccomp+icomp].lo(1) == REMORABCType::orlanski_rad_nudge);
93 const bool apply_south = (domain_bcs_type[bccomp+icomp].hi(1) == REMORABCType::clamped) ||
94 (domain_bcs_type[bccomp+icomp].hi(1) == REMORABCType::flather) ||
95 (domain_bcs_type[bccomp+icomp].hi(1) == REMORABCType::chapman) ||
96 (domain_bcs_type[bccomp+icomp].hi(1) == REMORABCType::orlanski_rad_nudge);
97
98 const bool cell_centered = (mf_index_type[0] == 0 and mf_index_type[1] == 0);
99
100 const Real obcfac = solverChoice.obcfac;
101
102#ifdef AMREX_USE_OMP
103#pragma omp parallel if (Gpu::notInLaunchRegion())
104#endif
105 // Currently no tiling in order to get the logic right
106 for (MFIter mfi(mf_to_fill,false); mfi.isValid(); ++mfi)
107 {
108 Box mf_box(mf_to_fill[mfi.index()].box());
109
110 // Compute intersections of the FAB to be filled and the bdry data boxes
111 Box xlo = bx_bdatxlo & mf_box;
112 Box xhi = bx_bdatxhi & mf_box;
113 Box ylo = bx_bdatylo & mf_box;
114 Box yhi = bx_bdatyhi & mf_box;
115
116 xlo.setSmall(0,lbound(mf_box).x);
117 xhi.setBig (0,ubound(mf_box).x);
118 ylo.setSmall(1,lbound(mf_box).y);
119 yhi.setBig (1,ubound(mf_box).y);
120
121 Box xlo_ylo = xlo & ylo;
122 Box xlo_yhi = xlo & yhi;
123 Box xhi_ylo = xhi & ylo;
124 Box xhi_yhi = xhi & yhi;
125
126 Box xlo_edge = xlo; xlo_edge.setSmall(0,ubound(xlo).x); xlo_edge.setBig(0,ubound(xlo).x);
127 Box xhi_edge = xhi; xhi_edge.setSmall(0,lbound(xhi).x); xhi_edge.setBig(0,lbound(xhi).x);
128 Box ylo_edge = ylo; ylo_edge.setSmall(1,ubound(ylo).y); ylo_edge.setBig(1,ubound(ylo).y);
129 Box yhi_edge = yhi; yhi_edge.setSmall(1,lbound(yhi).y); yhi_edge.setBig(1,lbound(yhi).y);
130
131 Box xlo_ghost = xlo; xlo_ghost.setBig(0,ubound(xlo).x-1);
132 Box xhi_ghost = xhi; xhi_ghost.setSmall(0,lbound(xhi).x+1);
133 Box ylo_ghost = ylo; ylo_ghost.setBig(1,ubound(ylo).y-1);
134 Box yhi_ghost = yhi; yhi_ghost.setSmall(1,lbound(yhi).y+1);
135
136 const Array4<Real>& dest_arr = mf_to_fill.array(mfi);
137 const Array4<const Real>& mask_arr = mf_mask.array(mfi);
138 const Array4<const Real>& calc_arr = (!null_mf_calc) ? mf_calc.array(mfi) : Array4<amrex::Real>();
139 const Array4<const Real>& h_arr = vec_h[lev]->const_array(mfi);
140 const Array4<const Real>& zeta_arr = vec_zeta[lev]->const_array(mfi);
141 const Array4<const Real>& pm = vec_pm[lev]->const_array(mfi);
142 const Array4<const Real>& pn = vec_pn[lev]->const_array(mfi);
143
144 const Array4<const Real>& msku = vec_msku[lev]->const_array(mfi);
145 const Array4<const Real>& mskv = vec_mskv[lev]->const_array(mfi);
146
147 const Array4<const Real> nudg_coeff_out = vec_nudg_coeff[bdy_var_type][lev]->const_array(mfi);
148
149 //
150 // We are inside a loop over components so we do one at a time here
151 //
152 Vector<BCRec> bcrs(1);
153 amrex::setBC(mf_box, domain, bccomp+icomp, 0, 1, domain_bcs_type, bcrs);
154
155 // xlo: ori = 0
156 // ylo: ori = 1
157 // zlo: ori = 2
158 // xhi: ori = 3
159 // yhi: ori = 4
160 // zhi: ori = 5
161
162 auto bcr = bcrs[0];
163
164 // Even though we don't loop over xlo itself, this is the right condition to check, since xlo_edge will always be the same for each grid,
165 // but if the grid doesn't include the low x-boundary, the xlo box will be invalid and the execution will be skipped.
166 if (!xlo.isEmpty() && apply_west) {
167 ParallelFor(grow(xlo_edge,IntVect(0,-1,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
168 {
169 Real bry_val = bdatxlo(ubound(xlo).x,j,k,0);
170 if (bcr.lo(0) == REMORABCType::clamped) {
171 dest_arr(i,j,k,icomp+icomp_to_fill) = bry_val * mask_arr(i,j,0);
172 } else if (bcr.lo(0) == REMORABCType::flather) {
173 Real bry_val_zeta = bdatxlo_zeta(ubound(xlo).x-1,j,k,0);
174 Real cff = 1.0_rt / (0.5_rt * (h_arr(dom_lo.x-1,j,0) + zeta_arr(dom_lo.x-1,j,0,icomp_calc)
175 + h_arr(dom_lo.x,j,0) + zeta_arr(dom_lo.x,j,0,icomp_calc)));
176 Real Cx = std::sqrt(g * cff);
177 dest_arr(i,j,k,icomp+icomp_to_fill) = (bry_val
178 - Cx * (0.5_rt * (zeta_arr(dom_lo.x-1,j,0,icomp_calc) + zeta_arr(dom_lo.x,j,0,icomp_calc))
179 - bry_val_zeta)) * mask_arr(i,j,0);
180 } else if (bcr.lo(0) == REMORABCType::chapman) {
181 Real cff = dt_calc * 0.5_rt * (pm(dom_lo.x,j-mf_index_type[1],0) + pm(dom_lo.x,j,0));
182 Real cff1 = std::sqrt(g * 0.5_rt * (h_arr(dom_lo.x,j-mf_index_type[1],0)
183 + zeta_arr(dom_lo.x,j-mf_index_type[1],0,icomp_calc) + h_arr(dom_lo.x,j,0)
184 + zeta_arr(dom_lo.x,j,0,icomp_calc)));
185 Real Cx = cff * cff1;
186 Real cff2 = 1.0_rt / (1.0_rt + Cx);
187 dest_arr(i,j,k,icomp+icomp_to_fill) = cff2 * (dest_arr(dom_lo.x-1,j,k,icomp_calc)
188 + Cx * dest_arr(dom_lo.x,j,k,icomp+icomp_to_fill)) * mask_arr(i,j,0);
189 } else if (bcr.lo(0) == REMORABCType::orlanski_rad_nudge) {
190 Real grad_lo_im1 = (calc_arr(dom_lo.x+mf_index_type[0]-1,j ,k,icomp+icomp_to_fill) - calc_arr(dom_lo.x-1+mf_index_type[0],j-1,k,icomp+icomp_to_fill));
191 Real grad_lo = (calc_arr(dom_lo.x+mf_index_type[0] ,j ,k,icomp+icomp_to_fill) - calc_arr(dom_lo.x +mf_index_type[0],j-1,k,icomp+icomp_to_fill));
192 Real grad_lo_imjp1 = (calc_arr(dom_lo.x+mf_index_type[0]-1,j+1,k,icomp+icomp_to_fill) - calc_arr(dom_lo.x-1+mf_index_type[0],j ,k,icomp+icomp_to_fill));
193 Real grad_lo_jp1 = (calc_arr(dom_lo.x+mf_index_type[0] ,j+1,k,icomp+icomp_to_fill) - calc_arr(dom_lo.x +mf_index_type[0],j ,k,icomp+icomp_to_fill));
194 if (cell_centered) {
195 grad_lo_im1 *= mskv(i,j,0);
196 grad_lo *= mskv(i,j,0);
197 grad_lo_imjp1 *= mskv(i,j,0);
198 grad_lo_jp1 *= mskv(i,j,0);
199 }
200 Real dTdt = calc_arr(dom_lo.x+mf_index_type[0],j,k,icomp+icomp_to_fill) - dest_arr(dom_lo.x+mf_index_type[0] ,j,k,icomp+icomp_to_fill);
201 Real dTdx = dest_arr(dom_lo.x+mf_index_type[0],j,k,icomp+icomp_to_fill) - dest_arr(dom_lo.x+mf_index_type[0]+1,j,k,icomp+icomp_to_fill);
202 Real tau;
203 Real nudg_coeff_out_local = (nudg_coeff_out(i-mf_index_type[0],j-mf_index_type[1],k) +
204 nudg_coeff_out(i,j,k)) * 0.5_rt;
205 if (dTdt*dTdx < 0.0_rt) {
206 tau = nudg_coeff_out_local * obcfac * dt_calc;
207 dTdt = 0.0_rt;
208 } else {
209 tau = nudg_coeff_out_local * dt_calc;
210 }
211 Real dTde = (dTdt * (grad_lo+grad_lo_jp1) > 0.0_rt) ? grad_lo : grad_lo_jp1;
212 Real cff = std::max(dTdx*dTdx+dTde*dTde,eps);
213 Real Cx = dTdt * dTdx;
214 dest_arr(i,j,k,icomp+icomp_to_fill) = (cff * calc_arr(dom_lo.x-1+mf_index_type[0],j,k,icomp+icomp_to_fill) + Cx * dest_arr(dom_lo.x+mf_index_type[0],j,k,icomp+icomp_to_fill)) / (cff+Cx);
215 dest_arr(i,j,k,icomp+icomp_to_fill) = mask_arr(i,j,0) * (dest_arr(dom_lo.x-1+mf_index_type[0],j,k,icomp+icomp_to_fill) + tau * (bry_val - calc_arr(dom_lo.x-1+mf_index_type[0],j,k,icomp+icomp_to_fill)));
216 }
217 });
218 ParallelFor(grow(xlo_ghost,IntVect(0,-1,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
219 {
220 dest_arr(i,j,k,icomp+icomp_to_fill) = dest_arr(ubound(xlo).x,j,k,icomp+icomp_to_fill);
221 });
222 }
223
224 // See comment on xlo
225 if (!xhi.isEmpty() && apply_east) {
226 ParallelFor(grow(xhi_edge,IntVect(0,-1,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
227 {
228 Real bry_val = bdatxhi(lbound(xhi).x,j,k,0);
229 if (bcr.hi(0) == REMORABCType::clamped) {
230 dest_arr(i,j,k,icomp+icomp_to_fill) = bry_val * mask_arr(i,j,0);
231 } else if (bcr.hi(0) == REMORABCType::flather) {
232 Real bry_val_zeta = bdatxhi_zeta(lbound(xhi).x,j,k,0);
233 Real cff = 1.0_rt / (0.5_rt * (h_arr(dom_hi.x-1,j,0) + zeta_arr(dom_hi.x-1,j,0,icomp_calc)
234 + h_arr(dom_hi.x,j,0) + zeta_arr(dom_hi.x,j,0,icomp_calc)));
235 Real Cx = std::sqrt(g * cff);
236 dest_arr(i,j,k,icomp+icomp_to_fill) = (bry_val
237 + Cx * (0.5_rt * (zeta_arr(dom_hi.x-1,j,0,icomp_calc) + zeta_arr(dom_hi.x,j,0,icomp_calc))
238 - bry_val_zeta)) * mask_arr(i,j,0);
239 } else if (bcr.hi(0) == REMORABCType::chapman) {
240 Real cff = dt_calc * 0.5_rt * (pm(dom_hi.x,j-mf_index_type[1],0) + pm(dom_hi.x,j,0));
241 Real cff1 = std::sqrt(g * 0.5_rt * (h_arr(dom_hi.x,j-mf_index_type[1],0)
242 + zeta_arr(dom_hi.x,j-mf_index_type[1],0,icomp_calc) + h_arr(dom_hi.x,j,0)
243 + zeta_arr(dom_hi.x,j,0,icomp_calc)));
244 Real Cx = cff * cff1;
245 Real cff2 = 1.0_rt / (1.0_rt + Cx);
246 dest_arr(i,j,k,icomp+icomp_to_fill) = cff2 * (dest_arr(dom_hi.x+1,j,k,icomp_calc)
247 + Cx * dest_arr(dom_hi.x,j,k,icomp+icomp_to_fill)) * mask_arr(i,j,0);
248 } else if (bcr.hi(0) == REMORABCType::orlanski_rad_nudge) {
249 Real grad_hi = (calc_arr(dom_hi.x-mf_index_type[0] ,j ,k,icomp+icomp_to_fill) - calc_arr(dom_hi.x-mf_index_type[0] ,j-1,k,icomp+icomp_to_fill));
250 Real grad_hi_ip1 = (calc_arr(dom_hi.x-mf_index_type[0]+1,j ,k,icomp+icomp_to_fill) - calc_arr(dom_hi.x-mf_index_type[0]+1,j-1,k,icomp+icomp_to_fill));
251 Real grad_hi_jp1 = (calc_arr(dom_hi.x-mf_index_type[0] ,j+1,k,icomp+icomp_to_fill) - calc_arr(dom_hi.x-mf_index_type[0] ,j ,k,icomp+icomp_to_fill));
252 Real grad_hi_ijp1 = (calc_arr(dom_hi.x-mf_index_type[0]+1,j+1,k,icomp+icomp_to_fill) - calc_arr(dom_hi.x-mf_index_type[0]+1,j ,k,icomp+icomp_to_fill));
253 if (cell_centered) {
254 grad_hi *= mskv(i,j,0);
255 grad_hi_ip1 *= mskv(i,j,0);
256 grad_hi_jp1 *= mskv(i,j,0);
257 grad_hi_ijp1 *= mskv(i,j,0);
258 }
259 Real dTdt = calc_arr(dom_hi.x-mf_index_type[0],j,k,icomp+icomp_to_fill) - dest_arr(dom_hi.x-mf_index_type[0] ,j,k,icomp+icomp_to_fill);
260 Real dTdx = dest_arr(dom_hi.x-mf_index_type[0],j,k,icomp+icomp_to_fill) - dest_arr(dom_hi.x-mf_index_type[0]-1,j,k,icomp+icomp_to_fill);
261 Real tau;
262 Real nudg_coeff_out_local = (nudg_coeff_out(i-mf_index_type[0],j-mf_index_type[1],k) +
263 nudg_coeff_out(i,j,k)) * 0.5_rt;
264 if (dTdt*dTdx < 0.0_rt) {
265 tau = nudg_coeff_out_local * obcfac * dt_calc;
266 dTdt = 0.0_rt;
267 } else {
268 tau = nudg_coeff_out_local * dt_calc;
269 }
270 if (dTdt * dTdx < 0.0_rt) dTdt = 0.0_rt;
271 Real dTde = (dTdt * (grad_hi + grad_hi_jp1) > 0.0_rt) ? grad_hi : grad_hi_jp1;
272 Real cff = std::max(dTdx*dTdx + dTde*dTde,eps);
273 Real Cx = dTdt * dTdx;
274 dest_arr(i,j,k,icomp+icomp_to_fill) = (cff * calc_arr(dom_hi.x+1-mf_index_type[0],j,k,icomp+icomp_to_fill) + Cx * dest_arr(dom_hi.x-mf_index_type[0],j,k,icomp+icomp_to_fill)) * mask_arr(i,j,0) / (cff+Cx);
275 dest_arr(i,j,k,icomp+icomp_to_fill) = mask_arr(i,j,0) * (dest_arr(dom_hi.x+1-mf_index_type[0],j,k,icomp+icomp_to_fill) + tau * (bry_val - calc_arr(dom_hi.x+1-mf_index_type[0],j,k,icomp+icomp_to_fill)));
276 }
277 });
278 ParallelFor(grow(xhi_ghost,IntVect(0,-1,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
279 {
280 dest_arr(i,j,k,icomp+icomp_to_fill) = dest_arr(lbound(xhi).x,j,k,icomp+icomp_to_fill);
281 });
282 }
283
284 // See comment on xlo
285 if (!ylo.isEmpty() && apply_south) {
286 ParallelFor(grow(ylo_edge,IntVect(-1,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
287 {
288 Real bry_val = bdatylo(i,ubound(ylo).y,k,0);
289 if (bcr.lo(1) == REMORABCType::clamped) {
290 dest_arr(i,j,k,icomp+icomp_to_fill) = bry_val * mask_arr(i,j,0);
291 } else if (bcr.lo(1) == REMORABCType::flather) {
292 Real bry_val_zeta = bdatylo_zeta(i,ubound(ylo).y-1,k,0);
293 Real cff = 1.0_rt / (0.5_rt * (h_arr(i,dom_lo.y-1,0) + zeta_arr(i,dom_lo.y-1,0,icomp_calc)
294 + h_arr(i,dom_lo.y,0) + zeta_arr(i,dom_lo.y,0,icomp_calc)));
295 Real Ce = std::sqrt(g * cff);
296 dest_arr(i,j,k,icomp+icomp_to_fill) = (bry_val
297 - Ce * (0.5_rt * (zeta_arr(i,dom_lo.y-1,0,icomp_calc) + zeta_arr(i,dom_lo.y,0,icomp_calc))
298 - bry_val_zeta)) * mask_arr(i,j,0);
299 } else if (bcr.lo(1) == REMORABCType::chapman) {
300 Real cff = dt_calc * 0.5_rt * (pn(i-mf_index_type[0],dom_lo.y,0) + pn(i,dom_lo.y,0));
301 Real cff1 = std::sqrt(g * 0.5_rt * (h_arr(i-mf_index_type[0],dom_lo.y,0) +
302 zeta_arr(i-mf_index_type[0],dom_lo.y,0,icomp_calc) + h_arr(i,dom_lo.y,0)
303 + zeta_arr(i,dom_lo.y,0,icomp_calc)));
304 Real Ce = cff * cff1;
305 Real cff2 = 1.0_rt / (1.0_rt + Ce);
306 dest_arr(i,j,k,icomp+icomp_to_fill) = cff2 * (dest_arr(i,dom_lo.y-1,k,icomp_calc)
307 + Ce * dest_arr(i,dom_lo.y,k,icomp+icomp_to_fill)) * mask_arr(i,j,0);
308 } else if (bcr.lo(1) == REMORABCType::orlanski_rad_nudge) {
309 Real grad_lo = (calc_arr(i ,dom_lo.y+mf_index_type[1], k,icomp+icomp_to_fill) - calc_arr(i-1,dom_lo.y+mf_index_type[1] ,k,icomp+icomp_to_fill));
310 Real grad_lo_jm1 = (calc_arr(i ,dom_lo.y+mf_index_type[1]-1,k,icomp+icomp_to_fill) - calc_arr(i-1,dom_lo.y+mf_index_type[1]-1,k,icomp+icomp_to_fill));
311 Real grad_lo_ip1 = (calc_arr(i+1,dom_lo.y+mf_index_type[1] ,k,icomp+icomp_to_fill) - calc_arr(i ,dom_lo.y+mf_index_type[1] ,k,icomp+icomp_to_fill));
312 Real grad_lo_ipjm1 = (calc_arr(i+1,dom_lo.y+mf_index_type[1]-1,k,icomp+icomp_to_fill) - calc_arr(i ,dom_lo.y+mf_index_type[1]-1,k,icomp+icomp_to_fill));
313 if (cell_centered) {
314 grad_lo *= msku(i,j,0);
315 grad_lo_jm1 *= msku(i,j,0);
316 grad_lo_ip1 *= msku(i,j,0);
317 grad_lo_ipjm1 *= msku(i,j,0);
318 }
319 Real dTdt = calc_arr(i,dom_lo.y+mf_index_type[1],k,icomp+icomp_to_fill) - dest_arr(i,dom_lo.y +mf_index_type[1],k,icomp+icomp_to_fill);
320 Real dTde = dest_arr(i,dom_lo.y+mf_index_type[1],k,icomp+icomp_to_fill) - dest_arr(i,dom_lo.y+1+mf_index_type[1],k,icomp+icomp_to_fill);
321 Real tau;
322 Real nudg_coeff_out_local = (nudg_coeff_out(i-mf_index_type[0],j-mf_index_type[1],k) +
323 nudg_coeff_out(i,j,k)) * 0.5_rt;
324 if (dTdt*dTde < 0.0_rt) {
325 tau = nudg_coeff_out_local * obcfac * dt_calc;
326 dTdt = 0.0_rt;
327 } else {
328 tau = nudg_coeff_out_local * dt_calc;
329 }
330 if (dTdt * dTde < 0.0_rt) dTdt = 0.0_rt;
331 Real dTdx = (dTdt * (grad_lo + grad_lo_ip1) > 0.0_rt) ? grad_lo : grad_lo_ip1;
332 Real cff = std::max(dTdx*dTdx + dTde*dTde, eps);
333 Real Ce = dTdt*dTde;
334 dest_arr(i,j,k,icomp+icomp_to_fill) = (cff * calc_arr(i,dom_lo.y-1+mf_index_type[1],k,icomp+icomp_to_fill) + Ce * dest_arr(i,dom_lo.y+mf_index_type[1],k,icomp+icomp_to_fill)) / (cff+Ce);
335 dest_arr(i,j,k,icomp+icomp_to_fill) = mask_arr(i,j,0) * (dest_arr(i,dom_lo.y-1+mf_index_type[1],k,icomp+icomp_to_fill) + tau * (bry_val - calc_arr(i,dom_lo.y-1+mf_index_type[1],k,icomp+icomp_to_fill)));
336 }
337 });
338 ParallelFor(grow(ylo_ghost,IntVect(-1,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
339 {
340 dest_arr(i,j,k,icomp+icomp_to_fill) = dest_arr(i,ubound(ylo).y,k,icomp+icomp_to_fill);
341 });
342 }
343
344 // See comment on xlo
345 if (!yhi.isEmpty() && apply_north) {
346 ParallelFor(grow(yhi_edge,IntVect(-1,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
347 {
348 Real bry_val = bdatyhi(i,lbound(yhi).y,k,0);
349 if (bcr.hi(1) == REMORABCType::clamped) {
350 dest_arr(i,j,k,icomp+icomp_to_fill) = bry_val * mask_arr(i,j,0);
351 } else if (bcr.hi(1) == REMORABCType::flather) {
352 Real bry_val_zeta = bdatyhi_zeta(i,lbound(yhi).y,k,0);
353 Real cff = 1.0_rt / (0.5_rt * (h_arr(i,dom_hi.y-1,0) + zeta_arr(i,dom_hi.y-1,0,icomp_calc)
354 + h_arr(i,dom_hi.y,0) + zeta_arr(i,dom_hi.y,0,icomp_calc)));
355 Real Ce = std::sqrt(g * cff);
356 dest_arr(i,j,k,icomp+icomp_to_fill) = (bry_val
357 + Ce * (0.5_rt * (zeta_arr(i,dom_hi.y-1,0,icomp_calc) + zeta_arr(i,dom_hi.y,0,icomp_calc))
358 - bry_val_zeta)) * mask_arr(i,j,0);
359 } else if (bcr.hi(1) == REMORABCType::chapman) {
360 Real cff = dt_calc * 0.5_rt * (pn(i-mf_index_type[0],dom_hi.y,0) + pn(i,dom_hi.y,0));
361 Real cff1 = std::sqrt(g * 0.5_rt * (h_arr(i-mf_index_type[0],dom_hi.y,0)
362 + zeta_arr(i-mf_index_type[0],dom_hi.y,0,icomp_calc) +
363 h_arr(i,dom_hi.y,0) + zeta_arr(i,dom_hi.y,0,icomp_calc)));
364 Real Ce = cff * cff1;
365 Real cff2 = 1.0_rt / (1.0_rt + Ce);
366 dest_arr(i,j,k,icomp+icomp_to_fill) = cff2 * (dest_arr(i,dom_hi.y+1,k,icomp_calc)
367 + Ce * dest_arr(i,dom_hi.y,k,icomp+icomp_to_fill)) * mask_arr(i,j,0);
368 } else if (bcr.hi(1) == REMORABCType::orlanski_rad_nudge) {
369 Real grad_hi = calc_arr(i ,dom_hi.y-mf_index_type[1] ,k,icomp+icomp_to_fill) - calc_arr(i-1,dom_hi.y-mf_index_type[1] ,k,icomp+icomp_to_fill);
370 Real grad_hi_jp1 = calc_arr(i ,dom_hi.y-mf_index_type[1]+1,k,icomp+icomp_to_fill) - calc_arr(i-1,dom_hi.y-mf_index_type[1]+1,k,icomp+icomp_to_fill);
371 Real grad_hi_ip1 = calc_arr(i+1,dom_hi.y-mf_index_type[1] ,k,icomp+icomp_to_fill) - calc_arr(i ,dom_hi.y-mf_index_type[1] ,k,icomp+icomp_to_fill);
372 Real grad_hi_ijp1 = calc_arr(i+1,dom_hi.y-mf_index_type[1]+1,k,icomp+icomp_to_fill) - calc_arr(i ,dom_hi.y-mf_index_type[1]+1,k,icomp+icomp_to_fill);
373 if (cell_centered) {
374 grad_hi *= msku(i,j,0);
375 grad_hi_jp1 *= msku(i,j,0);
376 grad_hi_ip1 *= msku(i,j,0);
377 grad_hi_ijp1 *= msku(i,j,0);
378 }
379 Real dTdt = calc_arr(i,dom_hi.y-mf_index_type[1],k,icomp+icomp_to_fill) - dest_arr(i,dom_hi.y -mf_index_type[1],k,icomp+icomp_to_fill);
380 Real dTde = dest_arr(i,dom_hi.y-mf_index_type[1],k,icomp+icomp_to_fill) - dest_arr(i,dom_hi.y-1-mf_index_type[1],k,icomp+icomp_to_fill);
381 Real tau;
382 Real nudg_coeff_out_local = (nudg_coeff_out(i-mf_index_type[0],j-mf_index_type[1],k) +
383 nudg_coeff_out(i,j,k)) * 0.5_rt;
384 if (dTdt*dTde < 0.0_rt) {
385 tau = nudg_coeff_out_local * obcfac * dt_calc;
386 dTdt = 0.0_rt;
387 } else {
388 tau = nudg_coeff_out_local * dt_calc;
389 }
390 if (dTdt * dTde < 0.0_rt) dTdt = 0.0_rt;
391 Real dTdx = (dTdt * (grad_hi + grad_hi_ip1) > 0.0_rt) ? grad_hi : grad_hi_ip1;
392 Real cff = std::max(dTdx*dTdx + dTde*dTde, eps);
393 Real Ce = dTdt*dTde;
394 dest_arr(i,j,k,icomp+icomp_to_fill) = (cff*calc_arr(i,dom_hi.y+1-mf_index_type[1],k,icomp+icomp_to_fill) + Ce*dest_arr(i,dom_hi.y-mf_index_type[1],k,icomp+icomp_to_fill)) * mask_arr(i,j,0) / (cff+Ce);
395 dest_arr(i,j,k,icomp+icomp_to_fill) = mask_arr(i,j,0) * (dest_arr(i,dom_hi.y+1-mf_index_type[1],k,icomp+icomp_to_fill) + tau * (bry_val - calc_arr(i,dom_hi.y+1-mf_index_type[1],k,icomp+icomp_to_fill)));
396 }
397 });
398 ParallelFor(grow(yhi_ghost,IntVect(-1,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
399 {
400 dest_arr(i,j,k,icomp+icomp_to_fill) = dest_arr(i,lbound(yhi).y,k,icomp+icomp_to_fill);
401 });
402
403 }
404 // If we've applied boundary conditions to either side, update the corner
405 if (!xlo_ylo.isEmpty() && (apply_west || apply_south)) {
406 ParallelFor(xlo_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k)
407 {
408 dest_arr(i,j,k,icomp+icomp_to_fill) = 0.5 * (dest_arr(i,dom_lo.y+mf_index_type[1],k,icomp+icomp_to_fill)
409 + dest_arr(dom_lo.x+mf_index_type[0],j,k,icomp+icomp_to_fill));
410 });
411 }
412 if (!xlo_yhi.isEmpty() && (apply_west || apply_north)) {
413 ParallelFor(xlo_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k)
414 {
415 dest_arr(i,j,k,icomp+icomp_to_fill) = 0.5 * (dest_arr(i,dom_hi.y-mf_index_type[1],k,icomp+icomp_to_fill)
416 + dest_arr(dom_lo.x+mf_index_type[0],j,k,icomp+icomp_to_fill));
417 });
418 }
419 if (!xhi_ylo.isEmpty() && (apply_east || apply_south)) {
420 ParallelFor(xhi_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k)
421 {
422 dest_arr(i,j,k,icomp+icomp_to_fill) = 0.5 * (dest_arr(i,dom_lo.y+mf_index_type[1],k,icomp+icomp_to_fill)
423 + dest_arr(dom_hi.x-mf_index_type[0],j,k,icomp+icomp_to_fill));
424 });
425 }
426 if (!xhi_yhi.isEmpty() && (apply_east || apply_north)) {
427 ParallelFor(xhi_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k)
428 {
429 dest_arr(i,j,k,icomp+icomp_to_fill) = 0.5 * (dest_arr(i,dom_hi.y-mf_index_type[1],k,icomp+icomp_to_fill)
430 + dest_arr(dom_hi.x-mf_index_type[0],j,k,icomp+icomp_to_fill));
431 });
432 }
433 } // mfi
434 } // icomp
435}
436#endif
g
constexpr amrex::Real g
Definition REMORA_Constants.H:16
Coord::y
@ y
Coord::x
@ x
Temp_comp
#define Temp_comp
Definition REMORA_IndexDefines.H:8
Salt_comp
#define Salt_comp
Definition REMORA_IndexDefines.H:9
REMORA::domain_bcs_type
amrex::Vector< amrex::BCRec > domain_bcs_type
vector (over BCVars) of BCRecs
Definition REMORA.H:1179
REMORA::vec_h
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_h
Bathymetry data (2D, positive valued, h in ROMS)
Definition REMORA.H:241
REMORA::vec_pm
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pm
horizontal scaling factor: 1 / dx (2D)
Definition REMORA.H:372
REMORA::vec_msku
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_msku
land/sea mask at x-faces (2D)
Definition REMORA.H:365
REMORA::boundary_series
amrex::Vector< NCTimeSeriesBoundary * > boundary_series
Vector over BdyVars of boundary series data containers.
Definition REMORA.H:1091
REMORA::vec_mskv
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskv
land/sea mask at y-faces (2D)
Definition REMORA.H:367
REMORA::fill_from_bdyfiles
void fill_from_bdyfiles(amrex::MultiFab &mf_to_fill, const amrex::MultiFab &mf_mask, const amrex::Real time, const int bccomp, const int bdy_var_type, const int icomp_to_fill, const int icomp_calc=0, const amrex::MultiFab &mf_calc=amrex::MultiFab(), const amrex::Real=amrex::Real(0.0))
Fill boundary data from netcdf file.
Definition REMORA_BoundaryConditions_netcdf.cpp:19
REMORA::solverChoice
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1280
REMORA::vec_zeta
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_zeta
free surface height (2D)
Definition REMORA.H:360
REMORA::vec_nudg_coeff
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > vec_nudg_coeff
Climatology nudging coefficients.
Definition REMORA.H:431
REMORA::vec_pn
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pn
horizontal scaling factor: 1 / dy (2D)
Definition REMORA.H:374
SolverChoice::obcfac
amrex::Real obcfac
Definition REMORA_DataStruct.H:720