ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Macros | Functions | Variables
ERF_HurricaneDiagnostics.H File Reference
#include <AMReX.H>
#include <AMReX_MultiFab.H>
#include <AMReX_ParallelReduce.H>
#include <limits>
#include "ERF_DataStruct.H"
Include dependency graph for ERF_HurricaneDiagnostics.H:
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Macros

#define M_PI   3.14159265358979323846
 

Functions

void ComputeGlobalMinLocation (const amrex::Geometry &geom, const amrex::Vector< amrex::MultiFab > &S_data, const amrex::Vector< amrex::MultiFab > *forecast_state_at_lev, amrex::Real *d_val_min_ptr, int *d_i_min_ptr, int *d_j_min_ptr, amrex::Real &global_val_min, int &global_i_min, int &global_j_min, amrex::Vector< std::array< amrex::Real, 2 >> &hurricane_eye_track_xy, amrex::Vector< std::array< amrex::Real, 2 >> &hurricane_eye_track_latlon)
 
void HurricaneTrackerCircle (const amrex::Vector< std::array< amrex::Real, 2 >> &hurricane_eye_track_xy, amrex::Vector< std::array< amrex::Real, 2 >> &hurricane_tracker_circle)
 
void HurricaneEyeTrackerInitial (const amrex::Geometry &geom, const amrex::Vector< amrex::MultiFab > &S_data, const amrex::Vector< amrex::MultiFab > *forecast_state_at_lev, amrex::Vector< std::array< amrex::Real, 2 >> &hurricane_eye_track_xy, amrex::Vector< std::array< amrex::Real, 2 >> &hurricane_eye_track_latlon, const amrex::Real &hurricane_eye_latitude, const amrex::Real &hurricane_eye_longitude)
 
void HurricaneEyeTrackerNotInitial (const amrex::Geometry &geom, const amrex::Vector< amrex::MultiFab > &S_data, MoistureType moisture_type, const amrex::Vector< amrex::MultiFab > *forecast_state_at_lev, amrex::Vector< std::array< amrex::Real, 2 >> &hurricane_eye_track_xy, amrex::Vector< std::array< amrex::Real, 2 >> &hurricane_eye_track_latlon)
 
void HurricaneEyeTracker (const amrex::Geometry &geom, const amrex::Vector< amrex::MultiFab > &S_data, MoistureType moisture_type, const amrex::Vector< amrex::MultiFab > *forecast_state_at_lev, const amrex::Real &hurricane_eye_latitude, const amrex::Real &hurricane_eye_longitude, amrex::Vector< std::array< amrex::Real, 2 >> &hurricane_eye_track_xy, amrex::Vector< std::array< amrex::Real, 2 >> &hurricane_eye_track_latlon, amrex::Vector< std::array< amrex::Real, 2 >> &hurricane_tracker_circle)
 
void HurricaneMaxVelTracker (const amrex::Geometry &geom, const amrex::MultiFab &mf_cc_vel, const amrex::Real &time, const amrex::Vector< std::array< amrex::Real, 2 >> &hurricane_eye_track_xy, amrex::Vector< std::array< amrex::Real, 2 >> &hurricane_maxvel_vs_time)
 

Variables

struct {
   amrex::Real   value
 
   int   rank
 
} in
 
struct {
   amrex::Real   value
 
   int   rank
 
} out
 

Macro Definition Documentation

◆ M_PI

#define M_PI   3.14159265358979323846

Routines to compute hurricane diagnostics

Function Documentation

◆ ComputeGlobalMinLocation()

