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(model_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(model_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 const Real l_g = solverChoice.g;
116
117#ifdef AMREX_USE_OMP
118#pragma omp parallel if (Gpu::notInLaunchRegion())
119#endif
120 // Currently no tiling in order to get the logic right
121 for (MFIter mfi(mf_to_fill,false); mfi.isValid(); ++mfi)
122 {
123 Box mf_box(mf_to_fill[mfi.index()].box());
124
125 // Compute intersections of the FAB to be filled and the bdry data boxes
126 Box xlo = bx_bdatxlo & mf_box;
127 Box xhi = bx_bdatxhi & mf_box;
128 Box ylo = bx_bdatylo & mf_box;
129 Box yhi = bx_bdatyhi & mf_box;
130
131 xlo.setSmall(0,lbound(mf_box).x);
132 xhi.setBig (0,ubound(mf_box).x);
133 ylo.setSmall(1,lbound(mf_box).y);
134 yhi.setBig (1,ubound(mf_box).y);
135
136 Box xlo_ylo = xlo & ylo;
137 Box xlo_yhi = xlo & yhi;
138 Box xhi_ylo = xhi & ylo;
139 Box xhi_yhi = xhi & yhi;
140
141 Box xlo_edge = xlo; xlo_edge.setSmall(0,ubound(xlo).x); xlo_edge.setBig(0,ubound(xlo).x);
142 Box xhi_edge = xhi; xhi_edge.setSmall(0,lbound(xhi).x); xhi_edge.setBig(0,lbound(xhi).x);
143 Box ylo_edge = ylo; ylo_edge.setSmall(1,ubound(ylo).y); ylo_edge.setBig(1,ubound(ylo).y);
144 Box yhi_edge = yhi; yhi_edge.setSmall(1,lbound(yhi).y); yhi_edge.setBig(1,lbound(yhi).y);
145
146 Box xlo_ghost = xlo; xlo_ghost.setBig(0,ubound(xlo).x-1);
147 Box xhi_ghost = xhi; xhi_ghost.setSmall(0,lbound(xhi).x+1);
148 Box ylo_ghost = ylo; ylo_ghost.setBig(1,ubound(ylo).y-1);
149 Box yhi_ghost = yhi; yhi_ghost.setSmall(1,lbound(yhi).y+1);
150
151 const Array4<Real>& dest_arr = mf_to_fill.array(mfi);
152 const Array4<const Real>& mask_arr = mf_mask.array(mfi);
153 const Array4<const Real>& calc_arr = (!null_mf_calc) ? mf_calc.array(mfi) : Array4<amrex::Real>();
154 const Array4<const Real>& h_arr = vec_h[lev]->const_array(mfi);
155 const Array4<const Real>& zeta_arr = vec_zeta[lev]->const_array(mfi);
156 const Array4<const Real>& pm = vec_pm[lev]->const_array(mfi);
157 const Array4<const Real>& pn = vec_pn[lev]->const_array(mfi);
158
159 const Array4<const Real>& msku = vec_msku[lev]->const_array(mfi);
160 const Array4<const Real>& mskv = vec_mskv[lev]->const_array(mfi);
161
162 // Same ivar+icomp mapping as boundary_series above: the BdyVars tracer slots
163 // are contiguous from BdyVars::t in cons component order, so this is the
164 // nudging coefficient for cons component icomp rather than temperature's.
165 const Array4<const Real> nudg_coeff_out = vec_nudg_coeff[ivar+icomp][lev]->const_array(mfi);
166
167 //
168 // We are inside a loop over components so we do one at a time here
169 //
170 Vector<BCRec> bcrs(1);
171 amrex::setBC(mf_box, domain, bccomp+icomp, 0, 1, domain_bcs_type, bcrs);
172
173 // xlo: ori = 0
174 // ylo: ori = 1
175 // zlo: ori = 2
176 // xhi: ori = 3
177 // yhi: ori = 4
178 // zhi: ori = 5
179
180 auto bcr = bcrs[0];
181
182 // 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,
183 // but if the grid doesn't include the low x-boundary, the xlo box will be invalid and the execution will be skipped.
184 if (!xlo.isEmpty() && apply_west) {
185 ParallelFor(grow(xlo_edge,IntVect(0,-1,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
186 {
187 Real bry_val = bdatxlo(ubound(xlo).x,j,k,0);
188 if (bcr.lo(0) == REMORABCType::clamped) {
189 dest_arr(i,j,k,icomp+icomp_to_fill) = bry_val * mask_arr(i,j,0);
190 } else if (bcr.lo(0) == REMORABCType::flather) {
191 Real bry_val_zeta = bdatxlo_zeta(ubound(xlo).x-1,j,k,0);
192 Real cff = one / (Real(0.5) * (h_arr(dom_lo.x-1,j,0) + zeta_arr(dom_lo.x-1,j,0,icomp_calc)
193 + h_arr(dom_lo.x,j,0) + zeta_arr(dom_lo.x,j,0,icomp_calc)));
194 Real Cx = std::sqrt(l_g * cff);
195 dest_arr(i,j,k,icomp+icomp_to_fill) = (bry_val
196 - Cx * (Real(0.5) * (zeta_arr(dom_lo.x-1,j,0,icomp_calc) + zeta_arr(dom_lo.x,j,0,icomp_calc))
197 - bry_val_zeta)) * mask_arr(i,j,0);
198 } else if (bcr.lo(0) == REMORABCType::chapman) {
199 Real cff = dt_calc * Real(0.5) * (pm(dom_lo.x,j-mf_index_type[1],0) + pm(dom_lo.x,j,0));
200 Real cff1 = std::sqrt(l_g * Real(0.5) * (h_arr(dom_lo.x,j-mf_index_type[1],0)
201 + zeta_arr(dom_lo.x,j-mf_index_type[1],0,icomp_calc) + h_arr(dom_lo.x,j,0)
202 + zeta_arr(dom_lo.x,j,0,icomp_calc)));
203 Real Cx = cff * cff1;
204 Real cff2 = one / (one + Cx);
205 dest_arr(i,j,k,icomp+icomp_to_fill) = cff2 * (dest_arr(dom_lo.x-1,j,k,icomp_calc)
206 + Cx * dest_arr(dom_lo.x,j,k,icomp+icomp_to_fill)) * mask_arr(i,j,0);
207 } else if (bcr.lo(0) == REMORABCType::orlanski_rad_nudge) {
208 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));
209 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));
210 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));
211 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));
212 if (cell_centered) {
213 grad_lo_im1 *= mskv(dom_lo.x+mf_index_type[0]-1,j ,0);
214 grad_lo *= mskv(dom_lo.x+mf_index_type[0] ,j ,0);
215 grad_lo_imjp1 *= mskv(dom_lo.x+mf_index_type[0]-1,j+1,0);
216 grad_lo_jp1 *= mskv(dom_lo.x+mf_index_type[0] ,j+1,0);
217 }
218 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);
219 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);
220 Real tau;
221 Real nudg_coeff_out_local = (nudg_coeff_out(i-mf_index_type[0],j-mf_index_type[1],k) +
222 nudg_coeff_out(i,j,k)) * Real(0.5);
223 if (dTdt*dTdx < zero) {
224 tau = nudg_coeff_out_local * obcfac * dt_calc;
225 dTdt = zero;
226 } else {
227 tau = nudg_coeff_out_local * dt_calc;
228 }
229 Real dTde = (dTdt * (grad_lo+grad_lo_jp1) > zero) ? grad_lo : grad_lo_jp1;
230 Real cff = std::max(dTdx*dTdx+dTde*dTde,eps);
231 Real Cx = dTdt * dTdx;
232 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);
233 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)));
234 }
235 });
236 ParallelFor(grow(xlo_ghost,IntVect(0,-1,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
237 {
238 dest_arr(i,j,k,icomp+icomp_to_fill) = dest_arr(ubound(xlo).x,j,k,icomp+icomp_to_fill);
239 });
240 }
241
242 // See comment on xlo
243 if (!xhi.isEmpty() && apply_east) {
244 ParallelFor(grow(xhi_edge,IntVect(0,-1,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
245 {
246 Real bry_val = bdatxhi(lbound(xhi).x,j,k,0);
247 if (bcr.hi(0) == REMORABCType::clamped) {
248 dest_arr(i,j,k,icomp+icomp_to_fill) = bry_val * mask_arr(i,j,0);
249 } else if (bcr.hi(0) == REMORABCType::flather) {
250 Real bry_val_zeta = bdatxhi_zeta(lbound(xhi).x,j,k,0);
251 Real cff = one / (Real(0.5) * (h_arr(dom_hi.x-1,j,0) + zeta_arr(dom_hi.x-1,j,0,icomp_calc)
252 + h_arr(dom_hi.x,j,0) + zeta_arr(dom_hi.x,j,0,icomp_calc)));
253 Real Cx = std::sqrt(l_g * cff);
254 dest_arr(i,j,k,icomp+icomp_to_fill) = (bry_val
255 + Cx * (Real(0.5) * (zeta_arr(dom_hi.x-1,j,0,icomp_calc) + zeta_arr(dom_hi.x,j,0,icomp_calc))
256 - bry_val_zeta)) * mask_arr(i,j,0);
257 } else if (bcr.hi(0) == REMORABCType::chapman) {
258 Real cff = dt_calc * Real(0.5) * (pm(dom_hi.x,j-mf_index_type[1],0) + pm(dom_hi.x,j,0));
259 Real cff1 = std::sqrt(l_g * Real(0.5) * (h_arr(dom_hi.x,j-mf_index_type[1],0)
260 + zeta_arr(dom_hi.x,j-mf_index_type[1],0,icomp_calc) + h_arr(dom_hi.x,j,0)
261 + zeta_arr(dom_hi.x,j,0,icomp_calc)));
262 Real Cx = cff * cff1;
263 Real cff2 = one / (one + Cx);
264 dest_arr(i,j,k,icomp+icomp_to_fill) = cff2 * (dest_arr(dom_hi.x+1,j,k,icomp_calc)
265 + Cx * dest_arr(dom_hi.x,j,k,icomp+icomp_to_fill)) * mask_arr(i,j,0);
266 } else if (bcr.hi(0) == REMORABCType::orlanski_rad_nudge) {
267 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));
268 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));
269 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));
270 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));
271 if (cell_centered) {
272 grad_hi *= mskv(dom_hi.x-mf_index_type[0] ,j ,0);
273 grad_hi_ip1 *= mskv(dom_hi.x-mf_index_type[0]+1,j ,0);
274 grad_hi_jp1 *= mskv(dom_hi.x-mf_index_type[0] ,j+1,0);
275 grad_hi_ijp1 *= mskv(dom_hi.x-mf_index_type[0]+1,j+1,0);
276 }
277 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);
278 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);
279 Real tau;
280 Real nudg_coeff_out_local = (nudg_coeff_out(i-mf_index_type[0],j-mf_index_type[1],k) +
281 nudg_coeff_out(i,j,k)) * Real(0.5);
282 if (dTdt*dTdx < zero) {
283 tau = nudg_coeff_out_local * obcfac * dt_calc;
284 dTdt = zero;
285 } else {
286 tau = nudg_coeff_out_local * dt_calc;
287 }
288 if (dTdt * dTdx < zero) dTdt = zero;
289 Real dTde = (dTdt * (grad_hi + grad_hi_jp1) > zero) ? grad_hi : grad_hi_jp1;
290 Real cff = std::max(dTdx*dTdx + dTde*dTde,eps);
291 Real Cx = dTdt * dTdx;
292 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);
293 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)));
294 }
295 });
296 ParallelFor(grow(xhi_ghost,IntVect(0,-1,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
297 {
298 dest_arr(i,j,k,icomp+icomp_to_fill) = dest_arr(lbound(xhi).x,j,k,icomp+icomp_to_fill);
299 });
300 }
301
302 // See comment on xlo
303 if (!ylo.isEmpty() && apply_south) {
304 ParallelFor(grow(ylo_edge,IntVect(-1,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
305 {
306 Real bry_val = bdatylo(i,ubound(ylo).y,k,0);
307 if (bcr.lo(1) == REMORABCType::clamped) {
308 dest_arr(i,j,k,icomp+icomp_to_fill) = bry_val * mask_arr(i,j,0);
309 } else if (bcr.lo(1) == REMORABCType::flather) {
310 Real bry_val_zeta = bdatylo_zeta(i,ubound(ylo).y-1,k,0);
311 Real cff = one / (Real(0.5) * (h_arr(i,dom_lo.y-1,0) + zeta_arr(i,dom_lo.y-1,0,icomp_calc)
312 + h_arr(i,dom_lo.y,0) + zeta_arr(i,dom_lo.y,0,icomp_calc)));
313 Real Ce = std::sqrt(l_g * cff);
314 dest_arr(i,j,k,icomp+icomp_to_fill) = (bry_val
315 - Ce * (Real(0.5) * (zeta_arr(i,dom_lo.y-1,0,icomp_calc) + zeta_arr(i,dom_lo.y,0,icomp_calc))
316 - bry_val_zeta)) * mask_arr(i,j,0);
317 } else if (bcr.lo(1) == REMORABCType::chapman) {
318 Real cff = dt_calc * Real(0.5) * (pn(i-mf_index_type[0],dom_lo.y,0) + pn(i,dom_lo.y,0));
319 Real cff1 = std::sqrt(l_g * Real(0.5) * (h_arr(i-mf_index_type[0],dom_lo.y,0) +
320 zeta_arr(i-mf_index_type[0],dom_lo.y,0,icomp_calc) + h_arr(i,dom_lo.y,0)
321 + zeta_arr(i,dom_lo.y,0,icomp_calc)));
322 Real Ce = cff * cff1;
323 Real cff2 = one / (one + Ce);
324 dest_arr(i,j,k,icomp+icomp_to_fill) = cff2 * (dest_arr(i,dom_lo.y-1,k,icomp_calc)
325 + Ce * dest_arr(i,dom_lo.y,k,icomp+icomp_to_fill)) * mask_arr(i,j,0);
326 } else if (bcr.lo(1) == REMORABCType::orlanski_rad_nudge) {
327 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));
328 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));
329 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));
330 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));
331 if (cell_centered) {
332 grad_lo *= msku(i ,dom_lo.y+mf_index_type[1] ,0);
333 grad_lo_jm1 *= msku(i ,dom_lo.y+mf_index_type[1]-1,0);
334 grad_lo_ip1 *= msku(i+1,dom_lo.y+mf_index_type[1] ,0);
335 grad_lo_ipjm1 *= msku(i+1,dom_lo.y+mf_index_type[1]-1,0);
336 }
337 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);
338 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);
339 Real tau;
340 Real nudg_coeff_out_local = (nudg_coeff_out(i-mf_index_type[0],j-mf_index_type[1],k) +
341 nudg_coeff_out(i,j,k)) * Real(0.5);
342 if (dTdt*dTde < zero) {
343 tau = nudg_coeff_out_local * obcfac * dt_calc;
344 dTdt = zero;
345 } else {
346 tau = nudg_coeff_out_local * dt_calc;
347 }
348 if (dTdt * dTde < zero) dTdt = zero;
349 Real dTdx = (dTdt * (grad_lo + grad_lo_ip1) > zero) ? grad_lo : grad_lo_ip1;
350 Real cff = std::max(dTdx*dTdx + dTde*dTde, eps);
351 Real Ce = dTdt*dTde;
352 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);
353 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)));
354 }
355 });
356 ParallelFor(grow(ylo_ghost,IntVect(-1,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
357 {
358 dest_arr(i,j,k,icomp+icomp_to_fill) = dest_arr(i,ubound(ylo).y,k,icomp+icomp_to_fill);
359 });
360 }
361
362 // See comment on xlo
363 if (!yhi.isEmpty() && apply_north) {
364 ParallelFor(grow(yhi_edge,IntVect(-1,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
365 {
366 Real bry_val = bdatyhi(i,lbound(yhi).y,k,0);
367 if (bcr.hi(1) == REMORABCType::clamped) {
368 dest_arr(i,j,k,icomp+icomp_to_fill) = bry_val * mask_arr(i,j,0);
369 } else if (bcr.hi(1) == REMORABCType::flather) {
370 Real bry_val_zeta = bdatyhi_zeta(i,lbound(yhi).y,k,0);
371 Real cff = one / (Real(0.5) * (h_arr(i,dom_hi.y-1,0) + zeta_arr(i,dom_hi.y-1,0,icomp_calc)
372 + h_arr(i,dom_hi.y,0) + zeta_arr(i,dom_hi.y,0,icomp_calc)));
373 Real Ce = std::sqrt(l_g * cff);
374 dest_arr(i,j,k,icomp+icomp_to_fill) = (bry_val
375 + Ce * (Real(0.5) * (zeta_arr(i,dom_hi.y-1,0,icomp_calc) + zeta_arr(i,dom_hi.y,0,icomp_calc))
376 - bry_val_zeta)) * mask_arr(i,j,0);
377 } else if (bcr.hi(1) == REMORABCType::chapman) {
378 Real cff = dt_calc * Real(0.5) * (pn(i-mf_index_type[0],dom_hi.y,0) + pn(i,dom_hi.y,0));
379 Real cff1 = std::sqrt(l_g * Real(0.5) * (h_arr(i-mf_index_type[0],dom_hi.y,0)
380 + zeta_arr(i-mf_index_type[0],dom_hi.y,0,icomp_calc) +
381 h_arr(i,dom_hi.y,0) + zeta_arr(i,dom_hi.y,0,icomp_calc)));
382 Real Ce = cff * cff1;
383 Real cff2 = one / (one + Ce);
384 dest_arr(i,j,k,icomp+icomp_to_fill) = cff2 * (dest_arr(i,dom_hi.y+1,k,icomp_calc)
385 + Ce * dest_arr(i,dom_hi.y,k,icomp+icomp_to_fill)) * mask_arr(i,j,0);
386 } else if (bcr.hi(1) == REMORABCType::orlanski_rad_nudge) {
387 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);
388 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);
389 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);
390 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);
391 if (cell_centered) {
392 grad_hi *= msku(i ,dom_hi.y-mf_index_type[1] ,0);
393 grad_hi_jp1 *= msku(i ,dom_hi.y-mf_index_type[1]+1,0);
394 grad_hi_ip1 *= msku(i+1,dom_hi.y-mf_index_type[1] ,0);
395 grad_hi_ijp1 *= msku(i+1,dom_hi.y-mf_index_type[1]+1,0);
396 }
397 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);
398 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);
399 Real tau;
400 Real nudg_coeff_out_local = (nudg_coeff_out(i-mf_index_type[0],j-mf_index_type[1],k) +
401 nudg_coeff_out(i,j,k)) * Real(0.5);
402 if (dTdt*dTde < zero) {
403 tau = nudg_coeff_out_local * obcfac * dt_calc;
404 dTdt = zero;
405 } else {
406 tau = nudg_coeff_out_local * dt_calc;
407 }
408 if (dTdt * dTde < zero) dTdt = zero;
409 Real dTdx = (dTdt * (grad_hi + grad_hi_ip1) > zero) ? grad_hi : grad_hi_ip1;
410 Real cff = std::max(dTdx*dTdx + dTde*dTde, eps);
411 Real Ce = dTdt*dTde;
412 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);
413 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)));
414 }
415 });
416 ParallelFor(grow(yhi_ghost,IntVect(-1,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
417 {
418 dest_arr(i,j,k,icomp+icomp_to_fill) = dest_arr(i,lbound(yhi).y,k,icomp+icomp_to_fill);
419 });
420
421 }
422 // If we've applied boundary conditions to either side, update the corner
423 if (!xlo_ylo.isEmpty() && (apply_west || apply_south)) {
424 ParallelFor(xlo_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k)
425 {
426 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)
427 + dest_arr(dom_lo.x+mf_index_type[0],j,k,icomp+icomp_to_fill));
428 });
429 }
430 if (!xlo_yhi.isEmpty() && (apply_west || apply_north)) {
431 ParallelFor(xlo_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k)
432 {
433 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)
434 + dest_arr(dom_lo.x+mf_index_type[0],j,k,icomp+icomp_to_fill));
435 });
436 }
437 if (!xhi_ylo.isEmpty() && (apply_east || apply_south)) {
438 ParallelFor(xhi_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k)
439 {
440 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)
441 + dest_arr(dom_hi.x-mf_index_type[0],j,k,icomp+icomp_to_fill));
442 });
443 }
444 if (!xhi_yhi.isEmpty() && (apply_east || apply_north)) {
445 ParallelFor(xhi_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k)
446 {
447 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)
448 + dest_arr(dom_hi.x-mf_index_type[0],j,k,icomp+icomp_to_fill));
449 });
450 }
451 } // mfi
452 } // icomp
453}
454#endif
one
constexpr amrex::Real one
Definition REMORA_Constants.H:7
zero
constexpr amrex::Real zero
Definition REMORA_Constants.H:6
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:1797
REMORA::domain_bcs_type
amrex::Vector< amrex::BCRec > domain_bcs_type
vector (over BCVars) of BCRecs
Definition REMORA.H:1728
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:413
REMORA::vec_pm
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pm
horizontal scaling factor: 1 / dx (2D)
Definition REMORA.H:603
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:586
REMORA::zeta_bc
int zeta_bc() const noexcept
Definition REMORA.H:1453
REMORA::bdy_zeta
int bdy_zeta() const noexcept
Definition REMORA.H:1466
REMORA::vec_mskv
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskv
land/sea mask at y-faces (2D)
Definition REMORA.H:588
REMORA::solverChoice
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1949
REMORA::boundary_series
amrex::Vector< amrex::Vector< std::unique_ptr< NCTimeSeriesBoundary > > > boundary_series
Vector over BdyVars of boundary series data containers.
Definition REMORA.H:1623
REMORA::vec_zeta
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_zeta
free surface height (2D)
Definition REMORA.H:578
REMORA::ubar_bc
int ubar_bc() const noexcept
Definition REMORA.H:1451
REMORA::vbar_bc
int vbar_bc() const noexcept
Definition REMORA.H:1452
REMORA::vec_nudg_coeff
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > vec_nudg_coeff
Climatology nudging coefficients.
Definition REMORA.H:678
REMORA::vec_pn
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pn
horizontal scaling factor: 1 / dy (2D)
Definition REMORA.H:605
REMORA::model_time
double model_time(amrex::Real elapsed) const noexcept
Time on the model clock, in seconds, of an elapsed time such as t_new.
Definition REMORA.H:2170
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