REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_DepthStretchTransform.H
Go to the documentation of this file.
1#ifndef _REMORA_STRETCH_H_
2#define _REMORA_STRETCH_H_
3
4#include <cmath>
5#include <REMORA_DataStruct.H>
6#include <REMORA.H>
8
9using namespace amrex;
10
11void
13{
14 BL_PROFILE("REMORA::calc_stretch_coeffs()");
15 int nz = geom[0].Domain().length(2);
16
17 auto N = nz; // Number of vertical "levels" aka, NZ
18 Real ds = 1.0_rt / Real(N);
19
20 Box bx(IntVect(0,0,0),IntVect(0,0,N));
21 const auto local_theta_s = solverChoice.theta_s;
22 const auto local_theta_b = solverChoice.theta_b;
23
24 //amrex::Vector<Real> s_w_vec(N+1); auto s_w_dat = s_w_vec.data();
25 //amrex::Vector<Real> Cs_w_vec(N+1); auto Cs_w_dat = Cs_w_vec.data();
26 //amrex::Vector<Real> s_r_vec(N); auto s_r_dat = s_r_vec.data();
27 //amrex::Vector<Real> Cs_r_vec(N); auto Cs_r_dat = Cs_r_vec.data();
28 auto s_w_dat = s_w.data();
29 auto Cs_w_dat = Cs_w.data();
30 auto s_r_dat = s_r.data();
31 auto Cs_r_dat = Cs_r.data();
32
33 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int , int , int k)
34 {
35 Real Csur,Cbot;
36
37 s_w_dat[k]=ds*(k-N);
38 if (local_theta_s > 0.0_rt) {
39 Csur=(1.0_rt-std::cosh(local_theta_s*s_w_dat[k]))/
40 (std::cosh(local_theta_s)-1.0_rt);
41 } else {
42 Csur=-s_w_dat[k]*s_w_dat[k];
43 }
44
45 if (local_theta_b > 0.0_rt) {
46 Cbot=(std::exp(local_theta_b*Csur)-1.0_rt)/
47 (1.0_rt-std::exp(-local_theta_b));
48 Cs_w_dat[k]=Cbot;
49 } else {
50 Cs_w_dat[k]=Csur;
51 }
52
53 if (k<N) {
54 s_r_dat[k]=ds*(k-N+0.5_rt);
55
56 if (local_theta_s > 0.0_rt) {
57 Csur=(1.0_rt-std::cosh(local_theta_s*s_r_dat[k]))/
58 (std::cosh(local_theta_s)-1.0_rt);
59 } else {
60 Csur=-s_r_dat[k]*s_r_dat[k];
61 }
62
63 if (local_theta_b > 0.0_rt) {
64 Cbot=(std::exp(local_theta_b*Csur)-1.0_rt)/
65 (1.0_rt-std::exp(-local_theta_b));
66 Cs_r_dat[k]=Cbot;
67 } else {
68 Cs_r_dat[k]=Csur;
69 }
70 }
71 });
72}
73
74/**
75 * @param[in] lev level to operate on
76 */
77void
79{
80 BL_PROFILE("REMORA::stretch_transform()");
81 std::unique_ptr<MultiFab>& mf_z_w = vec_z_w[lev];
82 std::unique_ptr<MultiFab>& mf_z_r = vec_z_r[lev];
83 std::unique_ptr<MultiFab>& mf_Hz = vec_Hz[lev];
84 std::unique_ptr<MultiFab>& mf_h = vec_hOfTheConfusingName[lev];
85 std::unique_ptr<MultiFab>& mf_Zt_avg1 = vec_Zt_avg1[lev];
86 std::unique_ptr<MultiFab>& mf_z_phys_nd = vec_z_phys_nd[lev];
87
88 auto s_w_dat = s_w.data();
89 auto Cs_w_dat = Cs_w.data();
90 auto s_r_dat = s_r.data();
91 auto Cs_r_dat = Cs_r.data();
92
93 for ( MFIter mfi(*cons_new[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi )
94 {
95 Array4<Real> const& z_w = (mf_z_w)->array(mfi);
96 Array4<Real> const& z_r = (mf_z_r)->array(mfi);
97 Array4<Real> const& Hz = (mf_Hz)->array(mfi);
98 Array4<Real> const& h = (mf_h)->array(mfi);
99 Array4<Real> const& Zt_avg1 = (mf_Zt_avg1)->array(mfi);
100 Box bx = mfi.tilebox();
101 Box gbx2 = bx;
102 gbx2.grow(IntVect(NGROW,NGROW,0));
103 Box gbx3 = bx;
104 gbx3.grow(IntVect(NGROW+1,NGROW+1,0));
105 Box gbx2D = gbx2;
106 gbx2D.makeSlab(2,0);
107 Box gbx3D = gbx3;
108 gbx3D.makeSlab(2,0);
109 Box wgbx3 = gbx3;
110 wgbx3.surroundingNodes(2);
111
112 const auto & geomdata = Geom(lev).data();
113
114 int nz = geom[lev].Domain().length(2);
115
116 auto N = nz; // Number of vertical "levels" aka, NZ
117 //forcing tcline to be the same as probhi for now, one in DataStruct.H other in inputs
118 Real hc=-min(geomdata.ProbHi(2),-solverChoice.tcline); // Do we need to enforce min here?
119 Real ds = 1.0_rt / Real(N);
120
121 amrex::ParallelFor(wgbx3, [=] AMREX_GPU_DEVICE (int i, int j, int k)
122 {
123 z_w(i,j,k) = h(i,j,0);
124 });
125
126 // ROMS Transform 2
127 Gpu::streamSynchronize();
128
129 amrex::ParallelFor(wgbx3, [=] AMREX_GPU_DEVICE (int i, int j, int k)
130 {
131 Real sc_r=ds*(k-N+0.5_rt);
132 Real sc_w=ds*(k-N);
133
134 Real cff_r=hc*sc_r;
135 Real cff1_r=Cs_r_dat[k];
136
137 Real cff_w=hc*sc_w;
138 Real cff1_w=Cs_w_dat[k];
139
140 Real hwater=h(i,j,0);
141
142 Real hinv=1.0_rt/(hc+hwater);
143 Real cff2_r=(cff_r+cff1_r*hwater)*hinv;
144 Real cff2_w=(cff_w+cff1_w*hwater)*hinv;
145
146 if(k==N) {
147 // HACK: should actually be the normal expression with coeffs evaluated at k=N-1
148 z_w(i,j,N)=Zt_avg1(i,j,0);
149
150 } else if (k==0) {
151 h(i,j,0,1) = Zt_avg1(i,j,0)+(Zt_avg1(i,j,0)+hwater)*cff2_w;
152 z_w(i,j,0) = h(i,j,0,1);
153
154 } else {
155 z_w(i,j,k)=Zt_avg1(i,j,0)+(Zt_avg1(i,j,0)+hwater)*cff2_w;
156
157 }
158
159 if(k!=N) {
160 z_r(i,j,k)=Zt_avg1(i,j,0)+(Zt_avg1(i,j,0)+hwater)*cff2_r;
161 }
162 });
163
164 Gpu::streamSynchronize();
165
166 amrex::ParallelFor(gbx3, [=] AMREX_GPU_DEVICE (int i, int j, int k)
167 {
168 Hz(i,j,k)=z_w(i,j,k+1)-z_w(i,j,k);
169 });
170 } // mfi
171
172 vec_z_w[lev]->FillBoundary(geom[lev].periodicity());
173 vec_z_r[lev]->FillBoundary(geom[lev].periodicity());
174 vec_Hz[lev]->FillBoundary(geom[lev].periodicity());
175
176 // Define nodal z as average of z on w-faces
177 for ( MFIter mfi(*cons_new[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi )
178 {
179 Array4<Real> const& z_w = (mf_z_w)->array(mfi);
180 Array4<Real> const& z_phys_nd = (mf_z_phys_nd)->array(mfi);
181
182 Box z_w_box = Box(z_w);
183 auto const lo = amrex::lbound(z_w_box);
184 auto const hi = amrex::ubound(z_w_box);
185
186 //
187 // NOTE: we assume that all boxes extend the full extent of the domain in the vertical direction
188 //
189 // NOTE: z_phys_nd(i,j,k) and z_w(i,j,k) both refer to the node on the LOW side of cell (i,j,k)
190 //
191
192 // We shrink the box in the vertical since z_phys_nd has a ghost cell in the vertical
193 // which we will fill by extrapolation, not from z_w
194 Box bx = mfi.growntilebox(IntVect(NGROW,NGROW,0));//Box(z_phys_nd); bx.grow(2,-1);
195
196 ParallelFor(convert(bx,IntVect(0,0,1)), [=] AMREX_GPU_DEVICE (int i, int j, int k)
197 {
198 // For now assume all boundaries are constant height --
199 // we will enforce periodicity below
200 if ( i >= lo.x && i <= hi.x-1 && j >= lo.y && j <= hi.y-1 )
201 {
202 z_phys_nd(i,j,k)=0.25_rt*( z_w(i,j ,k) + z_w(i+1,j ,k) +
203 z_w(i,j+1,k) + z_w(i+1,j+1,k) );
204 } else {
205 int ii = std::min(std::max(i, lo.x), hi.x);
206 int jj = std::min(std::max(j, lo.y), hi.y);
207 z_phys_nd(i,j,k) = z_w(ii,jj,k);
208 }
209 });
210
211 // Fill nodes below the surface (to avoid out of bounds errors in particle functions)
212 int klo = -1;
213 ParallelFor(makeSlab(bx,2,0), [=] AMREX_GPU_DEVICE (int i, int j, int)
214 {
215 z_phys_nd(i,j,klo) = 2.0_rt * z_phys_nd(i,j,klo+1) - z_phys_nd(i,j,klo+2);
216 });
217 } // mf
218
219 // Note that we do *not* want to do a multilevel fill here -- we have
220 // already filled z_phys_nd on the grown boxes, but we enforce periodicity just in case
221 vec_z_phys_nd[lev]->FillBoundary(geom[lev].periodicity());
222}
223
224#endif
#define NGROW
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_hOfTheConfusingName
Bathymetry data (2D, positive valued, h in ROMS)
Definition REMORA.H:235
amrex::Vector< amrex::MultiFab * > cons_new
multilevel data container for current step's scalar data: temperature, salinity, passive scalar
Definition REMORA.H:223
void stretch_transform(int lev)
Calculate vertical stretched coordinates.
amrex::Gpu::DeviceVector< amrex::Real > s_w
Scaled vertical coordinate (range [0,1]) that transforms to z, defined at w-points (cell faces)
Definition REMORA.H:274
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::Gpu::DeviceVector< amrex::Real > s_r
Scaled vertical coordinate (range [0,1]) that transforms to z, defined at rho points (cell centers)
Definition REMORA.H:272
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_z_r
z coordinates at rho points (cell centers)
Definition REMORA.H:267
void calc_stretch_coeffs()
calculate vertical stretch coefficients
amrex::Gpu::DeviceVector< amrex::Real > Cs_r
Stretching coefficients at rho points.
Definition REMORA.H:282
static SolverChoice solverChoice
Container for algorithmic choices.
Definition REMORA.H:1279
amrex::Gpu::DeviceVector< amrex::Real > Cs_w
Stretching coefficients at w points.
Definition REMORA.H:284
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_z_phys_nd
z coordinates at psi points (cell nodes)
Definition REMORA.H:287
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_Zt_avg1
Average of the free surface, zeta (2D)
Definition REMORA.H:290
amrex::Vector< std::unique_ptr< amrex::MultiFab > > vec_z_w
z coordinates at w points (faces between z-cells)
Definition REMORA.H:270
amrex::Real theta_b
amrex::Real theta_s
amrex::Real tcline