void ComputeGlobalMinLocation ( const amrex::Geometry &  geom,
const amrex::Vector< amrex::MultiFab > &  S_data,
const amrex::Vector< amrex::MultiFab > *  forecast_state_at_lev,
amrex::Real *  d_val_min_ptr,
int *  d_i_min_ptr,
int *  d_j_min_ptr,
amrex::Real &  global_val_min,
int &  global_i_min,
int &  global_j_min,
amrex::Vector< std::array< amrex::Real, 2 >> &  hurricane_eye_track_xy,
amrex::Vector< std::array< amrex::Real, 2 >> &  hurricane_eye_track_latlon 
)
35 {
36  amrex::Real h_val_min;
37  int h_i_min, h_j_min;
38 
39  amrex::Gpu::copy(amrex::Gpu::deviceToHost, d_val_min_ptr, d_val_min_ptr + 1, &h_val_min);
40  amrex::Gpu::copy(amrex::Gpu::deviceToHost, d_i_min_ptr, d_i_min_ptr + 1, &h_i_min);
41  amrex::Gpu::copy(amrex::Gpu::deviceToHost, d_j_min_ptr, d_j_min_ptr + 1, &h_j_min);
42  amrex::Gpu::synchronize();
43 
44  amrex::Real local_val_min = h_val_min;
45  int local_i_min = h_i_min;
46  int local_j_min = h_j_min;
47 
48  int rank = amrex::ParallelDescriptor::MyProc();
49 
50  in.value = local_val_min;
51  in.rank = rank;
52 
53  #ifdef AMREX_USE_MPI
54  MPI_Allreduce(&in, &out, 1, MPI_DOUBLE_INT, MPI_MINLOC, MPI_COMM_WORLD);
55  #else
56  out = in;
57  #endif
58 
59  global_val_min = out.value;
60  int owner_rank = out.rank;
61 
62  // Broadcast the indices from the rank that owns the minimum
63  global_i_min = local_i_min;
64  global_j_min = local_j_min;
65 
66  amrex::ParallelDescriptor::Bcast(&global_i_min, 1, owner_rank);
67  amrex::ParallelDescriptor::Bcast(&global_j_min, 1, owner_rank);
68 
69  if (rank == 0) {
70  amrex::Print() << "Global minimum distance to hurricane eye (k=0): "
71  << global_val_min << " at (i,j) = ("
72  << global_i_min << ", " << global_j_min << ")\n";
73  }
74 
75  amrex::Gpu::DeviceScalar<amrex::Real> d_eye_lat(0.0), d_eye_lon(0.0);
76 
77  amrex::Real* d_eye_lat_ptr = d_eye_lat.dataPtr();
78  amrex::Real* d_eye_lon_ptr = d_eye_lon.dataPtr();
79 
80  // On owner_rank, compute eye_lat and eye_lon
81  if (rank == owner_rank) {
82  for (amrex::MFIter mfi(S_data[IntVars::cons]); mfi.isValid(); ++mfi) {
83  const amrex::Box& box = mfi.validbox();
84  const auto latlon_arr = (*forecast_state_at_lev)[4].array(mfi);
85  amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
86  if (i == global_i_min && j == global_j_min && k == 0) {
87  *d_eye_lat_ptr = latlon_arr(i,j,k,0);
88  *d_eye_lon_ptr = latlon_arr(i,j,k,1);
89  }
90  });
91  }
92  }
93 
94  amrex::Real eye_lat = d_eye_lat.dataValue();
95  amrex::Real eye_lon = d_eye_lon.dataValue();
96 
97  // Synchronize to ensure the owner has computed values
98  amrex::Gpu::synchronize();
99 
100  amrex::ParallelDescriptor::Bcast(&eye_lat, 1, owner_rank);
101  amrex::ParallelDescriptor::Bcast(&eye_lon, 1, owner_rank);
102 
103  const auto dx = geom.CellSizeArray();
104  const auto prob_lo = geom.ProbLoArray();
105 
106  amrex::Real eye_x = prob_lo[0] + (global_i_min+0.5)*dx[0];
107  amrex::Real eye_y = prob_lo[1] + (global_j_min+0.5)*dx[1];
108 
109  hurricane_eye_track_xy.push_back({eye_x, eye_y});
110  hurricane_eye_track_latlon.push_back({eye_lon, eye_lat});
111 }
out
struct @19 out
rank
int rank
Definition: ERF_HurricaneDiagnostics.H:21
in
struct @19 in
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:19
IntVars::cons
@ cons
Definition: ERF_IndexDefines.H:158

