REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_init_from_netcdf.cpp
Go to the documentation of this file.
1/**
2 * \file REMORA_init_from_netcdf.cpp
3 */
4
5#include <REMORA.H>
6#include <REMORA_Constants.H>
8#include <REMORA_DataStruct.H>
9
10using namespace amrex;
11
12#ifdef REMORA_USE_NETCDF
13
14/** \brief helper function for reading in initial state data from netcdf */
15void
16read_data_from_netcdf (int /*lev*/, const Box& domain, const std::string& fname,
17 FArrayBox& NC_temp_fab, FArrayBox& NC_salt_fab,
18 FArrayBox& NC_xvel_fab, FArrayBox& NC_yvel_fab,
19 FArrayBox& NC_ubar_fab, FArrayBox& NC_vbar_fab);
20
21/** \brief helper function for reading in land-sea masks from netcdf */
22void
23read_masks_from_netcdf (int /*lev*/, const Box& domain, const std::string& fname,
24 FArrayBox& NC_mskr_fab, FArrayBox& NC_msku_fab,
25 FArrayBox& NC_mskv_fab);
26
27/** \brief helper function for reading boundary data from netcdf */
28Real
29read_bdry_from_netcdf (const Box& domain, const std::string& fname,
30 Vector<Vector<FArrayBox>>& bdy_data_xlo,
31 Vector<Vector<FArrayBox>>& bdy_data_xhi,
32 Vector<Vector<FArrayBox>>& bdy_data_ylo,
33 Vector<Vector<FArrayBox>>& bdy_data_yhi,
34 int& width, amrex::Real& start_bdy_time,
35 std::string bdry_time_varname,
36 amrex::GpuArray<amrex::GpuArray<bool, AMREX_SPACEDIM*2>,BdyVars::NumTypes+1>&);
37
38/** \brief helper function to initialize state from netcdf */
39void
41 FArrayBox& temp_fab, FArrayBox& salt_fab,
42 FArrayBox& x_vel_fab, FArrayBox& y_vel_fab,
43 FArrayBox& ubar_fab, FArrayBox& vbar_fab,
44 const Vector<FArrayBox>& NC_temp_fab,
45 const Vector<FArrayBox>& NC_salt_fab,
46 const Vector<FArrayBox>& NC_xvel_fab,
47 const Vector<FArrayBox>& NC_yvel_fab,
48 const Vector<FArrayBox>& NC_ubar_fab,
49 const Vector<FArrayBox>& NC_vbar_fab);
50
51/** \brief helper function to read bathymetry from netcdf */
52void
53read_bathymetry_from_netcdf (int lev, const Box& domain, const std::string& fname,
54 FArrayBox& NC_h_fab,
55 FArrayBox& NC_pm_fab, FArrayBox& NC_pn_fab,
56 FArrayBox& NC_xr_fab, FArrayBox& NC_yr_fab,
57 FArrayBox& NC_xu_fab, FArrayBox& NC_yu_fab,
58 FArrayBox& NC_xv_fab, FArrayBox& NC_yv_fab,
59 FArrayBox& NC_xp_fab, FArrayBox& NC_yp_fab);
60
61/** \brief helper function to read coriolis factor from netcdf */
62void
63read_coriolis_from_netcdf (int lev, const Box& domain, const std::string& fname, FArrayBox& NC_fcor_fab);
64
65/** \brief helper function to read sea surface height from netcdf */
66void
67read_zeta_from_netcdf (int lev, const Box& domain, const std::string& fname,
68 FArrayBox& NC_zeta_fab);
69
70/** \brief helper function to read climatology nudging from netcdf */
71void
72read_clim_nudg_coeff_from_netcdf (int lev, const Box& domain, const std::string& fname,
73 bool do_m2_clim_nudg,
74 bool do_m3_clim_nudg,
75 bool do_temp_clim_nudg,
76 bool do_salt_clim_nudg,
77 FArrayBox& NC_M2NC_fab,
78 FArrayBox& NC_M3NC_fab,
79 FArrayBox& NC_TempNC_fab,
80 FArrayBox& NC_SaltNC_fab);
81
82/** \brief helper function to read in vector of data from netcdf */
83//template <typename DType>
84void read_vec_from_netcdf (int lev, const std::string& fname, const std::string& field_name, amrex::Vector<int>& vec_dat);
85
86/**
87 * @param lev Integer specifying the current level
88 */
89void
91{
92 // *** FArrayBox's at this level for holding the INITIAL data
93 Vector<FArrayBox> NC_temp_fab ; NC_temp_fab.resize(num_boxes_at_level[lev]);
94 Vector<FArrayBox> NC_salt_fab ; NC_salt_fab.resize(num_boxes_at_level[lev]);
95 Vector<FArrayBox> NC_xvel_fab ; NC_xvel_fab.resize(num_boxes_at_level[lev]);
96 Vector<FArrayBox> NC_yvel_fab ; NC_yvel_fab.resize(num_boxes_at_level[lev]);
97 Vector<FArrayBox> NC_ubar_fab ; NC_ubar_fab.resize(num_boxes_at_level[lev]);
98 Vector<FArrayBox> NC_vbar_fab ; NC_vbar_fab.resize(num_boxes_at_level[lev]);
99
100 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
101 {
102 read_data_from_netcdf(lev, boxes_at_level[lev][idx], nc_init_file[lev][idx],
103 NC_temp_fab[idx], NC_salt_fab[idx],
104 NC_xvel_fab[idx], NC_yvel_fab[idx],
105 NC_ubar_fab[idx], NC_vbar_fab[idx]);
106 }
107
108
109 MultiFab mf_temp(*cons_new[lev], make_alias, Temp_comp, 1);
110 MultiFab mf_salt(*cons_new[lev], make_alias, Salt_comp, 1);
111
112#ifdef _OPENMP
113#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
114#endif
115 {
116 // Don't tile this since we are operating on full FABs in this routine
117 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
118 {
119 // Define fabs for holding the initial data
120 FArrayBox &temp_fab = mf_temp[mfi];
121 FArrayBox &salt_fab = mf_salt[mfi];
122 FArrayBox &xvel_fab = (*xvel_new[lev])[mfi];
123 FArrayBox &yvel_fab = (*yvel_new[lev])[mfi];
124 FArrayBox &ubar_fab = (*vec_ubar[lev])[mfi];
125 FArrayBox &vbar_fab = (*vec_vbar[lev])[mfi];
126
127 init_state_from_netcdf(lev, temp_fab, salt_fab,
128 xvel_fab, yvel_fab,
129 ubar_fab, vbar_fab,
130 NC_temp_fab, NC_salt_fab,
131 NC_xvel_fab, NC_yvel_fab,
132 NC_ubar_fab, NC_vbar_fab);
133 } // mf
134 } // omp
135}
136
137/**
138 * @param lev Integer specifying the current level
139 */
140void
142{
143 // *** FArrayBox's at this level for holding the INITIAL data
144 Vector<FArrayBox> NC_zeta_fab ; NC_zeta_fab.resize(num_boxes_at_level[lev]);
145
146 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
147 {
148 read_zeta_from_netcdf(lev,boxes_at_level[lev][idx], nc_init_file[lev][idx],
149 NC_zeta_fab[idx]);
150
151#ifdef _OPENMP
152#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
153#endif
154 {
155 // Don't tile this since we are operating on full FABs in this routine
156 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
157 {
158 FArrayBox &zeta_fab = (*vec_zeta[lev])[mfi];
159
160 //
161 // FArrayBox to FArrayBox copy does "copy on intersection"
162 // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks
163 //
164
165 zeta_fab.template copy<RunOn::Device>(NC_zeta_fab[idx],0,0,1);
166 } // mf
167 } // omp
168 } // idx
169 vec_zeta[lev]->FillBoundary(geom[lev].periodicity());
170 (*physbcs[lev])(*vec_zeta[lev],*vec_mskr[lev].get(),0,1,vec_zeta[lev]->nGrowVect(),t_new[lev],BCVars::zeta_bc);
171// (*physbcs[lev])(*vec_zeta[lev],*vec_mskr[lev].get(),1,1,vec_zeta[lev]->nGrowVect(),t_new[lev],BCVars::zeta_bc);
172// (*physbcs[lev])(*vec_zeta[lev],*vec_mskr[lev].get(),2,1,vec_zeta[lev]->nGrowVect(),t_new[lev],BCVars::zeta_bc);
173
174 Real told = t_new[lev];
176// fill_from_bdyfiles(*vec_zeta[lev], *vec_mskr[lev], told, BCVars::zeta_bc,BdyVars::zeta,1,1);
177// fill_from_bdyfiles(*vec_zeta[lev], *vec_mskr[lev], told, BCVars::zeta_bc,BdyVars::zeta,2,2);
178}
179/**
180 * @param lev Integer specifying the current level
181 */
182void
184{
185 // *** FArrayBox's at this level for holding the INITIAL data
186 Vector<FArrayBox> NC_h_fab ; NC_h_fab.resize(num_boxes_at_level[lev]);
187 Vector<FArrayBox> NC_pm_fab ; NC_pm_fab.resize(num_boxes_at_level[lev]);
188 Vector<FArrayBox> NC_pn_fab ; NC_pn_fab.resize(num_boxes_at_level[lev]);
189
190 Vector<FArrayBox> NC_xr_fab ; NC_xr_fab.resize(num_boxes_at_level[lev]);
191 Vector<FArrayBox> NC_yr_fab ; NC_yr_fab.resize(num_boxes_at_level[lev]);
192 Vector<FArrayBox> NC_xu_fab ; NC_xu_fab.resize(num_boxes_at_level[lev]);
193 Vector<FArrayBox> NC_yu_fab ; NC_yu_fab.resize(num_boxes_at_level[lev]);
194 Vector<FArrayBox> NC_xv_fab ; NC_xv_fab.resize(num_boxes_at_level[lev]);
195 Vector<FArrayBox> NC_yv_fab ; NC_yv_fab.resize(num_boxes_at_level[lev]);
196 Vector<FArrayBox> NC_xp_fab ; NC_xp_fab.resize(num_boxes_at_level[lev]);
197 Vector<FArrayBox> NC_yp_fab ; NC_yp_fab.resize(num_boxes_at_level[lev]);
198
199 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
200 {
202 NC_h_fab[idx],
203 NC_pm_fab[idx], NC_pn_fab[idx],
204 NC_xr_fab[idx], NC_yr_fab[idx],
205 NC_xu_fab[idx], NC_yu_fab[idx],
206 NC_xv_fab[idx], NC_yv_fab[idx],
207 NC_xp_fab[idx], NC_yp_fab[idx]);
208
209#ifdef _OPENMP
210#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
211#endif
212 {
213 // Don't tile this since we are operating on full FABs in this routine
214 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
215 {
216 FArrayBox &h_fab = (*vec_hOfTheConfusingName[lev])[mfi];
217 FArrayBox &pm_fab = (*vec_pm[lev])[mfi];
218 FArrayBox &pn_fab = (*vec_pn[lev])[mfi];
219 FArrayBox &xr_fab = (*vec_xr[lev])[mfi];
220 FArrayBox &yr_fab = (*vec_yr[lev])[mfi];
221 FArrayBox &xu_fab = (*vec_xu[lev])[mfi];
222 FArrayBox &yu_fab = (*vec_yu[lev])[mfi];
223 FArrayBox &xv_fab = (*vec_xv[lev])[mfi];
224 FArrayBox &yv_fab = (*vec_yv[lev])[mfi];
225 FArrayBox &xp_fab = (*vec_xp[lev])[mfi];
226 FArrayBox &yp_fab = (*vec_yp[lev])[mfi];
227
228 //
229 // FArrayBox to FArrayBox copy does "copy on intersection"
230 // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks
231 //
232
233 // Copy into both components of h
234 h_fab.template copy<RunOn::Device>(NC_h_fab[idx],0,0,1);
235 h_fab.template copy<RunOn::Device>(NC_h_fab[idx],0,1,1);
236
237 pm_fab.template copy<RunOn::Device>(NC_pm_fab[idx]);
238 pn_fab.template copy<RunOn::Device>(NC_pn_fab[idx]);
239
240 xr_fab.template copy<RunOn::Device>(NC_xr_fab[idx]);
241 yr_fab.template copy<RunOn::Device>(NC_yr_fab[idx]);
242 xu_fab.template copy<RunOn::Device>(NC_xu_fab[idx]);
243 yu_fab.template copy<RunOn::Device>(NC_yu_fab[idx]);
244 xv_fab.template copy<RunOn::Device>(NC_xv_fab[idx]);
245 yv_fab.template copy<RunOn::Device>(NC_yv_fab[idx]);
246 xp_fab.template copy<RunOn::Device>(NC_xp_fab[idx]);
247 yp_fab.template copy<RunOn::Device>(NC_yp_fab[idx]);
248 } // mf
249 } // omp
250 } // idx
251
252 const double dummy_time = 0.0_rt;
253 // Unconditional foextrap will overwrite periodicity, but EnforcePeriodicity will
254 // be called on h afterwards
255 FillPatch(lev,dummy_time,*vec_hOfTheConfusingName[lev],GetVecOfPtrs(vec_hOfTheConfusingName),
257 BdyVars::null,0,false,true,1);
258 FillPatch(lev,dummy_time,*vec_hOfTheConfusingName[lev],GetVecOfPtrs(vec_hOfTheConfusingName),
260 BdyVars::null,1,false,true,1);
261
262 int ng = vec_pm[lev]->nGrow();
263
264 const auto& dom_lo = amrex::lbound(geom[lev].Domain());
265 const auto& dom_hi = amrex::ubound(geom[lev].Domain());
266
267 //
268 // We need values of pm and pn outside the domain so we fill
269 // them here with foextrap
270 //
271 // We first fill interior ghost cells because we will need to extrapolate
272 // from ghost cells inside the domain to ghost cells outside the domain
273 //
274 vec_pm[lev]->FillBoundary(geom[lev].periodicity());
275 vec_pn[lev]->FillBoundary(geom[lev].periodicity());
276
277 vec_xr[lev]->FillBoundary(geom[lev].periodicity());
278 vec_yr[lev]->FillBoundary(geom[lev].periodicity());
279 vec_xu[lev]->FillBoundary(geom[lev].periodicity());
280 vec_yu[lev]->FillBoundary(geom[lev].periodicity());
281 vec_xv[lev]->FillBoundary(geom[lev].periodicity());
282 vec_yv[lev]->FillBoundary(geom[lev].periodicity());
283 vec_xp[lev]->FillBoundary(geom[lev].periodicity());
284 vec_yp[lev]->FillBoundary(geom[lev].periodicity());
285
286 for ( MFIter mfi(*vec_pm[lev]); mfi.isValid(); ++mfi )
287 {
288 Box bx = mfi.tilebox();
289
290 auto pm_fab = vec_pm[lev]->array(mfi);
291 auto pn_fab = vec_pn[lev]->array(mfi);
292
293 Box gbx_lox = adjCellLo(bx,0,ng); gbx_lox.grow(1,ng); gbx_lox.setBig (0,dom_lo.x-2);
294 Box gbx_hix = adjCellHi(bx,0,ng); gbx_hix.grow(1,ng); gbx_hix.setSmall(0,dom_hi.x+2);
295 Box gbx_loy = adjCellLo(bx,1,ng); gbx_loy.grow(0,ng); gbx_loy.setBig (1,dom_lo.y-2);
296 Box gbx_hiy = adjCellHi(bx,1,ng); gbx_hiy.grow(0,ng); gbx_hiy.setSmall(1,dom_hi.y+2);
297
298 // if (gbx_lox.ok()) amrex::AllPrint() << "GBX_XLO " << gbx_lox << std::endl;
299 // if (gbx_hix.ok()) amrex::AllPrint() << "GBX_XHI " << gbx_hix << std::endl;
300 // if (gbx_loy.ok()) amrex::AllPrint() << "GBX_YLO " << gbx_loy << std::endl;
301 // if (gbx_hiy.ok()) amrex::AllPrint() << "GBX_YHI " << gbx_hiy << std::endl;
302
303 if (gbx_lox.ok()) {
304 ParallelFor(gbx_lox, [=] AMREX_GPU_DEVICE (int i, int j, int k)
305 {
306 pm_fab(i,j,k,0) = pm_fab(dom_lo.x-1,j,k,0);
307 pn_fab(i,j,k,0) = pn_fab(dom_lo.x-1,j,k,0);
308 });
309 }
310 if (gbx_hix.ok()) {
311 ParallelFor(gbx_hix, [=] AMREX_GPU_DEVICE (int i, int j, int k)
312 {
313 pm_fab(i,j,k,0) = pm_fab(dom_hi.x+1,j,k,0);
314 pn_fab(i,j,k,0) = pn_fab(dom_hi.x+1,j,k,0);
315 });
316 }
317 if (gbx_loy.ok()) {
318 ParallelFor(gbx_loy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
319 {
320 pm_fab(i,j,k,0) = pm_fab(i,dom_lo.y-1,k,0);
321 pn_fab(i,j,k,0) = pn_fab(i,dom_lo.y-1,k,0);
322 });
323 }
324 if (gbx_hiy.ok()) {
325 ParallelFor(gbx_hiy, [=] AMREX_GPU_DEVICE (int i, int j, int k)
326 {
327 pm_fab(i,j,k,0) = pm_fab(i,dom_hi.y+1,k,0);
328 pn_fab(i,j,k,0) = pn_fab(i,dom_hi.y+1,k,0);
329 });
330 }
331 } // mfi
332
333 vec_hOfTheConfusingName[lev]->FillBoundary(geom[lev].periodicity());
334}
335
336/**
337 * @param lev Integer specifying the current level
338 */
339void
341{
342 // *** FArrayBox's at this level for holding the INITIAL data
343 Vector<FArrayBox> NC_fcor_fab ; NC_fcor_fab.resize(num_boxes_at_level[lev]);
344
345 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
346 {
348 NC_fcor_fab[idx]);
349
350#ifdef _OPENMP
351#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
352#endif
353 {
354 // Don't tile this since we are operating on full FABs in this routine
355 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
356 {
357 FArrayBox &fcor_fab = (*vec_fcor[lev])[mfi];
358
359 //
360 // FArrayBox to FArrayBox copy does "copy on intersection"
361 // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks
362 //
363
364 fcor_fab.template copy<RunOn::Device>(NC_fcor_fab[idx]);
365 } // mf
366 } // omp
367 } // idx
368 vec_fcor[lev]->FillBoundary(geom[lev].periodicity());
369}
370
371/**
372 * @param lev Integer specifying the current level
373 */
374void
376{
377 // *** FArrayBox's at this level for holding the INITIAL data
378 Vector<FArrayBox> NC_mskr_fab ; NC_mskr_fab.resize(num_boxes_at_level[lev]);
379 Vector<FArrayBox> NC_msku_fab ; NC_msku_fab.resize(num_boxes_at_level[lev]);
380 Vector<FArrayBox> NC_mskv_fab ; NC_mskv_fab.resize(num_boxes_at_level[lev]);
381
382 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
383 {
384 read_masks_from_netcdf(lev,boxes_at_level[lev][idx], nc_grid_file[lev][idx],
385 NC_mskr_fab[idx],NC_msku_fab[idx],
386 NC_mskv_fab[idx]);
387
388#ifdef _OPENMP
389#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
390#endif
391 {
392 // Don't tile this since we are operating on full FABs in this routine
393 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
394 {
395 FArrayBox &mskr_fab = (*vec_mskr[lev])[mfi];
396 FArrayBox &msku_fab = (*vec_msku[lev])[mfi];
397 FArrayBox &mskv_fab = (*vec_mskv[lev])[mfi];
398
399 //
400 // FArrayBox to FArrayBox copy does "copy on intersection"
401 // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks
402 //
403
404 mskr_fab.template copy<RunOn::Device>(NC_mskr_fab[idx]);
405 msku_fab.template copy<RunOn::Device>(NC_msku_fab[idx]);
406 mskv_fab.template copy<RunOn::Device>(NC_mskv_fab[idx]);
407 } // mf
408 } // omp
409 } // idx
410
411 update_mskp(lev);
412 vec_mskr[lev]->FillBoundary(geom[lev].periodicity());
413 vec_msku[lev]->FillBoundary(geom[lev].periodicity());
414 vec_mskv[lev]->FillBoundary(geom[lev].periodicity());
415 vec_mskp[lev]->FillBoundary(geom[lev].periodicity());
416}
417
418/**
419 * @param lev Integer specifying the current level
420 */
421void
423{
424 if (nc_bdry_file.empty()) {
425 amrex::Error("NetCDF boundary file name must be provided via input");
426 }
427
428 int width = 1;
429 for (int ifile = 0; ifile < nc_bdry_file.size(); ifile++) {
430 amrex::Real bdy_time_interval_this = read_bdry_from_netcdf(geom[0].Domain(), nc_bdry_file[ifile],
432 width, start_bdy_time,
435 if (ifile == 0) {
436 bdy_time_interval = bdy_time_interval_this;
437 } else {
438 AMREX_ALWAYS_ASSERT(bdy_time_interval == bdy_time_interval_this);
439 }
440 }
441
442 amrex::Print() << "Read in boundary data with width " << width;
443}
444
445/**
446 * \brief Helper function to initialize state and velocity data in a Fab from a REMORAdataset.
447 *
448 * @param lev Integer specifying current level
449 * @param state_fab FArrayBox object holding the state data we initialize
450 * @param temp_fab FArrayBox object holding the temperature data we initialize
451 * @param salt_fab FArrayBox object holding the salt data we initialize
452 * @param x_vel_fab FArrayBox object holding the x-velocity data we initialize
453 * @param y_vel_fab FArrayBox object holding the y-velocity data we initialize
454 * @param ubar_fab FArrayBox object holding the ubar data we initialize
455 * @param vbar_fab FArrayBox object holding the vbar data we initialize
456 * @param zeta_fab FArrayBox object holding the zeta data we initialize
457 * @param NC_temp_fab Vector of FArrayBox objects with the REMORA dataset specifying temperature
458 * @param NC_salt_fab Vector of FArrayBox objects with the REMORA dataset specifying salinity
459 * @param NC_xvel_fab Vector of FArrayBox objects with the REMORA dataset specifying x-velocity
460 * @param NC_yvel_fab Vector of FArrayBox objects with the REMORA dataset specifying y-velocity
461 * @param NC_ubar_fab Vector of FArrayBox objects with the REMORA dataset specifying ubar
462 * @param NC_vbar_fab Vector of FArrayBox objects with the REMORA dataset specifying vbar
463 * @param NC_zeta_fab Vector of FArrayBox objects with the REMORA dataset specifying zeta
464 */
465void
467 FArrayBox& temp_fab, FArrayBox& salt_fab,
468 FArrayBox& x_vel_fab, FArrayBox& y_vel_fab,
469 FArrayBox& ubar_fab, FArrayBox& vbar_fab,
470 const Vector<FArrayBox>& NC_temp_fab,
471 const Vector<FArrayBox>& NC_salt_fab,
472 const Vector<FArrayBox>& NC_xvel_fab,
473 const Vector<FArrayBox>& NC_yvel_fab,
474 const Vector<FArrayBox>& NC_ubar_fab,
475 const Vector<FArrayBox>& NC_vbar_fab)
476{
477 int nboxes = NC_xvel_fab.size();
478 for (int idx = 0; idx < nboxes; idx++)
479 {
480 //
481 // FArrayBox to FArrayBox copy does "copy on intersection"
482 // This only works here because we have broadcast the FArrayBox of data from the netcdf file to all ranks
483 //
484 temp_fab.template copy<RunOn::Device>(NC_temp_fab[idx]);
485 salt_fab.template copy<RunOn::Device>(NC_salt_fab[idx]);
486 x_vel_fab.template copy<RunOn::Device>(NC_xvel_fab[idx]);
487 y_vel_fab.template copy<RunOn::Device>(NC_yvel_fab[idx]);
488 ubar_fab.template copy<RunOn::Device>(NC_ubar_fab[idx],0,0,1);
489 vbar_fab.template copy<RunOn::Device>(NC_vbar_fab[idx],0,0,1);
490 } // idx
491}
492
493/**
494 * @param lev Integer specifying the current level
495 */
496void
498{
499 // *** FArrayBox's at this level for holding the INITIAL data
500 Vector<FArrayBox> NC_M2NC_fab ; NC_M2NC_fab.resize(num_boxes_at_level[lev]);
501 Vector<FArrayBox> NC_M3NC_fab ; NC_M3NC_fab.resize(num_boxes_at_level[lev]);
502 Vector<FArrayBox> NC_TempNC_fab ; NC_TempNC_fab.resize(num_boxes_at_level[lev]);
503 Vector<FArrayBox> NC_SaltNC_fab ; NC_SaltNC_fab.resize(num_boxes_at_level[lev]);
504
505 for (int idx = 0; idx < num_boxes_at_level[lev]; idx++)
506 {
512 NC_M2NC_fab[idx],NC_M3NC_fab[idx],
513 NC_TempNC_fab[idx],NC_SaltNC_fab[idx]);
514
515#ifdef _OPENMP
516#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
517#endif
518 {
519 // Don't tile this since we are operating on full FABs in this routine
520 for ( MFIter mfi(*cons_new[lev], false); mfi.isValid(); ++mfi )
521 {
523 FArrayBox &ubarNC_fab = (*vec_nudg_coeff[BdyVars::ubar][lev])[mfi];
524 ubarNC_fab.template copy<RunOn::Device>(NC_M2NC_fab[idx]);
525 FArrayBox &vbarNC_fab = (*vec_nudg_coeff[BdyVars::vbar][lev])[mfi];
526 vbarNC_fab.template copy<RunOn::Device>(NC_M2NC_fab[idx]);
527 }
529 FArrayBox &uNC_fab = (*vec_nudg_coeff[BdyVars::u][lev])[mfi];
530 uNC_fab.template copy<RunOn::Device>(NC_M3NC_fab[idx]);
531 FArrayBox &vNC_fab = (*vec_nudg_coeff[BdyVars::v][lev])[mfi];
532 vNC_fab.template copy<RunOn::Device>(NC_M3NC_fab[idx]);
533 }
535 FArrayBox &TempNC_fab = (*vec_nudg_coeff[BdyVars::t][lev])[mfi];
536 TempNC_fab.template copy<RunOn::Device>(NC_TempNC_fab[idx]);
537 }
539 FArrayBox &SaltNC_fab = (*vec_nudg_coeff[BdyVars::s][lev])[mfi];
540 SaltNC_fab.template copy<RunOn::Device>(NC_SaltNC_fab[idx]);
541 }
542
543 } // mf
544 } // omp
545 } // idx
546
548 vec_nudg_coeff[BdyVars::ubar][lev]->FillBoundary(geom[lev].periodicity());
549 vec_nudg_coeff[BdyVars::vbar][lev]->FillBoundary(geom[lev].periodicity());
552 }
554 vec_nudg_coeff[BdyVars::u][lev]->FillBoundary(geom[lev].periodicity());
555 vec_nudg_coeff[BdyVars::v][lev]->FillBoundary(geom[lev].periodicity());
558 }
560 vec_nudg_coeff[BdyVars::t][lev]->FillBoundary(geom[lev].periodicity());
562 }
564 vec_nudg_coeff[BdyVars::s][lev]->FillBoundary(geom[lev].periodicity());
566 }
567}
568
569/**
570 * @param[in ] lev level to read in river data
571 */
572void
574{
575 amrex::Vector<int> river_pos_x;
576 amrex::Vector<int> river_pos_y;
577 amrex::Vector<int> river_direction_tmp;
578
579 std::string river_x_name = "river_Xposition";
580 std::string river_y_name = "river_Eposition";
581 std::string river_dir_name = "river_direction";
582
583 read_vec_from_netcdf(lev, nc_riv_file, river_x_name, river_pos_x);
584 read_vec_from_netcdf(lev, nc_riv_file, river_y_name, river_pos_y);
585 read_vec_from_netcdf(lev, nc_riv_file, river_dir_name, river_direction_tmp);
586
587 int nriv = river_pos_x.size();
588 amrex::Gpu::DeviceVector<int> xpos_d(nriv);
589 amrex::Gpu::DeviceVector<int> ypos_d(nriv);
590 river_direction.resize(nriv);
591#ifdef AMREX_USE_GPU
592 Gpu::htod_memcpy(xpos_d.data(), river_pos_x.data(), sizeof(int)*nriv);
593 Gpu::htod_memcpy(ypos_d.data(), river_pos_y.data(), sizeof(int)*nriv);
594 Gpu::htod_memcpy(river_direction.data(), river_direction_tmp.data(), sizeof(int)*nriv);
595#else
596 std::memcpy(xpos_d.data(), river_pos_x.data(), sizeof(int)*nriv);
597 std::memcpy(ypos_d.data(), river_pos_y.data(), sizeof(int)*nriv);
598 std::memcpy(river_direction.data(), river_direction_tmp.data(), sizeof(int)*nriv);
599#endif
600 const int* xpos_ptr = xpos_d.data();
601 const int* ypos_ptr = ypos_d.data();
602
603 for (amrex::MFIter mfi(*(vec_river_position[lev]).get(),true); mfi.isValid(); ++mfi) {
604 amrex::Box bx = mfi.growntilebox(amrex::IntVect(NGROW,NGROW,0));
605 auto river_pos = vec_river_position[lev]->array(mfi);
606 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int ) {
607 for (int iriv=0; iriv < nriv; iriv++) {
608 int xriv = xpos_ptr[iriv]-1;
609 int yriv = ypos_ptr[iriv]-1;
610 if (i==xriv && j==yriv) {
611 river_pos(i,j,0) = iriv;
612 }
613 }
614 });
615 }
616}
617
618/**
619 * @param[inout] mf multifab of data to convert
620 */
621void
623 Real inv_days_to_inv_s = 1.0_rt / (3600._rt * 24._rt);
624
625 for ( MFIter mfi(*mf, TilingIfNotGPU()); mfi.isValid(); ++mfi )
626 {
627 Array4<Real> const& arr = mf->array(mfi);
628 Box bx = mfi.growntilebox(IntVect(NGROW,NGROW,0));
629 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) {
630 arr(i,j,k) *= inv_days_to_inv_s;
631 });
632 }
633
634}
635
636#endif // REMORA_USE_NETCDF
#define NGROW
#define Temp_comp
#define Salt_comp
void init_state_from_netcdf(int lev, FArrayBox &temp_fab, FArrayBox &salt_fab, FArrayBox &x_vel_fab, FArrayBox &y_vel_fab, FArrayBox &ubar_fab, FArrayBox &vbar_fab, const Vector< FArrayBox > &NC_temp_fab, const Vector< FArrayBox > &NC_salt_fab, const Vector< FArrayBox > &NC_xvel_fab, const Vector< FArrayBox > &NC_yvel_fab, const Vector< FArrayBox > &NC_ubar_fab, const Vector< FArrayBox > &NC_vbar_fab)
helper function to initialize state from netcdf
void read_coriolis_from_netcdf(int lev, const Box &domain, const std::string &fname, FArrayBox &NC_fcor_fab)
helper function to read coriolis factor from netcdf
void read_vec_from_netcdf(int lev, const std::string &fname, const std::string &field_name, amrex::Vector< int > &vec_dat)
helper function to read in vector of data from netcdf
void read_data_from_netcdf(int, const Box &domain, const std::string &fname, FArrayBox &NC_temp_fab, FArrayBox &NC_salt_fab, FArrayBox &NC_xvel_fab, FArrayBox &NC_yvel_fab, FArrayBox &NC_ubar_fab, FArrayBox &NC_vbar_fab)
helper function for reading in initial state data from netcdf
void read_zeta_from_netcdf(int lev, const Box &domain, const std::string &fname, FArrayBox &NC_zeta_fab)
helper function to read sea surface height from netcdf
Real read_bdry_from_netcdf(const Box &domain, const std::string &fname, Vector< Vector< FArrayBox > > &bdy_data_xlo, Vector< Vector< FArrayBox > > &bdy_data_xhi, Vector< Vector< FArrayBox > > &bdy_data_ylo, Vector< Vector< FArrayBox > > &bdy_data_yhi, int &width, amrex::Real &start_bdy_time, std::string bdry_time_varname, amrex::GpuArray< amrex::GpuArray< bool, AMREX_SPACEDIM *2 >, BdyVars::NumTypes+1 > &)
helper function for reading boundary data from netcdf
void read_masks_from_netcdf(int, const Box &domain, const std::string &fname, FArrayBox &NC_mskr_fab, FArrayBox &NC_msku_fab, FArrayBox &NC_mskv_fab)
helper function for reading in land-sea masks from netcdf
void read_bathymetry_from_netcdf(int lev, const Box &domain, const std::string &fname, FArrayBox &NC_h_fab, FArrayBox &NC_pm_fab, FArrayBox &NC_pn_fab, FArrayBox &NC_xr_fab, FArrayBox &NC_yr_fab, FArrayBox &NC_xu_fab, FArrayBox &NC_yu_fab, FArrayBox &NC_xv_fab, FArrayBox &NC_yv_fab, FArrayBox &NC_xp_fab, FArrayBox &NC_yp_fab)
helper function to read bathymetry from netcdf
void read_clim_nudg_coeff_from_netcdf(int lev, const Box &domain, const std::string &fname, bool do_m2_clim_nudg, bool do_m3_clim_nudg, bool do_temp_clim_nudg, bool do_salt_clim_nudg, FArrayBox &NC_M2NC_fab, FArrayBox &NC_M3NC_fab, FArrayBox &NC_TempNC_fab, FArrayBox &NC_SaltNC_fab)
helper function to read climatology nudging from netcdf
void init_bathymetry_from_netcdf(int lev)
Bathymetry data initialization from NetCDF file.
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_fcor
coriolis factor (2D)
Definition REMORA.H:360
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_hOfTheConfusingName
Bathymetry data (2D, positive valued, h in ROMS)
Definition REMORA.H:229
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pm
horizontal scaling factor: 1 / dx (2D)
Definition REMORA.H:355
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_yv
y_grid on v-points (2D)
Definition REMORA.H:375
amrex::Vector< amrex::MultiFab * > cons_new
multilevel data container for current step's scalar data: temperature, salinity, passive scalar
Definition REMORA.H:217
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskr
land/sea mask at cell centers (2D)
Definition REMORA.H:346
std::string bdry_time_varname
Name of time field for boundary data.
Definition REMORA.H:1278
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_yp
y_grid on psi-points (2D)
Definition REMORA.H:380
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_xr
x_grid on rho points (2D)
Definition REMORA.H:363
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_xv
x_grid on v-points (2D)
Definition REMORA.H:373
amrex::Vector< amrex::Vector< amrex::Box > > boxes_at_level
the boxes specified at each level by tagging criteria
Definition REMORA.H:1091
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_msku
land/sea mask at x-faces (2D)
Definition REMORA.H:348
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
void init_data_from_netcdf(int lev)
Problem initialization from NetCDF file.
void init_masks_from_netcdf(int lev)
Mask data initialization from NetCDF file.
amrex::Vector< amrex::MultiFab * > yvel_new
multilevel data container for current step's y velocities (v in ROMS)
Definition REMORA.H:221
amrex::Real bdy_time_interval
Interval between boundary data times.
Definition REMORA.H:1001
void init_bdry_from_netcdf()
Boundary data initialization from NetCDF file.
amrex::Real start_bdy_time
Start time in the time series of boundary data.
Definition REMORA.H:999
amrex::Vector< int > num_boxes_at_level
how many boxes specified at each level by tagging criteria
Definition REMORA.H:1087
amrex::Vector< amrex::MultiFab * > xvel_new
multilevel data container for current step's x velocities (u in ROMS)
Definition REMORA.H:219
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskp
land/sea mask at cell corners (2D)
Definition REMORA.H:352
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskv
land/sea mask at y-faces (2D)
Definition REMORA.H:350
amrex::Vector< std::unique_ptr< REMORAPhysBCFunct > > physbcs
Vector (over level) of functors to apply physical boundary conditions.
Definition REMORA.H:1108
void FillPatch(int lev, amrex::Real time, amrex::MultiFab &mf_to_be_filled, amrex::Vector< amrex::MultiFab * > const &mfs, const int bccomp, const int bdy_var_type=BdyVars::null, const int icomp=0, const bool fill_all=true, const bool fill_set=true, const int n_not_fill=0, const int icomp_calc=0, const amrex::Real dt=amrex::Real(0.0), const amrex::MultiFab &mf_calc=amrex::MultiFab())
Fill a new MultiFab by copying in phi from valid region and filling ghost cells.
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
std::string nc_clim_coeff_file
NetCDF climatology coefficient file.
Definition REMORA.H:1275
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.
void init_riv_pos_from_netcdf(int lev)
static amrex::Vector< std::string > nc_bdry_file
NetCDF boundary data.
Definition REMORA.H:48
std::string nc_riv_file
Definition REMORA.H:1270
void update_mskp(int lev)
Set psi-point mask to be consistent with rho-point mask.
void init_zeta_from_netcdf(int lev)
Sea-surface height data initialization from NetCDF file.
void init_coriolis_from_netcdf(int lev)
Coriolis parameter data initialization from NetCDF file.
amrex::Vector< amrex::Real > t_new
new time at each level
Definition REMORA.H:1098
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1221
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > vec_river_position
iMultiFab for river positions; contents are indices of rivers
Definition REMORA.H:1036
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_xp
x_grid on psi-points (2D)
Definition REMORA.H:378
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
static amrex::Vector< amrex::Vector< std::string > > nc_grid_file
NetCDF grid file.
Definition REMORA.H:50
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_zeta
free surface height (2D)
Definition REMORA.H:343
amrex::Gpu::DeviceVector< int > river_direction
Vector over rivers of river direction: 0: u-face; 1: v-face; 2: w-face.
Definition REMORA.H:1038
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_vbar
barotropic y velocity (2D)
Definition REMORA.H:341
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_ubar
barotropic x velocity (2D)
Definition REMORA.H:339
void init_clim_nudg_coeff_from_netcdf(int lev)
Climatology nudging coefficient initialization from NetCDF file.
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_yu
y_grid on u-points (2D)
Definition REMORA.H:370
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_xu
x_grid on u-points (2D)
Definition REMORA.H:368
amrex::GpuArray< amrex::GpuArray< bool, AMREX_SPACEDIM *2 >, BdyVars::NumTypes+1 > phys_bc_need_data
These are flags that indicate whether we need to read in boundary data from file.
Definition REMORA.H:1134
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > vec_nudg_coeff
Climatology nudging coefficients.
Definition REMORA.H:414
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pn
horizontal scaling factor: 1 / dy (2D)
Definition REMORA.H:357
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
static amrex::Vector< amrex::Vector< std::string > > nc_init_file
NetCDF initialization file.
Definition REMORA.H:49
void convert_inv_days_to_inv_s(amrex::MultiFab *)
Convert data in a multifab from inverse days to inverse seconds.
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_yr
y_grid on rho points (2D)
Definition REMORA.H:365