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