Referenced by HurricaneEyeTrackerInitial(), and HurricaneEyeTrackerNotInitial().

Here is the caller graph for this function:

◆ HurricaneEyeTracker()

void HurricaneEyeTracker ( const amrex::Geometry &  geom,
const amrex::Vector< amrex::MultiFab > &  S_data,
MoistureType  moisture_type,
const amrex::Vector< amrex::MultiFab > *  forecast_state_at_lev,
const amrex::Real &  hurricane_eye_latitude,
const amrex::Real &  hurricane_eye_longitude,
amrex::Vector< std::array< amrex::Real, 2 >> &  hurricane_eye_track_xy,
amrex::Vector< std::array< amrex::Real, 2 >> &  hurricane_eye_track_latlon,
amrex::Vector< std::array< amrex::Real, 2 >> &  hurricane_tracker_circle 
)
264 {
265  static bool is_start = true;
266  if(is_start){
267  HurricaneEyeTrackerInitial(geom, S_data,
268  forecast_state_at_lev, hurricane_eye_track_xy,
269  hurricane_eye_track_latlon,
270  hurricane_eye_latitude, hurricane_eye_longitude);
271  is_start = false;
272  } else {
273  HurricaneEyeTrackerNotInitial(geom, S_data, moisture_type,
274  forecast_state_at_lev, hurricane_eye_track_xy,
275  hurricane_eye_track_latlon);
276  }
277  HurricaneTrackerCircle(hurricane_eye_track_xy, hurricane_tracker_circle);
278 }
HurricaneTrackerCircle
void HurricaneTrackerCircle(const amrex::Vector< std::array< amrex::Real, 2 >> &hurricane_eye_track_xy, amrex::Vector< std::array< amrex::Real, 2 >> &hurricane_tracker_circle)
Definition: ERF_HurricaneDiagnostics.H:114
HurricaneEyeTrackerNotInitial
void HurricaneEyeTrackerNotInitial(const amrex::Geometry &geom, const amrex::Vector< amrex::MultiFab > &S_data, MoistureType moisture_type, const amrex::Vector< amrex::MultiFab > *forecast_state_at_lev, amrex::Vector< std::array< amrex::Real, 2 >> &hurricane_eye_track_xy, amrex::Vector< std::array< amrex::Real, 2 >> &hurricane_eye_track_latlon)
Definition: ERF_HurricaneDiagnostics.H:187
HurricaneEyeTrackerInitial
void HurricaneEyeTrackerInitial(const amrex::Geometry &geom, const amrex::Vector< amrex::MultiFab > &S_data, const amrex::Vector< amrex::MultiFab > *forecast_state_at_lev, amrex::Vector< std::array< amrex::Real, 2 >> &hurricane_eye_track_xy, amrex::Vector< std::array< amrex::Real, 2 >> &hurricane_eye_track_latlon, const amrex::Real &hurricane_eye_latitude, const amrex::Real &hurricane_eye_longitude)
Definition: ERF_HurricaneDiagnostics.H:140

Referenced by ERF::post_timestep().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ HurricaneEyeTrackerInitial()

