REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_NCTimeSeriesBoundary.cpp
Go to the documentation of this file.
1#include "REMORA_NCTimeSeriesBoundary.H"
2#include "REMORA_NCFile.H"
3
4#include "AMReX_ParallelDescriptor.H"
5
6#include <string>
7
8#ifdef REMORA_USE_NETCDF
9/**
10 * @param[in ] a_lev level at which we will store the data
11 * @param[in ] a_geom Vector of Geometry objects at all levels
12 * @param[in ] a_file_name file name to read from
13 * @param[in ] a_field_name name of field to read in
14 * @param[in ] a_time_name name of time variable in NetCDF file
15 * @param[in ] a_index_type nodality of data field
16 * @param[in ] a_var_need_data array over boundaries of flags that indicate whether we need data for the variable
17 * @param[in ] a_is2d Whether the variable we're working with is 2D
18 * @param[in ] a_rx refinement ratio in x relative to level 0
19 * @param[in ] a_ry refinement ratio in y relative to level 0
20 */
21NCTimeSeriesBoundary::NCTimeSeriesBoundary (int a_lev, const amrex::Vector<amrex::Geometry> a_geom,
22 const amrex::Vector<std::string>& a_file_names, const std::string a_field_name,
23 const std::string a_time_name,
24 const amrex::IntVect a_index_type,
25 const amrex::GpuArray<bool, AMREX_SPACEDIM*2>* a_var_need_data,
26 bool a_is2d, int a_rx, int a_ry)
27{
28 m_lev = a_lev;
29 m_geom = a_geom;
30 file_names.assign(a_file_names.begin(), a_file_names.end());
31 time_name = a_time_name;
32 field_name = a_field_name;
33 domain = a_geom[a_lev].Domain();
34 index_type = a_index_type;
35 var_need_data = *a_var_need_data;
36 is2d = a_is2d;
37 m_rx = a_rx;
38 m_ry = a_ry;
39}
40
41void NCTimeSeriesBoundary::Initialize()
42{
43 // Initialize Fabs
44 amrex::Arena* Arena_Used = amrex::The_Arena();
45#ifdef AMREX_USE_GPU
46 Arena_Used = amrex::The_Pinned_Arena();
47#endif
48
49 // open file
50 amrex::Print() << "Setting up boundary data for " << field_name << " coming from NetCDF file " << std::endl;
51
52 // The time field can have any number of names, depending on the field.
53 // If not specified in input file (time_name.empty()) then set it by default
54 if (time_name.empty())
55 {
56 time_name = "ocean_time";
57 }
58
59 for (int ifile = 0; ifile < file_names.size(); ifile++) {
60 std::string file_name = file_names[ifile];
61 // Check units of time stamps; should be days
62 std::string unit_str = ReadNetCDFVarAttrStr(file_name, time_name, "units"); // works on proc 0
63 if (amrex::ParallelDescriptor::IOProcessor())
64 {
65 if (unit_str.find("days") == std::string::npos) {
66 amrex::Print() << "Units of ocean_time given as: " << unit_str << std::endl;
67 amrex::Abort("Units must be in days.");
68 }
69 }
70 // get times and put in array
71 using RARRAY = NDArray<amrex::Real>;
72 amrex::Vector<RARRAY> array_ts(1);
73 ReadNetCDFFile(file_name, {time_name}, array_ts); // filled only on proc 0
74 if (amrex::ParallelDescriptor::IOProcessor())
75 {
76 int ntimes_io = array_ts[0].get_vshape()[0];
77 for (int nt(0); nt < ntimes_io; nt++)
78 {
79 // Convert ocean time from days to seconds
80 bry_times.push_back((*(array_ts[0].get_data() + nt)) * amrex::Real(60.0) * amrex::Real(60.0) * amrex::Real(24.0));
81 file_for_time.push_back(ifile);
82 file_itime_offset.push_back(nt);
83 // amrex::Print() << "TIMES " << bry_times[nt] << std::endl;
84 }
85 }
86 }
87
88 AMREX_ASSERT(std::is_sorted(bry_times.begin(), bry_times.end()));
89
90 int ntimes = bry_times.size();
91 int ioproc = amrex::ParallelDescriptor::IOProcessorNumber();
92 amrex::ParallelDescriptor::Bcast(&ntimes,1,ioproc);
93 if (!(amrex::ParallelDescriptor::IOProcessor())) {
94 bry_times.resize(ntimes);
95 file_for_time.resize(ntimes);
96 file_itime_offset.resize(ntimes);
97 }
98 amrex::ParallelDescriptor::Bcast(bry_times.data(), bry_times.size(), ioproc);
99 amrex::ParallelDescriptor::Bcast(file_for_time.data(), file_for_time.size(), ioproc);
100 amrex::ParallelDescriptor::Bcast(file_itime_offset.data(), file_itime_offset.size(), ioproc);
101
102 const auto& lo = domain.loVect();
103 const auto& hi = domain.hiVect();
104
105 amrex::Box xlo_bx(amrex::IntVect(lo[0]+index_type[0]-1, lo[1]+index_type[1]-1, lo[2]),
106 amrex::IntVect(lo[0]+index_type[0]-1, hi[1]+1 , hi[2]), index_type);
107 amrex::Box xhi_bx(amrex::IntVect(hi[0]+1 , lo[1]+index_type[1]-1, lo[2]),
108 amrex::IntVect(hi[0]+1 , hi[1]+1 , hi[2]), index_type);
109 amrex::Box ylo_bx(amrex::IntVect(lo[0]+index_type[0]-1, lo[1]+index_type[1]-1, lo[2]),
110 amrex::IntVect(hi[0]+1 , lo[1]+index_type[1]-1, hi[2]), index_type);
111 amrex::Box yhi_bx(amrex::IntVect(lo[0]+index_type[0]-1, hi[1]+1 , lo[2]),
112 amrex::IntVect(hi[0]+1 , hi[1]+1 , hi[2]), index_type);
113 if (is2d) {
114 xlo_bx.makeSlab(2,0);
115 xhi_bx.makeSlab(2,0);
116 ylo_bx.makeSlab(2,0);
117 yhi_bx.makeSlab(2,0);
118 }
119
120 // amrex::Print() << xlo_bx << " " << xhi_bx << " " << ylo_bx << " " << yhi_bx << std::endl;
121
122 xlo_dat_before = amrex::FArrayBox(xlo_bx, 1, Arena_Used);
123 xhi_dat_before = amrex::FArrayBox(xhi_bx, 1, Arena_Used);
124 ylo_dat_before = amrex::FArrayBox(ylo_bx, 1, Arena_Used);
125 yhi_dat_before = amrex::FArrayBox(yhi_bx, 1, Arena_Used);
126
127 xlo_dat_after = amrex::FArrayBox(xlo_bx, 1, Arena_Used);
128 xhi_dat_after = amrex::FArrayBox(xhi_bx, 1, Arena_Used);
129 ylo_dat_after = amrex::FArrayBox(ylo_bx, 1, Arena_Used);
130 yhi_dat_after = amrex::FArrayBox(yhi_bx, 1, Arena_Used);
131
132 xlo_dat_interp = amrex::FArrayBox(xlo_bx, 1, Arena_Used);
133 xhi_dat_interp = amrex::FArrayBox(xhi_bx, 1, Arena_Used);
134 ylo_dat_interp = amrex::FArrayBox(ylo_bx, 1, Arena_Used);
135 yhi_dat_interp = amrex::FArrayBox(yhi_bx, 1, Arena_Used);
136
137 // dummy initialization
138 i_time_before = -100;
139
140 if (var_need_data[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)] == true) {
141 nc_var_names.push_back(field_name + "_west");
142 }
143 if (var_need_data[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)] == true) {
144 nc_var_names.push_back(field_name + "_east");
145 }
146 if (var_need_data[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)] == true) {
147 nc_var_names.push_back(field_name + "_south");
148 }
149 if (var_need_data[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)] == true) {
150 nc_var_names.push_back(field_name + "_north");
151 }
152}
153
154/**
155 * @param time time to interpolate to
156 */
157void NCTimeSeriesBoundary::update_interpolated_to_time (amrex::Real time)
158{
159 // Initialize Fabs
160 amrex::Arena* Arena_Used = amrex::The_Arena();
161#ifdef AMREX_USE_GPU
162 Arena_Used = amrex::The_Pinned_Arena();
163#endif
164
165 // Figure out time index:
166 AMREX_ASSERT(time >= bry_times[0]);
167 AMREX_ASSERT(time <= bry_times[bry_times.size()-1]);
168 int i_time_before_old = i_time_before;
169 for (int nt=0; nt < bry_times.size()-1; nt++) {
170 if ((bry_times[nt] <= time) and (bry_times[nt+1] >= time)) {
171 i_time_before = nt;
172 time_before = bry_times[nt];
173 time_after = bry_times[nt+1];
174 break;
175 }
176 }
177
178 int i_time_after = i_time_before + 1;
179
180 amrex::FArrayBox crse_xlo_dat;
181 amrex::FArrayBox crse_xhi_dat;
182 amrex::FArrayBox crse_ylo_dat;
183 amrex::FArrayBox crse_yhi_dat;
184
185 if (m_lev > 0) {
186 const auto& crse_lo = m_geom[0].Domain().loVect();
187 const auto& crse_hi = m_geom[0].Domain().hiVect();
188
189 amrex::Box crse_xlo_bx(amrex::IntVect(crse_lo[0]+index_type[0]-1, crse_lo[1]+index_type[1]-1, crse_lo[2]),
190 amrex::IntVect(crse_lo[0]+index_type[0]-1, crse_hi[1]+1 , crse_hi[2]), index_type);
191 amrex::Box crse_xhi_bx(amrex::IntVect(crse_hi[0]+1 , crse_lo[1]+index_type[1]-1, crse_lo[2]),
192 amrex::IntVect(crse_hi[0]+1 , crse_hi[1]+1 , crse_hi[2]), index_type);
193 amrex::Box crse_ylo_bx(amrex::IntVect(crse_lo[0]+index_type[0]-1, crse_lo[1]+index_type[1]-1, crse_lo[2]),
194 amrex::IntVect(crse_hi[0]+1 , crse_lo[1]+index_type[1]-1, crse_hi[2]), index_type);
195 amrex::Box crse_yhi_bx(amrex::IntVect(crse_lo[0]+index_type[0]-1, crse_hi[1]+1 , crse_lo[2]),
196 amrex::IntVect(crse_hi[0]+1 , crse_hi[1]+1 , crse_hi[2]), index_type);
197 if (is2d) {
198 crse_xlo_bx.makeSlab(2,0);
199 crse_xhi_bx.makeSlab(2,0);
200 crse_ylo_bx.makeSlab(2,0);
201 crse_yhi_bx.makeSlab(2,0);
202 }
203
204 crse_xlo_dat = amrex::FArrayBox(crse_xlo_bx, 1, Arena_Used);
205 crse_xhi_dat = amrex::FArrayBox(crse_xhi_bx, 1, Arena_Used);
206 crse_ylo_dat = amrex::FArrayBox(crse_ylo_bx, 1, Arena_Used);
207 crse_yhi_dat = amrex::FArrayBox(crse_yhi_bx, 1, Arena_Used);
208 }
209
210 if (i_time_before_old + 1 == i_time_before) {
211 // swap multifabs so we only have to read in one MultiFab
212 std::swap(xlo_dat_before, xlo_dat_after);
213 std::swap(xhi_dat_before, xhi_dat_after);
214 std::swap(ylo_dat_before, ylo_dat_after);
215 std::swap(yhi_dat_before, yhi_dat_after);
216
217 amrex::Print() << "Reading in " << field_name << " at (after) time " << i_time_after << std::endl;
218
219 if (m_lev == 0) {
220 read_in_at_time(xlo_dat_after, xhi_dat_after, ylo_dat_after, yhi_dat_after, i_time_after);
221 } else {
222
223 read_in_at_time(crse_xlo_dat, crse_xhi_dat, crse_ylo_dat, crse_yhi_dat, i_time_after);
224
225 interp_fab(crse_xlo_dat, xlo_dat_after);
226 interp_fab(crse_xhi_dat, xhi_dat_after);
227 interp_fab(crse_ylo_dat, ylo_dat_after);
228 interp_fab(crse_yhi_dat, yhi_dat_after);
229
230 }
231
232 } else if (i_time_before_old != i_time_before) {
233
234
235 if (m_lev == 0) {
236 amrex::Print() << "Reading in " << field_name << " at (before) time " << i_time_before << std::endl;
237 read_in_at_time(xlo_dat_before,xhi_dat_before, ylo_dat_before,yhi_dat_before, i_time_before);
238 amrex::Print() << "Reading in " << field_name << " at (after) time " << i_time_after << std::endl;
239 read_in_at_time(xlo_dat_after ,xhi_dat_after , ylo_dat_after , yhi_dat_after, i_time_after);
240 } else {
241
242 amrex::Print() << "Reading in " << field_name << " at (before) time " << i_time_before << std::endl;
243
244 read_in_at_time(crse_xlo_dat, crse_xhi_dat, crse_ylo_dat,crse_yhi_dat, i_time_before);
245
246 interp_fab(crse_xlo_dat, xlo_dat_before);
247 interp_fab(crse_xhi_dat, xhi_dat_before);
248 interp_fab(crse_ylo_dat, ylo_dat_before);
249 interp_fab(crse_yhi_dat, yhi_dat_before);
250
251 amrex::Print() << "Reading in " << field_name << " at time " << i_time_after << std::endl;
252
253 read_in_at_time(crse_xlo_dat, crse_xhi_dat, crse_ylo_dat, crse_yhi_dat, i_time_after);
254
255 interp_fab(crse_xlo_dat, xlo_dat_after);
256 interp_fab(crse_xhi_dat, xhi_dat_after);
257 interp_fab(crse_ylo_dat, ylo_dat_after);
258 interp_fab(crse_yhi_dat, yhi_dat_after);
259 } // lev
260 } // i_time
261
262 amrex::Real dt = time_after - time_before;
263 amrex::Real time_before_copy = time_before;
264
265 amrex::Array4<amrex::Real> xlo_interp_arr = xlo_dat_interp.array();
266 amrex::Array4<amrex::Real> xhi_interp_arr = xhi_dat_interp.array();
267 amrex::Array4<amrex::Real> ylo_interp_arr = ylo_dat_interp.array();
268 amrex::Array4<amrex::Real> yhi_interp_arr = yhi_dat_interp.array();
269
270 amrex::Array4<amrex::Real> xlo_before_arr = xlo_dat_before.array();
271 amrex::Array4<amrex::Real> xhi_before_arr = xhi_dat_before.array();
272 amrex::Array4<amrex::Real> ylo_before_arr = ylo_dat_before.array();
273 amrex::Array4<amrex::Real> yhi_before_arr = yhi_dat_before.array();
274
275 amrex::Array4<amrex::Real> xlo_after_arr = xlo_dat_after.array();
276 amrex::Array4<amrex::Real> xhi_after_arr = xhi_dat_after.array();
277 amrex::Array4<amrex::Real> ylo_after_arr = ylo_dat_after.array();
278 amrex::Array4<amrex::Real> yhi_after_arr = yhi_dat_after.array();
279
280 if (var_need_data[amrex::Orientation(amrex::Direction::x,amrex::Orientation::low)] == true) {
281 amrex::ParallelFor(xlo_dat_interp.box(), [=] AMREX_GPU_DEVICE (int i, int j, int k)
282 {
283 xlo_interp_arr(i,j,k) = xlo_before_arr(i,j,k) + (time - time_before_copy) * (xlo_after_arr(i,j,k) - xlo_before_arr(i,j,k)) / dt;
284 });
285 }
286
287 if (var_need_data[amrex::Orientation(amrex::Direction::x,amrex::Orientation::high)] == true) {
288 amrex::ParallelFor(xhi_dat_interp.box(), [=] AMREX_GPU_DEVICE (int i, int j, int k)
289 {
290 xhi_interp_arr(i,j,k) = xhi_before_arr(i,j,k) + (time - time_before_copy) * (xhi_after_arr(i,j,k) - xhi_before_arr(i,j,k)) / dt;
291 });
292 }
293
294 if (var_need_data[amrex::Orientation(amrex::Direction::y,amrex::Orientation::low)] == true) {
295 amrex::ParallelFor(ylo_dat_interp.box(), [=] AMREX_GPU_DEVICE (int i, int j, int k)
296 {
297 ylo_interp_arr(i,j,k) = ylo_before_arr(i,j,k) + (time - time_before_copy) * (ylo_after_arr(i,j,k) - ylo_before_arr(i,j,k)) / dt;
298 });
299 }
300
301 if (var_need_data[amrex::Orientation(amrex::Direction::y,amrex::Orientation::high)] == true) {
302 amrex::ParallelFor(yhi_dat_interp.box(), [=] AMREX_GPU_DEVICE (int i, int j, int k)
303 {
304 yhi_interp_arr(i,j,k) = yhi_before_arr(i,j,k) + (time - time_before_copy) * (yhi_after_arr(i,j,k) - yhi_before_arr(i,j,k)) / dt;
305 });
306 }
307}
308
309/**
310 * @param[inout] fab_xlo fab to store xlo boundary data into
311 * @param[inout] fab_xhi fab to store xhi boundary data into
312 * @param[inout] fab_ylo fab to store ylo boundary data into
313 * @param[inout] fab_yhi fab to store yhi boundary data into
314 * @param[in ] itime index of time step to read from file
315 */
316void NCTimeSeriesBoundary::read_in_at_time (amrex::FArrayBox& fab_xlo,
317 amrex::FArrayBox& fab_xhi,
318 amrex::FArrayBox& fab_ylo,
319 amrex::FArrayBox& fab_yhi,
320 int itime) {
321 using RARRAY = NDArray<amrex::Real>;
322 amrex::Vector<RARRAY> arrays(nc_var_names.size());
323
324 amrex::Print() << "Actually reading in " << field_name << " at time " << bry_times[itime] << std::endl;
325 std::string nc_bdry_file = file_names[file_for_time[itime]];
326 int itime_offset = file_itime_offset[itime];
327 ReadNetCDFFile(nc_bdry_file, nc_var_names, arrays, true, itime_offset); // does work on proc 0 only
328
329 for (int iv=0; iv < nc_var_names.size(); iv++) {
330 if (amrex::ParallelDescriptor::IOProcessor())
331 {
332 std::string last4 = nc_var_names[iv].substr(nc_var_names[iv].size()-4, 4);
333 std::string last5 = nc_var_names[iv].substr(nc_var_names[iv].size()-5, 5);
334 int nx, ny, nz, n_plane;
335 int i, j, k, ioff, joff;
336 if (last4 == "west") {
337 amrex::Box my_box = fab_xlo.box();
338 if (is2d) {
339 nz = 1;
340 ny = arrays[iv].get_vshape()[1];
341 } else {
342 nz = arrays[iv].get_vshape()[1];
343 ny = arrays[iv].get_vshape()[2];
344 }
345 n_plane = ny * nz;
346 AMREX_ALWAYS_ASSERT(my_box.numPts() == n_plane);
347
348 i = my_box.smallEnd()[0];
349 joff = my_box.smallEnd()[1];
350
351 amrex::Array4<amrex::Real> fab_arr = fab_xlo.array();
352
353 for (int n(0); n < n_plane; n++) {
354 k = n / ny;
355 j = n - (k * ny);
356 fab_arr(i, j+joff, k, 0) = static_cast<amrex::Real>(*(arrays[iv].get_data() + n));
357 }
358 } else if (last4 == "east") {
359 amrex::Box my_box = fab_xhi.box();
360 if (is2d) {
361 nz = 1;
362 ny = arrays[iv].get_vshape()[1];
363 } else {
364 nz = arrays[iv].get_vshape()[1];
365 ny = arrays[iv].get_vshape()[2];
366 }
367 n_plane = ny * nz;
368 AMREX_ALWAYS_ASSERT(my_box.numPts() == n_plane);
369
370 i = my_box.smallEnd()[0];
371 joff = my_box.smallEnd()[1];
372
373 amrex::Array4<amrex::Real> fab_arr = fab_xhi.array();
374 for (int n(0); n < n_plane; n++) {
375 k = n / ny;
376 j = n - (k * ny);
377 fab_arr(i, j+joff, k, 0) = static_cast<amrex::Real>(*(arrays[iv].get_data() + n));
378 }
379 } else if (last5 == "south") {
380 amrex::Box my_box = fab_ylo.box();
381 if (is2d) {
382 nz = 1;
383 nx = arrays[iv].get_vshape()[1];
384 } else {
385 nz = arrays[iv].get_vshape()[1];
386 nx = arrays[iv].get_vshape()[2];
387 }
388 n_plane = nx * nz;
389 AMREX_ALWAYS_ASSERT(my_box.numPts() == n_plane);
390
391 j = my_box.smallEnd()[1];
392 ioff = my_box.smallEnd()[0];
393
394 amrex::Array4<amrex::Real> fab_arr = fab_ylo.array();
395 for (int n(0); n < n_plane; n++) {
396 k = n / nx;
397 i = n - (k * nx);
398 fab_arr(i+ioff, j, k, 0) = static_cast<amrex::Real>(*(arrays[iv].get_data() + n));
399 }
400 } else if (last5 == "north") {
401 amrex::Box my_box = fab_yhi.box();
402 if (is2d) {
403 nz = 1;
404 nx = arrays[iv].get_vshape()[1];
405 } else {
406 nz = arrays[iv].get_vshape()[1];
407 nx = arrays[iv].get_vshape()[2];
408 }
409 n_plane = nx * nz;
410 AMREX_ALWAYS_ASSERT(my_box.numPts() == n_plane);
411
412 j = my_box.smallEnd()[1];
413 ioff = my_box.smallEnd()[0];
414
415 amrex::Array4<amrex::Real> fab_arr = fab_yhi.array();
416 for (int n(0); n < n_plane; n++) {
417 k = n / nx;
418 i = n - (k * nx);
419 fab_arr(i+ioff, j, k, 0) = static_cast<amrex::Real>(*(arrays[iv].get_data() + n));
420 }
421 }
422 }
423 }
424
425 amrex::ParallelDescriptor::Barrier();
426
427 // When an FArrayBox is built, space is allocated on every rank. However, we only
428 // filled the data in these FABs on the IOProcessor. So here we broadcast
429 // the data to every rank.
430 int ioproc = amrex::ParallelDescriptor::IOProcessorNumber(); // I/O rank
431 amrex::ParallelDescriptor::Bcast(fab_xlo.dataPtr(),fab_xlo.box().numPts(),ioproc);
432 amrex::ParallelDescriptor::Bcast(fab_xhi.dataPtr(),fab_xhi.box().numPts(),ioproc);
433 amrex::ParallelDescriptor::Bcast(fab_ylo.dataPtr(),fab_ylo.box().numPts(),ioproc);
434 amrex::ParallelDescriptor::Bcast(fab_yhi.dataPtr(),fab_yhi.box().numPts(),ioproc);
435
436 amrex::Print() << "DONE reading in " << field_name << " at time " << bry_times[itime] << std::endl;
437}
438
439/**
440 * @param[in ] dat_crse fab of coarse data to interpolate from
441 * @param[ out] dat_fine fab of fine data to interpoalte to
442 */
443void NCTimeSeriesBoundary::interp_fab(amrex::FArrayBox& dat_crse, amrex::FArrayBox& dat_fine)//, int rx, int ry)
444{
445 amrex::Array4<amrex::Real> crse_arr = dat_crse.array();
446 amrex::Array4<amrex::Real> fine_arr = dat_fine.array();
447
448 const auto& bhi = ubound(dat_crse.box());
449
450 amrex::Real xfac = 1.0 / static_cast<amrex::Real>(m_rx);
451 amrex::Real yfac = 1.0 / static_cast<amrex::Real>(m_ry);
452
453 amrex::Real rx = m_rx;
454 amrex::Real ry = m_ry;
455
456 // Doing box on x-face
457 if (dat_crse.box().length(0) == 1) {
458 // amrex::Print() << "DOING INTERP ON XFACE " << dat_crse.box() << " " << dat_fine.box() << std::endl;
459 if (dat_crse.box().ixType()[1] == 0) {
460 amrex::ParallelFor(dat_crse.box(), [=] AMREX_GPU_DEVICE (int i, int j, int k)
461 {
462 int i_f = (i == -1) ? -1 : rx*i;
463 if (j == -1) {
464 fine_arr(i_f,j,k) = crse_arr(i,j,k);
465 } else {
466 fine_arr(i_f,ry*j ,k) = crse_arr(i,j,k);
467 if (j < bhi.y) {
468 for (int n = 1; n < ry; n++) {
469 fine_arr(i_f,ry*j+n,k) = crse_arr(i,j,k) + n * yfac * (crse_arr(i,j+1,k) - crse_arr(i,j,k));
470 }
471 }
472 }
473 });
474 } else {
475 amrex::ParallelFor(dat_crse.box(), [=] AMREX_GPU_DEVICE (int i, int j, int k)
476 {
477 int i_f = (i == -1) ? -1 : rx*i;
478 if (j == -1) {
479 fine_arr(i_f,j,k) = crse_arr(i,j,k);
480 } else {
481 fine_arr(i_f,ry*j,k) = crse_arr(i,j,k);
482 if (j < bhi.y) {
483 for (int n = 1; n < ry; n++) {
484 fine_arr(i_f,ry*j+n,k) = crse_arr(i,j,k) + n * yfac * (crse_arr(i,j+1,k) - crse_arr(i,j,k));
485 }
486 }
487 }
488 });
489 }
490 } else if (dat_crse.box().length(1) == 1) {
491 // Doing box on y-face
492 // amrex::Print() << "DOING INTERP ON YFACE " << dat_crse.box() << " " << dat_fine.box() << std::endl;
493 if (dat_crse.box().ixType()[0] == 0) {
494 amrex::ParallelFor(dat_crse.box(), [=] AMREX_GPU_DEVICE (int i, int j, int k)
495 {
496 int j_f = (j == -1) ? -1 : ry*j;
497 if (i == -1) {
498 fine_arr(i,j_f,k) = crse_arr(i,j,k);
499 } else {
500 fine_arr(rx*i ,j_f,k) = crse_arr(i,j,k);
501 if (i < bhi.x) {
502 for (int n = 1; n < rx; n++) {
503 fine_arr(rx*i+n,j_f,k) = crse_arr(i,j,k) + n * xfac * (crse_arr(i+1,j,k) - crse_arr(i,j,k));
504 }
505 }
506 }
507 });
508 } else {
509 amrex::ParallelFor(dat_crse.box(), [=] AMREX_GPU_DEVICE (int i, int j, int k)
510 {
511 int j_f = (j == -1) ? -1 : ry*j;
512 if (i == -1) {
513 fine_arr(i,j_f,k) = crse_arr(i,j,k);
514 } else {
515 fine_arr(rx*i,j_f,k) = crse_arr(i,j,k);
516 if (i < bhi.x) {
517 for (int n = 1; n < rx; n++) {
518 fine_arr(rx*i+n,j_f,k) = crse_arr(i,j,k) + n * xfac * (crse_arr(i+1,j,k) - crse_arr(i,j,k));
519 }
520 }
521 }
522 });
523 }
524 } else {
525 amrex::Abort(" What am I doing here??");
526 }
527}
528
529#endif // REMORA_USE_NETCDF
REMORA_NCFile.H
ReadNetCDFFile
void ReadNetCDFFile(const std::string &fname, amrex::Vector< std::string > names, amrex::Vector< NDArray< DType > > &arrays, bool one_time=false, int fill_time=0)
Read in data from netcdf file and save to data arrays.
Definition REMORA_NCFile.H:119
ReadNetCDFVarAttrStr
std::string ReadNetCDFVarAttrStr(const std::string &fname, const std::string &var_name, const std::string &attr_name)
Helper function for reading a single variable attribute.
Definition REMORA_NCFile.cpp:13
REMORA_NCTimeSeriesBoundary.H
NCTimeSeriesBoundary::yhi_dat_before
amrex::FArrayBox yhi_dat_before
FArrayBox to store data at time_before at yhi boundary.
Definition REMORA_NCTimeSeriesBoundary.H:114
NCTimeSeriesBoundary::yhi_dat_after
amrex::FArrayBox yhi_dat_after
FArrayBox to store data at time_after at yhi boundary.
Definition REMORA_NCTimeSeriesBoundary.H:116
NCTimeSeriesBoundary::var_need_data
amrex::GpuArray< bool, AMREX_SPACEDIM *2 > var_need_data
Array over boundaries indicating whether we need physical data for this variable.
Definition REMORA_NCTimeSeriesBoundary.H:119
NCTimeSeriesBoundary::interp_fab
void interp_fab(amrex::FArrayBox &dat_crse, amrex::FArrayBox &dat_fine)
spatially interpolate boundary fabs from level 0 data (read in) to current level
Definition REMORA_NCTimeSeriesBoundary.cpp:443
NCTimeSeriesBoundary::ylo_dat_before
amrex::FArrayBox ylo_dat_before
FArrayBox to store data at time_before at ylo boundary.
Definition REMORA_NCTimeSeriesBoundary.H:104
NCTimeSeriesBoundary::xlo_dat_after
amrex::FArrayBox xlo_dat_after
FArrayBox to store data at time_after at xlo boundary.
Definition REMORA_NCTimeSeriesBoundary.H:101
NCTimeSeriesBoundary::nc_var_names
amrex::Vector< std::string > nc_var_names
Variable names that will be read from file.
Definition REMORA_NCTimeSeriesBoundary.H:68
NCTimeSeriesBoundary::ylo_dat_after
amrex::FArrayBox ylo_dat_after
FArrayBox to store data at time_after at ylo boundary.
Definition REMORA_NCTimeSeriesBoundary.H:106
NCTimeSeriesBoundary::time_before
amrex::Real time_before
Time in ocean_times immediately before the last time interpolated to.
Definition REMORA_NCTimeSeriesBoundary.H:94
NCTimeSeriesBoundary::Initialize
void Initialize()
Read in time array from file and allocate data arrays.
Definition REMORA_NCTimeSeriesBoundary.cpp:41
NCTimeSeriesBoundary::yhi_dat_interp
amrex::FArrayBox yhi_dat_interp
FArrayBox to store data at inteprolated time at yhi boundary *‍/.
Definition REMORA_NCTimeSeriesBoundary.H:46
NCTimeSeriesBoundary::xlo_dat_interp
amrex::FArrayBox xlo_dat_interp
FArrayBox to store data at inteprolated time at xlo boundary *‍/.
Definition REMORA_NCTimeSeriesBoundary.H:37
NCTimeSeriesBoundary::bry_times
amrex::Vector< amrex::Real > bry_times
Time points in netcdf file.
Definition REMORA_NCTimeSeriesBoundary.H:86
NCTimeSeriesBoundary::time_name
std::string time_name
Field name for time series in netcdf file.
Definition REMORA_NCTimeSeriesBoundary.H:70
NCTimeSeriesBoundary::xhi_dat_interp
amrex::FArrayBox xhi_dat_interp
FArrayBox to store data at inteprolated time at xhi boundary *‍/.
Definition REMORA_NCTimeSeriesBoundary.H:43
NCTimeSeriesBoundary::index_type
amrex::IntVect index_type
Index type for field to fill.
Definition REMORA_NCTimeSeriesBoundary.H:74
NCTimeSeriesBoundary::i_time_before
int i_time_before
Time index immediately before the last time interpolated to.
Definition REMORA_NCTimeSeriesBoundary.H:92
NCTimeSeriesBoundary::field_name
std::string field_name
Field name in netcdf file.
Definition REMORA_NCTimeSeriesBoundary.H:66
NCTimeSeriesBoundary::ylo_dat_interp
amrex::FArrayBox ylo_dat_interp
FArrayBox to store data at inteprolated time at ylo boundary *‍/.
Definition REMORA_NCTimeSeriesBoundary.H:40
NCTimeSeriesBoundary::m_geom
amrex::Vector< amrex::Geometry > m_geom
Geometry at all levels.
Definition REMORA_NCTimeSeriesBoundary.H:80
NCTimeSeriesBoundary::xhi_dat_after
amrex::FArrayBox xhi_dat_after
FArrayBox to store data at time_after at xhi boundary.
Definition REMORA_NCTimeSeriesBoundary.H:111
NCTimeSeriesBoundary::read_in_at_time
void read_in_at_time(amrex::FArrayBox &fab_xlo, amrex::FArrayBox &fab_xhi, amrex::FArrayBox &fab_ylo, amrex::FArrayBox &fab_yhi, int itime)
Read in data from file at time index itime and fill into mf.
Definition REMORA_NCTimeSeriesBoundary.cpp:316
NCTimeSeriesBoundary::file_names
amrex::Vector< std::string > file_names
File name to read from.
Definition REMORA_NCTimeSeriesBoundary.H:64
NCTimeSeriesBoundary::file_itime_offset
amrex::Vector< int > file_itime_offset
Offset to access a particular time within its file.
Definition REMORA_NCTimeSeriesBoundary.H:90
NCTimeSeriesBoundary::xhi_dat_before
amrex::FArrayBox xhi_dat_before
FArrayBox to store data at time_before at xhi boundary.
Definition REMORA_NCTimeSeriesBoundary.H:109
NCTimeSeriesBoundary::is2d
bool is2d
Whether the field we're reading in is 2d.
Definition REMORA_NCTimeSeriesBoundary.H:83
NCTimeSeriesBoundary::update_interpolated_to_time
void update_interpolated_to_time(amrex::Real time)
Calculate interpolated values at time, reading in data as necessary.
Definition REMORA_NCTimeSeriesBoundary.cpp:157
NCTimeSeriesBoundary::m_lev
int m_lev
Level at which we are holding the boundary data.
Definition REMORA_NCTimeSeriesBoundary.H:77
NCTimeSeriesBoundary::m_rx
int m_rx
Refinement ratios relative to level 0.
Definition REMORA_NCTimeSeriesBoundary.H:49
NCTimeSeriesBoundary::NCTimeSeriesBoundary
NCTimeSeriesBoundary(int a_lev, const amrex::Vector< amrex::Geometry > a_geom, const amrex::Vector< std::string > &a_file_name, const std::string a_field_name, const std::string a_time_name, const amrex::IntVect a_index_type, const amrex::GpuArray< bool, AMREX_SPACEDIM *2 > *a_var_need_data, bool a_is2d, int rx, int ry)
Constructor.
Definition REMORA_NCTimeSeriesBoundary.cpp:21
NCTimeSeriesBoundary::xlo_dat_before
amrex::FArrayBox xlo_dat_before
FArrayBox to store data at time_before at xlo boundary.
Definition REMORA_NCTimeSeriesBoundary.H:99
NCTimeSeriesBoundary::time_after
amrex::Real time_after
Time in ocean_times immediately after the last time interpolated to.
Definition REMORA_NCTimeSeriesBoundary.H:96
NCTimeSeriesBoundary::file_for_time
amrex::Vector< int > file_for_time
File index to access a particular time.
Definition REMORA_NCTimeSeriesBoundary.H:88
NDArray
NDArray is the datatype designed to hold any data, including scalars, multidimensional arrays,...
Definition REMORA_NCFile.H:36