ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
Macros | Enumerations | Functions
ERF_WeatherDataInterpolation.cpp File Reference
#include <filesystem>
#include <stdexcept>
#include "ERF.H"
#include "ERF_ReadCustomBinaryIC.H"
#include "ERF_Interpolation_Bilinear.H"
Include dependency graph for ERF_WeatherDataInterpolation.cpp:

Macros

#define ERF_WEATHERDATAINTERPOLATION_H_
 

Enumerations

enum class  BoundType { Lo , Hi }
 
enum class  MultiFabType { CC , NC }
 

Functions

void fill_weather_data_multifab (MultiFab &mf, const Geometry &geom_weather, const int nx, const int ny, const int nz, const Vector< Real > &latvec_h, const Vector< Real > &lonvec_h, const Vector< Real > &zvec_h, const Vector< Real > &rho_h, const Vector< Real > &uvel_h, const Vector< Real > &vvel_h, const Vector< Real > &wvel_h, const Vector< Real > &theta_h, const Vector< Real > &qv_h, const Vector< Real > &qc_h, const Vector< Real > &qr_h)
 
void PlotMultiFab (const MultiFab &mf, const Geometry &geom_mf, const std::string plotfilename, MultiFabType mftype)
 
IntVect find_bound_idx (const Real &x, const Real &y, const Real &z, const BoxList &bl_weather, const Geometry &geom_weather, BoundType bound_type)
 

Macro Definition Documentation

◆ ERF_WEATHERDATAINTERPOLATION_H_

#define ERF_WEATHERDATAINTERPOLATION_H_

Enumeration Type Documentation

◆ BoundType

enum BoundType
strong
Enumerator
Lo 
Hi 
107 { Lo, Hi };
BoundType::Hi
@ Hi
BoundType::Lo
@ Lo

◆ MultiFabType

enum MultiFabType
strong
Enumerator
CC 
NC 
108 { CC, NC };
MultiFabType::NC
@ NC
MultiFabType::CC
@ CC

Function Documentation

◆ fill_weather_data_multifab()

void fill_weather_data_multifab ( MultiFab &  mf,
const Geometry &  geom_weather,
const int  nx,
const int  ny,
const int  nz,
const Vector< Real > &  latvec_h,
const Vector< Real > &  lonvec_h,
const Vector< Real > &  zvec_h,
const Vector< Real > &  rho_h,
const Vector< Real > &  uvel_h,
const Vector< Real > &  vvel_h,
const Vector< Real > &  wvel_h,
const Vector< Real > &  theta_h,
const Vector< Real > &  qv_h,
const Vector< Real > &  qc_h,
const Vector< Real > &  qr_h 
)
34 {
35  const int nx_d = nx;
36  const int ny_d = ny;
37  const int nz_d = nz;
38  const int tot_size = nx*ny*nz;
39 
40  amrex::Gpu::DeviceVector<Real> latvec_d(nx*ny), lonvec_d(nx*ny), zvec_d(nz);
41  amrex::Gpu::DeviceVector<Real> rho_d(tot_size), uvel_d(tot_size),
42  vvel_d(tot_size), wvel_d(tot_size), theta_d(tot_size),
43  qv_d(tot_size), qc_d(tot_size), qr_d(tot_size);
44 
45  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, latvec_h.begin(), latvec_h.end(), latvec_d.begin());
46  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, lonvec_h.begin(), lonvec_h.end(), lonvec_d.begin());
47  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, zvec_h.begin(), zvec_h.end(), zvec_d.begin());
48  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, rho_h.begin(), rho_h.end(), rho_d.begin());
49  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, theta_h.begin(), theta_h.end(), theta_d.begin());
50  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, uvel_h.begin(), uvel_h.end(), uvel_d.begin());
51  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, vvel_h.begin(), vvel_h.end(), vvel_d.begin());
52  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, wvel_h.begin(), wvel_h.end(), wvel_d.begin());
53  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, qv_h.begin(), qv_h.end(), qv_d.begin());
54  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, qc_h.begin(), qc_h.end(), qc_d.begin());
55  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, qr_h.begin(), qr_h.end(), qr_d.begin());
56 
57  Real* latvec_d_ptr = latvec_d.data();
58  Real* lonvec_d_ptr = lonvec_d.data();
59  Real* zvec_d_ptr = zvec_d.data();
60  Real* rho_d_ptr = rho_d.data();
61  Real* uvel_d_ptr = uvel_d.data();
62  Real* vvel_d_ptr = vvel_d.data();
63  Real* wvel_d_ptr = wvel_d.data();
64  Real* theta_d_ptr = theta_d.data();
65  Real* qv_d_ptr = qv_d.data();
66  Real* qc_d_ptr = qc_d.data();
67  Real* qr_d_ptr = qr_d.data();
68 
69  const auto prob_lo = geom_weather.ProbLoArray();
70  const auto dx = geom_weather.CellSizeArray();
71 
72  for (MFIter mfi(mf, TilingIfNotGPU()); mfi.isValid(); ++mfi)
73  {
74  const Box& bx = mfi.nodaltilebox();
75  Array4<Real> const& arr = mf.array(mfi);
76 
77  ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
78  {
79  const Real z = prob_lo[2] + k * dx[2];
80  int kloc = -1;
81  for (int kk = 0; kk < nz_d; kk++) {
82  if (zvec_d_ptr[kk] > z) {
83  kloc = kk;
84  break;
85  }
86  }
87  // Replace this with logic using your coarse input data
88  int idx1 = get_single_index(i,j,kloc-1,nx_d,ny_d);
89  int idx2 = get_single_index(i,j,kloc,nx_d,ny_d);
90  Real dz = zvec_d_ptr[kloc] - zvec_d_ptr[kloc-1];
91  Real fac = (z-zvec_d_ptr[kloc-1])/dz;
92  arr(i,j,k,0) = rho_d_ptr[idx1] + fac*(rho_d_ptr[idx2]-rho_d_ptr[idx1]);
93  arr(i,j,k,1) = uvel_d_ptr[idx1] + fac*(uvel_d_ptr[idx2]-uvel_d_ptr[idx1]);
94  arr(i,j,k,2) = vvel_d_ptr[idx1] + fac*(vvel_d_ptr[idx2]-vvel_d_ptr[idx1]);
95  arr(i,j,k,3) = wvel_d_ptr[idx1] + fac*(wvel_d_ptr[idx2]-wvel_d_ptr[idx1]);
96  arr(i,j,k,4) = theta_d_ptr[idx1] + fac*(theta_d_ptr[idx2]-theta_d_ptr[idx1]);
97  arr(i,j,k,5) = qv_d_ptr[idx1] + fac*(qv_d_ptr[idx2]-qv_d_ptr[idx1]);
98  arr(i,j,k,6) = qc_d_ptr[idx1] + fac*(qc_d_ptr[idx2]-qc_d_ptr[idx1]);
99  arr(i,j,k,7) = qr_d_ptr[idx1] + fac*(qr_d_ptr[idx2]-qr_d_ptr[idx1]);
100  idx1 = get_single_index(i,j,0,nx_d,ny_d);
101  arr(i,j,k,8) = latvec_d_ptr[idx1];
102  arr(i,j,k,9) = lonvec_d_ptr[idx1];
103  });
104  }
105 }
Coord::z
@ z
get_single_index
AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE int get_single_index(int i, int j, int k, int nx, int ny)
Definition: ERF_Interpolation_Bilinear.H:13
Real
amrex::Real Real
Definition: ERF_ShocInterface.H:16