void HurricaneEyeTrackerInitial ( const amrex::Geometry &  geom,
const amrex::Vector< amrex::MultiFab > &  S_data,
const amrex::Vector< amrex::MultiFab > *  forecast_state_at_lev,
amrex::Vector< std::array< amrex::Real, 2 >> &  hurricane_eye_track_xy,
amrex::Vector< std::array< amrex::Real, 2 >> &  hurricane_eye_track_latlon,
const amrex::Real &  hurricane_eye_latitude,
const amrex::Real &  hurricane_eye_longitude 
)
147 {
148  amrex::Gpu::DeviceScalar<amrex::Real> d_val_min(1e10);
149  amrex::Gpu::DeviceScalar<int> d_i_min(-1), d_j_min(-1);
150 
151  amrex::Real* d_val_min_ptr = d_val_min.dataPtr();
152  int* d_i_min_ptr = d_i_min.dataPtr();
153  int* d_j_min_ptr = d_j_min.dataPtr();
154 
155  for (amrex::MFIter mfi(S_data[IntVars::cons]); mfi.isValid(); ++mfi) {
156  const amrex::Box& box = mfi.validbox();
157  const auto latlon_arr = (*forecast_state_at_lev)[4].array(mfi);
158 
159  amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
160  if (k==0) {
161 
162  amrex::Real dlat = latlon_arr(i,j,k,0) - hurricane_eye_latitude;
163  amrex::Real dlon = latlon_arr(i,j,k,1) - hurricane_eye_longitude;
164  amrex::Real dist = std::sqrt(dlat*dlat + dlon*dlon);
165  // Atomic min using device pointer from DeviceVector
166  amrex::Real old = amrex::Gpu::Atomic::Min(&d_val_min_ptr[0], dist);
167  //amrex::Gpu::Atomic::Min(&d_val_min_ptr[0], dist);
168  if (dist < old) {
169  // We are the new minimum; record indices
170  d_i_min_ptr[0] = i;
171  d_j_min_ptr[0] = j;
172  }
173  }
174  });
175  }
176 
177  amrex::Real global_val_min;
178  int global_i_min, global_j_min;
179 
180  ComputeGlobalMinLocation(geom, S_data, forecast_state_at_lev,
181  d_val_min_ptr, d_i_min_ptr, d_j_min_ptr,
182  global_val_min, global_i_min, global_j_min,
183  hurricane_eye_track_xy,
184  hurricane_eye_track_latlon);
185 }
ComputeGlobalMinLocation
void ComputeGlobalMinLocation(const amrex::Geometry &geom, const amrex::Vector< amrex::MultiFab > &S_data, const amrex::Vector< amrex::MultiFab > *forecast_state_at_lev, amrex::Real *d_val_min_ptr, int *d_i_min_ptr, int *d_j_min_ptr, amrex::Real &global_val_min, int &global_i_min, int &global_j_min, amrex::Vector< std::array< amrex::Real, 2 >> &hurricane_eye_track_xy, amrex::Vector< std::array< amrex::Real, 2 >> &hurricane_eye_track_latlon)
Definition: ERF_HurricaneDiagnostics.H:24

Referenced by HurricaneEyeTracker().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ HurricaneEyeTrackerNotInitial()

