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 // amrex::Print() << "TIME " << time << std::endl;
24
25 // Time interpolation
26 Real dT = bdy_time_interval;
27 // amrex::Print() << "DT " << dT << std::endl;
28 // amrex::Print() << "START " << start_bdy_time << std::endl;
29
30 Real time_since_start = time - start_bdy_time;
31 int n_time = static_cast<int>( time_since_start / dT);
32
33 amrex::Real alpha = (time_since_start - n_time * dT) / dT;
34 AMREX_ALWAYS_ASSERT( alpha >= 0. && alpha <= 1.0_rt);
35 amrex::Real oma = 1.0_rt - alpha;
36
37 // Which variable are we filling
38 int ivar = bdy_var_type;
39
40 //
41 // Note that "domain" is mapped onto the type of box the data is in
42 //
43 Box domain = geom[lev].Domain();
44
45 const auto& mf_index_type = mf_to_fill.boxArray().ixType();
46 domain.convert(mf_index_type);
47
48 const auto& dom_lo = amrex::lbound(domain);
49 const auto& dom_hi = amrex::ubound(domain);
50
51 int ncomp;
52
53 // If we are doing the scalars then do salt as well as temp
54 if (ivar == BdyVars::t) {
55 ncomp = 2;
56 } else {
57 ncomp = 1;
58 }
59
60 // This must be true for the logic below to work
61 AMREX_ALWAYS_ASSERT(Temp_comp == 0);
62 AMREX_ALWAYS_ASSERT(Salt_comp == 1);
63
64 // Make sure we can interpolate in time
65 AMREX_ALWAYS_ASSERT(n_time + 1 < bdy_data_xlo.size());
66
67 const Real eps= 1.0e-20_rt;
68 const bool null_mf_calc = (!mf_calc.ok());
69
70 for (int icomp = 0; icomp < ncomp; icomp++) // This is to do both temp and salt if doing scalars
71 {
72 // We have data at fixed time intervals we will call dT
73 // Then to interpolate, given time, we can define n = (time/dT)
74 // and alpha = (time - n*dT) / dT, then we define the data at time
75 // as alpha * (data at time n+1) + (1 - alpha) * (data at time n)
76 const auto& bdatxlo_n = bdy_data_xlo[n_time ][ivar+icomp].const_array();
77 const auto& bdatxlo_np1 = bdy_data_xlo[n_time+1][ivar+icomp].const_array();
78 const auto& bdatxhi_n = bdy_data_xhi[n_time ][ivar+icomp].const_array();
79 const auto& bdatxhi_np1 = bdy_data_xhi[n_time+1][ivar+icomp].const_array();
80 const auto& bdatylo_n = bdy_data_ylo[n_time ][ivar+icomp].const_array();
81 const auto& bdatylo_np1 = bdy_data_ylo[n_time+1][ivar+icomp].const_array();
82 const auto& bdatyhi_n = bdy_data_yhi[n_time ][ivar+icomp].const_array();
83 const auto& bdatyhi_np1 = bdy_data_yhi[n_time+1][ivar+icomp].const_array();
84
85 const auto& bdatzetaxlo_n = bdy_data_xlo[n_time ][BdyVars::zeta].const_array();
86 const auto& bdatzetaxlo_np1 = bdy_data_xlo[n_time+1][BdyVars::zeta].const_array();
87 const auto& bdatzetaxhi_n = bdy_data_xhi[n_time ][BdyVars::zeta].const_array();
88 const auto& bdatzetaxhi_np1 = bdy_data_xhi[n_time+1][BdyVars::zeta].const_array();
89 const auto& bdatzetaylo_n = bdy_data_ylo[n_time ][BdyVars::zeta].const_array();
90 const auto& bdatzetaylo_np1 = bdy_data_ylo[n_time+1][BdyVars::zeta].const_array();
91 const auto& bdatzetayhi_n = bdy_data_yhi[n_time ][BdyVars::zeta].const_array();
92 const auto& bdatzetayhi_np1 = bdy_data_yhi[n_time+1][BdyVars::zeta].const_array();
93
94 const bool apply_west = (domain_bcs_type[bccomp+icomp].lo(0) == REMORABCType::clamped) ||
95 (domain_bcs_type[bccomp+icomp].lo(0) == REMORABCType::flather) ||
96 (domain_bcs_type[bccomp+icomp].lo(0) == REMORABCType::chapman) ||
97 (domain_bcs_type[bccomp+icomp].lo(0) == REMORABCType::orlanski_rad_nudge);
98 const bool apply_east = (domain_bcs_type[bccomp+icomp].hi(0) == REMORABCType::clamped) ||
99 (domain_bcs_type[bccomp+icomp].hi(0) == REMORABCType::flather) ||
100 (domain_bcs_type[bccomp+icomp].hi(0) == REMORABCType::chapman) ||
101 (domain_bcs_type[bccomp+icomp].hi(0) == REMORABCType::orlanski_rad_nudge);
102 const bool apply_north = (domain_bcs_type[bccomp+icomp].lo(1) == REMORABCType::clamped) ||
103 (domain_bcs_type[bccomp+icomp].lo(1) == REMORABCType::flather) ||
104 (domain_bcs_type[bccomp+icomp].lo(1) == REMORABCType::chapman) ||
105 (domain_bcs_type[bccomp+icomp].lo(1) == REMORABCType::orlanski_rad_nudge);
106 const bool apply_south = (domain_bcs_type[bccomp+icomp].hi(1) == REMORABCType::clamped) ||
107 (domain_bcs_type[bccomp+icomp].hi(1) == REMORABCType::flather) ||
108 (domain_bcs_type[bccomp+icomp].hi(1) == REMORABCType::chapman) ||
109 (domain_bcs_type[bccomp+icomp].hi(1) == REMORABCType::orlanski_rad_nudge);
110
111 const bool cell_centered = (mf_index_type[0] == 0 and mf_index_type[1] == 0);
112
113 const Real obcfac = solverChoice.obcfac;
114
115#ifdef AMREX_USE_OMP
116#pragma omp parallel if (Gpu::notInLaunchRegion())
117#endif
118 // Currently no tiling in order to get the logic right
119 for (MFIter mfi(mf_to_fill,false); mfi.isValid(); ++mfi)
120 {
121 Box mf_box(mf_to_fill[mfi.index()].box());
122
123 // Compute intersections of the FAB to be filled and the bdry data boxes
124 Box xlo_tmp(bdy_data_xlo[n_time][ivar].box());
125 Box xlo = xlo_tmp & mf_box;
126 Box xhi_tmp(bdy_data_xhi[n_time][ivar].box()); Box xhi = xhi_tmp & mf_box;
127 Box ylo_tmp(bdy_data_ylo[n_time][ivar].box()); Box ylo = ylo_tmp & mf_box;
128 Box yhi_tmp(bdy_data_yhi[n_time][ivar].box()); Box yhi = yhi_tmp & 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,dom_lo.x-1+mf_index_type[0]); xlo_edge.setBig(0,dom_lo.x-1+mf_index_type[0]);
141// Box xhi_edge = xhi; xhi_edge.setSmall(0,dom_hi.x+1-mf_index_type[0]); xhi_edge.setBig(0,dom_hi.x+1-mf_index_type[0]);
142// Box ylo_edge = ylo; ylo_edge.setSmall(1,dom_lo.y-1+mf_index_type[1]); ylo_edge.setBig(1,dom_lo.x-1+mf_index_type[1]);
143// Box yhi_edge = yhi; yhi_edge.setSmall(1,dom_hi.y+1-mf_index_type[1]); yhi_edge.setBig(1,dom_hi.x+1-mf_index_type[1]);
144 Box xlo_edge = xlo; xlo_edge.setSmall(0,ubound(xlo).x); xlo_edge.setBig(0,ubound(xlo).x);
145 Box xhi_edge = xhi; xhi_edge.setSmall(0,lbound(xhi).x); xhi_edge.setBig(0,lbound(xhi).x);
146 Box ylo_edge = ylo; ylo_edge.setSmall(1,ubound(ylo).y); ylo_edge.setBig(1,ubound(ylo).y);
147 Box yhi_edge = yhi; yhi_edge.setSmall(1,lbound(yhi).y); yhi_edge.setBig(1,lbound(yhi).y);
148
149 Box xlo_ghost = xlo; xlo_ghost.setBig(0,ubound(xlo).x-1);
150 Box xhi_ghost = xhi; xhi_ghost.setSmall(0,lbound(xhi).x+1);
151 Box ylo_ghost = ylo; ylo_ghost.setBig(1,ubound(ylo).y-1);
152 Box yhi_ghost = yhi; yhi_ghost.setSmall(1,lbound(yhi).y+1);
153
154 const Array4<Real>& dest_arr = mf_to_fill.array(mfi);
155 const Array4<const Real>& mask_arr = mf_mask.array(mfi);
156 const Array4<const Real>& calc_arr = (!null_mf_calc) ? mf_calc.array(mfi) : Array4<amrex::Real>();
157 const Array4<const Real>& h_arr = vec_hOfTheConfusingName[lev]->const_array(mfi);
158 const Array4<const Real>& zeta_arr = vec_zeta[lev]->const_array(mfi);
159 const Array4<const Real>& pm = vec_pm[lev]->const_array(mfi);
160 const Array4<const Real>& pn = vec_pn[lev]->const_array(mfi);
161
162 const Array4<const Real>& msku = vec_msku[lev]->const_array(mfi);
163 const Array4<const Real>& mskv = vec_mskv[lev]->const_array(mfi);
164
165 const Array4<const Real> nudg_coeff_out = vec_nudg_coeff[bdy_var_type][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 = (oma * bdatxlo_n (ubound(xlo).x,j,k,0)
188 + alpha * bdatxlo_np1(ubound(xlo).x,j,k,0));
189 if (bcr.lo(0) == REMORABCType::clamped) {
190 dest_arr(i,j,k,icomp+icomp_to_fill) = bry_val * mask_arr(i,j,0);
191 } else if (bcr.lo(0) == REMORABCType::flather) {
192 Real bry_val_zeta = (oma * bdatzetaxlo_n(ubound(xlo).x-1,j,k,0) + alpha * bdatzetaxlo_np1(ubound(xlo).x-1,j,k,0));
193 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)
194 + h_arr(dom_lo.x,j,0) + zeta_arr(dom_lo.x,j,0,icomp_calc)));
195 Real Cx = std::sqrt(g * cff);
196 dest_arr(i,j,k,icomp+icomp_to_fill) = (bry_val
197 - Cx * (0.5_rt * (zeta_arr(dom_lo.x-1,j,0,icomp_calc) + zeta_arr(dom_lo.x,j,0,icomp_calc))
198 - bry_val_zeta)) * mask_arr(i,j,0);
199 } else if (bcr.lo(0) == REMORABCType::chapman) {
200 Real cff = dt_calc * 0.5_rt * (pm(dom_lo.x,j-mf_index_type[1],0) + pm(dom_lo.x,j,0));
201 Real cff1 = std::sqrt(g * 0.5_rt * (h_arr(dom_lo.x,j-mf_index_type[1],0)
202 + zeta_arr(dom_lo.x,j-mf_index_type[1],0,icomp_calc) + h_arr(dom_lo.x,j,0)
203 + zeta_arr(dom_lo.x,j,0,icomp_calc)));
204 Real Cx = cff * cff1;
205 Real cff2 = 1.0_rt / (1.0_rt + Cx);
206 dest_arr(i,j,k,icomp+icomp_to_fill) = cff2 * (dest_arr(dom_lo.x-1,j,k,icomp_calc)
207 + Cx * dest_arr(dom_lo.x,j,k,icomp+icomp_to_fill)) * mask_arr(i,j,0);
208 } else if (bcr.lo(0) == REMORABCType::orlanski_rad_nudge) {
209 Real grad_lo_im1 = (calc_arr(dom_lo.x+mf_index_type[0]-1,j ,k,icomp+icomp_to_fill) - calc_arr(dom_lo.x-1+mf_index_type[0],j-1,k,icomp+icomp_to_fill));
210 Real grad_lo = (calc_arr(dom_lo.x+mf_index_type[0] ,j ,k,icomp+icomp_to_fill) - calc_arr(dom_lo.x +mf_index_type[0],j-1,k,icomp+icomp_to_fill));
211 Real grad_lo_imjp1 = (calc_arr(dom_lo.x+mf_index_type[0]-1,j+1,k,icomp+icomp_to_fill) - calc_arr(dom_lo.x-1+mf_index_type[0],j ,k,icomp+icomp_to_fill));
212 Real grad_lo_jp1 = (calc_arr(dom_lo.x+mf_index_type[0] ,j+1,k,icomp+icomp_to_fill) - calc_arr(dom_lo.x +mf_index_type[0],j ,k,icomp+icomp_to_fill));
213 if (cell_centered) {
214 grad_lo_im1 *= mskv(i,j,0);
215 grad_lo *= mskv(i,j,0);
216 grad_lo_imjp1 *= mskv(i,j,0);
217 grad_lo_jp1 *= mskv(i,j,0);
218 }
219 Real dTdt = calc_arr(dom_lo.x+mf_index_type[0],j,k,icomp+icomp_to_fill) - dest_arr(dom_lo.x+mf_index_type[0] ,j,k,icomp+icomp_to_fill);
220 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);
221 Real tau;
222 Real nudg_coeff_out_local = (nudg_coeff_out(i-mf_index_type[0],j-mf_index_type[1],k) +
223 nudg_coeff_out(i,j,k)) * 0.5_rt;
224 if (dTdt*dTdx < 0.0_rt) {
225 tau = nudg_coeff_out_local * obcfac * dt_calc;
226 dTdt = 0.0_rt;
227 } else {
228 tau = nudg_coeff_out_local * dt_calc;
229 }
230 Real dTde = (dTdt * (grad_lo+grad_lo_jp1) > 0.0_rt) ? grad_lo : grad_lo_jp1;
231 Real cff = std::max(dTdx*dTdx+dTde*dTde,eps);
232 Real Cx = dTdt * dTdx;
233 dest_arr(i,j,k,icomp+icomp_to_fill) = (cff * calc_arr(dom_lo.x-1+mf_index_type[0],j,k,icomp+icomp_to_fill) + Cx * dest_arr(dom_lo.x+mf_index_type[0],j,k,icomp+icomp_to_fill)) / (cff+Cx);
234 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)));
235 }
236 });
237 ParallelFor(grow(xlo_ghost,IntVect(0,-1,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
238 {
239 dest_arr(i,j,k,icomp+icomp_to_fill) = dest_arr(ubound(xlo).x,j,k,icomp+icomp_to_fill);
240 });
241 }
242
243 // See comment on xlo
244 if (!xhi.isEmpty() && apply_east) {
245 ParallelFor(grow(xhi_edge,IntVect(0,-1,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
246 {
247 Real bry_val = (oma * bdatxhi_n (lbound(xhi).x,j,k,0)
248 + alpha * bdatxhi_np1(lbound(xhi).x,j,k,0));
249 if (bcr.hi(0) == REMORABCType::clamped) {
250 dest_arr(i,j,k,icomp+icomp_to_fill) = bry_val * mask_arr(i,j,0);
251 } else if (bcr.hi(0) == REMORABCType::flather) {
252 Real bry_val_zeta = (oma * bdatzetaxhi_n(lbound(xhi).x,j,k,0) + alpha * bdatzetaxhi_np1(lbound(xhi).x,j,k,0));
253 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)
254 + h_arr(dom_hi.x,j,0) + zeta_arr(dom_hi.x,j,0,icomp_calc)));
255 Real Cx = std::sqrt(g * cff);
256 dest_arr(i,j,k,icomp+icomp_to_fill) = (bry_val
257 + Cx * (0.5_rt * (zeta_arr(dom_hi.x-1,j,0,icomp_calc) + zeta_arr(dom_hi.x,j,0,icomp_calc))
258 - bry_val_zeta)) * mask_arr(i,j,0);
259 } else if (bcr.hi(0) == REMORABCType::chapman) {
260 Real cff = dt_calc * 0.5_rt * (pm(dom_hi.x,j-mf_index_type[1],0) + pm(dom_hi.x,j,0));
261 Real cff1 = std::sqrt(g * 0.5_rt * (h_arr(dom_hi.x,j-mf_index_type[1],0)
262 + zeta_arr(dom_hi.x,j-mf_index_type[1],0,icomp_calc) + h_arr(dom_hi.x,j,0)
263 + zeta_arr(dom_hi.x,j,0,icomp_calc)));
264 Real Cx = cff * cff1;
265 Real cff2 = 1.0_rt / (1.0_rt + Cx);
266 dest_arr(i,j,k,icomp+icomp_to_fill) = cff2 * (dest_arr(dom_hi.x+1,j,k,icomp_calc)
267 + Cx * dest_arr(dom_hi.x,j,k,icomp+icomp_to_fill)) * mask_arr(i,j,0);
268 } else if (bcr.hi(0) == REMORABCType::orlanski_rad_nudge) {
269 Real grad_hi = (calc_arr(dom_hi.x-mf_index_type[0] ,j ,k,icomp+icomp_to_fill) - calc_arr(dom_hi.x-mf_index_type[0] ,j-1,k,icomp+icomp_to_fill));
270 Real grad_hi_ip1 = (calc_arr(dom_hi.x-mf_index_type[0]+1,j ,k,icomp+icomp_to_fill) - calc_arr(dom_hi.x-mf_index_type[0]+1,j-1,k,icomp+icomp_to_fill));
271 Real grad_hi_jp1 = (calc_arr(dom_hi.x-mf_index_type[0] ,j+1,k,icomp+icomp_to_fill) - calc_arr(dom_hi.x-mf_index_type[0] ,j ,k,icomp+icomp_to_fill));
272 Real grad_hi_ijp1 = (calc_arr(dom_hi.x-mf_index_type[0]+1,j+1,k,icomp+icomp_to_fill) - calc_arr(dom_hi.x-mf_index_type[0]+1,j ,k,icomp+icomp_to_fill));
273 if (cell_centered) {
274 grad_hi *= mskv(i,j,0);
275 grad_hi_ip1 *= mskv(i,j,0);
276 grad_hi_jp1 *= mskv(i,j,0);
277 grad_hi_ijp1 *= mskv(i,j,0);
278 }
279 Real dTdt = calc_arr(dom_hi.x-mf_index_type[0],j,k,icomp+icomp_to_fill) - dest_arr(dom_hi.x-mf_index_type[0] ,j,k,icomp+icomp_to_fill);
280 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);
281 Real tau;
282 Real nudg_coeff_out_local = (nudg_coeff_out(i-mf_index_type[0],j-mf_index_type[1],k) +
283 nudg_coeff_out(i,j,k)) * 0.5_rt;
284 if (dTdt*dTdx < 0.0_rt) {
285 tau = nudg_coeff_out_local * obcfac * dt_calc;
286 dTdt = 0.0_rt;
287 } else {
288 tau = nudg_coeff_out_local * dt_calc;
289 }
290 if (dTdt * dTdx < 0.0_rt) dTdt = 0.0_rt;
291 Real dTde = (dTdt * (grad_hi + grad_hi_jp1) > 0.0_rt) ? grad_hi : grad_hi_jp1;
292 Real cff = std::max(dTdx*dTdx + dTde*dTde,eps);
293 Real Cx = dTdt * dTdx;
294 dest_arr(i,j,k,icomp+icomp_to_fill) = (cff * calc_arr(dom_hi.x+1-mf_index_type[0],j,k,icomp+icomp_to_fill) + Cx * dest_arr(dom_hi.x-mf_index_type[0],j,k,icomp+icomp_to_fill)) * mask_arr(i,j,0) / (cff+Cx);
295 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)));
296 }
297 });
298 ParallelFor(grow(xhi_ghost,IntVect(0,-1,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
299 {
300 dest_arr(i,j,k,icomp+icomp_to_fill) = dest_arr(lbound(xhi).x,j,k,icomp+icomp_to_fill);
301 });
302 }
303
304 // See comment on xlo
305 if (!ylo.isEmpty() && apply_south) {
306 ParallelFor(grow(ylo_edge,IntVect(-1,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
307 {
308 Real bry_val = (oma * bdatylo_n (i,ubound(ylo).y,k,0)
309 + alpha * bdatylo_np1(i,ubound(ylo).y,k,0));
310 if (bcr.lo(1) == REMORABCType::clamped) {
311 dest_arr(i,j,k,icomp+icomp_to_fill) = bry_val * mask_arr(i,j,0);
312 } else if (bcr.lo(1) == REMORABCType::flather) {
313 Real bry_val_zeta = (oma * bdatzetaylo_n (i,ubound(ylo).y-1,k,0)
314 + alpha * bdatzetaylo_np1(i,ubound(ylo).y-1,k,0));
315 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)
316 + h_arr(i,dom_lo.y,0) + zeta_arr(i,dom_lo.y,0,icomp_calc)));
317 Real Ce = std::sqrt(g * cff);
318 dest_arr(i,j,k,icomp+icomp_to_fill) = (bry_val
319 - Ce * (0.5_rt * (zeta_arr(i,dom_lo.y-1,0,icomp_calc) + zeta_arr(i,dom_lo.y,0,icomp_calc))
320 - bry_val_zeta)) * mask_arr(i,j,0);
321 } else if (bcr.lo(1) == REMORABCType::chapman) {
322 Real cff = dt_calc * 0.5_rt * (pn(i-mf_index_type[0],dom_lo.y,0) + pn(i,dom_lo.y,0));
323 Real cff1 = std::sqrt(g * 0.5_rt * (h_arr(i-mf_index_type[0],dom_lo.y,0) +
324 zeta_arr(i-mf_index_type[0],dom_lo.y,0,icomp_calc) + h_arr(i,dom_lo.y,0)
325 + zeta_arr(i,dom_lo.y,0,icomp_calc)));
326 Real Ce = cff * cff1;
327 Real cff2 = 1.0_rt / (1.0_rt + Ce);
328 dest_arr(i,j,k,icomp+icomp_to_fill) = cff2 * (dest_arr(i,dom_lo.y-1,k,icomp_calc)
329 + Ce * dest_arr(i,dom_lo.y,k,icomp+icomp_to_fill)) * mask_arr(i,j,0);
330 } else if (bcr.lo(1) == REMORABCType::orlanski_rad_nudge) {
331 Real grad_lo = (calc_arr(i ,dom_lo.y+mf_index_type[1], k,icomp+icomp_to_fill) - calc_arr(i-1,dom_lo.y+mf_index_type[1] ,k,icomp+icomp_to_fill));
332 Real grad_lo_jm1 = (calc_arr(i ,dom_lo.y+mf_index_type[1]-1,k,icomp+icomp_to_fill) - calc_arr(i-1,dom_lo.y+mf_index_type[1]-1,k,icomp+icomp_to_fill));
333 Real grad_lo_ip1 = (calc_arr(i+1,dom_lo.y+mf_index_type[1] ,k,icomp+icomp_to_fill) - calc_arr(i ,dom_lo.y+mf_index_type[1] ,k,icomp+icomp_to_fill));
334 Real grad_lo_ipjm1 = (calc_arr(i+1,dom_lo.y+mf_index_type[1]-1,k,icomp+icomp_to_fill) - calc_arr(i ,dom_lo.y+mf_index_type[1]-1,k,icomp+icomp_to_fill));
335 if (cell_centered) {
336 grad_lo *= msku(i,j,0);
337 grad_lo_jm1 *= msku(i,j,0);
338 grad_lo_ip1 *= msku(i,j,0);
339 grad_lo_ipjm1 *= msku(i,j,0);
340 }
341 Real dTdt = calc_arr(i,dom_lo.y+mf_index_type[1],k,icomp+icomp_to_fill) - dest_arr(i,dom_lo.y +mf_index_type[1],k,icomp+icomp_to_fill);
342 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);
343 Real tau;
344 Real nudg_coeff_out_local = (nudg_coeff_out(i-mf_index_type[0],j-mf_index_type[1],k) +
345 nudg_coeff_out(i,j,k)) * 0.5_rt;
346 if (dTdt*dTde < 0.0_rt) {
347 tau = nudg_coeff_out_local * obcfac * dt_calc;
348 dTdt = 0.0_rt;
349 } else {
350 tau = nudg_coeff_out_local * dt_calc;
351 }
352 if (dTdt * dTde < 0.0_rt) dTdt = 0.0_rt;
353 Real dTdx = (dTdt * (grad_lo + grad_lo_ip1) > 0.0_rt) ? grad_lo : grad_lo_ip1;
354 Real cff = std::max(dTdx*dTdx + dTde*dTde, eps);
355 Real Ce = dTdt*dTde;
356 dest_arr(i,j,k,icomp+icomp_to_fill) = (cff * calc_arr(i,dom_lo.y-1+mf_index_type[1],k,icomp+icomp_to_fill) + Ce * dest_arr(i,dom_lo.y+mf_index_type[1],k,icomp+icomp_to_fill)) / (cff+Ce);
357 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)));
358 }
359 });
360 ParallelFor(grow(ylo_ghost,IntVect(-1,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
361 {
362 dest_arr(i,j,k,icomp+icomp_to_fill) = dest_arr(i,ubound(ylo).y,k,icomp+icomp_to_fill);
363 });
364 }
365
366 // See comment on xlo
367 if (!yhi.isEmpty() && apply_north) {
368 ParallelFor(grow(yhi_edge,IntVect(-1,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
369 {
370 Real bry_val = (oma * bdatyhi_n (i,lbound(yhi).y,k,0)
371 + alpha * bdatyhi_np1(i,lbound(yhi).y,k,0)) * mask_arr(i,j,0);
372 if (bcr.hi(1) == REMORABCType::clamped) {
373 dest_arr(i,j,k,icomp+icomp_to_fill) = bry_val * mask_arr(i,j,0);
374 } else if (bcr.hi(1) == REMORABCType::flather) {
375 Real bry_val_zeta = (oma * bdatzetayhi_n (i,lbound(yhi).y,k,0)
376 + alpha * bdatzetayhi_np1(i,lbound(yhi).y,k,0));
377 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)
378 + h_arr(i,dom_hi.y,0) + zeta_arr(i,dom_hi.y,0,icomp_calc)));
379 Real Ce = std::sqrt(g * cff);
380 dest_arr(i,j,k,icomp+icomp_to_fill) = (bry_val
381 + Ce * (0.5_rt * (zeta_arr(i,dom_hi.y-1,0,icomp_calc) + zeta_arr(i,dom_hi.y,0,icomp_calc))
382 - bry_val_zeta)) * mask_arr(i,j,0);
383 } else if (bcr.hi(1) == REMORABCType::chapman) {
384 Real cff = dt_calc * 0.5_rt * (pn(i-mf_index_type[0],dom_hi.y,0) + pn(i,dom_hi.y,0));
385 Real cff1 = std::sqrt(g * 0.5_rt * (h_arr(i-mf_index_type[0],dom_hi.y,0)
386 + zeta_arr(i-mf_index_type[0],dom_hi.y,0,icomp_calc) +
387 h_arr(i,dom_hi.y,0) + zeta_arr(i,dom_hi.y,0,icomp_calc)));
388 Real Ce = cff * cff1;
389 Real cff2 = 1.0_rt / (1.0_rt + Ce);
390 dest_arr(i,j,k,icomp+icomp_to_fill) = cff2 * (dest_arr(i,dom_hi.y+1,k,icomp_calc)
391 + Ce * dest_arr(i,dom_hi.y,k,icomp+icomp_to_fill)) * mask_arr(i,j,0);
392 } else if (bcr.hi(1) == REMORABCType::orlanski_rad_nudge) {
393 Real grad_hi = calc_arr(i ,dom_hi.y-mf_index_type[1] ,k,icomp+icomp_to_fill) - calc_arr(i-1,dom_hi.y-mf_index_type[1] ,k,icomp+icomp_to_fill);
394 Real grad_hi_jp1 = calc_arr(i ,dom_hi.y-mf_index_type[1]+1,k,icomp+icomp_to_fill) - calc_arr(i-1,dom_hi.y-mf_index_type[1]+1,k,icomp+icomp_to_fill);
395 Real grad_hi_ip1 = calc_arr(i+1,dom_hi.y-mf_index_type[1] ,k,icomp+icomp_to_fill) - calc_arr(i ,dom_hi.y-mf_index_type[1] ,k,icomp+icomp_to_fill);
396 Real grad_hi_ijp1 = calc_arr(i+1,dom_hi.y-mf_index_type[1]+1,k,icomp+icomp_to_fill) - calc_arr(i ,dom_hi.y-mf_index_type[1]+1,k,icomp+icomp_to_fill);
397 if (cell_centered) {
398 grad_hi *= msku(i,j,0);
399 grad_hi_jp1 *= msku(i,j,0);
400 grad_hi_ip1 *= msku(i,j,0);
401 grad_hi_ijp1 *= msku(i,j,0);
402 }
403 Real dTdt = calc_arr(i,dom_hi.y-mf_index_type[1],k,icomp+icomp_to_fill) - dest_arr(i,dom_hi.y -mf_index_type[1],k,icomp+icomp_to_fill);
404 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);
405 Real tau;
406 Real nudg_coeff_out_local = (nudg_coeff_out(i-mf_index_type[0],j-mf_index_type[1],k) +
407 nudg_coeff_out(i,j,k)) * 0.5_rt;
408 if (dTdt*dTde < 0.0_rt) {
409 tau = nudg_coeff_out_local * obcfac * dt_calc;
410 dTdt = 0.0_rt;
411 } else {
412 tau = nudg_coeff_out_local * dt_calc;
413 }
414 if (dTdt * dTde < 0.0_rt) dTdt = 0.0_rt;
415 Real dTdx = (dTdt * (grad_hi + grad_hi_ip1) > 0.0_rt) ? grad_hi : grad_hi_ip1;
416 Real cff = std::max(dTdx*dTdx + dTde*dTde, eps);
417 Real Ce = dTdt*dTde;
418 dest_arr(i,j,k,icomp+icomp_to_fill) = (cff*calc_arr(i,dom_hi.y+1-mf_index_type[1],k,icomp+icomp_to_fill) + Ce*dest_arr(i,dom_hi.y-mf_index_type[1],k,icomp+icomp_to_fill)) * mask_arr(i,j,0) / (cff+Ce);
419 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)));
420 }
421 });
422 ParallelFor(grow(yhi_ghost,IntVect(-1,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
423 {
424 dest_arr(i,j,k,icomp+icomp_to_fill) = dest_arr(i,lbound(yhi).y,k,icomp+icomp_to_fill);
425 });
426
427 }
428 // If we've applied boundary conditions to either side, update the corner
429 if (!xlo_ylo.isEmpty() && (apply_west || apply_south)) {
430 ParallelFor(xlo_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k)
431 {
432 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)
433 + dest_arr(dom_lo.x+mf_index_type[0],j,k,icomp+icomp_to_fill));
434 });
435 }
436 if (!xlo_yhi.isEmpty() && (apply_west || apply_north)) {
437 ParallelFor(xlo_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k)
438 {
439 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)
440 + dest_arr(dom_lo.x+mf_index_type[0],j,k,icomp+icomp_to_fill));
441 });
442 }
443 if (!xhi_ylo.isEmpty() && (apply_east || apply_south)) {
444 ParallelFor(xhi_ylo, [=] AMREX_GPU_DEVICE (int i, int j, int k)
445 {
446 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)
447 + dest_arr(dom_hi.x-mf_index_type[0],j,k,icomp+icomp_to_fill));
448 });
449 }
450 if (!xhi_yhi.isEmpty() && (apply_east || apply_north)) {
451 ParallelFor(xhi_yhi, [=] AMREX_GPU_DEVICE (int i, int j, int k)
452 {
453 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)
454 + dest_arr(dom_hi.x-mf_index_type[0],j,k,icomp+icomp_to_fill));
455 });
456 }
457 } // mfi
458 } // icomp
459}
460#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:1120
REMORA::vec_hOfTheConfusingName
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_hOfTheConfusingName
Bathymetry data (2D, positive valued, h in ROMS)
Definition REMORA.H:229
REMORA::vec_pm
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pm
horizontal scaling factor: 1 / dx (2D)
Definition REMORA.H:355
REMORA::vec_msku
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_msku
land/sea mask at x-faces (2D)
Definition REMORA.H:348
REMORA::bdy_data_yhi
amrex::Vector< amrex::Vector< amrex::FArrayBox > > bdy_data_yhi
Vectors (over time) of Vector (over variables) of FArrayBoxs for holding Northern boundary data from ...
Definition REMORA.H:996
REMORA::bdy_time_interval
amrex::Real bdy_time_interval
Interval between boundary data times.
Definition REMORA.H:1001
REMORA::start_bdy_time
amrex::Real start_bdy_time
Start time in the time series of boundary data.
Definition REMORA.H:999
REMORA::vec_mskv
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskv
land/sea mask at y-faces (2D)
Definition REMORA.H:350
REMORA::bdy_data_ylo
amrex::Vector< amrex::Vector< amrex::FArrayBox > > bdy_data_ylo
Vectors (over time) of Vector (over variables) of FArrayBoxs for holding Southern boundary data from ...
Definition REMORA.H:994
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:1221
REMORA::bdy_data_xhi
amrex::Vector< amrex::Vector< amrex::FArrayBox > > bdy_data_xhi
Vectors (over time) of Vector (over variables) of FArrayBoxs for holding Eastern boundary data from f...
Definition REMORA.H:992
REMORA::vec_zeta
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_zeta
free surface height (2D)
Definition REMORA.H:343
REMORA::vec_nudg_coeff
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > vec_nudg_coeff
Climatology nudging coefficients.
Definition REMORA.H:414
REMORA::vec_pn
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pn
horizontal scaling factor: 1 / dy (2D)
Definition REMORA.H:357
REMORA::bdy_data_xlo
amrex::Vector< amrex::Vector< amrex::FArrayBox > > bdy_data_xlo
Vectors (over time) of Vector (over variables) of FArrayBoxs for holding Western boundary data from f...
Definition REMORA.H:990
SolverChoice::obcfac
amrex::Real obcfac
Definition REMORA_DataStruct.H:691