REMORA
Energy Research and Forecasting: An Atmospheric Modeling Code
DepthStretchTransform.H
Go to the documentation of this file.
1 #ifndef STRETCH_H_
2 #define STRETCH_H_
3 
4 #include <cmath>
5 #include <DataStruct.H>
6 #include <REMORA.H>
7 #include <prob_common.H>
8 
9 using namespace amrex;
10 
11 void
13 {
14  std::unique_ptr<MultiFab>& mf_z_w = vec_z_w[lev];
15  std::unique_ptr<MultiFab>& mf_z_r = vec_z_r[lev];
16  std::unique_ptr<MultiFab>& mf_s_r = vec_s_r[lev];
17  std::unique_ptr<MultiFab>& mf_Hz = vec_Hz[lev];
18  std::unique_ptr<MultiFab>& mf_h = vec_hOfTheConfusingName[lev];
19  std::unique_ptr<MultiFab>& mf_Zt_avg1 = vec_Zt_avg1[lev];
20  std::unique_ptr<MultiFab>& mf_z_phys_nd = vec_z_phys_nd[lev];
21  auto N_loop = Geom(lev).Domain().size()[2]-1; // Number of vertical "levs" aka, NZ
22 
23  for ( MFIter mfi(*cons_new[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi )
24  {
25  Array4<Real> const& z_w = (mf_z_w)->array(mfi);
26  Array4<Real> const& z_r = (mf_z_r)->array(mfi);
27  Array4<Real> const& s_r = (mf_s_r)->array(mfi);
28  Array4<Real> const& Hz = (mf_Hz)->array(mfi);
29  Array4<Real> const& h = (mf_h)->array(mfi);
30  Array4<Real> const& Zt_avg1 = (mf_Zt_avg1)->array(mfi);
31  Box bx = mfi.tilebox();
32  Box gbx2 = bx;
33  gbx2.grow(IntVect(NGROW,NGROW,0));
34  Box gbx3 = bx;
35  gbx3.grow(IntVect(NGROW+1,NGROW+1,0));
36  Box gbx2D = gbx2;
37  gbx2D.makeSlab(2,0);
38  Box gbx3D = gbx3;
39  gbx3D.makeSlab(2,0);
40  Box wgbx3 = gbx3.surroundingNodes(2);
41 
42  const auto & geomdata = Geom(lev).data();
43 
44  int nz = geom[lev].Domain().length(2);
45 
46  auto N = nz; // Number of vertical "levels" aka, NZ
47  //forcing tcline to be the same as probhi for now, one in DataStruct.H other in inputs
48  Real hc=-min(geomdata.ProbHi(2),-solverChoice.tcline); // Do we need to enforce min here?
49  const auto local_theta_s = solverChoice.theta_s;
50  const auto local_theta_b = solverChoice.theta_b;
51 
52  amrex::ParallelFor(wgbx3, [=] AMREX_GPU_DEVICE (int i, int j, int k)
53  {
54  z_w(i,j,k) = h(i,j,0);
55  });
56 
57  // ROMS Transform 2
58  Gpu::streamSynchronize();
59 
60  amrex::ParallelFor(gbx3D, [=] AMREX_GPU_DEVICE (int i, int j, int )
61  {
62  for (int k=-1; k<=N_loop; k++) {
63  // const Real z = prob_lo[2] + (k + 0.5) * dx[2];
64  // const auto prob_lo = geomdata.ProbLo();
65  // const auto dx = geomdata.CellSize();
66  // This is the z for the bottom of the cell this k corresponds to
67  // if we weren't stretching and transforming
68  // const Real z = prob_lo[2] + (k) * dx[2];
69  // h(i,j,0) = -prob_lo[2]; // conceptually
70 
71  ////////////////////////////////////////////////////////////////////
72  //ROMS Stretching 4
73  // Move this block to it's own function for maintainability if needed
74  // Information about the problem dimension would need to be added
75  // This file would need a k dependent function to return the
76  // stretching scalars, or access to 4 vectors of length prob_length(2)
77  /////////////////////////////////////////////////////////////////////
78  Real ds = 1.0_rt / Real(N);
79 
80  Real cff_r, cff_w, cff1_r, cff1_w, cff2_r, cff2_w, Csur, Cbot;
81  Real sc_r,sc_w,Cs_r,Cs_w;
82 
83  if (k==N) // end of array // pretend we're storing 0?
84  {
85  sc_w=0.0; //sc_w / hc
86  Cs_w=0.0; //Cs_w
87  }
88  else if (k==0) // beginning of array
89  {
90  sc_w=-1.0; //sc_w / hc
91  Cs_w=-1.0; //Cs_w
92  }
93  else
94  {
95  sc_w=ds*(k-N);
96 
97  if (local_theta_s > 0.0_rt) {
98  Csur=(1.0_rt-std::cosh(local_theta_s*sc_w))/
99  (std::cosh(local_theta_s)-1.0_rt);
100  } else {
101  Csur=-sc_w*sc_w;
102  }
103 
104  if (local_theta_b > 0.0_rt) {
105  Cbot=(std::exp(local_theta_b*Csur)-1.0_rt)/
106  (1.0_rt-std::exp(-local_theta_b));
107  Cs_w=Cbot;
108  } else {
109  Cs_w=Csur;
110  }
111  } // k test
112 
113  cff_w=hc*sc_w;
114  cff1_w=Cs_w;
115 
116  //cff_r => sc_r *hc
117  //cff1_r => Cs_r
118  //Don't do anything special for first/last index
119  {
120  sc_r=ds*(k-N+0.5_rt);
121 
122  if (local_theta_s > 0.0_rt) {
123  Csur=(1.0_rt-std::cosh(local_theta_s*sc_r))/
124  (std::cosh(local_theta_s)-1.0_rt);
125  } else {
126  Csur=-sc_r*sc_r;
127  }
128 
129  if (local_theta_b > 0.0_rt) {
130  Cbot=(std::exp(local_theta_b*Csur)-1.0_rt)/
131  (1.0_rt-std::exp(-local_theta_b));
132  Cs_r=Cbot;
133  } else {
134  Cs_r=Csur;
135  }
136  }
137 
138  if (i==0&&j==0&&k<N&&k>=0) {
139  s_r(0,0,k) = sc_r;
140  }
141 
142  cff_r=hc*sc_r;
143  cff1_r=Cs_r;
144 
145  ////////////////////////////////////////////////////////////////////
146  Real hwater=h(i,j,0);
147  //
148  // if (k==0) //extra guess added (maybe not actually defined in ROMS)
149  // {
150  // Real hinv=1.0_rt/(hc+hwater);
151  // cff2_r=(cff_r+cff1_r*hwater)*hinv;
152  // // z_w(i,j,k-2) = hwater;
153  // // z_w(i,j,k-1)= -hwater;
154 
155  // z_r(i,j,k) = Zt_avg1(i,j,0)+(Zt_avg1(i,j,0)+hwater)*cff2_r;
156  // Hz(i,j,k)=z_w(i,j,k)+hwater;//-z_w(i,j,k-1);
157  // } else
158 
159  //Note, we are not supporting ICESHELF flag
160 
161  Real hinv=1.0_rt/(hc+hwater);
162  cff2_r=(cff_r+cff1_r*hwater)*hinv;
163  cff2_w=(cff_w+cff1_w*hwater)*hinv;
164 
165  if(k==0) {
166  // HACK: should actually be the normal expression with coeffs evaluated at k=N-1
167  z_w(i,j,N-1)=Zt_avg1(i,j,0);
168 
169  } else if (k==-1) {
170  h(i,j,0,1) = Zt_avg1(i,j,0)+(Zt_avg1(i,j,0)+hwater)*cff2_w;
171 
172  } else {
173  z_w(i,j,k-1)=Zt_avg1(i,j,0)+(Zt_avg1(i,j,0)+hwater)*cff2_w;
174 
175  }
176 
177  if(k!=-1) {
178  z_r(i,j,k)=Zt_avg1(i,j,0)+(Zt_avg1(i,j,0)+hwater)*cff2_r;
179  }
180  } // k
181  });
182 
183  Gpu::streamSynchronize();
184 
185  amrex::ParallelFor(gbx3, [=] AMREX_GPU_DEVICE (int i, int j, int k)
186  {
187  if (k==0) {
188  Hz(i,j,k)=z_w(i,j,k)+h(i,j,0);
189  } else {
190  Hz(i,j,k)=z_w(i,j,k)-z_w(i,j,k-1);
191  }
192  });
193  } // mfi
194 
195  Real time = t_new[lev];
196  FillPatch(lev, time, *vec_z_w[lev], GetVecOfPtrs(vec_z_w));
197  FillPatch(lev, time, *vec_z_r[lev], GetVecOfPtrs(vec_z_r));
198  FillPatch(lev, time, *vec_s_r[lev], GetVecOfPtrs(vec_s_r));
199  FillPatch(lev, time, *vec_Hz[lev] , GetVecOfPtrs(vec_Hz));
200 
201  // Define nodal z as average of z on w-faces
202  for ( MFIter mfi(*cons_new[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi )
203  {
204  Array4<Real> const& z_w = (mf_z_w)->array(mfi);
205  Array4<Real> const& z_phys_nd = (mf_z_phys_nd)->array(mfi);
206 
207  Box z_w_box = Box(z_w);
208  auto const lo = amrex::lbound(z_w_box);
209  auto const hi = amrex::ubound(z_w_box);
210 
211  //
212  // WARNING: z_w(i,j,k) refers to the face on the HIGH side of cell (i,j,k)
213  // WARNING: z_phys_nd(i,j,k) refers to the node on the LOW side of cell (i,j,k)
214  //
215  ParallelFor(Box(z_phys_nd), [=] AMREX_GPU_DEVICE (int i, int j, int k)
216  {
217  // For now assume all boundaries are constant height --
218  // we will enforce periodicity below
219  int kk = (k == 0) ? hi.z : k-1;
220  if ( i >= lo.x && i <= hi.x-1 && j >= lo.y && j <= hi.y-1 )
221  {
222  z_phys_nd(i,j,k)=0.25*( z_w(i,j ,kk) + z_w(i+1,j ,kk) +
223  z_w(i,j+1,kk) + z_w(i+1,j+1,kk) );
224  } else {
225  int ii = std::min(std::max(i, lo.x), hi.x);
226  int jj = std::min(std::max(j, lo.y), hi.y);
227  z_phys_nd(i,j,k) = z_w(ii,jj,kk);
228  }
229  if (k == 0) z_phys_nd(i,j,k) *= -1.;
230  });
231  } // mf
232 
233  // Note that we do *not* want to do a multilevel fill here -- we have
234  // already filled z_phys_nd on the grown boxes, but we enforce periodicity just in case
235  vec_z_phys_nd[lev]->FillBoundary(geom[lev].periodicity());
236 }
237 
238 #endif
#define NGROW
Definition: IndexDefines.H:13
void stretch_transform(int lev)
Definition: DepthStretchTransform.H:12