void HurricaneEyeTrackerNotInitial ( const amrex::Geometry &  geom,
const amrex::Vector< amrex::MultiFab > &  S_data,
MoistureType  moisture_type,
const amrex::Vector< amrex::MultiFab > *  forecast_state_at_lev,
amrex::Vector< std::array< amrex::Real, 2 >> &  hurricane_eye_track_xy,
amrex::Vector< std::array< amrex::Real, 2 >> &  hurricane_eye_track_latlon 
)
193 {
194 
195  if (hurricane_eye_track_xy.empty()) {
196  amrex::Print() << "Error: hurricane_eye_track_xy is empty!\n";
197  amrex::Abort("Attempted to access hurricane_eye_track_xy[0]");
198  }
199 
200  amrex::Real tmp_x_eye = hurricane_eye_track_xy.back()[0];
201  amrex::Real tmp_y_eye = hurricane_eye_track_xy.back()[1];
202 
203  if(amrex::ParallelDescriptor::IOProcessor()){
204  std::cout << "The value of x y are " << tmp_x_eye << " " << tmp_y_eye << std::endl;
205  }
206 
207  amrex::Gpu::DeviceScalar<amrex::Real> d_val_min(1e10);
208  amrex::Gpu::DeviceScalar<int> d_i_min(-1), d_j_min(-1);
209 
210  amrex::Real* d_val_min_ptr = d_val_min.dataPtr();
211  int* d_i_min_ptr = d_i_min.dataPtr();
212  int* d_j_min_ptr = d_j_min.dataPtr();
213 
214  bool use_moisture = (moisture_type != MoistureType::None);
215  const int ncomp = S_data[IntVars::cons].nComp();
216 
217  const auto dx = geom.CellSizeArray();
218  const auto prob_lo = geom.ProbLoArray();
219 
220  for (amrex::MFIter mfi(S_data[IntVars::cons]); mfi.isValid(); ++mfi) {
221  const amrex::Box& box = mfi.validbox();
222  const amrex::Array4<amrex::Real const>& S_arr = S_data[IntVars::cons].const_array(mfi);
223 
224  amrex::ParallelFor(box,[=] AMREX_GPU_DEVICE(int i, int j, int k) {
225  if(k==0) {
226  amrex::Real x = prob_lo[0] + (i+0.5)*dx[0];
227  amrex::Real y = prob_lo[1] + (j+0.5)*dx[1];
228  amrex::Real dist = std::sqrt((x-tmp_x_eye)*(x-tmp_x_eye) + (y-tmp_y_eye)*(y-tmp_y_eye));
229  if(dist < 200e3) {
230  amrex::Real qv_for_p = (use_moisture && (ncomp > RhoQ1_comp)) ? S_arr(i,j,k,RhoQ1_comp)/S_arr(i,j,k,Rho_comp) : 0;
231  const amrex::Real rhotheta = S_arr(i,j,k,RhoTheta_comp);
232  amrex::Real pressure = getPgivenRTh(rhotheta,qv_for_p);
233  amrex::Real old = amrex::Gpu::Atomic::Min(&d_val_min_ptr[0], pressure);
234  //amrex::Gpu::Atomic::Min(&d_val_min_ptr[0], dist);
235  if (old > pressure) {
236  // We are the new minimum; record indices
237  d_i_min_ptr[0] = i;
238  d_j_min_ptr[0] = j;
239  }
240  }
241  }
242  });
243  }
244 
245  amrex::Real global_val_min;
246  int global_i_min, global_j_min;
247 
248  ComputeGlobalMinLocation(geom, S_data, forecast_state_at_lev,
249  d_val_min_ptr, d_i_min_ptr, d_j_min_ptr,
250  global_val_min, global_i_min, global_j_min,
251  hurricane_eye_track_xy,
252  hurricane_eye_track_latlon);
253 }
Coord::y
@ y
Coord::x
@ x
getPgivenRTh
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv=0.)
Definition: ERF_EOS.H:81
Rho_comp
#define Rho_comp
Definition: ERF_IndexDefines.H:36
RhoTheta_comp
#define RhoTheta_comp
Definition: ERF_IndexDefines.H:37
RhoQ1_comp
#define RhoQ1_comp
Definition: ERF_IndexDefines.H:42

Referenced by HurricaneEyeTracker().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ HurricaneMaxVelTracker()

void HurricaneMaxVelTracker ( const amrex::Geometry &  geom,
const amrex::MultiFab &  mf_cc_vel,
const amrex::Real &  time,
const amrex::Vector< std::array< amrex::Real, 2 >> &  hurricane_eye_track_xy,
amrex::Vector< std::array< amrex::Real, 2 >> &  hurricane_maxvel_vs_time 
)
286 {
287  const int ncomp = AMREX_SPACEDIM;
288 
289  amrex::Real* d_val_max_ptr;
290  amrex::Gpu::DeviceVector<amrex::Real> d_val_max(1, -1.0e30);
291  d_val_max_ptr = d_val_max.data();
292 
293  const auto [x_last, y_last] = hurricane_eye_track_xy.back();
294  const auto dx = geom.CellSizeArray();
295  const auto prob_lo = geom.ProbLoArray();
296 
297  amrex::Real x_eye = x_last;
298  amrex::Real y_eye = y_last;
299 
300  for (amrex::MFIter mfi(mf_cc_vel); mfi.isValid(); ++mfi) {
301  const amrex::Box& box = mfi.validbox();
302  const auto& vel_arr = mf_cc_vel.const_array(mfi);
303 
304  amrex::ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
305  amrex::Real x = prob_lo[0] + (i+0.5)*dx[0];
306  amrex::Real y = prob_lo[1] + (j+0.5)*dx[1];
307  amrex::Real dist = std::sqrt((x-x_eye)*(x-x_eye) +
308  (y-y_eye)*(y-y_eye));
309  if(k==1 && dist < 200e3) {
310  amrex::Real velmag = 0.0;
311  for (int comp = 0; comp < ncomp; ++comp) {
312  amrex::Real vel = vel_arr(i, j, k, comp);
313  velmag += vel * vel;
314  }
315  velmag = std::sqrt(velmag)*3.6; // km/hr
316  amrex::Gpu::Atomic::Max(&d_val_max_ptr[0], velmag);
317  }
318  });
319  }
320 
321  amrex::Gpu::synchronize();
322 
323  amrex::Real h_val_max_local = -1.0e30;
324  amrex::Gpu::copy(amrex::Gpu::deviceToHost, d_val_max.begin(), d_val_max.end(), &h_val_max_local);
325 
326  amrex::Real h_val_max_global = -1.0e30;
327  #ifdef AMREX_USE_MPI
328  MPI_Allreduce(&h_val_max_local, &h_val_max_global, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
329  h_val_max_global = h_val_max_local;
330  #else
331  h_val_max_global = h_val_max_local;
332  #endif
333 
334  amrex::Real time_in_hrs = time / 3600.0;
335  hurricane_maxvel_vs_time.push_back({time_in_hrs, h_val_max_global});
336 }

