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 auto phi_bxD=phi_bx;
41 phi_bxD.makeSlab(2,0);
42 auto phi_gbxD=phi_gbx & phi_bx;
43 phi_gbxD.makeSlab(2,0);
44 Box phi_ubx = surroundingNodes(phi_bx,0);
45 Box phi_vbx = surroundingNodes(phi_bx,1);
46 auto utbxD = utbx;
47 auto vtbxD = vtbx;
48 utbxD.makeSlab(2,0);
49 vtbxD.makeSlab(2,0);
50
51 const Real OneFifth = 0.2_rt;
52 const Real OneTwelfth = 1.0_rt/12.0_rt;
53 const Real eps = 1.0E-10_rt;
54 Real GRho = g/solverChoice.rho0;
55 Real GRho0 = 1000.0_rt * GRho;
56 Real HalfGRho = 0.5_rt * GRho;
57
58 FArrayBox fab_P(phi_bx,1,The_Async_Arena());
59 FArrayBox fab_aux(Box(z_r),1,The_Async_Arena());
60 FArrayBox fab_dR(phi_bx,1,The_Async_Arena());
61 FArrayBox fab_dZ(phi_bx,1,The_Async_Arena());
62 FArrayBox fab_dRx(phi_bx,1,The_Async_Arena());
63 FArrayBox fab_dZx(phi_bx,1,The_Async_Arena());
64
65 auto P=fab_P.array();
66 auto aux=fab_aux.array();
67 auto dR=fab_dR.array();
68 auto dZ=fab_dZ.array();
69 auto dRx=fab_dRx.array();
70 auto dZx=fab_dZx.array();
71
72 // Derivatives in the z direction
73 ParallelFor(phi_bx,
74 [=] AMREX_GPU_DEVICE (int i, int j, int k)
75 {
76 if(k>=0&&k<N) {
77 dR(i,j,k)=rho(i,j,k+1)-rho(i,j,k);
78 dZ(i,j,k)=z_r(i,j,k+1)-z_r(i,j,k);
79 } else {
80 dR(i,j,N)=rho(i,j,N)-rho(i,j,N-1);
81 dZ(i,j,N)=z_r(i,j,N)-z_r(i,j,N-1);
82 }
83 });
84
85 ParallelFor(phi_bxD,
86 [=] AMREX_GPU_DEVICE (int i, int j, int )
87 {
88 for(int k=N;k>=0;k--) {
89 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);
90 if (cff>eps) {
91 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));
92 } else {
93 dR(i,j,k)=0.0_rt;
94 }
95 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)) :
96 2.0_rt*dZ(i,j,k)*dZ(i,j,k)/(dZ(i,j,k)+dZ(i,j,k));
97 }
98 });
99
100 ParallelFor(phi_bxD,
101 [=] AMREX_GPU_DEVICE (int i, int j, int )
102 {
103 Real cff1=1.0_rt/(z_r(i,j,N)-z_r(i,j,N-1));
104 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;
105
106 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));
107
108 for (int k=N-1;k>=0;k--)
109 {
110 Real rho_diff = rho(i,j,k+1)-rho(i,j,k) - OneTwelfth* (dR(i,j,k+1)+dR(i,j,k));
111 Real z_diff = z_r(i,j,k+1)-z_r(i,j,k) - OneTwelfth* (dZ(i,j,k+1)+dZ(i,j,k));
112 Real rz_avg = (rho(i,j,k+1)+rho(i,j,k)) * (z_r(i,j,k+1)-z_r(i,j,k));
113
114 P(i,j,k) = P(i,j,k+1) + HalfGRho * ( rz_avg -
115 OneFifth* ( (dR(i,j,k+1)-dR(i,j,k)) * z_diff -
116 (dZ(i,j,k+1)-dZ(i,j,k)) * rho_diff ) );
117 }
118 });
119
120 //This should be nodal
121 // Derivatives in the x direction
122 ParallelFor(phi_ubx,
123 [=] AMREX_GPU_DEVICE (int i, int j, int k)
124 {
125 FC(i,j,k)=(rho(i,j,k)-rho(i-1,j,k)) * msku(i,j,0);
126 aux(i,j,k)=(z_r(i,j,k)-z_r(i-1,j,k)) * msku(i,j,0);
127 });
128
129 //This should be nodal aux and FC need wider boxes above
130 //dZx and dRx may have index mismatch issues at k=2 and k=N
131 ParallelFor(phi_bxD,
132 [=] AMREX_GPU_DEVICE (int i, int j, int )
133 {
134 for(int k=N;k>=0;k--) {
135 Real cff= 2.0_rt*aux(i,j,k)*aux(i+1,j,k);
136 if (cff>eps) {
137 Real cff1= 1.0_rt/(aux(i+1,j,k)+aux(i,j,k));
138 dZx(i,j,k)=cff*cff1;
139 } else {
140 dZx(i,j,k)=0.0_rt;
141 }
142 Real cff1= 2.0_rt*FC(i,j,k)*FC(i+1,j,k);
143 if (cff1>eps) {
144 Real cff2= 1.0_rt/(FC(i,j,k)+FC(i+1,j,k));
145 dRx(i,j,k)=cff1*cff2;
146 } else {
147 dRx(i,j,k)=0.0_rt;
148 }
149 }
150 });
151
152 //This should be nodal aux and FC need wider boxes above
153 ParallelFor(utbxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
154 {
155 for(int k=N;k>=0;k--)
156 {
157 Real rho_diff = rho(i,j,k)-rho(i-1,j,k)- OneTwelfth* (dRx(i,j,k)+dRx(i-1,j,k));
158 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));
159 Real Hz_avg = 0.5_rt * (Hz(i,j,k)+Hz(i-1,j,k));
160
161 Real on_u = 2.0_rt / (pn(i-1,j,0)+pn(i,j,0));
162 ru(i,j,k,nrhs) = on_u * Hz_avg * (
163 P(i-1,j,k) - P(i,j,k) - HalfGRho *
164 ( (rho(i,j,k)+rho(i-1,j,k))*(z_r(i,j,k)-z_r(i-1,j,k))-
165 OneFifth * ( (dRx(i,j,k)-dRx(i-1,j,k)) * z_r_diff -
166 (dZx(i,j,k)-dZx(i-1,j,k)) * rho_diff ) ) );
167 }
168 });
169
170 //This should be nodal
171 ParallelFor(phi_vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
172 {
173 FC(i,j,k)= (rho(i,j,k)-rho(i,j-1,k)) * mskv(i,j,0);
174 aux(i,j,k)= (z_r(i,j,k)-z_r(i,j-1,k)) * mskv(i,j,0);
175 });
176
177 //This should be nodal aux and FC need wider boxes above
178 //dZx and dRx may have index mismatch issues at k=2 and k=N
179 ParallelFor(phi_bxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
180 {
181 for(int k=N;k>=0;k--) {
182 Real cff= 2.0_rt*aux(i,j,k)*aux(i,j+1,k);
183 if (cff>eps) {
184 Real cff1= 1.0_rt/(aux(i,j+1,k)+aux(i,j,k));
185 dZx(i,j,k)=cff*cff1;
186 } else {
187 dZx(i,j,k)=0.0_rt;
188 }
189 Real cff1= 2.0_rt*FC(i,j,k)*FC(i,j+1,k);
190 if (cff1>eps) {
191 Real cff2= 1.0_rt/(FC(i,j,k)+FC(i,j+1,k));
192 dRx(i,j,k)=cff1*cff2;
193 } else {
194 dRx(i,j,k)=0.0_rt;
195 }
196 }
197 });
198
199 //This should be nodal aux and FC need wider boxes above
200 ParallelFor(vtbxD, [=] AMREX_GPU_DEVICE (int i, int j, int )
201 {
202 for (int k=N;k>=0;k--)
203 {
204 Real rho_diff = rho(i,j,k)-rho(i,j-1,k)- OneTwelfth* (dRx(i,j,k)+dRx(i,j-1,k));
205 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));
206 Real Hz_avg = 0.5_rt * (Hz(i,j,k)+Hz(i,j-1,k));
207
208 Real om_v = 2.0_rt / (pm(i,j-1,0)+pm(i,j,0));
209 rv(i,j,k,nrhs) = om_v * Hz_avg * (
210 P(i,j-1,k) - P(i,j,k) - HalfGRho *
211 ( (rho(i,j,k)+rho(i,j-1,k))*(z_r(i,j,k)-z_r(i,j-1,k))-
212 OneFifth * ( (dRx(i,j,k)-dRx(i,j-1,k)) * z_r_diff -
213 (dZx(i,j,k)-dZx(i,j-1,k)) * rho_diff ) ) );
214 } // k
215 });
216}
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:1237