Referenced by ERF::FillWeatherDataMultiFab().

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

◆ find_bound_idx()

IntVect find_bound_idx ( const Real &  x,
const Real &  y,
const Real &  z,
const BoxList &  bl_weather,
const Geometry &  geom_weather,
BoundType  bound_type 
)
220 {
221  const auto prob_lo_weather = geom_weather.ProbLoArray();
222  const auto dx_weather = geom_weather.CellSizeArray();
223 
224  int i, j, k;
225 
226  if (bound_type == BoundType::Lo) {
227  i = static_cast<int>(std::floor((x - prob_lo_weather[0]) / dx_weather[0]));
228  j = static_cast<int>(std::floor((y - prob_lo_weather[1]) / dx_weather[1]));
229  k = static_cast<int>(std::floor((z - prob_lo_weather[2]) / dx_weather[2]));
230  } else { // BoundType::Hi
231  i = static_cast<int>(std::ceil((x - prob_lo_weather[0]) / dx_weather[0]));
232  j = static_cast<int>(std::ceil((y - prob_lo_weather[1]) / dx_weather[1]));
233  k = static_cast<int>(std::ceil((z - prob_lo_weather[2]) / dx_weather[2]));
234  }
235 
236  IntVect idx(i, j, k);
237 
238  for (const auto& b : bl_weather) {
239  if (b.contains(idx)) {
240  return idx;
241  }
242  }
243 
244  amrex::Abort("Bound index not found in any box in BoxList!");
245  return IntVect::TheZeroVector(); // unreachable if Abort
246 }

Referenced by ERF::InterpWeatherDataOntoMesh().

Here is the caller graph for this function:

◆ PlotMultiFab()

void PlotMultiFab ( const MultiFab &  mf,
const Geometry &  geom_mf,
const std::string  plotfilename,
MultiFabType  mftype 
)
115 {
116 
117  Vector<std::string> varnames = {
118  "rho", "uvel", "vvel", "wvel", "theta", "qv", "qc", "qr", "latitude", "longitude"
119  }; // Customize variable names
120 
121  const Real time = 0.0;
122 
123 
124  // Assume weather_mf is nodal in all directions
125  if(mftype == MultiFabType::NC) {
126  BoxArray cba = mf.boxArray();
127  cba = amrex::convert(mf.boxArray(), IntVect::TheCellVector());
128 
129  MultiFab cc_mf(cba, mf.DistributionMap(),
130  mf.nComp(), 0);
131 
132  amrex::average_node_to_cellcenter(cc_mf, 0, mf, 0, mf.nComp());
133 
134  WriteSingleLevelPlotfile(
135  plotfilename,
136  cc_mf,
137  varnames,
138  geom_mf,
139  time,
140  0 // level
141  );
142  } else {
143  WriteSingleLevelPlotfile(
144  plotfilename,
145  mf,
146  varnames,
147  geom_mf,
148  time,
149  0 // level
150  );
151  }
152 }

Referenced by ERF::FillWeatherDataMultiFab(), and ERF::InterpWeatherDataOntoMesh().

Here is the caller graph for this function: