REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_t3dmix.cpp
Go to the documentation of this file.
1#include <REMORA.H>
2
3using namespace amrex;
4
5/**
6 * @param[in ] bx box to operate on
7 * @param[inout] state scalar data
8 * @param[ out] state_rhs scalar data RHS
9 * @param[in ] diff2 diffusivity
10 * @param[in ] Hz vertical cell height
11 * @param[in ] z_r rho-point z coordinate
12 * @param[in ] pm 1/dx
13 * @param[in ] pn 1/dy
14 * @param[in ] msku land-sea mask on u-points
15 * @param[in ] mskv land-sea mask on v-points
16 * @param[in ] dt_lev time step at level
17 * @param[in ] ncomp number of components to do this calculation on
18 * @param[in ] N number of vertical levels
19 */
20void
21REMORA::t3dmix2 (const Box& bx,
22 const Array4<Real >& state,
23 const Array4<Real >& state_rhs,
24 const Array4<Real const>& diff2,
25 const Array4<Real const>& Hz,
26 const Array4<Real const>& z_r,
27 const Array4<Real const>& pm,
28 const Array4<Real const>& pn,
29 const Array4<Real const>& msku,
30 const Array4<Real const>& mskv,
31 const Real dt_lev, const int ncomp,
32 const int N) {
33 if (solverChoice.harmonic_mixing_type == HarmonicMixingType::s) {
34 t3dmix2_s(bx, state, state_rhs, diff2, Hz, pm, pn, msku, mskv, dt_lev, ncomp);
35 } else if (solverChoice.harmonic_mixing_type == HarmonicMixingType::geopotential) {
36 t3dmix2_geo(bx, state, state_rhs, diff2, Hz, z_r, pm, pn, msku, mskv, dt_lev, ncomp, N);
37 }
38}
39
40
41/**
42 * @param[in ] bx box to operate on
43 * @param[inout] state scalar data
44 * @param[ out] state_rhs scalar data RHS
45 * @param[in ] diff2 diffusivity
46 * @param[in ] Hz vertical cell height
47 * @param[in ] pm 1/dx
48 * @param[in ] pn 1/dy
49 * @param[in ] msku land-sea mask on u-points
50 * @param[in ] mskv land-sea mask on v-points
51 * @param[in ] dt_lev time step at level
52 * @param[in ] ncomp number of components to do this calculation on
53 */
54void
55REMORA::t3dmix2_s (const Box& bx,
56 const Array4<Real >& state,
57 const Array4<Real >& state_rhs,
58 const Array4<Real const>& diff2,
59 const Array4<Real const>& Hz,
60 const Array4<Real const>& pm,
61 const Array4<Real const>& pn,
62 const Array4<Real const>& msku,
63 const Array4<Real const>& mskv,
64 const Real dt_lev, const int ncomp)
65{
66 BL_PROFILE("REMORA::t3dmix2_s()");
67 //-----------------------------------------------------------------------
68 // Add in harmonic diffusivity s terms.
69 //-----------------------------------------------------------------------
70
71 Box xbx(bx); xbx.surroundingNodes(0);
72 Box ybx(bx); ybx.surroundingNodes(1);
73
74 FArrayBox fab_FX(xbx,ncomp,The_Async_Arena());
75 FArrayBox fab_FE(ybx,ncomp,The_Async_Arena());
76
77 auto FX=fab_FX.array();
78 auto FE=fab_FE.array();
79
80 ParallelFor(xbx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
81 {
82 const Real pmon_u = (pm(i-1,j,0)+pm(i,j,0))/(pn(i-1,j,0)+pn(i,j,0));
83
84 const Real cff = 0.25_rt * (diff2(i,j,0,n) + diff2(i-1,j,0,n)) * pmon_u;
85 FX(i,j,k,n) = cff * (Hz(i,j,k) + Hz(i-1,j,k)) * (state_rhs(i,j,k,n)-state_rhs(i-1,j,k,n));
86 FX(i,j,k,n) *= msku(i,j,0);
87 });
88
89 ParallelFor(ybx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
90 {
91 const Real pnom_v = (pn(i,j-1,0)+pn(i,j,0))/(pm(i,j-1,0)+pm(i,j,0));
92
93 const Real cff = 0.25_rt*(diff2(i,j,0,n)+diff2(i,j-1,0,n)) * pnom_v;
94 FE(i,j,k,n) = cff * (Hz(i,j,k) + Hz(i,j-1,k)) * (state_rhs(i,j,k,n) - state_rhs(i,j-1,k,n));
95 FE(i,j,k,n) *= mskv(i,j,0);
96 });
97
98 /*
99 Time-step harmonic, S-surfaces diffusion term.
100 */
101 ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
102 {
103 const Real cff = dt_lev*pm(i,j,0)*pn(i,j,0);
104
105 state(i,j,k,n) += cff * ( (FX(i+1,j ,k,n)-FX(i,j,k,n))
106 +(FE(i ,j+1,k,n)-FE(i,j,k,n)) );
107 });
108}
109
110/**
111 * @param[in ] bx box to operate on
112 * @param[inout] state scalar data
113 * @param[ out] state_rhs scalar data RHS
114 * @param[in ] diff2 diffusivity
115 * @param[in ] Hz vertical cell height
116 * @param[in ] z_r rho-point z coordinates
117 * @param[in ] pm 1/dx
118 * @param[in ] pn 1/dy
119 * @param[in ] msku land-sea mask on u-points
120 * @param[in ] mskv land-sea mask on v-points
121 * @param[in ] dt_lev time step at level
122 * @param[in ] ncomp number of components to do this calculation on
123 * @param[in ] N number of vertical levels
124 */
125void
126REMORA::t3dmix2_geo(const Box& bx,
127 const Array4<Real >& state,
128 const Array4<Real >& state_rhs,
129 const Array4<Real const>& diff2,
130 const Array4<Real const>& Hz,
131 const Array4<Real const>& z_r,
132 const Array4<Real const>& pm,
133 const Array4<Real const>& pn,
134 const Array4<Real const>& msku,
135 const Array4<Real const>& mskv,
136 const Real dt_lev, const int ncomp, const int N)
137{
138 BL_PROFILE("REMORA::t3dmix2_geo()");
139 //-----------------------------------------------------------------------
140 // Add in harmonic diffusivity s terms.
141 //-----------------------------------------------------------------------
142
143 Box xbx(bx); xbx.surroundingNodes(0);
144 Box ybx(bx); ybx.surroundingNodes(1);
145 Box zbx(bx); zbx.surroundingNodes(2);
146
147 FArrayBox fab_dZdx(xbx,ncomp,The_Async_Arena());
148 FArrayBox fab_dTdx(xbx,ncomp,The_Async_Arena());
149 FArrayBox fab_dZde(ybx,ncomp,The_Async_Arena());
150 FArrayBox fab_dTde(ybx,ncomp,The_Async_Arena());
151 FArrayBox fab_dTdz(grow(zbx,IntVect(1,1,0)),ncomp,The_Async_Arena());
152 FArrayBox fab_FS(grow(zbx,IntVect(1,1,0)),ncomp,The_Async_Arena());
153 FArrayBox fab_FX(xbx,ncomp,The_Async_Arena());
154 FArrayBox fab_FE(ybx,ncomp,The_Async_Arena());
155
156 auto dZdx = fab_dZdx.array();
157 auto dTdx = fab_dTdx.array();
158 auto dZde = fab_dZde.array();
159 auto dTde = fab_dTde.array();
160 auto dTdz = fab_dTdz.array();
161 auto FS = fab_FS.array();
162 auto FX = fab_FX.array();
163 auto FE = fab_FE.array();
164
165 ParallelFor(xbx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
166 {
167 Real cff = 0.5_rt * (pm(i,j,0) + pm(i-1,j,0)) * msku(i,j,0);
168 dZdx(i,j,k,n) = cff * (z_r(i,j,k) - z_r(i-1,j,k));
169 dTdx(i,j,k,n)=cff*(state_rhs(i ,j,k,n)-state_rhs(i-1,j,k,n));
170 });
171 ParallelFor(ybx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
172 {
173 Real cff = 0.5_rt * (pn(i,j,0) + pn(i,j-1,0)) * mskv(i,j,0);
174 dZde(i,j,k,n) = cff * (z_r(i,j,k) - z_r(i,j-1,k));
175 dTde(i,j,k,n)=cff*(state_rhs(i,j,k,n)-state_rhs(i,j-1,k,n));
176 });
177 ParallelFor(makeSlab(grow(bx,IntVect(1,1,0)),2,0), ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int , int n)
178 {
179 dTdz(i,j,0,n) = 0.0_rt;
180 dTdz(i,j,N+1,n) = 0.0_rt;
181 FS(i,j,0,n) = 0.0_rt;
182 FS(i,j,N+1,n) = 0.0_rt;
183 });
184 ParallelFor(grow(zbx,IntVect(1,1,-1)), ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
185 {
186 Real cff = 1.0_rt / (z_r(i,j,k)-z_r(i,j,k-1));
187 dTdz(i,j,k,n) = cff * (state_rhs(i,j,k,n) - state_rhs(i,j,k-1,n));
188 });
189
190 // Compute components of the rotated tracer flux (T m3/s) along
191 // geopotential surfaces.
192 ParallelFor(xbx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
193 {
194 Real on_u = 2.0_rt / (pn(i,j,0) + pn(i-1,j,0));
195 Real cff = 0.25_rt * (diff2(i,j,0,n)+diff2(i-1,j,0,n)) * on_u;
196 FX(i,j,k,n) = cff *
197 (Hz(i,j,k)+Hz(i-1,j,k))*
198 (dTdx(i,j,k,n)-
199 0.5_rt*(std::min(dZdx(i,j,k,n),0.0_rt)*
200 (dTdz(i-1,j,k ,n)+
201 dTdz(i ,j,k+1,n))+
202 std::max(dZdx(i,j,k,n),0.0_rt)*
203 (dTdz(i-1,j,k+1,n)+
204 dTdz(i ,j,k ,n))));
205 });
206 ParallelFor(ybx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
207 {
208 Real om_v = 2.0_rt / (pm(i,j,0) + pm(i,j-1,0));
209 Real cff = 0.25_rt * (diff2(i,j,0,n)+diff2(i,j-1,0,n)) * om_v;
210 FE(i,j,k,n) = cff *
211 (Hz(i,j,k)+Hz(i,j-1,k))*
212 (dTde(i,j,k,n)-
213 0.5_rt*(std::min(dZde(i,j,k,n),0.0_rt)*
214 (dTdz(i,j-1,k ,n)+
215 dTdz(i,j ,k+1,n))+
216 std::max(dZde(i,j,k,n),0.0_rt)*
217 (dTdz(i,j-1,k+1,n)+
218 dTdz(i,j ,k ,n))));
219 });
220 ParallelFor(grow(zbx,IntVect(0,0,-1)), ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
221 {
222 Real cff = 0.5_rt * diff2(i,j,0,n);
223 Real cff1=std::min(dZdx(i ,j,k-1,n),0.0_rt);
224 Real cff2=std::min(dZdx(i+1,j,k ,n),0.0_rt);
225 Real cff3=std::max(dZdx(i ,j,k ,n),0.0_rt);
226 Real cff4=std::max(dZdx(i+1,j,k-1,n),0.0_rt);
227 FS(i,j,k,n) = cff *
228 (cff1*(cff1*dTdz(i,j,k,n)-dTdx(i ,j,k-1,n))+
229 cff2*(cff2*dTdz(i,j,k,n)-dTdx(i+1,j,k ,n))+
230 cff3*(cff3*dTdz(i,j,k,n)-dTdx(i ,j,k ,n))+
231 cff4*(cff4*dTdz(i,j,k,n)-dTdx(i+1,j,k-1,n)));
232 cff1=std::min(dZde(i,j ,k-1,n),0.0_rt);
233 cff2=std::min(dZde(i,j+1,k ,n),0.0_rt);
234 cff3=std::max(dZde(i,j ,k ,n),0.0_rt);
235 cff4=std::max(dZde(i,j+1,k-1,n),0.0_rt);
236 FS(i,j,k,n)=FS(i,j,k,n)+
237 cff*
238 (cff1*(cff1*dTdz(i,j,k,n)-dTde(i,j ,k-1,n))+
239 cff2*(cff2*dTdz(i,j,k,n)-dTde(i,j+1,k ,n))+
240 cff3*(cff3*dTdz(i,j,k,n)-dTde(i,j ,k ,n))+
241 cff4*(cff4*dTdz(i,j,k,n)-dTde(i,j+1,k-1,n)));
242 });
243
244 // Time-step harmonic, geopotential diffusion term (m Tunits).
245
246 ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
247 {
248 Real cff=dt_lev*pm(i,j,0)*pn(i,j,0);
249 Real cff1=cff*(FX(i+1,j ,k,n)-FX(i,j,k,n));
250 Real cff2=cff*(FE(i ,j+1,k,n)-FE(i,j,k,n));
251 Real cff3=dt_lev*(FS(i,j,k+1,n)-FS(i,j,k,n));
252 Real cff4=cff1+cff2+cff3;
253 state(i,j,k,n)=state(i,j,k,n)+cff4;
254 });
255}
HarmonicMixingType::s
@ s
HarmonicMixingType::geopotential
@ geopotential
REMORA::t3dmix2
void t3dmix2(const amrex::Box &bx, const amrex::Array4< amrex::Real > &state, const amrex::Array4< amrex::Real > &state_rhs, const amrex::Array4< amrex::Real const > &diff2, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real const > &z_r, 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::Real dt_lev, const int ncomp, const int N)
Wrapper for harmonic diffusivity for tracers.
Definition REMORA_t3dmix.cpp:21
REMORA::solverChoice
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1403
REMORA::t3dmix2_geo
void t3dmix2_geo(const amrex::Box &bx, const amrex::Array4< amrex::Real > &state, const amrex::Array4< amrex::Real > &state_rhs, const amrex::Array4< amrex::Real const > &diff2, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real const > &z_r, 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::Real dt_lev, const int ncomp, const int N)
Harmonic diffusivity for tracers along geopotential surfaces.
Definition REMORA_t3dmix.cpp:126
REMORA::t3dmix2_s
void t3dmix2_s(const amrex::Box &bx, const amrex::Array4< amrex::Real > &state, const amrex::Array4< amrex::Real > &state_rhs, const amrex::Array4< amrex::Real const > &diff2, const amrex::Array4< amrex::Real const > &Hz, 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::Real dt_lev, const int ncomp)
Harmonic diffusivity for tracers along S-coordinate level surfaces.
Definition REMORA_t3dmix.cpp:55
SolverChoice::harmonic_mixing_type
HarmonicMixingType harmonic_mixing_type
Definition REMORA_DataStruct.H:675