Referenced by ERF::post_timestep().

Here is the caller graph for this function:

◆ HurricaneTrackerCircle()

void HurricaneTrackerCircle ( const amrex::Vector< std::array< amrex::Real, 2 >> &  hurricane_eye_track_xy,
amrex::Vector< std::array< amrex::Real, 2 >> &  hurricane_tracker_circle 
)
116 {
117  // Check that there is at least one eye position
118  if (hurricane_eye_track_xy.empty()) return;
119 
120  // Get the last known (x, y) position of the eye
121  const auto [x_last, y_last] = hurricane_eye_track_xy.back();
122 
123  // Define circle properties
124  const int n_points = 100; // number of points on the circle
125  const amrex::Real radius = 200e3; // radius in meters (example: 50 km)
126 
127  // Clear previous points and reserve space
128  hurricane_tracker_circle.clear();
129  hurricane_tracker_circle.reserve(n_points);
130 
131  // Fill the circle points
132  for (int i = 0; i < n_points; ++i) {
133  amrex::Real theta = 2.0 * M_PI * i / n_points;
134  amrex::Real x = x_last + radius * std::cos(theta);
135  amrex::Real y = y_last + radius * std::sin(theta);
136  hurricane_tracker_circle.push_back({x, y});
137  }
138 }
M_PI
#define M_PI
Definition: ERF_HurricaneDiagnostics.H:16
LsmVar_MM5::theta
@ theta
Definition: ERF_MM5.H:20

Referenced by HurricaneEyeTracker().

Here is the caller graph for this function:

Variable Documentation

◆ 

struct { ... } in

Referenced by ComputeGlobalMinLocation(), ForestDrag::ForestDrag(), ERF::ReadCheckpointFile(), and redistribute_term().

◆ 

struct { ... } out

Referenced by ComputeGlobalMinLocation(), LineSampler::LineSampler(), ERF::print_usage(), redistribute_term(), IRadiation::setRecordDataInfo(), ERF::setRecordDataInfo(), ERF::setRecordDerDataInfo(), ERF::setRecordEnergyDataInfo(), ERF::setRecordSampleLineInfo(), ERF::setRecordSamplePointInfo(), WriteBndryPlanes::write_planes(), ERF::WriteCheckpointFile(), ERF::writeJobInfo(), and ERF::WriteMultiLevelPlotfileWithTerrain().

◆ rank

int rank

Referenced by ComputeGlobalMinLocation().

◆ value

amrex::Real value

Referenced by closest_index(), ParallelBox(), PerpendicularBox(), ncutils::NCVar::put_attr(), ncutils::NCGroup::put_attr(), ReadCustomBinaryIC(), and ERF::refinement_criteria_setup().