REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_make_new_level.cpp
Go to the documentation of this file.
1/**
2 * \file REMORA_make_new_level.cpp
3 */
4
5#include <REMORA.H>
7
8#include <AMReX_buildInfo.H>
9
10using namespace amrex;
11
12/**
13 * Make a new level using provided BoxArray and DistributionMapping and
14 * fill with interpolated coarse level data (overrides the pure virtual function in AmrCore)
15 * regrid --> RemakeLevel (if level already existed)
16 * regrid --> MakeNewLevelFromCoarse (if adding new level)
17 *
18 * @param[in ] lev level to make
19 * @param[in ] time current time
20 * @param[in ] ba BoxArray for the level
21 * @param[in ] dm DistributionMapping for the level
22 */
23void
24REMORA::MakeNewLevelFromCoarse (int lev, Real time, const BoxArray& ba,
25 const DistributionMapping& dm)
26{
27 BoxList bl2d = ba.boxList();
28 for (auto& b : bl2d) {
29 b.setRange(2,0);
30 }
31 BoxArray ba2d(std::move(bl2d));
32
33 amrex::Print() << "Making level " << lev << " from coarse" << std::endl;
34 amrex::Print() << "GRIDS AT LEVEL " << lev << " ARE " << ba << std::endl;
35
36 cons_new[lev] = new MultiFab(ba, dm, NCONS, cons_new[lev-1]->nGrowVect());
37 cons_old[lev] = new MultiFab(ba, dm, NCONS, cons_new[lev-1]->nGrowVect());
38
39 xvel_new[lev] = new MultiFab(convert(ba, IntVect(1,0,0)), dm, 1, xvel_new[lev-1]->nGrowVect());
40 xvel_old[lev] = new MultiFab(convert(ba, IntVect(1,0,0)), dm, 1, xvel_new[lev-1]->nGrowVect());
41
42 yvel_new[lev] = new MultiFab(convert(ba, IntVect(0,1,0)), dm, 1, yvel_new[lev-1]->nGrowVect());
43 yvel_old[lev] = new MultiFab(convert(ba, IntVect(0,1,0)), dm, 1, yvel_new[lev-1]->nGrowVect());
44
45 zvel_new[lev] = new MultiFab(convert(ba, IntVect(0,0,1)), dm, 1, zvel_new[lev-1]->nGrowVect());
46 zvel_old[lev] = new MultiFab(convert(ba, IntVect(0,0,1)), dm, 1, zvel_new[lev-1]->nGrowVect());
47
48 resize_stuff(lev);
49
50 vec_Zt_avg1[lev].reset(new MultiFab(ba2d ,dm,1,IntVect(NGROW+1,NGROW+1,0))); //2d, average of the free surface (zeta)
51 vec_h[lev].reset(new MultiFab(ba2d ,dm,2,IntVect(NGROW+1,NGROW+1,0))); //2d, average of the free surface (zeta)
52 vec_ubar[lev].reset(new MultiFab(convert(ba2d,IntVect(1,0,0)),dm,3,IntVect(NGROW,NGROW,0)));
53 vec_vbar[lev].reset(new MultiFab(convert(ba2d,IntVect(0,1,0)),dm,3,IntVect(NGROW,NGROW,0)));
54
55 vec_ru[lev].reset(new MultiFab(convert(ba,IntVect(1,0,0)),dm,2,IntVect(NGROW,NGROW,0))); // RHS u (incl horizontal and vertical advection)
56 vec_rv[lev].reset(new MultiFab(convert(ba,IntVect(0,1,0)),dm,2,IntVect(NGROW,NGROW,0))); // RHS v
57
58 vec_ru2d[lev].reset(new MultiFab(convert(ba2d,IntVect(1,0,0)),dm,2,IntVect(NGROW,NGROW,0))); // RHS u for 2d
59 vec_rv2d[lev].reset(new MultiFab(convert(ba2d,IntVect(0,1,0)),dm,2,IntVect(NGROW,NGROW,0))); // RHS v for 2d
60
61 t_new[lev] = time;
62 t_old[lev] = time - 1.e200_rt;
63
64 init_masks(lev, ba, dm);
65
66 init_stuff(lev, ba, dm);
67
68 cons_new[lev]->setVal(0.0_rt);
69 xvel_new[lev]->setVal(0.0_rt);
70 yvel_new[lev]->setVal(0.0_rt);
71 zvel_new[lev]->setVal(0.0_rt);
72
73 cons_old[lev]->setVal(0.0_rt);
74 xvel_old[lev]->setVal(0.0_rt);
75 yvel_old[lev]->setVal(0.0_rt);
76 zvel_old[lev]->setVal(0.0_rt);
77
78 vec_ru[lev]->setVal(0.0_rt);
79 vec_rv[lev]->setVal(0.0_rt);
80
81 vec_ru2d[lev]->setVal(0.0_rt);
82 vec_rv2d[lev]->setVal(0.0_rt);
83
84 vec_ubar[lev]->setVal(0.0_rt);
85 vec_vbar[lev]->setVal(0.0_rt);
86
87
92
93 FillCoarsePatch(lev, time, vec_h[lev].get(), vec_h[lev-1].get(),
95 FillCoarsePatch(lev, time, vec_Zt_avg1[lev].get(), vec_Zt_avg1[lev-1].get(),BCVars::cons_bc);
96 for (int icomp=0; icomp<3; icomp++) {
97 FillCoarsePatch(lev, time, vec_ubar[lev].get(), vec_ubar[lev-1].get(),BCVars::ubar_bc,
98 BdyVars::ubar,icomp,false);
99 FillCoarsePatch(lev, time, vec_vbar[lev].get(), vec_vbar[lev-1].get(),BCVars::vbar_bc,
100 BdyVars::vbar,icomp,false);
101 }
102 for (int icomp=0; icomp<2; icomp++) {
103 FillCoarsePatch(lev, time, vec_ru[lev].get(), vec_ru[lev-1].get(),BCVars::xvel_bc,
104 BdyVars::null,icomp,false);
105 FillCoarsePatch(lev, time, vec_rv[lev].get(), vec_rv[lev-1].get(),BCVars::yvel_bc,
106 BdyVars::null,icomp,false);
107 FillCoarsePatch(lev, time, vec_ru2d[lev].get(), vec_ru2d[lev-1].get(),BCVars::xvel_bc,
108 BdyVars::null,icomp,false);
109 FillCoarsePatch(lev, time, vec_rv2d[lev].get(), vec_rv2d[lev-1].get(),BCVars::yvel_bc,
110 BdyVars::null,icomp,false);
111 }
112
113 // Not totally sure foextrap is right here
114 FillCoarsePatchPC(lev, time, vec_mskr[lev].get(), vec_mskr[lev-1].get(),
116
118
119
120 set_grid_scale(lev);
122
123 init_set_vmix(lev);
124 set_hmixcoef(lev);
125 set_coriolis(lev);
127 // Previously set smflux
128
129 // ********************************************************************************************
130 // If we are making a new level then the FillPatcher for this level hasn't been allocated yet
131 // ********************************************************************************************
132 if (cf_width >= 0) {
135 }
136
137#ifdef REMORA_USE_PARTICLES
138 // particleData.Redistribute();
139#endif
140}
141
142/**
143 * Remake an existing level using provided BoxArray and DistributionMapping and
144 * fill with existing fine and coarse data.
145 * overrides the pure virtual function in AmrCore
146 * @param[in ] lev level to make
147 * @param[in ] time current time
148 * @param[in ] ba BoxArray for the level
149 * @param[in ] dm DistributionMapping for the level
150 */
151void
152REMORA::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapping& dm)
153{
154 BoxArray ba_old(cons_new[lev]->boxArray());
155 DistributionMapping dm_old(cons_new[lev]->DistributionMap());
156
157 BoxList bl2d = ba.boxList();
158 for (auto& b : bl2d) {
159 b.setRange(2,0);
160 }
161 BoxArray ba2d(std::move(bl2d));
162
163 amrex::Print() << "Remaking level " << lev << std::endl;
164 amrex::Print() << "GRIDS AT LEVEL " << lev << " ARE " << ba << std::endl;
165
166#if (NGROW==2)
167 int ngrow_state = ComputeGhostCells(solverChoice.spatial_order)+1;
168 int ngrow_vels = ComputeGhostCells(solverChoice.spatial_order)+1;
169 int ngrow_zeta = ComputeGhostCells(solverChoice.spatial_order)+1;
171 int ngrow_velbar = ComputeGhostCells(solverChoice.spatial_order);
172#else
173 int ngrow_state = ComputeGhostCells(solverChoice.spatial_order)+2;
174 int ngrow_vels = ComputeGhostCells(solverChoice.spatial_order)+2;
175 int ngrow_zeta = ComputeGhostCells(solverChoice.spatial_order)+2;
177 int ngrow_velbar = ComputeGhostCells(solverChoice.spatial_order)+1;
178#endif
179
180 MultiFab tmp_cons_new(ba, dm, NCONS, ngrow_state);
181 MultiFab tmp_cons_old(ba, dm, NCONS, ngrow_state);
182
183 MultiFab tmp_xvel_new(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels);
184 MultiFab tmp_xvel_old(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels);
185
186 MultiFab tmp_yvel_new(convert(ba, IntVect(0,1,0)), dm, 1, ngrow_vels);
187 MultiFab tmp_yvel_old(convert(ba, IntVect(0,1,0)), dm, 1, ngrow_vels);
188
189 MultiFab tmp_zvel_new(convert(ba, IntVect(0,0,1)), dm, 1, IntVect(ngrow_vels,ngrow_vels,0));
190 MultiFab tmp_zvel_old(convert(ba, IntVect(0,0,1)), dm, 1, IntVect(ngrow_vels,ngrow_vels,0));
191
192 MultiFab tmp_Zt_avg1_new(ba2d, dm, 1, IntVect(ngrow_zeta,ngrow_zeta,0));
193 MultiFab tmp_h(ba2d, dm, 2, IntVect(ngrow_h,ngrow_h,0));
194
195 MultiFab tmp_ubar_new(convert(ba2d, IntVect(1,0,0)), dm, 3, IntVect(ngrow_velbar,ngrow_velbar,0));
196
197 MultiFab tmp_vbar_new(convert(ba2d, IntVect(0,1,0)), dm, 3, IntVect(ngrow_velbar,ngrow_velbar,0));
198
199 MultiFab tmp_ru_new(convert(ba, IntVect(1,0,0)),dm,2,IntVect(NGROW,NGROW,0));
200 MultiFab tmp_rv_new(convert(ba, IntVect(0,1,0)),dm,2,IntVect(NGROW,NGROW,0));
201
202 MultiFab tmp_ru2d_new(convert(ba2d, IntVect(1,0,0)),dm,2,IntVect(NGROW,NGROW,0));
203 MultiFab tmp_rv2d_new(convert(ba2d, IntVect(0,1,0)),dm,2,IntVect(NGROW,NGROW,0));
204
205 init_masks(lev, ba, dm);
206
207 tmp_cons_new.setVal(0.0_rt);
208 tmp_xvel_new.setVal(0.0_rt);
209 tmp_yvel_new.setVal(0.0_rt);
210 tmp_zvel_new.setVal(0.0_rt);
211
212 tmp_cons_old.setVal(0.0_rt);
213 tmp_xvel_old.setVal(0.0_rt);
214 tmp_yvel_old.setVal(0.0_rt);
215 tmp_zvel_old.setVal(0.0_rt);
216
217 tmp_ru_new.setVal(0.0_rt);
218 tmp_rv_new.setVal(0.0_rt);
219
220 tmp_ru2d_new.setVal(0.0_rt);
221 tmp_rv2d_new.setVal(0.0_rt);
222
223 tmp_ubar_new.setVal(0.0_rt);
224 tmp_vbar_new.setVal(0.0_rt);
225
226
227 // This will fill the temporary MultiFabs with data from previous fine data as well as coarse where needed
228 FillPatch(lev, time, tmp_cons_new, cons_new, BCVars::cons_bc, BdyVars::t,0,true,false);
229 FillPatch(lev, time, tmp_xvel_new, xvel_new, BCVars::xvel_bc, BdyVars::u,0,true,false,0,0,0.0,tmp_xvel_new);
230 FillPatch(lev, time, tmp_yvel_new, yvel_new, BCVars::yvel_bc, BdyVars::v,0,true,false,0,0,0.0,tmp_yvel_new);
231 FillPatch(lev, time, tmp_zvel_new, zvel_new, BCVars::zvel_bc, BdyVars::null,0,true,false);
232
233 FillPatch(lev, time, tmp_h, GetVecOfPtrs(vec_h), BCVars::cons_bc, BdyVars::null,0,false,false);
234 FillPatch(lev, time, tmp_h, GetVecOfPtrs(vec_h), BCVars::cons_bc, BdyVars::null,1,false,false);
235 FillPatch(lev, time, tmp_Zt_avg1_new, GetVecOfPtrs(vec_Zt_avg1), BCVars::zeta_bc, BdyVars::null,0,true,false);
236 for (int icomp=0; icomp<3; icomp++) {
237 FillPatch(lev, time, tmp_ubar_new, GetVecOfPtrs(vec_ubar), BCVars::ubar_bc, BdyVars::ubar, icomp,false,false);
238 FillPatch(lev, time, tmp_vbar_new, GetVecOfPtrs(vec_vbar), BCVars::vbar_bc, BdyVars::vbar, icomp,false,false);
239 }
240 for (int icomp=0; icomp<2; icomp++) {
241 FillPatch(lev, time, tmp_ru_new, GetVecOfPtrs(vec_ru),BCVars::xvel_bc, BdyVars::null, icomp,false,false);
242 FillPatch(lev, time, tmp_rv_new, GetVecOfPtrs(vec_rv),BCVars::yvel_bc, BdyVars::null, icomp,false,false);
243 // These might want to have BCVars::ubar_bc and vbar_bc
244 FillPatch(lev, time, tmp_ru2d_new, GetVecOfPtrs(vec_ru2d),BCVars::xvel_bc, BdyVars::null, icomp,false,false);
245 FillPatch(lev, time, tmp_rv2d_new, GetVecOfPtrs(vec_rv2d),BCVars::yvel_bc, BdyVars::null, icomp,false,false);
246 }
247
248 MultiFab::Copy(tmp_cons_old,tmp_cons_new,0,0,NCONS,tmp_cons_new.nGrowVect());
249 MultiFab::Copy(tmp_xvel_old,tmp_xvel_new,0,0, 1,tmp_xvel_new.nGrowVect());
250 MultiFab::Copy(tmp_yvel_old,tmp_yvel_new,0,0, 1,tmp_yvel_new.nGrowVect());
251 MultiFab::Copy(tmp_zvel_old,tmp_zvel_new,0,0, 1,tmp_zvel_new.nGrowVect());
252
253 std::swap(tmp_cons_new, *cons_new[lev]);
254 std::swap(tmp_cons_old, *cons_old[lev]);
255 std::swap(tmp_xvel_new, *xvel_new[lev]);
256 std::swap(tmp_xvel_old, *xvel_old[lev]);
257 std::swap(tmp_yvel_new, *yvel_new[lev]);
258 std::swap(tmp_yvel_old, *yvel_old[lev]);
259 std::swap(tmp_zvel_new, *zvel_new[lev]);
260 std::swap(tmp_zvel_old, *zvel_old[lev]);
261 std::swap(tmp_Zt_avg1_new, *vec_Zt_avg1[lev]);
262 std::swap(tmp_h, *vec_h[lev]);
263 std::swap(tmp_ubar_new, *vec_ubar[lev]);
264 std::swap(tmp_vbar_new, *vec_vbar[lev]);
265 std::swap(tmp_ru_new, *vec_ru[lev]);
266 std::swap(tmp_rv_new, *vec_rv[lev]);
267 std::swap(tmp_ru2d_new, *vec_ru2d[lev]);
268 std::swap(tmp_rv2d_new, *vec_rv2d[lev]);
269
270 t_new[lev] = time;
271 t_old[lev] = time - 1.e200_rt;
272
273 init_masks(lev, ba, dm);
274 FillCoarsePatchPC(lev, time, vec_mskr[lev].get(), vec_mskr[lev-1].get(),
277
278 init_stuff(lev, ba, dm);
279
280 set_grid_scale(lev);
282
283 init_set_vmix(lev);
284 set_hmixcoef(lev);
285 set_coriolis(lev);
287 // Previously set smflux here
288
289 // We need to re-define the FillPatcher if the grids have changed
290 if (lev > 0 && cf_width >= 0) {
291 bool ba_changed = (ba != ba_old);
292 bool dm_changed = (dm != dm_old);
293 if (ba_changed || dm_changed) {
295 }
296 }
297
298#ifdef REMORA_USE_PARTICLES
299 particleData.Redistribute();
300#endif
301}
302
303/**
304 * Make a new level from scratch using provided BoxArray and DistributionMapping.
305 * This is called both for initialization and for restart
306 * (overrides the pure virtual function in AmrCore)
307 * main.cpp --> REMORA::InitData --> InitFromScratch --> MakeNewGrids --> MakeNewLevelFromScratch
308 * restart --> MakeNewGrids --> MakeNewLevelFromScratch
309 *
310 * @param[in ] lev level to make
311 * @param[in ] time current time
312 * @param[in ] ba BoxArray for the level
313 * @param[in ] dm DistributionMapping for the level
314 */
315void REMORA::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba,
316 const DistributionMapping& dm)
317{
318 // Set BoxArray grids and DistributionMapping dmap in AMReX_AmrMesh.H class
319 SetBoxArray(lev, ba);
320 SetDistributionMap(lev, dm);
321
322 BoxList bl2d = ba.boxList();
323 for (auto& b : bl2d) {
324 b.setRange(2,0);
325 }
326 BoxArray ba2d(std::move(bl2d));
327
328 amrex::Print() << "Making level " << lev << " from scratch" << std::endl;
329 amrex::Print() << "GRIDS AT LEVEL " << lev << " ARE " << ba << std::endl;
330
331 // The number of ghost cells for density must be 1 greater than that for velocity
332 // so that we can go back in forth between velocity and momentum on all faces
333#if NGROW==2
334 int ngrow_state = ComputeGhostCells(solverChoice.spatial_order)+1;
335 int ngrow_vels = ComputeGhostCells(solverChoice.spatial_order)+1;
336#else
337 int ngrow_state = ComputeGhostCells(solverChoice.spatial_order)+2;
338 int ngrow_vels = ComputeGhostCells(solverChoice.spatial_order)+2;
339#endif
340
341 cons_old[lev] = new MultiFab(ba, dm, NCONS, ngrow_state);
342 cons_new[lev] = new MultiFab(ba, dm, NCONS, ngrow_state);
343
344 xvel_new[lev] = new MultiFab(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels);
345 xvel_old[lev] = new MultiFab(convert(ba, IntVect(1,0,0)), dm, 1, ngrow_vels);
346
347 yvel_new[lev] = new MultiFab(convert(ba, IntVect(0,1,0)), dm, 1, ngrow_vels);
348 yvel_old[lev] = new MultiFab(convert(ba, IntVect(0,1,0)), dm, 1, ngrow_vels);
349
350 zvel_new[lev] = new MultiFab(convert(ba, IntVect(0,0,1)), dm, 1, IntVect(ngrow_vels,ngrow_vels,0));
351 zvel_old[lev] = new MultiFab(convert(ba, IntVect(0,0,1)), dm, 1, IntVect(ngrow_vels,ngrow_vels,0));
352
353 resize_stuff(lev);
354
355 vec_Zt_avg1[lev].reset(new MultiFab(ba2d ,dm,1,IntVect(NGROW+1,NGROW+1,0))); //2d, average of the free surface (zeta)
356 vec_h[lev].reset(new MultiFab(ba2d ,dm,2,IntVect(NGROW+1,NGROW+1,0))); //2d, bathymetry
357 vec_ubar[lev].reset(new MultiFab(convert(ba2d,IntVect(1,0,0)),dm,3,IntVect(NGROW,NGROW,0)));
358 vec_vbar[lev].reset(new MultiFab(convert(ba2d,IntVect(0,1,0)),dm,3,IntVect(NGROW,NGROW,0)));
359
360 vec_ru[lev].reset(new MultiFab(convert(ba,IntVect(1,0,0)),dm,2,IntVect(NGROW,NGROW,0))); // RHS u (incl horizontal and vertical advection)
361 vec_rv[lev].reset(new MultiFab(convert(ba,IntVect(0,1,0)),dm,2,IntVect(NGROW,NGROW,0))); // RHS v
362
363 vec_ru2d[lev].reset(new MultiFab(convert(ba2d,IntVect(1,0,0)),dm,2,IntVect(NGROW,NGROW,0))); // RHS u (incl horizontal and vertical advection)
364 vec_rv2d[lev].reset(new MultiFab(convert(ba2d,IntVect(0,1,0)),dm,2,IntVect(NGROW,NGROW,0))); // RHS v
365
366 init_masks(lev, ba, dm);
367 init_stuff(lev, ba, dm);
368
369 init_only(lev, time);
370
371#ifdef REMORA_USE_PARTICLES
372 if (restart_chkfile.empty()) {
373 if (lev == 0) {
374 initializeTracers((ParGDBBase*)GetParGDB(),vec_z_phys_nd);
375 } else {
376 particleData.Redistribute();
377 }
378 }
379#endif
380}
381
382/**
383 * @param[in ] lev level do operate on
384 */
386{
387 vec_z_phys_nd.resize(lev+1);
388
389 vec_h.resize(lev+1);
390 vec_Zt_avg1.resize(lev+1);
391 vec_z_w.resize(lev+1);
392 vec_z_r.resize(lev+1);
393 vec_Hz.resize(lev+1);
394 vec_Huon.resize(lev+1);
395 vec_Hvom.resize(lev+1);
396 vec_Akv.resize(lev+1);
397 vec_Akt.resize(lev+1);
398 vec_visc2_p.resize(lev+1);
399 vec_visc2_r.resize(lev+1);
400 vec_diff2.resize(lev+1);
401 vec_ru.resize(lev+1);
402 vec_rv.resize(lev+1);
403 vec_ru2d.resize(lev+1);
404 vec_rv2d.resize(lev+1);
405 vec_rufrc.resize(lev+1);
406 vec_rvfrc.resize(lev+1);
407 vec_sustr.resize(lev+1);
408 vec_svstr.resize(lev+1);
409 vec_btflx.resize(lev+1);
410 vec_stflx.resize(lev+1);
411 vec_btflux.resize(lev+1);
412 vec_stflux.resize(lev+1);
413 vec_lrflx.resize(lev+1);
414 vec_lhflx.resize(lev+1);
415 vec_shflx.resize(lev+1);
416 vec_rain.resize(lev+1);
417 vec_evap.resize(lev+1);
418 vec_rdrag.resize(lev+1);
419 vec_rdrag2.resize(lev+1);
420 vec_ZoBot.resize(lev+1);
421 vec_bustr.resize(lev+1);
422 vec_bvstr.resize(lev+1);
423 vec_uwind.resize(lev+1);
424 vec_vwind.resize(lev+1);
425 vec_alpha.resize(lev+1);
426 vec_beta.resize(lev+1);
427
428 vec_DU_avg1.resize(lev+1);
429 vec_DU_avg2.resize(lev+1);
430 vec_DV_avg1.resize(lev+1);
431 vec_DV_avg2.resize(lev+1);
432 vec_rubar.resize(lev+1);
433 vec_rvbar.resize(lev+1);
434 vec_rzeta.resize(lev+1);
435 vec_ubar.resize(lev+1);
436 vec_vbar.resize(lev+1);
437 vec_zeta.resize(lev+1);
438 vec_mskr.resize(lev+1);
439 vec_msku.resize(lev+1);
440 vec_mskv.resize(lev+1);
441 vec_mskp.resize(lev+1);
442 vec_sstore.resize(lev+1);
443
444 vec_pm.resize(lev+1);
445 vec_pn.resize(lev+1);
446 vec_fcor.resize(lev+1);
447
448 vec_xr.resize(lev+1);
449 vec_yr.resize(lev+1);
450 vec_xu.resize(lev+1);
451 vec_yu.resize(lev+1);
452 vec_xv.resize(lev+1);
453 vec_yv.resize(lev+1);
454 vec_xp.resize(lev+1);
455 vec_yp.resize(lev+1);
456
457 vec_rhoS.resize(lev+1);
458 vec_rhoA.resize(lev+1);
459 vec_bvf.resize(lev+1);
460
461 vec_tke.resize(lev+1);
462 vec_gls.resize(lev+1);
463 vec_Lscale.resize(lev+1);
464 vec_Akk.resize(lev+1);
465 vec_Akp.resize(lev+1);
466
467 vec_river_position.resize(lev+1);
468
469 if (lev==0) vec_nudg_coeff.resize(BdyVars::NumTypes);
470
471 vec_nudg_coeff[BdyVars::u].resize(lev+1);
472 vec_nudg_coeff[BdyVars::v].resize(lev+1);
473 vec_nudg_coeff[BdyVars::t].resize(lev+1);
474 vec_nudg_coeff[BdyVars::s].resize(lev+1);
475 vec_nudg_coeff[BdyVars::ubar].resize(lev+1);
476 vec_nudg_coeff[BdyVars::vbar].resize(lev+1);
477 vec_nudg_coeff[BdyVars::zeta].resize(lev+1);
478}
479
480/**
481 * @param[in ] lev level to operate on
482 * @param[in ] ba BoxArray for the level
483 * @param[in ] dm DistributionMapping for the level
484 */
485void REMORA::init_masks (int lev, const BoxArray& ba, const DistributionMapping& dm)
486{
487 BoxList bl2d = ba.boxList();
488 for (auto& b : bl2d) {
489 b.setRange(2,0);
490 }
491
492 BoxArray ba2d(std::move(bl2d));
493 vec_mskr[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW+1,NGROW+1,0)));
494 vec_msku[lev].reset(new MultiFab(convert(ba2d,IntVect(1,0,0)),dm,1,IntVect(NGROW+1,NGROW+1,0)));
495 vec_mskv[lev].reset(new MultiFab(convert(ba2d,IntVect(0,1,0)),dm,1,IntVect(NGROW+1,NGROW+1,0)));
496 vec_mskp[lev].reset(new MultiFab(convert(ba2d,IntVect(1,1,0)),dm,1,IntVect(NGROW+1,NGROW+1,0)));
497
498 vec_mskr[lev]->setVal(1.0_rt);
499 vec_msku[lev]->setVal(1.0_rt);
500 vec_mskv[lev]->setVal(1.0_rt);
501 vec_mskp[lev]->setVal(1.0_rt);
502}
503
504/**
505 * @param[in ] lev level to operate on
506 * @param[in ] ba BoxArray for the level
507 * @param[in ] dm DistributionMapping for the level
508 */
509void REMORA::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm)
510{
511 // ********************************************************************************************
512 // Initialize the boundary conditions
513 // ********************************************************************************************
514 physbcs[lev] = std::make_unique<REMORAPhysBCFunct> (lev, geom[lev], domain_bcs_type, domain_bcs_type_d,
516
517 BoxList bl2d = ba.boxList();
518 for (auto& b : bl2d) {
519 b.setRange(2,0);
520 }
521 BoxArray ba2d(std::move(bl2d));
522
523 BoxList bl1d = ba.boxList();
524 for (auto& b : bl1d) {
525 b.setRange(0,0);
526 b.setRange(1,0);
527 }
528 BoxArray ba1d(std::move(bl1d));
529
530 BoxArray ba_nd(ba);
531 ba_nd.surroundingNodes();
532 BoxArray ba_w(ba);
533 ba_w.surroundingNodes(2);
534
535 vec_z_phys_nd[lev].reset (new MultiFab(ba_nd,dm,1,IntVect(NGROW,NGROW,1))); // z at psi points (nodes) MIGHT NEED NGROW+1
536 vec_z_w[lev].reset (new MultiFab(convert(ba,IntVect(0,0,1)),dm,1,IntVect(NGROW+1,NGROW+1,0))); // z at w points (cell faces)
537 vec_z_r[lev].reset (new MultiFab(ba,dm,1,IntVect(NGROW+1,NGROW+1,0))); // z at r points (cell center)
538 vec_Hz[lev].reset (new MultiFab(ba,dm,1,IntVect(NGROW+1,NGROW+1,NGROW+1))); // like in ROMS, thickness of cell in z
539
540 vec_Huon[lev].reset (new MultiFab(convert(ba,IntVect(1,0,0)),dm,1,IntVect(NGROW,NGROW,0))); // mass flux for u component
541 vec_Hvom[lev].reset (new MultiFab(convert(ba,IntVect(0,1,0)),dm,1,IntVect(NGROW,NGROW,0))); // mass flux for v component
542
543 vec_Akv[lev].reset (new MultiFab(convert(ba,IntVect(0,0,1)),dm,1,IntVect(NGROW,NGROW,0))); // vertical mixing coefficient (.in)
544 vec_Akt[lev].reset (new MultiFab(convert(ba,IntVect(0,0,1)),dm,NCONS,IntVect(NGROW,NGROW,0))); // vertical mixing coefficient (.in)
545
546 // check dimensionality
547 vec_visc2_p[lev].reset(new MultiFab(ba,dm,1,IntVect(NGROW,NGROW,0))); // harmonic viscosity at psi points -- difference to 3d?
548 vec_visc2_r[lev].reset(new MultiFab(ba,dm,1,IntVect(NGROW,NGROW,0))); // harmonic viscosity at rho points
549 vec_diff2[lev].reset(new MultiFab(ba,dm,NCONS,IntVect(NGROW,NGROW,0))); // harmonic diffusivity temperature/salt
550
551 //2d, (incl advection terms and surface/bottom stresses, integral over the whole column, k=0)
552 vec_rufrc[lev].reset(new MultiFab(convert(ba2d,IntVect(1,0,0)),dm,2,IntVect(NGROW,NGROW,0)));
553 vec_rvfrc[lev].reset(new MultiFab(convert(ba2d,IntVect(0,1,0)),dm,2,IntVect(NGROW,NGROW,0))); //2d, same as above but v
554
555 vec_sustr[lev].reset(new MultiFab(convert(ba2d,IntVect(1,0,0)),dm,1,IntVect(NGROW,NGROW,0))); //2d, surface stress
556 vec_svstr[lev].reset(new MultiFab(convert(ba2d,IntVect(0,1,0)),dm,1,IntVect(NGROW,NGROW,0))); //2d
557
559 //2d, linear drag coefficient [m/s], defined at rho, somehow related to rdrg
560 vec_rdrag[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW,NGROW,0)));
562 vec_rdrag2[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW,NGROW,0)));
563 }
564
567 vec_ZoBot[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW,NGROW,0)));
568 }
569
570 vec_bustr[lev].reset(new MultiFab(convert(ba2d,IntVect(1,0,0)),dm,1,IntVect(NGROW,NGROW,0))); //2d, bottom stress
571 vec_bvstr[lev].reset(new MultiFab(convert(ba2d,IntVect(0,1,0)),dm,1,IntVect(NGROW,NGROW,0)));
572
573 //all 2d -- all associated with the 2D advance
574 //2d DU: sum(height[incl free surface?] * u)
575 vec_DU_avg1[lev].reset(new MultiFab(convert(ba2d,IntVect(1,0,0)),dm,1,IntVect(NGROW,NGROW,0)));
576
577 //2d like above, but correct(or)?
578 vec_DU_avg2[lev].reset(new MultiFab(convert(ba2d,IntVect(1,0,0)),dm,1,IntVect(NGROW,NGROW,0)));
579
580 vec_DV_avg1[lev].reset(new MultiFab(convert(ba2d,IntVect(0,1,0)),dm,1,IntVect(NGROW,NGROW,0)));
581 vec_DV_avg2[lev].reset(new MultiFab(convert(ba2d,IntVect(0,1,0)),dm,1,IntVect(NGROW,NGROW,0)));
582
583 vec_rubar[lev].reset(new MultiFab(convert(ba2d,IntVect(1,0,0)),dm,4,IntVect(NGROW,NGROW,0))); // 2d RHS ubar
584 vec_rvbar[lev].reset(new MultiFab(convert(ba2d,IntVect(0,1,0)),dm,4,IntVect(NGROW,NGROW,0)));
585 vec_rzeta[lev].reset(new MultiFab(ba2d,dm,4,IntVect(NGROW,NGROW,0))); // 2d RHS zeta
586
587 // starts off kind of like a depth-averaged u, but exists at more points and more timesteps (b/c fast 2D update) than full u
588 vec_zeta[lev].reset(new MultiFab(ba2d,dm,3,IntVect(NGROW+1,NGROW+1,0))); // 2d free surface
589
590 vec_pm[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW+1,NGROW+2,0)));
591 vec_pn[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW+2,NGROW+1,0)));
592 vec_fcor[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW+1,NGROW+1,0)));
593
594 vec_xr[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW+1,NGROW+1,0)));
595 vec_yr[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW+1,NGROW+1,0)));
596
597 vec_xu[lev].reset(new MultiFab(convert(ba2d,IntVect(1,0,0)),dm,1,IntVect(NGROW,NGROW,0)));
598 vec_yu[lev].reset(new MultiFab(convert(ba2d,IntVect(1,0,0)),dm,1,IntVect(NGROW,NGROW,0)));
599 vec_xv[lev].reset(new MultiFab(convert(ba2d,IntVect(0,1,0)),dm,1,IntVect(NGROW,NGROW,0)));
600 vec_yv[lev].reset(new MultiFab(convert(ba2d,IntVect(0,1,0)),dm,1,IntVect(NGROW,NGROW,0)));
601 vec_xp[lev].reset(new MultiFab(convert(ba2d,IntVect(1,1,0)),dm,1,IntVect(NGROW,NGROW,0)));
602 vec_yp[lev].reset(new MultiFab(convert(ba2d,IntVect(1,1,0)),dm,1,IntVect(NGROW,NGROW,0)));
603
604
605 // tempstore, saltstore, etc
606 vec_sstore[lev].reset(new MultiFab(ba,dm,NCONS,IntVect(NGROW,NGROW,0)));
607
608 vec_rhoS[lev].reset(new MultiFab(ba,dm,1,IntVect(NGROW,NGROW,0)));
609 vec_rhoA[lev].reset(new MultiFab(ba,dm,1,IntVect(NGROW,NGROW,0)));
610 vec_bvf[lev].reset(new MultiFab(convert(ba,IntVect(0,0,1)),dm,1,IntVect(NGROW,NGROW,0)));
611
612 vec_tke[lev].reset(new MultiFab(convert(ba,IntVect(0,0,1)),dm,3,IntVect(NGROW,NGROW,0)));
613 vec_gls[lev].reset(new MultiFab(convert(ba,IntVect(0,0,1)),dm,3,IntVect(NGROW,NGROW,0)));
614 vec_Lscale[lev].reset(new MultiFab(convert(ba,IntVect(0,0,1)),dm,1,IntVect(NGROW,NGROW,0)));
615 vec_Akk[lev].reset(new MultiFab(convert(ba,IntVect(0,0,1)),dm,1,IntVect(NGROW,NGROW,0)));
616 vec_Akp[lev].reset(new MultiFab(convert(ba,IntVect(0,0,1)),dm,1,IntVect(NGROW,NGROW,0)));
617
618 // surface/bottom tracer fluxes for update
619 vec_stflx[lev].reset(new MultiFab(ba2d,dm,NCONS,IntVect(NGROW,NGROW,0)));
620 vec_btflx[lev].reset(new MultiFab(ba2d,dm,NCONS,IntVect(NGROW,NGROW,0)));
621 // surface/bottom tracer fluxes to be filled by inputs
622 vec_stflux[lev].reset(new MultiFab(ba2d,dm,NCONS,IntVect(NGROW,NGROW,0)));
623 vec_btflux[lev].reset(new MultiFab(ba2d,dm,NCONS,IntVect(NGROW,NGROW,0)));
624
626 vec_uwind[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW,NGROW,0))); //2d, surface wind u
627 vec_vwind[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW,NGROW,0))); //2d, surface wind v
628 vec_alpha[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW,NGROW,0)));
629 vec_beta[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW,NGROW,0)));
630 vec_lrflx[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW,NGROW,0)));
631 vec_lhflx[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW,NGROW,0)));
632 vec_shflx[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW,NGROW,0)));
633 vec_rain[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW,NGROW,0)));
634 vec_evap[lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW,NGROW,0)));
635 vec_lhflx[lev]->setVal(0.0_rt);
636 vec_shflx[lev]->setVal(0.0_rt);
637 vec_rain[lev]->setVal(solverChoice.rain);
638 }
639
641 vec_river_position[lev].reset(new iMultiFab(ba2d,dm,1,IntVect(NGROW,NGROW,0)));
642 vec_river_position[lev]->setVal(-1);
643 }
644
645 vec_nudg_coeff[BdyVars::u][lev].reset(new MultiFab(ba,dm,1,IntVect(NGROW,NGROW,0)));
646 vec_nudg_coeff[BdyVars::v][lev].reset(new MultiFab(ba,dm,1,IntVect(NGROW,NGROW,0)));
647 vec_nudg_coeff[BdyVars::t][lev].reset(new MultiFab(ba,dm,1,IntVect(NGROW,NGROW,0)));
648 vec_nudg_coeff[BdyVars::s][lev].reset(new MultiFab(ba,dm,1,IntVect(NGROW,NGROW,0)));
649 vec_nudg_coeff[BdyVars::ubar][lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW,NGROW,0)));
650 vec_nudg_coeff[BdyVars::vbar][lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW,NGROW,0)));
651 vec_nudg_coeff[BdyVars::zeta][lev].reset(new MultiFab(ba2d,dm,1,IntVect(NGROW,NGROW,0)));
652
653 set_weights(lev);
654
655 vec_DU_avg1[lev]->setVal(0.0_rt);
656 vec_DU_avg2[lev]->setVal(0.0_rt);
657 vec_DV_avg1[lev]->setVal(0.0_rt);
658 vec_DV_avg2[lev]->setVal(0.0_rt);
659 vec_rubar[lev]->setVal(0.0_rt);
660 vec_rvbar[lev]->setVal(0.0_rt);
661 vec_rzeta[lev]->setVal(0.0_rt);
662
663 // Initialize these vars even if we aren't using GLS to
664 // avoid issues on e.g. checkpoint
665 vec_tke[lev]->setVal(solverChoice.gls_Kmin);
666 vec_gls[lev]->setVal(solverChoice.gls_Pmin);
667 vec_Lscale[lev]->setVal(0.0_rt);
668 vec_Akk[lev]->setVal(solverChoice.Akk_bak);
669 vec_Akp[lev]->setVal(solverChoice.Akp_bak);
670
671 vec_stflx[lev]->setVal(0.0_rt);
672 vec_btflx[lev]->setVal(0.0_rt);
673 vec_stflux[lev]->setVal(0.0_rt);
674 vec_btflux[lev]->setVal(0.0_rt);
675
676 // NOTE: Used to set vec_pm and vec_pn to 1e34 here to make foextrap work
677 // when init_type = real. However, this does not appear to be necessary so removing
678
679 // Set initial linear drag coefficient
681 vec_rdrag[lev]->setVal(solverChoice.rdrag);
683 vec_rdrag2[lev]->setVal(solverChoice.rdrag2);
684 }
685
688 vec_ZoBot[lev]->setVal(solverChoice.Zob);
689 }
690
691
692 // ********************************************************************************************
693 // Create the REMORAFillPatcher object
694 // ********************************************************************************************
695 if (lev > 0 && cf_width >= 0) {
698 }
699}
700
701/**
702 * Delete level data. Overrides the pure virtual function in AmrCore
703 *
704 * @param[in ] lev level to operate on
705 */
706void
708{
709 delete cons_new[lev]; delete xvel_new[lev]; delete yvel_new[lev]; delete zvel_new[lev];
710 delete cons_old[lev]; delete xvel_old[lev]; delete yvel_old[lev]; delete zvel_old[lev];
711}
712
713/**
714 * @param[in ] lev level to operate on
715 */
716void
718{
719 AMREX_ASSERT(solverChoice.ic_type == IC_Type::analytic);
721 const auto dxi = Geom(lev).InvCellSize();
722 vec_pm[lev]->setVal(dxi[0]); vec_pm[lev]->FillBoundary(geom[lev].periodicity());
723 vec_pn[lev]->setVal(dxi[1]); vec_pn[lev]->FillBoundary(geom[lev].periodicity());
725 prob->init_analytic_grid_scale(lev, Geom(lev), solverChoice, *this, *vec_pm[lev].get(), *vec_pn[lev].get());
726 vec_pm[lev]->FillBoundary(geom[lev].periodicity());
727 vec_pn[lev]->FillBoundary(geom[lev].periodicity());
728 }
729
730 for ( MFIter mfi(*vec_xr[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi )
731 {
732 Array4<const Real> const& pm = vec_pm[lev]->const_array(mfi);
733 Array4<const Real> const& pn = vec_pn[lev]->const_array(mfi);
734 Array4<Real> const& xr = vec_xr[lev]->array(mfi);
735 Array4<Real> const& yr = vec_yr[lev]->array(mfi);
736 Array4<Real> const& xu = vec_xu[lev]->array(mfi);
737 Array4<Real> const& yu = vec_yu[lev]->array(mfi);
738 Array4<Real> const& xv = vec_xv[lev]->array(mfi);
739 Array4<Real> const& yv = vec_yv[lev]->array(mfi);
740 Array4<Real> const& xp = vec_xp[lev]->array(mfi);
741 Array4<Real> const& yp = vec_yp[lev]->array(mfi);
742
743 Box bx = mfi.growntilebox(IntVect(NGROW,NGROW,0));
744 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int)
745 {
746 xr(i,j,0) = (i + 0.5_rt) / pm(i,j,0);
747 yr(i,j,0) = (j + 0.5_rt) / pn(i,j,0);
748 });
749
750 ParallelFor(grow(convert(bx,IntVect(1,0,0)),IntVect(-1,0,0)), [=] AMREX_GPU_DEVICE (int i, int j, int)
751 {
752 xu(i,j,0) = i / pm(i,j,0);
753 yu(i,j,0) = (j + 0.5_rt) / pn(i,j,0);
754 });
755
756 ParallelFor(grow(convert(bx,IntVect(0,1,0)),IntVect(0,-1,0)), [=] AMREX_GPU_DEVICE (int i, int j, int)
757 {
758 xv(i,j,0) = (i + 0.5_rt) / pm(i,j,0);
759 yv(i,j,0) = j / pn(i,j,0);
760 });
761
762 ParallelFor(grow(convert(bx,IntVect(1,1,0)),IntVect(-1,-1,0)), [=] AMREX_GPU_DEVICE (int i, int j, int)
763 {
764 xp(i,j,0) = i / pm(i,j,0);
765 yp(i,j,0) = j / pn(i,j,0);
766 });
767 }
768}
769
770/**
771 * @param[in ] lev level to operate on
772 */
773void
775{
776 BL_PROFILE("REMORA::set_zeta_to_Ztavg()");
777 std::unique_ptr<MultiFab>& mf_zeta = vec_zeta[lev];
778 std::unique_ptr<MultiFab>& mf_Zt_avg1 = vec_Zt_avg1[lev];
780 for ( MFIter mfi(*vec_zeta[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi )
781 {
782 Array4<Real> const& Zt_avg1 = (mf_Zt_avg1)->array(mfi);
783 Array4<const Real> const& evap = vec_evap[lev]->const_array(mfi);
784 Array4<const Real> const& rain = vec_rain[lev]->const_array(mfi);
785
786 Box bx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));// bx2.grow(IntVect(NGROW,NGROW,0));
787
788 Real cff = dt[lev] / rhow;
789
790 ParallelFor(bx2, [=] AMREX_GPU_DEVICE (int i, int j, int )
791 {
792 Zt_avg1(i,j,0) = Zt_avg1(i,j,0) - (evap(i,j,0) - rain(i,j,0)) * cff;
793 });
794 }
795 }
796 Gpu::streamSynchronize();
797
798 vec_Zt_avg1[lev]->FillBoundary(geom[lev].periodicity());
799
800 for ( MFIter mfi(*vec_zeta[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi )
801 {
802 Box bx3 = mfi.tilebox(); bx3.grow(IntVect(NGROW+1,NGROW+1,0));
803 Array4<Real> const& zeta = mf_zeta->array(mfi);
804 Array4<Real> const& Zt_avg1 = (mf_Zt_avg1)->array(mfi);
805
806 ParallelFor(bx3, 3, [=] AMREX_GPU_DEVICE (int i, int j, int , int n)
807 {
808 zeta(i,j,0,n) = Zt_avg1(i,j,0);
809 });
810 }
811}
812
813/**
814 * @param[in ] lev level to operate on
815 */
816void
818{
819 for ( MFIter mfi(*vec_mskr[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi )
820 {
821 Array4<const Real> const& mskr = vec_mskr[lev]->const_array(mfi);
822 Array4< Real> const& mskp = vec_mskp[lev]->array(mfi);
823
824 Box bx = mfi.tilebox(); bx.grow(IntVect(1,1,0)); bx.makeSlab(2,0);
825
826 Real cff1 = 1.0_rt;
827 Real cff2 = 2.0_rt;
828
829 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int)
830 {
831 if ((mskr(i-1,j,0) > 0.5) and (mskr(i,j,0) > 0.5) and (mskr(i-1,j-1,0) > 0.5) and (mskr(i,j-1,0) > 0.5)) {
832 mskp(i,j,0) = 1.0_rt;
833 } else if ((mskr(i-1,j,0) < 0.5) and (mskr(i,j,0) > 0.5) and (mskr(i-1,j-1,0) > 0.5) and (mskr(i,j-1,0) > 0.5)) {
834 mskp(i,j,0) = cff1;
835 } else if ((mskr(i-1,j,0) > 0.5) and (mskr(i,j,0) < 0.5) and (mskr(i-1,j-1,0) > 0.5) and (mskr(i,j-1,0) > 0.5)) {
836 mskp(i,j,0) = cff1;
837 } else if ((mskr(i-1,j,0) > 0.5) and (mskr(i,j,0) > 0.5) and (mskr(i-1,j-1,0) < 0.5) and (mskr(i,j-1,0) > 0.5)) {
838 mskp(i,j,0) = cff1;
839 } else if ((mskr(i-1,j,0) > 0.5) and (mskr(i,j,0) > 0.5) and (mskr(i-1,j-1,0) > 0.5) and (mskr(i,j-1,0) < 0.5)) {
840 mskp(i,j,0) = cff1;
841 } else if ((mskr(i-1,j,0) > 0.5) and (mskr(i,j,0) < 0.5) and (mskr(i-1,j-1,0) > 0.5) and (mskr(i,j-1,0) < 0.5)) {
842 mskp(i,j,0) = cff2;
843 } else if ((mskr(i-1,j,0) < 0.5) and (mskr(i,j,0) > 0.5) and (mskr(i-1,j-1,0) < 0.5) and (mskr(i,j-1,0) > 0.5)) {
844 mskp(i,j,0) = cff2;
845 } else if ((mskr(i-1,j,0) > 0.5) and (mskr(i,j,0) > 0.5) and (mskr(i-1,j-1,0) < 0.5) and (mskr(i,j-1,0) < 0.5)) {
846 mskp(i,j,0) = cff2;
847 } else if ((mskr(i-1,j,0) < 0.5) and (mskr(i,j,0) < 0.5) and (mskr(i-1,j-1,0) > 0.5) and (mskr(i,j-1,0) > 0.5)) {
848 mskp(i,j,0) = cff2;
849 } else {
850 mskp(i,j,0) = 0.0_rt;
851 }
852
853 });
854 }
855}
856/**
857 * @param[in ] lev level to operate on
858 */
859void
861{
862 for ( MFIter mfi(*vec_mskr[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi )
863 {
864 Array4<const Real> const& mskr = vec_mskr[lev]->const_array(mfi);
865 Array4< Real> const& msku = vec_msku[lev]->array(mfi);
866 Array4< Real> const& mskv = vec_mskv[lev]->array(mfi);
867 Array4< Real> const& mskp = vec_mskp[lev]->array(mfi);
868
869 Box bx = mfi.tilebox(); bx.grow(IntVect(1,1,0)); bx.makeSlab(2,0);
870
871 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int)
872 {
873 msku(i,j,0) = mskr(i-1,j ,0) * mskr(i,j,0);
874 mskv(i,j,0) = mskr(i ,j-1,0) * mskr(i,j,0);
875 mskp(i,j,0) = mskr(i-1,j-1,0) * mskr(i,j,0) * mskr(i-1,j,0) * mskr(i,j-1,0);
876 });
877 }
878}
constexpr amrex::Real rhow
#define NGROW
#define NCONS
void FillCoarsePatchPC(int lev, amrex::Real time, amrex::MultiFab *mf_fine, amrex::MultiFab *mf_crse, const int bccomp, const int bdy_var_type=BdyVars::null, const int icomp=0, const bool fill_all=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 an entire multifab by interpolating from the coarser level using the piecewise constant interpol...
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_rv2d
v velocity RHS (2D, includes horizontal and vertical advection)
Definition REMORA.H:256
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_evap
evaporation rate [kg/m^2/s]
Definition REMORA.H:327
amrex::Vector< amrex::BCRec > domain_bcs_type
vector (over BCVars) of BCRecs
Definition REMORA.H:1179
void calculate_nodal_masks(int lev)
Calculate u-, v-, and psi-point masks based on rho-point masks after analytic initialization.
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_fcor
coriolis factor (2D)
Definition REMORA.H:377
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_btflux
Bottom tracer flux; input arrays.
Definition REMORA.H:322
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_rubar
barotropic x velocity for the RHS (2D)
Definition REMORA.H:350
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_h
Bathymetry data (2D, positive valued, h in ROMS)
Definition REMORA.H:241
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pm
horizontal scaling factor: 1 / dx (2D)
Definition REMORA.H:372
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_ZoBot
Bottom roughness length [m], defined at rho points.
Definition REMORA.H:334
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_DU_avg2
correct time average of barotropic x velocity flux for coupling (2D)
Definition REMORA.H:344
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_lrflx
longwave radiation
Definition REMORA.H:309
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_yv
y_grid on v-points (2D)
Definition REMORA.H:392
amrex::Vector< amrex::MultiFab * > cons_new
multilevel data container for current step's scalar data: temperature, salinity, passive scalar
Definition REMORA.H:229
virtual void MakeNewLevelFromCoarse(int lev, amrex::Real time, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm) override
Make a new level using provided BoxArray and DistributionMapping and fill with interpolated coarse le...
void FillCoarsePatch(int lev, amrex::Real time, amrex::MultiFab *mf_fine, amrex::MultiFab *mf_crse, const int bccomp, const int bdy_var_type=BdyVars::null, const int icomp=0, const bool fill_all=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 an entire multifab by interpolating from the coarser level
void stretch_transform(int lev)
Calculate vertical stretched coordinates.
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_vwind
Wind in the v direction, defined at rho-points.
Definition REMORA.H:306
std::unique_ptr< ProblemBase > prob
Pointer to container of analytical functions for problem definition.
Definition REMORA.H:1143
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskr
land/sea mask at cell centers (2D)
Definition REMORA.H:363
void Construct_REMORAFillPatchers(int lev)
Construct FillPatchers.
Definition REMORA.cpp:373
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_rain
precipitation rate [kg/m^2/s]
Definition REMORA.H:325
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_tke
Turbulent kinetic energy.
Definition REMORA.H:420
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_stflx
Surface tracer flux; working arrays.
Definition REMORA.H:316
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_gls
Turbulent generic length scale.
Definition REMORA.H:422
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_sustr
Surface stress in the u direction.
Definition REMORA.H:299
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_yp
y_grid on psi-points (2D)
Definition REMORA.H:397
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_xr
x_grid on rho points (2D)
Definition REMORA.H:380
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_ru2d
u velocity RHS (2D, includes horizontal and vertical advection)
Definition REMORA.H:254
virtual void ClearLevel(int lev) override
Delete level data Overrides the pure virtual function in AmrCore.
amrex::Vector< amrex::MultiFab * > zvel_new
multilevel data container for current step's z velocities (largely unused; W stored separately)
Definition REMORA.H:235
AMREX_FORCE_INLINE int ComputeGhostCells(const int &spatial_order)
Helper function to determine number of ghost cells.
Definition REMORA.H:1377
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_sstore
additional scratch space for calculations on temp, salt, etc
Definition REMORA.H:400
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_xv
x_grid on v-points (2D)
Definition REMORA.H:390
void init_only(int lev, amrex::Real time)
Init (NOT restart or regrid)
Definition REMORA.cpp:739
void init_set_vmix(int lev)
Initialize vertical mixing coefficients from file or analytic.
Definition REMORA.cpp:615
void set_grid_scale(int lev)
Set pm and pn arrays on level lev. Only works if using analytic initialization.
void set_coriolis(int lev)
Initialize Coriolis factor from file or analytic.
Definition REMORA.cpp:591
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Lscale
Vertical mixing turbulent length scale.
Definition REMORA.H:424
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Hz
Width of cells in the vertical (z-) direction (3D, Hz in ROMS)
Definition REMORA.H:244
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akt
Vertical diffusion coefficient (3D)
Definition REMORA.H:264
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_msku
land/sea mask at x-faces (2D)
Definition REMORA.H:365
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_rvfrc
v velocity RHS, integrated, including advection and bottom/surface stresses (2D)
Definition REMORA.H:260
amrex::Vector< amrex::MultiFab * > xvel_old
multilevel data container for last step's x velocities (u in ROMS)
Definition REMORA.H:222
amrex::Vector< amrex::MultiFab * > yvel_new
multilevel data container for current step's y velocities (v in ROMS)
Definition REMORA.H:233
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_uwind
Wind in the u direction, defined at rho-points.
Definition REMORA.H:304
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_rufrc
u velocity RHS, integrated, including advection and bottom/surface stresses (2D)
Definition REMORA.H:258
void Define_REMORAFillPatchers(int lev)
Define FillPatchers.
Definition REMORA.cpp:422
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_shflx
sensible heat flux
Definition REMORA.H:313
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_visc2_p
Harmonic viscosity defined on the psi points (corners of horizontal grid cells)
Definition REMORA.H:266
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_rvbar
barotropic y velocity for the RHS (2D)
Definition REMORA.H:352
amrex::Vector< amrex::MultiFab * > zvel_old
multilevel data container for last step's z velocities (largely unused; W stored separately)
Definition REMORA.H:226
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_z_r
z coordinates at rho points (cell centers)
Definition REMORA.H:273
amrex::Vector< amrex::MultiFab * > xvel_new
multilevel data container for current step's x velocities (u in ROMS)
Definition REMORA.H:231
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_lhflx
latent heat flux
Definition REMORA.H:311
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskp
land/sea mask at cell corners (2D)
Definition REMORA.H:369
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_bvf
Brunt-Vaisala frequency (3D)
Definition REMORA.H:407
virtual void RemakeLevel(int lev, amrex::Real time, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm) override
Remake an existing level using provided BoxArray and DistributionMapping and fill with existing fine ...
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_mskv
land/sea mask at y-faces (2D)
Definition REMORA.H:367
void init_masks(int lev, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm)
Allocate MultiFabs for masks.
amrex::Vector< std::unique_ptr< REMORAPhysBCFunct > > physbcs
Vector (over level) of functors to apply physical boundary conditions.
Definition REMORA.H:1167
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< std::unique_ptr< amrex::MultiFab > > vec_rhoS
density perturbation
Definition REMORA.H:403
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_visc2_r
Harmonic viscosity defined on the rho points (centers)
Definition REMORA.H:268
void update_mskp(int lev)
Set psi-point mask to be consistent with rho-point mask.
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_svstr
Surface stress in the v direction.
Definition REMORA.H:301
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Huon
u-volume flux (3D)
Definition REMORA.H:246
amrex::Vector< amrex::MultiFab * > yvel_old
multilevel data container for last step's y velocities (v in ROMS)
Definition REMORA.H:224
void init_stuff(int lev, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm)
Allocate MultiFabs for state and evolution variables.
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_rhoA
vertically-averaged density
Definition REMORA.H:405
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_DV_avg1
time average of barotropic y velocity flux
Definition REMORA.H:346
amrex::Vector< amrex::Real > t_new
new time at each level
Definition REMORA.H:1157
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1280
void resize_stuff(int lev)
Resize variable containers to accommodate data on levels 0 to max_lev.
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > vec_river_position
iMultiFab for river positions; contents are indices of rivers
Definition REMORA.H:1095
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akk
Turbulent kinetic energy vertical diffusion coefficient.
Definition REMORA.H:426
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_rdrag2
Quadratic drag coefficient [unitless], defined at rho points.
Definition REMORA.H:332
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_ru
u velocity RHS (3D, includes horizontal and vertical advection)
Definition REMORA.H:250
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_xp
x_grid on psi-points (2D)
Definition REMORA.H:395
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_zeta
free surface height (2D)
Definition REMORA.H:360
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_vbar
barotropic y velocity (2D)
Definition REMORA.H:358
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_DU_avg1
time average of barotropic x velocity flux (2D)
Definition REMORA.H:342
void set_zeta_to_Ztavg(int lev)
Set zeta components to be equal to time-averaged Zt_avg1.
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_alpha
Thermal expansion coefficient (3D)
Definition REMORA.H:409
amrex::Array< amrex::Array< amrex::Real, AMREX_SPACEDIM *2 >, AMREX_SPACEDIM+NCONS+8 > m_bc_extdir_vals
Array holding the Dirichlet values at walls which need them.
Definition REMORA.H:1187
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_ubar
barotropic x velocity (2D)
Definition REMORA.H:356
amrex::Vector< amrex::MultiFab * > cons_old
multilevel data container for last step's scalar data: temperature, salinity, passive scalar
Definition REMORA.H:220
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_bustr
Bottom stress in the u direction.
Definition REMORA.H:337
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_DV_avg2
correct time average of barotropic y velocity flux for coupling (2D)
Definition REMORA.H:348
void set_hmixcoef(int lev)
Initialize horizontal mixing coefficients.
Definition REMORA.cpp:646
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_yu
y_grid on u-points (2D)
Definition REMORA.H:387
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_bvstr
Bottom stress in the v direction.
Definition REMORA.H:339
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_rzeta
free surface height for the RHS (2D)
Definition REMORA.H:354
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_z_phys_nd
z coordinates at psi points (cell nodes)
Definition REMORA.H:293
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_xu
x_grid on u-points (2D)
Definition REMORA.H:385
amrex::Vector< amrex::Vector< std::unique_ptr< amrex::MultiFab > > > vec_nudg_coeff
Climatology nudging coefficients.
Definition REMORA.H:431
void set_weights(int lev)
Set weights for averaging 3D variables to 2D.
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_pn
horizontal scaling factor: 1 / dy (2D)
Definition REMORA.H:374
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akv
Vertical viscosity coefficient (3D)
Definition REMORA.H:262
virtual void MakeNewLevelFromScratch(int lev, amrex::Real time, const amrex::BoxArray &ba, const amrex::DistributionMapping &dm) override
Make a new level from scratch using provided BoxArray and DistributionMapping. Only used during initi...
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_stflux
Surface tracer flux; input arrays.
Definition REMORA.H:318
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_rdrag
Linear drag coefficient [m/s], defined at rho points.
Definition REMORA.H:330
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Zt_avg1
Average of the free surface, zeta (2D)
Definition REMORA.H:296
std::string restart_chkfile
If set, restart from this checkpoint file.
Definition REMORA.H:1225
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_rv
v velocity RHS (3D, includes horizontal and vertical advection)
Definition REMORA.H:252
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_btflx
Bottom tracer flux; working arrays.
Definition REMORA.H:320
int cf_width
Nudging width at coarse-fine interface.
Definition REMORA.H:1105
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_beta
Saline contraction coefficient (3D)
Definition REMORA.H:411
amrex::Vector< amrex::Real > t_old
old time at each level
Definition REMORA.H:1159
amrex::Vector< amrex::Real > dt
time step at each level
Definition REMORA.H:1161
amrex::Gpu::DeviceVector< amrex::BCRec > domain_bcs_type_d
GPU vector (over BCVars) of BCRecs.
Definition REMORA.H:1181
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_z_w
z coordinates at w points (faces between z-cells)
Definition REMORA.H:276
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_yr
y_grid on rho points (2D)
Definition REMORA.H:382
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Hvom
v-volume flux (3D)
Definition REMORA.H:248
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akp
Turbulent length scale vertical diffusion coefficient.
Definition REMORA.H:428
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_diff2
Harmonic diffusivity for temperature / salinity.
Definition REMORA.H:270
amrex::Real rdrag2
amrex::Real Akk_bak
BottomStressType bottom_stress_type
amrex::Real gls_Kmin
amrex::Real rdrag
VertMixingType vert_mixing_type
GridScaleType grid_scale_type
amrex::Real gls_Pmin
amrex::Real Akp_bak