REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_prestep.cpp
Go to the documentation of this file.
1#include <REMORA.H>
2
3using namespace amrex;
4
5/**
6 * @param[in ] lev level to operate on
7 * @param[in ] mf_uold u-velocity from last time step
8 * @param[in ] mf_vold v-velocity from last time step
9 * @param[ out] mf_u u-velocity at current time step
10 * @param[ out] mf_v v-velocity at current time step
11 * @param[inout] mf_ru u-velocity RHS at current time step
12 * @param[inout] mf_rv v-velocity RHS at current time step
13 * @param[in ] S_old scalar variables at last time step
14 * @param[inout] S_new scalar variables at current time step
15 * @param[in ] mf_W vertical velocity
16 * @param mf_DC temporary variable container
17 * @param[in ] mf_z_r z coordinates at rho points (cell centers)
18 * @param[in ] mf_z_w z coordinates at w points
19 * @param[in ] mf_h bathymetry
20 * @param[in ] mf_pm 1/dx
21 * @param[in ] mf_pn 1/dy
22 * @param[in ] mf_sustr u-direction surface momentum flux
23 * @param[in ] mf_svstr v-direction surface momentum flux
24 * @param[in ] mf_bustr u-direction bottom stress
25 * @param[in ] mf_bvstr v-direction bottom stress
26 * @param[in ] mf_msku land-sea mask on u-points
27 * @param[in ] mf_mskv land-sea mask on v-points
28 * @param[in ] iic which time step we're on
29 * @param[in ] ntfirst what is the first time step?
30 * @param[in ] nnew index of time step to update
31 * @param[in ] nstp index of last time step
32 * @param[in ] nrhs index of RHS component
33 * @param[in ] N number of vertical levels
34 * @param[in ] dt_lev time step at this level
35 */
36
37void
39 MultiFab& mf_uold, MultiFab& mf_vold,
40 MultiFab& mf_u, MultiFab& mf_v,
41 MultiFab* mf_ru,
42 MultiFab* mf_rv,
43 MultiFab& S_old, MultiFab& S_new,
44 MultiFab& mf_W, MultiFab& mf_DC,
45 const MultiFab* mf_z_r,
46 const MultiFab* mf_z_w,
47 const MultiFab* mf_h,
48 const MultiFab* mf_pm,
49 const MultiFab* mf_pn,
50 const MultiFab* mf_sustr,
51 const MultiFab* mf_svstr,
52 const MultiFab* mf_bustr,
53 const MultiFab* mf_bvstr,
54 const MultiFab* mf_msku,
55 const MultiFab* mf_mskv,
56 const int iic, const int ntfirst,
57 const int nnew, int nstp, int nrhs,
58 int N, const Real dt_lev)
59{
60 BL_PROFILE("REMORA::prestep()");
61 const BoxArray& ba = S_old.boxArray();
62 const DistributionMapping& dm = S_old.DistributionMap();
63
64 // Maybe not the best way to do this, but need to cache salt and temp since
65 // they get rewritten by prestep_t
66 MultiFab mf_saltcache(ba,dm,1,IntVect(NGROW,NGROW,0));
67 MultiFab mf_tempcache(ba,dm,1,IntVect(NGROW,NGROW,0));
68
69 MultiFab mf_scalarcache(ba,dm,NCONS,IntVect(NGROW,NGROW,0));
70 MultiFab::Copy(mf_scalarcache,S_new,0,0,NCONS,IntVect(NGROW,NGROW,0));
71
72 MultiFab::Copy(mf_saltcache,S_new,Salt_comp,0,1,IntVect(NGROW,NGROW,0));
73 MultiFab::Copy(mf_tempcache,S_new,Temp_comp,0,1,IntVect(NGROW,NGROW,0));
74
75 for ( MFIter mfi(S_new, TilingIfNotGPU()); mfi.isValid(); ++mfi )
76 {
77 Array4<Real> const& DC = mf_DC.array(mfi);
78 Array4<Real> const& Akv = vec_Akv[lev]->array(mfi);
79 Array4<Real> const& Akt = vec_Akt[lev]->array(mfi);
80 Array4<Real> const& Hz = vec_Hz[lev]->array(mfi);
81 Array4<Real> const& Huon = vec_Huon[lev]->array(mfi);
82 Array4<Real> const& Hvom = vec_Hvom[lev]->array(mfi);
83 Array4<Real const> const& uold = mf_uold.const_array(mfi);
84 Array4<Real const> const& vold = mf_vold.const_array(mfi);
85 Array4<Real> const& u = (mf_u).array(mfi);
86 Array4<Real> const& v = (mf_v).array(mfi);
87
88 Array4<Real const> const& z_r = mf_z_r->const_array(mfi);
89 Array4<Real const> const& z_w = mf_z_w->const_array(mfi);
90 Array4<Real const> const& h = mf_h->const_array(mfi);
91 Array4<Real const> const& pm = mf_pm->const_array(mfi);
92 Array4<Real const> const& pn = mf_pn->const_array(mfi);
93
94 Array4<Real > const& ru = mf_ru->array(mfi);
95 Array4<Real > const& rv = mf_rv->array(mfi);
96 Array4<Real > const& W = (mf_W).array(mfi);
97 Array4<Real const> const& sustr = mf_sustr->const_array(mfi);
98 Array4<Real const> const& svstr = mf_svstr->const_array(mfi);
99 Array4<Real const> const& bustr = mf_bustr->const_array(mfi);
100 Array4<Real const> const& bvstr = mf_bvstr->const_array(mfi);
101 Array4<Real const> const& msku = mf_msku->const_array(mfi);
102 Array4<Real const> const& mskv = mf_mskv->const_array(mfi);
103
104 Real lambda = 1.0_rt;
105
106 Box bx = mfi.tilebox();
107 Box gbx = mfi.growntilebox();
108 Box gbx1 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,0));
109 Box gbx2 = mfi.growntilebox(IntVect(NGROW,NGROW,0));
110 //Box gbx11 = mfi.growntilebox(IntVect(NGROW-1,NGROW-1,NGROW-1));
111
112 //copy the tilebox
113 Box tbxp1 = bx;
114 Box tbxp11 = bx;
115 Box tbxp2 = bx;
116
117 //TODO: adjust for tiling
118 //Box gbx3uneven(IntVect(AMREX_D_DECL(bx.smallEnd(0)-3,bx.smallEnd(1)-3,bx.smallEnd(2))),
119 // IntVect(AMREX_D_DECL(bx.bigEnd(0)+2,bx.bigEnd(1)+2,bx.bigEnd(2))));
120 //Box gbx2uneven(IntVect(AMREX_D_DECL(bx.smallEnd(0)-2,bx.smallEnd(1)-2,bx.smallEnd(2))),
121 // IntVect(AMREX_D_DECL(bx.bigEnd(0)+1,bx.bigEnd(1)+1,bx.bigEnd(2))));
122 //make only gbx be grown to match multifabs
123 tbxp2.grow(IntVect(NGROW,NGROW,0));
124 tbxp1.grow(IntVect(NGROW-1,NGROW-1,0));
125 tbxp11.grow(IntVect(NGROW-1,NGROW-1,NGROW-1));
126
127 Box bxD = bx;
128 bxD.makeSlab(2,0);
129 Box gbx1D = gbx1;
130 gbx1D.makeSlab(2,0);
131 Box gbx2D = gbx2;
132 gbx2D.makeSlab(2,0);
133
134 Box tbxp1D = tbxp1;
135 tbxp1D.makeSlab(2,0);
136 Box tbxp2D = tbxp2;
137 tbxp2D.makeSlab(2,0);
138
139 FArrayBox fab_FC(convert(tbxp2,IntVect(0,0,1)),1,amrex::The_Async_Arena()); //3D
140
141 auto FC=fab_FC.array();
142
143
144 for (int i_comp=0; i_comp < NCONS; i_comp++) {
145 Array4<Real> const& sstore = (vec_sstore[lev])->array(mfi,i_comp);
146#ifdef REMORA_USE_NETCDF
147 FArrayBox* fab_river_source;
148 if (solverChoice.do_rivers_cons[i_comp]) {
149 river_source_cons[i_comp]->update_interpolated_to_time(t_old[lev]);
150 fab_river_source = river_source_cons[i_comp]->fab_interp;
151 }
152 const Array4<const int>& river_pos = (solverChoice.do_rivers) ? vec_river_position[lev]->const_array(mfi) : Array4<const int>();
153 const Array4<const Real>& river_source = (solverChoice.do_rivers_cons[i_comp]) ? fab_river_source->const_array() : Array4<const Real>();
154#else
155 const Array4<const int >& river_pos = Array4<const int>();
156 const Array4<const Real>& river_source = Array4<const Real>();
157#endif
158 prestep_t_advection(lev, bx, gbx, S_old.array(mfi,i_comp),
159 mf_scalarcache.array(mfi,i_comp), Hz, Huon, Hvom,
160 W, DC, FC, sstore, z_w, h, pm, pn, msku, mskv, river_pos,
161 river_source, iic, ntfirst,
162 nrhs, N, dt_lev);
163 }
164
165 // Only do diffusion for salt and temperature, not other tracer(s)
166 for (int i_comp=0; i_comp < NCONS; i_comp++) {
167 const Array4<Real const>& stflx = vec_stflx[lev]->const_array(mfi,i_comp);
168 const Array4<Real const>& btflx = vec_btflx[lev]->const_array(mfi,i_comp);
169 prestep_diffusion(bx,gbx,0,0,S_new.array(mfi,i_comp), S_old.array(mfi,i_comp), ru,
170 Hz, Akt, FC, stflx, btflx, z_r, pm, pn, iic, iic, nnew, nstp,
171 nrhs, N, lambda, dt_lev);
172 }
173
174 //
175 //-----------------------------------------------------------------------
176 // prestep_uv_3d
177 //-----------------------------------------------------------------------
178 //
179 //updates u,v,ru,rv (ru and rv have multiple components)
180
181 prestep_diffusion(tbxp1, gbx, 1, 0, u, uold, ru, Hz, Akv, FC,
182 sustr, bustr, z_r, pm, pn, iic, ntfirst, nnew, nstp, nrhs, N, lambda, dt_lev);
183
184 prestep_diffusion(tbxp1, gbx, 0, 1, v, vold, rv, Hz, Akv, FC,
185 svstr, bvstr, z_r, pm, pn, iic, ntfirst, nnew, nstp, nrhs, N, lambda, dt_lev);
186 }
187}
#define NGROW
#define Temp_comp
#define Salt_comp
#define NCONS
amrex::Vector< NCTimeSeriesRiver * > river_source_cons
Vector of data containers for scalar data in rivers.
Definition REMORA.H:1087
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_stflx
Surface tracer flux; working arrays.
Definition REMORA.H:310
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_sstore
additional scratch space for calculations on temp, salt, etc
Definition REMORA.H:394
void prestep_diffusion(const amrex::Box &bx, const amrex::Box &gbx, const int ioff, const int joff, const amrex::Array4< amrex::Real > &vel, const amrex::Array4< amrex::Real const > &vel_old, const amrex::Array4< amrex::Real > &rvel, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real const > &Akv, const amrex::Array4< amrex::Real > &FC, const amrex::Array4< amrex::Real const > &sstr, const amrex::Array4< amrex::Real const > &bstr, const amrex::Array4< amrex::Real const > &z_r, const amrex::Array4< amrex::Real const > &pm, const amrex::Array4< amrex::Real const > &pn, const int iic, const int ntfirst, const int nnew, int nstp, int nrhs, int N, const amrex::Real lambda, const amrex::Real dt_lev)
Update velocities or tracers with diffusion/viscosity as the last part of the prestep.
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Hz
Width of cells in the vertical (z-) direction (3D, Hz in ROMS)
Definition REMORA.H:238
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akt
Vertical diffusion coefficient (3D)
Definition REMORA.H:258
void prestep(int lev, amrex::MultiFab &mf_uold, amrex::MultiFab &mf_vold, amrex::MultiFab &mf_u, amrex::MultiFab &mf_v, amrex::MultiFab *mf_ru, amrex::MultiFab *mf_rv, amrex::MultiFab &S_old, amrex::MultiFab &S_new, amrex::MultiFab &mf_W, amrex::MultiFab &mf_DC, const amrex::MultiFab *mf_z_r, const amrex::MultiFab *mf_z_w, const amrex::MultiFab *mf_h, const amrex::MultiFab *mf_pm, const amrex::MultiFab *mf_pn, const amrex::MultiFab *mf_sustr, const amrex::MultiFab *mf_svstr, const amrex::MultiFab *mf_bustr, const amrex::MultiFab *mf_bvstr, const amrex::MultiFab *mf_msku, const amrex::MultiFab *mf_mskv, const int iic, const int nfirst, const int nnew, int nstp, int nrhs, int N, const amrex::Real dt_lev)
Wrapper function for prestep.
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Huon
u-volume flux (3D)
Definition REMORA.H:240
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1279
amrex::Vector< std::unique_ptr< amrex::iMultiFab > > vec_river_position
iMultiFab for river positions; contents are indices of rivers
Definition REMORA.H:1094
void prestep_t_advection(int lev, const amrex::Box &tbx, const amrex::Box &gbx, const amrex::Array4< amrex::Real > &tempold, const amrex::Array4< amrex::Real > &tempcache, const amrex::Array4< amrex::Real > &Hz, const amrex::Array4< amrex::Real > &Huon, const amrex::Array4< amrex::Real > &Hvom, const amrex::Array4< amrex::Real > &W, const amrex::Array4< amrex::Real > &DC, const amrex::Array4< amrex::Real > &FC, const amrex::Array4< amrex::Real > &sstore, const amrex::Array4< amrex::Real const > &z_w, const amrex::Array4< amrex::Real const > &h, const amrex::Array4< amrex::Real const > &pm, const amrex::Array4< amrex::Real const > &pn, const amrex::Array4< amrex::Real const > &msku, const amrex::Array4< amrex::Real const > &mskv, const amrex::Array4< int const > &river_pos, const amrex::Array4< amrex::Real const > &river_source, int iic, int ntfirst, int nrhs, int N, const amrex::Real dt_lev)
Prestep advection calculations for the tracers.
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Akv
Vertical viscosity coefficient (3D)
Definition REMORA.H:256
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_btflx
Bottom tracer flux; working arrays.
Definition REMORA.H:314
amrex::Vector< amrex::Real > t_old
old time at each level
Definition REMORA.H:1158
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Hvom
v-volume flux (3D)
Definition REMORA.H:242
amrex::Vector< int > do_rivers_cons