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 // If we're doing zeta, ubar, or vbar, then calc_arr only has a single component
56 // corresponding to the component to be used in calculating the boundary
57 // value. Since we access icomp + icomp_to_fill_calc, we need icomp_to_fill_calc to be zero.
58 // If it's another variable, either we aren't using calc_arr
59 // or the components correspond to salt, temp, etc so we leave it as is.
60 int icomp_to_fill_calc = (bccomp == BCVars::zeta_bc || bccomp == BCVars::ubar_bc ||
61 bccomp == BCVars::vbar_bc) ? 0 : icomp_to_fill;
62
63 boundary_series[ivar+icomp]->update_interpolated_to_time(time);
64 const auto& bdatxlo = boundary_series[ivar+icomp]->xlo_dat_interp.const_array();
65 const auto& bdatxhi = boundary_series[ivar+icomp]->xhi_dat_interp.const_array();
66 const auto& bdatylo = boundary_series[ivar+icomp]->ylo_dat_interp.const_array();
67 const auto& bdatyhi = boundary_series[ivar+icomp]->yhi_dat_interp.const_array();
68
69 const auto& bx_bdatxlo = boundary_series[ivar+icomp]->xlo_dat_interp.box();
70 const auto& bx_bdatxhi = boundary_series[ivar+icomp]->xhi_dat_interp.box();
71 const auto& bx_bdatylo = boundary_series[ivar+icomp]->ylo_dat_interp.box();
72 const auto& bx_bdatyhi = boundary_series[ivar+icomp]->yhi_dat_interp.box();
73
74 if (domain_bcs_type[bccomp+icomp].lo(0) == REMORABCType::flather ||
75 domain_bcs_type[bccomp+icomp].hi(0) == REMORABCType::flather ||
76 domain_bcs_type[bccomp+icomp].lo(1) == REMORABCType::flather ||
77 domain_bcs_type[bccomp+icomp].hi(1) == REMORABCType::flather) {
78 boundary_series[BdyVars::zeta]->update_interpolated_to_time(time);
79 }
80 const auto& bdatxlo_zeta = domain_bcs_type[bccomp+icomp].lo(0) == REMORABCType::flather ?
81 boundary_series[BdyVars::zeta]->xlo_dat_interp.const_array() : Array4<Real>();
82 const auto& bdatxhi_zeta = domain_bcs_type[bccomp+icomp].hi(0) == REMORABCType::flather ?
83 boundary_series[BdyVars::zeta]->xhi_dat_interp.const_array() : Array4<Real>();
84 const auto& bdatylo_zeta = domain_bcs_type[bccomp+icomp].lo(1) == REMORABCType::flather ?
85 boundary_series[BdyVars::zeta]->ylo_dat_interp.const_array() : Array4<Real>();
86 const auto& bdatyhi_zeta = domain_bcs_type[bccomp+icomp].hi(1) == REMORABCType::flather ?
87 boundary_series[BdyVars::zeta]->yhi_dat_interp.const_array() : Array4<Real>();
88
89 const bool apply_west = (domain_bcs_type[bccomp+icomp].lo(0) == REMORABCType::clamped) ||
90 (domain_bcs_type[bccomp+icomp].lo(0) == REMORABCType::flather) ||
91 (domain_bcs_type[bccomp+icomp].lo(0) == REMORABCType::chapman) ||
92 (domain_bcs_type[bccomp+icomp].lo(0) == REMORABCType::orlanski_rad_nudge);
93 const bool apply_east = (domain_bcs_type[bccomp+icomp].hi(0) == REMORABCType::clamped) ||
94 (domain_bcs_type[bccomp+icomp].hi(0) == REMORABCType::flather) ||
95 (domain_bcs_type[bccomp+icomp].hi(0) == REMORABCType::chapman) ||
96 (domain_bcs_type[bccomp+icomp].hi(0) == REMORABCType::orlanski_rad_nudge);
97 const bool apply_north = (domain_bcs_type[bccomp+icomp].lo(1) == REMORABCType::clamped) ||
98 (domain_bcs_type[bccomp+icomp].lo(1) == REMORABCType::flather) ||
99 (domain_bcs_type[bccomp+icomp].lo(1) == REMORABCType::chapman) ||
100 (domain_bcs_type[bccomp+icomp].lo(1) == REMORABCType::orlanski_rad_nudge);
101 const bool apply_south = (domain_bcs_type[bccomp+icomp].hi(1) == REMORABCType::clamped) ||
102 (domain_bcs_type[bccomp+icomp].hi(1) == REMORABCType::flather) ||
103 (domain_bcs_type[bccomp+icomp].hi(1) == REMORABCType::chapman) ||
104 (domain_bcs_type[bccomp+icomp].hi(1) == REMORABCType::orlanski_rad_nudge);
105
106 const bool cell_centered = (mf_index_type[0] == 0 and mf_index_type[1] == 0);
107
108 const Real obcfac = solverChoice.obcfac;
109
110#ifdef AMREX_USE_OMP
111#pragma omp parallel if (Gpu::notInLaunchRegion())
112#endif
113 // Currently no tiling in order to get the logic right
114 for (MFIter mfi(mf_to_fill,false); mfi.isValid(); ++mfi)
115 {
116 Box mf_box(mf_to_fill[mfi.index()].box());
117
118 // Compute intersections of the FAB to be filled and the bdry data boxes
119 Box xlo = bx_bdatxlo & mf_box;
120 Box xhi = bx_bdatxhi & mf_box;
121 Box ylo = bx_bdatylo & mf_box;
122 Box yhi = bx_bdatyhi & mf_box;
123
124 xlo.setSmall(0,lbound(mf_box).x);
125 xhi.setBig (0,ubound(mf_box).x);
126 ylo.setSmall(1,lbound(mf_box).y);
127 yhi.setBig (1,ubound(mf_box).y);
128
129 Box xlo_ylo = xlo & ylo;
130 Box xlo_yhi = xlo & yhi;
131 Box xhi_ylo = xhi & ylo;
132 Box xhi_yhi = xhi & yhi;
133
134 Box xlo_edge = xlo; xlo_edge.setSmall(0,ubound(xlo).x); xlo_edge.setBig(0,ubound(xlo).x);
135 Box xhi_edge = xhi; xhi_edge.setSmall(0,lbound(xhi).x); xhi_edge.setBig(0,lbound(xhi).x);
136 Box ylo_edge = ylo; ylo_edge.setSmall(1,ubound(ylo).y); ylo_edge.setBig(1,ubound(ylo).y);
137 Box yhi_edge = yhi; yhi_edge.setSmall(1,lbound(yhi).y); yhi_edge.setBig(1,lbound(yhi).y);
138
139 Box xlo_ghost = xlo; xlo_ghost.setBig(0,ubound(xlo).x-1);
140 Box xhi_ghost = xhi; xhi_ghost.setSmall(0,lbound(xhi).x+1);
141 Box ylo_ghost = ylo; ylo_ghost.setBig(1,ubound(ylo).y-1);
142 Box yhi_ghost = yhi; yhi_ghost.setSmall(1,lbound(yhi).y+1);
143
144 const Array4<Real>& dest_arr = mf_to_fill.array(mfi);
145 const Array4<const Real>& mask_arr = mf_mask.array(mfi);
146 const Array4<const Real>& calc_arr = (!null_mf_calc) ? mf_calc.array(mfi) : Array4<amrex::Real>();
147 const Array4<const Real>& h_arr = vec_h[lev]->const_array(mfi);
148 const Array4<const Real>& zeta_arr = vec_zeta[lev]->const_array(mfi);
149 const Array4<const Real>& pm = vec_pm[lev]->const_array(mfi);
150 const Array4<const Real>& pn = vec_pn[lev]->const_array(mfi);
151
152 const Array4<const Real>& msku = vec_msku[lev]->const_array(mfi);
153 const Array4<const Real>& mskv = vec_mskv[lev]->const_array(mfi);
154
155 const Array4<const Real> nudg_coeff_out = vec_nudg_coeff[bdy_var_type][lev]->const_array(mfi);
156
157 //
158 // We are inside a loop over components so we do one at a time here
159 //
160 Vector<BCRec> bcrs(1);
161 amrex::setBC(mf_box, domain, bccomp+icomp, 0, 1, domain_bcs_type, bcrs);
162
163 // xlo: ori = 0
164 // ylo: ori = 1
165 // zlo: ori = 2
166 // xhi: ori = 3
167 // yhi: ori = 4
168 // zhi: ori = 5
169
170 auto bcr = bcrs[0];
171
172 // 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,
173 // but if the grid doesn't include the low x-boundary, the xlo box will be invalid and the execution will be skipped.
174 if (!xlo.isEmpty() && apply_west) {
175 ParallelFor(grow(xlo_edge,IntVect(0,-1,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
176 {
177 Real bry_val = bdatxlo(ubound(xlo).x,j,k,0);
178 if (bcr.lo(0) == REMORABCType::clamped) {
179 dest_arr(i,j,k,icomp+icomp_to_fill) = bry_val * mask_arr(i,j,0);
180 } else if (bcr.lo(0) == REMORABCType::flather) {
181 Real bry_val_zeta = bdatxlo_zeta(ubound(xlo).x-1,j,k,0);
182 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)
183 + h_arr(dom_lo.x,j,0) + zeta_arr(dom_lo.x,j,0,icomp_calc)));
184 Real Cx = std::sqrt(g * cff);
185 dest_arr(i,j,k,icomp+icomp_to_fill) = (bry_val
186 - Cx * (0.5_rt * (zeta_arr(dom_lo.x-1,j,0,icomp_calc) + zeta_arr(dom_lo.x,j,0,icomp_calc))
187 - bry_val_zeta)) * mask_arr(i,j,0);
188 } else if (bcr.lo(0) == REMORABCType::chapman) {
189 Real cff = dt_calc * 0.5_rt * (pm(dom_lo.x,j-mf_index_type[1],0) + pm(dom_lo.x,j,0));
190 Real cff1 = std::sqrt(g * 0.5_rt * (h_arr(dom_lo.x,j-mf_index_type[1],0)
191 + zeta_arr(dom_lo.x,j-mf_index_type[1],0,icomp_calc) + h_arr(dom_lo.x,j,0)
192 + zeta_arr(dom_lo.x,j,0,icomp_calc)));
193 Real Cx = cff * cff1;
194 Real cff2 = 1.0_rt / (1.0_rt + Cx);
195 dest_arr(i,j,k,icomp+icomp_to_fill) = cff2 * (dest_arr(dom_lo.x-1,j,k,icomp_calc)
196 + Cx * dest_arr(dom_lo.x,j,k,icomp+icomp_to_fill)) * mask_arr(i,j,0);
197 } else if (bcr.lo(0) == REMORABCType::orlanski_rad_nudge) {
198 Real grad_lo_im1 = (calc_arr(dom_lo.x+mf_index_type[0]-1,j ,k,icomp+icomp_to_fill_calc) - calc_arr(dom_lo.x-1+mf_index_type[0],j-1,k,icomp+icomp_to_fill_calc));
199 Real grad_lo = (calc_arr(dom_lo.x+mf_index_type[0] ,j ,k,icomp+icomp_to_fill_calc) - calc_arr(dom_lo.x +mf_index_type[0],j-1,k,icomp+icomp_to_fill_calc));
200 Real grad_lo_imjp1 = (calc_arr(dom_lo.x+mf_index_type[0]-1,j+1,k,icomp+icomp_to_fill_calc) - calc_arr(dom_lo.x-1+mf_index_type[0],j ,k,icomp+icomp_to_fill_calc));
201 Real grad_lo_jp1 = (calc_arr(dom_lo.x+mf_index_type[0] ,j+1,k,icomp+icomp_to_fill_calc) - calc_arr(dom_lo.x +mf_index_type[0],j ,k,icomp+icomp_to_fill_calc));
202 if (cell_centered) {
203 grad_lo_im1 *= mskv(i,j,0);
204 grad_lo *= mskv(i,j,0);
205 grad_lo_imjp1 *= mskv(i,j,0);
206 grad_lo_jp1 *= mskv(i,j,0);
207 }
208 Real dTdt = calc_arr(dom_lo.x+mf_index_type[0],j,k,icomp+icomp_to_fill_calc) - dest_arr(dom_lo.x+mf_index_type[0] ,j,k,icomp+icomp_to_fill);
209 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);
210 Real tau;
211 Real nudg_coeff_out_local = (nudg_coeff_out(i-mf_index_type[0],j-mf_index_type[1],k) +
212 nudg_coeff_out(i,j,k)) * 0.5_rt;
213 if (dTdt*dTdx < 0.0_rt) {
214 tau = nudg_coeff_out_local * obcfac * dt_calc;
215 dTdt = 0.0_rt;
216 } else {
217 tau = nudg_coeff_out_local * dt_calc;
218 }
219 Real dTde = (dTdt * (grad_lo+grad_lo_jp1) > 0.0_rt) ? grad_lo : grad_lo_jp1;
220 Real cff = std::max(dTdx*dTdx+dTde*dTde,eps);
221 Real Cx = dTdt * dTdx;
222 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_calc) + Cx * dest_arr(dom_lo.x+mf_index_type[0],j,k,icomp+icomp_to_fill)) / (cff+Cx);
223 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_calc)));
224 }
225 });
226 ParallelFor(grow(xlo_ghost,IntVect(0,-1,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
227 {
228 dest_arr(i,j,k,icomp+icomp_to_fill) = dest_arr(ubound(xlo).x,j,k,icomp+icomp_to_fill);
229 });
230 }
231
232 // See comment on xlo
233 if (!xhi.isEmpty() && apply_east) {
234 ParallelFor(grow(xhi_edge,IntVect(0,-1,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
235 {
236 Real bry_val = bdatxhi(lbound(xhi).x,j,k,0);
237 if (bcr.hi(0) == REMORABCType::clamped) {
238 dest_arr(i,j,k,icomp+icomp_to_fill) = bry_val * mask_arr(i,j,0);
239 } else if (bcr.hi(0) == REMORABCType::flather) {
240 Real bry_val_zeta = bdatxhi_zeta(lbound(xhi).x,j,k,0);
241 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)
242 + h_arr(dom_hi.x,j,0) + zeta_arr(dom_hi.x,j,0,icomp_calc)));
243 Real Cx = std::sqrt(g * cff);
244 dest_arr(i,j,k,icomp+icomp_to_fill) = (bry_val
245 + Cx * (0.5_rt * (zeta_arr(dom_hi.x-1,j,0,icomp_calc) + zeta_arr(dom_hi.x,j,0,icomp_calc))
246 - bry_val_zeta)) * mask_arr(i,j,0);
247 } else if (bcr.hi(0) == REMORABCType::chapman) {
248 Real cff = dt_calc * 0.5_rt * (pm(dom_hi.x,j-mf_index_type[1],0) + pm(dom_hi.x,j,0));
249 Real cff1 = std::sqrt(g * 0.5_rt * (h_arr(dom_hi.x,j-mf_index_type[1],0)
250 + zeta_arr(dom_hi.x,j-mf_index_type[1],0,icomp_calc) + h_arr(dom_hi.x,j,0)
251 + zeta_arr(dom_hi.x,j,0,icomp_calc)));
252 Real Cx = cff * cff1;
253 Real cff2 = 1.0_rt / (1.0_rt + Cx);
254 dest_arr(i,j,k,icomp+icomp_to_fill) = cff2 * (dest_arr(dom_hi.x+1,j,k,icomp_calc)
255 + Cx * dest_arr(dom_hi.x,j,k,icomp+icomp_to_fill)) * mask_arr(i,j,0);
256 } else if (bcr.hi(0) == REMORABCType::orlanski_rad_nudge) {
257 Real grad_hi = (calc_arr(dom_hi.x-mf_index_type[0] ,j ,k,icomp+icomp_to_fill_calc) - calc_arr(dom_hi.x-mf_index_type[0] ,j-1,k,icomp+icomp_to_fill_calc));
258 Real grad_hi_ip1 = (calc_arr(dom_hi.x-mf_index_type[0]+1,j ,k,icomp+icomp_to_fill_calc) - calc_arr(dom_hi.x-mf_index_type[0]+1,j-1,k,icomp+icomp_to_fill_calc));
259 Real grad_hi_jp1 = (calc_arr(dom_hi.x-mf_index_type[0] ,j+1,k,icomp+icomp_to_fill_calc) - calc_arr(dom_hi.x-mf_index_type[0] ,j ,k,icomp+icomp_to_fill_calc));
260 Real grad_hi_ijp1 = (calc_arr(dom_hi.x-mf_index_type[0]+1,j+1,k,icomp+icomp_to_fill_calc) - calc_arr(dom_hi.x-mf_index_type[0]+1,j ,k,icomp+icomp_to_fill_calc));
261 if (cell_centered) {
262 grad_hi *= mskv(i,j,0);
263 grad_hi_ip1 *= mskv(i,j,0);
264 grad_hi_jp1 *= mskv(i,j,0);
265 grad_hi_ijp1 *= mskv(i,j,0);
266 }
267 Real dTdt = calc_arr(dom_hi.x-mf_index_type[0],j,k,icomp+icomp_to_fill_calc) - dest_arr(dom_hi.x-mf_index_type[0] ,j,k,icomp+icomp_to_fill);
268 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);
269 Real tau;
270 Real nudg_coeff_out_local = (nudg_coeff_out(i-mf_index_type[0],j-mf_index_type[1],k) +
271 nudg_coeff_out(i,j,k)) * 0.5_rt;
272 if (dTdt*dTdx < 0.0_rt) {
273 tau = nudg_coeff_out_local * obcfac * dt_calc;
274 dTdt = 0.0_rt;
275 } else {
276 tau = nudg_coeff_out_local * dt_calc;
277 }
278 if (dTdt * dTdx < 0.0_rt) dTdt = 0.0_rt;
279 Real dTde = (dTdt * (grad_hi + grad_hi_jp1) > 0.0_rt) ? grad_hi : grad_hi_jp1;
280 Real cff = std::max(dTdx*dTdx + dTde*dTde,eps);
281 Real Cx = dTdt * dTdx;
282 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_calc) + Cx * dest_arr(dom_hi.x-mf_index_type[0],j,k,icomp+icomp_to_fill)) * mask_arr(i,j,0) / (cff+Cx);
283 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_calc)));
284 }
285 });
286 ParallelFor(grow(xhi_ghost,IntVect(0,-1,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
287 {
288 dest_arr(i,j,k,icomp+icomp_to_fill) = dest_arr(lbound(xhi).x,j,k,icomp+icomp_to_fill);
289 });
290 }
291
292 // See comment on xlo
293 if (!ylo.isEmpty() && apply_south) {
294 ParallelFor(grow(ylo_edge,IntVect(-1,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
295 {
296 Real bry_val = bdatylo(i,ubound(ylo).y,k,0);
297 if (bcr.lo(1) == REMORABCType::clamped) {
298 dest_arr(i,j,k,icomp+icomp_to_fill) = bry_val * mask_arr(i,j,0);
299 } else if (bcr.lo(1) == REMORABCType::flather) {
300 Real bry_val_zeta = bdatylo_zeta(i,ubound(ylo).y-1,k,0);
301 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)
302 + h_arr(i,dom_lo.y,0) + zeta_arr(i,dom_lo.y,0,icomp_calc)));
303 Real Ce = std::sqrt(g * cff);
304 dest_arr(i,j,k,icomp+icomp_to_fill) = (bry_val
305 - Ce * (0.5_rt * (zeta_arr(i,dom_lo.y-1,0,icomp_calc) + zeta_arr(i,dom_lo.y,0,icomp_calc))
306 - bry_val_zeta)) * mask_arr(i,j,0);
307 } else if (bcr.lo(1) == REMORABCType::chapman) {
308 Real cff = dt_calc * 0.5_rt * (pn(i-mf_index_type[0],dom_lo.y,0) + pn(i,dom_lo.y,0));
309 Real cff1 = std::sqrt(g * 0.5_rt * (h_arr(i-mf_index_type[0],dom_lo.y,0) +
310 zeta_arr(i-mf_index_type[0],dom_lo.y,0,icomp_calc) + h_arr(i,dom_lo.y,0)
311 + zeta_arr(i,dom_lo.y,0,icomp_calc)));
312 Real Ce = cff * cff1;
313 Real cff2 = 1.0_rt / (1.0_rt + Ce);
314 dest_arr(i,j,k,icomp+icomp_to_fill) = cff2 * (dest_arr(i,dom_lo.y-1,k,icomp_calc)
315 + Ce * dest_arr(i,dom_lo.y,k,icomp+icomp_to_fill)) * mask_arr(i,j,0);
316 } else if (bcr.lo(1) == REMORABCType::orlanski_rad_nudge) {
317 Real grad_lo = (calc_arr(i ,dom_lo.y+mf_index_type[1], k,icomp+icomp_to_fill_calc) - calc_arr(i-1,dom_lo.y+mf_index_type[1] ,k,icomp+icomp_to_fill_calc));
318 Real grad_lo_jm1 = (calc_arr(i ,dom_lo.y+mf_index_type[1]-1,k,icomp+icomp_to_fill_calc) - calc_arr(i-1,dom_lo.y+mf_index_type[1]-1,k,icomp+icomp_to_fill_calc));
319 Real grad_lo_ip1 = (calc_arr(i+1,dom_lo.y+mf_index_type[1] ,k,icomp+icomp_to_fill_calc) - calc_arr(i ,dom_lo.y+mf_index_type[1] ,k,icomp+icomp_to_fill_calc));
320 Real grad_lo_ipjm1 = (calc_arr(i+1,dom_lo.y+mf_index_type[1]-1,k,icomp+icomp_to_fill_calc) - calc_arr(i ,dom_lo.y+mf_index_type[1]-1,k,icomp+icomp_to_fill_calc));
321 if (cell_centered) {
322 grad_lo *= msku(i,j,0);
323 grad_lo_jm1 *= msku(i,j,0);
324 grad_lo_ip1 *= msku(i,j,0);
325 grad_lo_ipjm1 *= msku(i,j,0);
326 }
327 Real dTdt = calc_arr(i,dom_lo.y+mf_index_type[1],k,icomp+icomp_to_fill_calc) - dest_arr(i,dom_lo.y +mf_index_type[1],k,icomp+icomp_to_fill);
328 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);
329 Real tau;
330 Real nudg_coeff_out_local = (nudg_coeff_out(i-mf_index_type[0],j-mf_index_type[1],k) +
331 nudg_coeff_out(i,j,k)) * 0.5_rt;
332 if (dTdt*dTde < 0.0_rt) {
333 tau = nudg_coeff_out_local * obcfac * dt_calc;
334 dTdt = 0.0_rt;
335 } else {
336 tau = nudg_coeff_out_local * dt_calc;
337 }
338 if (dTdt * dTde < 0.0_rt) dTdt = 0.0_rt;
339 Real dTdx = (dTdt * (grad_lo + grad_lo_ip1) > 0.0_rt) ? grad_lo : grad_lo_ip1;
340 Real cff = std::max(dTdx*dTdx + dTde*dTde, eps);
341 Real Ce = dTdt*dTde;
342 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_calc) + Ce * dest_arr(i,dom_lo.y+mf_index_type[1],k,icomp+icomp_to_fill)) / (cff+Ce);
343 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_calc)));
344 }
345 });
346 ParallelFor(grow(ylo_ghost,IntVect(-1,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
347 {
348 dest_arr(i,j,k,icomp+icomp_to_fill) = dest_arr(i,ubound(ylo).y,k,icomp+icomp_to_fill);
349 });
350 }
351
352 // See comment on xlo
353 if (!yhi.isEmpty() && apply_north) {
354 ParallelFor(grow(yhi_edge,IntVect(-1,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
355 {
356 Real bry_val = bdatyhi(i,lbound(yhi).y,k,0);
357 if (bcr.hi(1) == REMORABCType::clamped) {
358 dest_arr(i,j,k,icomp+icomp_to_fill) = bry_val * mask_arr(i,j,0);
359 } else if (bcr.hi(1) == REMORABCType::flather) {
360 Real bry_val_zeta = bdatyhi_zeta(i,lbound(yhi).y,k,0);
361 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)
362 + h_arr(i,dom_hi.y,0) + zeta_arr(i,dom_hi.y,0,icomp_calc)));
363 Real Ce = std::sqrt(g * cff);
364 dest_arr(i,j,k,icomp+icomp_to_fill) = (bry_val
365 + Ce * (0.5_rt * (zeta_arr(i,dom_hi.y-1,0,icomp_calc) + zeta_arr(i,dom_hi.y,0,icomp_calc))
366 - bry_val_zeta)) * mask_arr(i,j,0);
367 } else if (bcr.hi(1) == REMORABCType::chapman) {
368 Real cff = dt_calc * 0.5_rt * (pn(i-mf_index_type[0],dom_hi.y,0) + pn(i,dom_hi.y,0));
369 Real cff1 = std::sqrt(g * 0.5_rt * (h_arr(i-mf_index_type[0],dom_hi.y,0)
370 + zeta_arr(i-mf_index_type[0],dom_hi.y,0,icomp_calc) +
371 h_arr(i,dom_hi.y,0) + zeta_arr(i,dom_hi.y,0,icomp_calc)));
372 Real Ce = cff * cff1;
373 Real cff2 = 1.0_rt / (1.0_rt + Ce);
374 dest_arr(i,j,k,icomp+icomp_to_fill) = cff2 * (dest_arr(i,dom_hi.y+1,k,icomp_calc)
375 + Ce * dest_arr(i,dom_hi.y,k,icomp+icomp_to_fill)) * mask_arr(i,j,0);
376 } else if (bcr.hi(1) == REMORABCType::orlanski_rad_nudge) {
377 Real grad_hi = calc_arr(i ,dom_hi.y-mf_index_type[1] ,k,icomp+icomp_to_fill_calc) - calc_arr(i-1,dom_hi.y-mf_index_type[1] ,k,icomp+icomp_to_fill_calc);
378 Real grad_hi_jp1 = calc_arr(i ,dom_hi.y-mf_index_type[1]+1,k,icomp+icomp_to_fill_calc) - calc_arr(i-1,dom_hi.y-mf_index_type[1]+1,k,icomp+icomp_to_fill_calc);
379 Real grad_hi_ip1 = calc_arr(i+1,dom_hi.y-mf_index_type[1] ,k,icomp+icomp_to_fill_calc) - calc_arr(i ,dom_hi.y-mf_index_type[1] ,k,icomp+icomp_to_fill_calc);
380 Real grad_hi_ijp1 = calc_arr(i+1,dom_hi.y-mf_index_type[1]+1,k,icomp+icomp_to_fill_calc) - calc_arr(i ,dom_hi.y-mf_index_type[1]+1,k,icomp+icomp_to_fill_calc);
381 if (cell_centered) {
382 grad_hi *= msku(i,j,0);
383 grad_hi_jp1 *= msku(i,j,0);
384 grad_hi_ip1 *= msku(i,j,0);
385 grad_hi_ijp1 *= msku(i,j,0);
386 }
387 Real dTdt = calc_arr(i,dom_hi.y-mf_index_type[1],k,icomp+icomp_to_fill_calc) - dest_arr(i,dom_hi.y -mf_index_type[1],k,icomp+icomp_to_fill);
388 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);
389 Real tau;
390 Real nudg_coeff_out_local = (nudg_coeff_out(i-mf_index_type[0],j-mf_index_type[1],k) +
391 nudg_coeff_out(i,j,k)) * 0.5_rt;
392 if (dTdt*dTde < 0.0_rt) {
393 tau = nudg_coeff_out_local * obcfac * dt_calc;
394 dTdt = 0.0_rt;
395 } else {
396 tau = nudg_coeff_out_local * dt_calc;
397 }
398 if (dTdt * dTde < 0.0_rt) dTdt = 0.0_rt;
399 Real dTdx = (dTdt * (grad_hi + grad_hi_ip1) > 0.0_rt) ? grad_hi : grad_hi_ip1;
400 Real cff = std::max(dTdx*dTdx + dTde*dTde, eps);
401 Real Ce = dTdt*dTde;
402 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_calc) + Ce*dest_arr(i,dom_hi.y-mf_index_type[1],k,icomp+icomp_to_fill)) * mask_arr(i,j,0) / (cff+Ce);
403 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_calc)));
404 }
405 });
406 ParallelFor(grow(yhi_ghost,IntVect(-1,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
407 {
408 dest_arr(i,j,k,icomp+icomp_to_fill) = dest_arr(i,lbound(yhi).y,k,icomp+icomp_to_fill);
409 });
410
411 }
412 // If we've applied boundary conditions to either side, update the corner
413 if (!xlo_ylo.isEmpty() && (apply_west || apply_south)) {
414 ParallelFor(xlo_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k)
415 {
416 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)
417 + dest_arr(dom_lo.x+mf_index_type[0],j,k,icomp+icomp_to_fill));
418 });
419 }
420 if (!xlo_yhi.isEmpty() && (apply_west || apply_north)) {
421 ParallelFor(xlo_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k)
422 {
423 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)
424 + dest_arr(dom_lo.x+mf_index_type[0],j,k,icomp+icomp_to_fill));
425 });
426 }
427 if (!xhi_ylo.isEmpty() && (apply_east || apply_south)) {
428 ParallelFor(xhi_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k)
429 {
430 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)
431 + dest_arr(dom_hi.x-mf_index_type[0],j,k,icomp+icomp_to_fill));
432 });
433 }
434 if (!xhi_yhi.isEmpty() && (apply_east || apply_north)) {
435 ParallelFor(xhi_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k)
436 {
437 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)
438 + dest_arr(dom_hi.x-mf_index_type[0],j,k,icomp+icomp_to_fill));
439 });
440 }
441 } // mfi
442 } // icomp
443}
444#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