REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_prsgrd.cpp
Go to the documentation of this file.
1#include <REMORA.H>
2
3using namespace amrex;
4
5/**
6 * @param[in ] phi_bx tilebox
7 * @param[in ] phi_gbx grownbox
8 * @param[in ] utbx u-nodal tilebox
9 * @param[in ] vtbx v-nodal tilebox
10 * @param[ out] ru u-velocity RHS
11 * @param[ out] rv v-velocity RHS
12 * @param[in ] pm 1/dx
13 * @param[in ] pn 1/dy
14 * @param[in ] rho density
15 * @param FC temporary
16 * @param[in ] Hz vertical cell height
17 * @param[in ] z_r z coordinates at rho points
18 * @param[in ] z_w z coordinates at w points
19 * @param[in ] msku land-sea mask on u-points
20 * @param[in ] mskv land-sea mask on v-points
21 * @param[in ] nrhs index of RHS component
22 * @param[in ] N number of vertical levels
23 */
24void
25REMORA::prsgrd (const Box& phi_bx, const Box& phi_gbx,
26 const Box& utbx, const Box& vtbx,
27 const Array4<Real >& ru,
28 const Array4<Real >& rv,
29 const Array4<Real const>& pn,
30 const Array4<Real const>& pm,
31 const Array4<Real const>& rho,
32 const Array4<Real >& FC,
33 const Array4<Real const>& Hz,
34 const Array4<Real const>& z_r,
35 const Array4<Real const>& z_w,
36 const Array4<Real const>& msku,
37 const Array4<Real const>& mskv,
38 const int nrhs, const int N)
39{
40 BL_PROFILE("REMORA::prsgrd()");
41 auto phi_bxD=phi_bx;
42 phi_bxD.makeSlab(2,0);
43 auto phi_gbxD=phi_gbx & phi_bx;
44 phi_gbxD.makeSlab(2,0);
45 Box phi_ubx = surroundingNodes(phi_bx,0);
46 Box phi_vbx = surroundingNodes(phi_bx,1);
47 auto utbxD = utbx;
48 auto vtbxD = vtbx;
49 utbxD.makeSlab(2,0);
50 vtbxD.makeSlab(2,0);
51
52 const Real OneFifth = 0.2_rt;
53 const Real OneTwelfth = 1.0_rt/12.0_rt;
54 const Real eps = 1.0E-10_rt;
55 Real GRho = g/solverChoice.rho0;
56 Real GRho0 = 1000.0_rt * GRho;
57 Real HalfGRho = 0.5_rt * GRho;
58
59 int ncomp = 0;
60 int P_comp = ncomp++;
61 int aux_comp = ncomp++;
62 int dR_comp = ncomp++;
63 int dZ_comp = ncomp++;
64 int dRx_comp = ncomp++;
65 int dZx_comp = ncomp++;
66
67 FArrayBox fab(grow(phi_bx,IntVect(1,1,0)),ncomp,The_Async_Arena());
68
69 auto P= fab.array(P_comp);
70 auto aux=fab.array(aux_comp);
71 auto dR= fab.array(dR_comp);
72 auto dZ= fab.array(dZ_comp);
73 auto dRx=fab.array(dRx_comp);
74 auto dZx=fab.array(dZx_comp);
75
76 // Derivatives in the z direction
77 ParallelFor(phi_bx,
78 [=] AMREX_GPU_DEVICE (int i, int j, int k)
79 {
80 if(k>=0&&k<N) {
81 dR(i,j,k)=rho(i,j,k+1)-rho(i,j,k);
82 dZ(i,j,k)=z_r(i,j,k+1)-z_r(i,j,k);
83 } else {
84 dR(i,j,N)=rho(i,j,N)-rho(i,j,N-1);
85 dZ(i,j,N)=z_r(i,j,N)-z_r(i,j,N-1);
86 }
87 });
88
89 ParallelFor(phi_bxD,
90 [=] AMREX_GPU_DEVICE (int i, int j, int )
91 {
92 for(int k=N;k>=0;k--) {
93 Real cff= k>0 ? 2.0_rt*dR(i,j,k)*dR(i,j,k-1) : 2.0_rt*dR(i,j,k)*dR(i,j,k);
94 if (cff>eps) {
95 dR(i,j,k)= k>0 ? cff/(dR(i,j,k)+dR(i,j,k-1)) : cff/(dR(i,j,k)+dR(i,j,k));
96 } else {
97 dR(i,j,k)=0.0_rt;
98 }
99 dZ(i,j,k)= k>0 ? 2.0_rt*dZ(i,j,k)*dZ(i,j,k-1)/(dZ(i,j,k)+dZ(i,j,k-1)) :
100 2.0_rt*dZ(i,j,k)*dZ(i,j,k)/(dZ(i,j,k)+dZ(i,j,k));
101 }
102 });
103
104 ParallelFor(phi_bxD,
105 [=] AMREX_GPU_DEVICE (int i, int j, int )
106 {
107 Real cff1=1.0_rt/(z_r(i,j,N)-z_r(i,j,N-1));
108 Real cff2=0.5_rt*(rho(i,j,N)-rho(i,j,N-1))*(z_w(i,j,N+1)-z_r(i,j,N))*cff1;
109
110 P(i,j,N)=GRho0*z_w(i,j,N+1)+GRho*(rho(i,j,N)+cff2)*(z_w(i,j,N+1)-z_r(i,j,N));
111
112 for (int k=N-1;k>=0;k--)
113 {
114 Real rho_diff = rho(i,j,k+1)-rho(i,j,k) - OneTwelfth* (dR(i,j,k+1)+dR(i,j,k));
115 Real z_diff = z_r(i,j,k+1)-z_r(i,j,k) - OneTwelfth* (dZ(i,j,k+1)+dZ(i,j,k));
116 Real rz_avg = (rho(i,j,k+1)+rho(i,j,k)) * (z_r(i,j,k+1)-z_r(i,j,k));
117
118 P(i,j,k) = P(i,j,k+1) + HalfGRho * ( rz_avg -
119 OneFifth* ( (dR(i,j,k+1)-dR(i,j,k)) * z_diff -
120 (dZ(i,j,k+1)-dZ(i,j,k)) * rho_diff ) );
121 }
122 });
123
124 //This should be nodal
125 // Derivatives in the x direction
126 ParallelFor(phi_ubx,
127 [=] AMREX_GPU_DEVICE (int i, int j, int k)
128 {
129 FC(i,j,k)=(rho(i,j,k)-rho(i-1,j,k)) * msku(i,j,0);
130 aux(i,j,k)=(z_r(i,j,k)-z_r(i-1,j,k)) * msku(i,j,0);
131 });
132
133 //This should be nodal aux and FC need wider boxes above
134 //dZx and dRx may have index mismatch issues at k=2 and k=N
135 ParallelFor(phi_bxD,
136 [=] AMREX_GPU_DEVICE (int i, int j, int )
137 {
138 for(int k=N;k>=0;k--) {
139 Real cff= 2.0_rt*aux(i,j,k)*aux(i+1,j,k);
140 if (cff>eps) {
141 Real cff1= 1.0_rt/(aux(i+1,j,k)+aux(i,j,k));
142 dZx(i,j,k)=cff*cff1;
143 } else {
144 dZx(i,j,k)=0.0_rt;
145 }
146 Real cff1= 2.0_rt*FC(i,j,k)*FC(i+1,j,k);
147 if (cff1>eps) {
148 Real cff2= 1.0_rt/(FC(i,j,k)+FC(i+1,j,k));
149 dRx(i,j,k)=cff1*cff2;
150 } else {
151 dRx(i,j,k)=0.0_rt;
152 }
153 }
154 });
155
156 //This should be nodal aux and FC need wider boxes above
157 ParallelFor(utbxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
158 {
159 for(int k=N;k>=0;k--)
160 {
161 Real rho_diff = rho(i,j,k)-rho(i-1,j,k)- OneTwelfth* (dRx(i,j,k)+dRx(i-1,j,k));
162 Real z_r_diff = z_r(i,j,k)-z_r(i-1,j,k)- OneTwelfth* (dZx(i,j,k)+dZx(i-1,j,k));
163 Real Hz_avg = 0.5_rt * (Hz(i,j,k)+Hz(i-1,j,k));
164
165 Real on_u = 2.0_rt / (pn(i-1,j,0)+pn(i,j,0));
166 ru(i,j,k,nrhs) = on_u * Hz_avg * (
167 P(i-1,j,k) - P(i,j,k) - HalfGRho *
168 ( (rho(i,j,k)+rho(i-1,j,k))*(z_r(i,j,k)-z_r(i-1,j,k))-
169 OneFifth * ( (dRx(i,j,k)-dRx(i-1,j,k)) * z_r_diff -
170 (dZx(i,j,k)-dZx(i-1,j,k)) * rho_diff ) ) );
171 }
172 });
173
174 //This should be nodal
175 ParallelFor(phi_vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
176 {
177 FC(i,j,k)= (rho(i,j,k)-rho(i,j-1,k)) * mskv(i,j,0);
178 aux(i,j,k)= (z_r(i,j,k)-z_r(i,j-1,k)) * mskv(i,j,0);
179 });
180
181 //This should be nodal aux and FC need wider boxes above
182 //dZx and dRx may have index mismatch issues at k=2 and k=N
183 ParallelFor(phi_bxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
184 {
185 for(int k=N;k>=0;k--) {
186 Real cff= 2.0_rt*aux(i,j,k)*aux(i,j+1,k);
187 if (cff>eps) {
188 Real cff1= 1.0_rt/(aux(i,j+1,k)+aux(i,j,k));
189 dZx(i,j,k)=cff*cff1;
190 } else {
191 dZx(i,j,k)=0.0_rt;
192 }
193 Real cff1= 2.0_rt*FC(i,j,k)*FC(i,j+1,k);
194 if (cff1>eps) {
195 Real cff2= 1.0_rt/(FC(i,j,k)+FC(i,j+1,k));
196 dRx(i,j,k)=cff1*cff2;
197 } else {
198 dRx(i,j,k)=0.0_rt;
199 }
200 }
201 });
202
203 //This should be nodal aux and FC need wider boxes above
204 ParallelFor(vtbxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
205 {
206 for (int k=N;k>=0;k--)
207 {
208 Real rho_diff = rho(i,j,k)-rho(i,j-1,k)- OneTwelfth* (dRx(i,j,k)+dRx(i,j-1,k));
209 Real z_r_diff = z_r(i,j,k)-z_r(i,j-1,k)- OneTwelfth* (dZx(i,j,k)+dZx(i,j-1,k));
210 Real Hz_avg = 0.5_rt * (Hz(i,j,k)+Hz(i,j-1,k));
211
212 Real om_v = 2.0_rt / (pm(i,j-1,0)+pm(i,j,0));
213 rv(i,j,k,nrhs) = om_v * Hz_avg * (
214 P(i,j-1,k) - P(i,j,k) - HalfGRho *
215 ( (rho(i,j,k)+rho(i,j-1,k))*(z_r(i,j,k)-z_r(i,j-1,k))-
216 OneFifth * ( (dRx(i,j,k)-dRx(i,j-1,k)) * z_r_diff -
217 (dZx(i,j,k)-dZx(i,j-1,k)) * rho_diff ) ) );
218 } // k
219 });
220}
g
constexpr amrex::Real g
Definition REMORA_Constants.H:16
REMORA::prsgrd
void prsgrd(const amrex::Box &bx, const amrex::Box &gbx, const amrex::Box &utbx, const amrex::Box &vtbx, const amrex::Array4< amrex::Real > &ru, const amrex::Array4< amrex::Real > &rv, const amrex::Array4< amrex::Real const > &pn, const amrex::Array4< amrex::Real const > &pm, const amrex::Array4< amrex::Real const > &rho, const amrex::Array4< amrex::Real > &FC, const amrex::Array4< amrex::Real const > &Hz, const amrex::Array4< amrex::Real const > &z_r, const amrex::Array4< amrex::Real const > &z_w, const amrex::Array4< amrex::Real const > &msku, const amrex::Array4< amrex::Real const > &mskv, const int nrhs, const int N)
Calculate pressure gradient.
Definition REMORA_prsgrd.cpp:25
REMORA::solverChoice
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1279