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