Skip to content

Commit dc89817

Browse files
ewquonewquon
authored
Implement ReadTimeSliceFromNetCDFFile, update time output precision (#2192)
* Add fractional seconds to timestamp output * Increase output time precision throughout, cleanup verbose output * Tell user that we've set the start time from the wrfbdy * Check that first dimension is Time * Implement ReadTimeSliceFromNetCDFFile * Replace sprintf with snprintf
1 parent 3d88325 commit dc89817

17 files changed

+195
-54
lines changed

‎Exec/MoistRegTests/WPS_Test/inputs_real_ChisholmView

-1
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,5 @@
11
# ------------------ INPUTS TO MAIN PROGRAM -------------------
22
max_step = 1
3-
start_time = 1385042400.
43

54
amrex.fpe_trap_invalid = 1
65

‎Source/ERF.H

+7
Original file line numberDiff line numberDiff line change
@@ -916,13 +916,20 @@ private:
916916
int last_check_file_step;
917917
int plot_file_on_restart = 1;
918918

919+
// Output control
920+
const int datwidth = 14;
921+
const int datprecision = 6;
922+
const int timeprecision = 13; // e.g., 1-yr LES: 31,536,000 s with dt ~ 0.01 ==> min prec = 10
923+
// epoch time for 2030-01-01 00:00:00 is 1,893,456,000 s with dt ~ 0.01 ==> min prec = 12
924+
919925
////////////////
920926
// runtime parameters
921927

922928
// maximum number of steps and stop time
923929
int max_step = std::numeric_limits<int>::max();
924930
amrex::Real start_time = 0.0;
925931
amrex::Real stop_time = std::numeric_limits<amrex::Real>::max();
932+
926933
bool use_datetime = false;
927934
const std::string datetime_format = "%Y-%m-%d %H:%M:%S";
928935

‎Source/ERF.cpp

+11-4
Original file line numberDiff line numberDiff line change
@@ -412,9 +412,8 @@ ERF::Evolve ()
412412
for (int step = istep[0]; step < max_step && cur_time < stop_time; ++step)
413413
{
414414
if (use_datetime) {
415-
Print() << "\n" << getTimestamp(static_cast<std::time_t>(cur_time),
416-
datetime_format)
417-
<< " (" << cur_time-start_time << " s elapsed)"
415+
Print() << "\n" << getTimestamp(cur_time, datetime_format)
416+
<< " (" << cur_time-start_time << " s elapsed)"
418417
<< std::endl;
419418
}
420419
Print() << "\nCoarse STEP " << step+1 << " starts ..." << std::endl;
@@ -1419,6 +1418,15 @@ ERF::init_only (int lev, Real time)
14191418
// The base state is initialized from WRF wrfinput data, output by
14201419
// ideal.exe or real.exe
14211420
init_from_wrfinput(lev);
1421+
if (lev==0) {
1422+
if ((start_time > 0) && (start_time != t_new[lev])) {
1423+
Print() << "Ignoring specified start_time="
1424+
<< std::setprecision(timeprecision) << start_time
1425+
<< std::endl;
1426+
}
1427+
start_time = t_new[lev];
1428+
}
1429+
use_datetime = true;
14221430

14231431
// The physbc's need the terrain but are needed for initHSE
14241432
if (!solverChoice.use_real_bcs) {
@@ -1492,7 +1500,6 @@ ERF::ReadParameters ()
14921500

14931501
std::string start_datetime, stop_datetime;
14941502
if (pp.query("start_datetime", start_datetime)) {
1495-
// Both start and stop datetimes must be provided
14961503
start_time = getEpochTime(start_datetime, datetime_format);
14971504
if (pp.query("stop_datetime", stop_datetime)) {
14981505
stop_time = getEpochTime(stop_datetime, datetime_format);

‎Source/IO/ERF_NCInterface.H

+3
Original file line numberDiff line numberDiff line change
@@ -66,6 +66,9 @@ struct NCVar
6666
//! Shape of the array (size in each array dimension)
6767
[[nodiscard]] std::vector<size_t> shape () const;
6868

69+
//! Dimension names of the array (str in each array dimension)
70+
[[nodiscard]] std::vector<std::string> dimnames () const;
71+
6972
//! Write out the entire variable
7073
void put (const double*) const;
7174
void put (const float*) const;

‎Source/IO/ERF_NCInterface.cpp

+21
Original file line numberDiff line numberDiff line change
@@ -64,6 +64,27 @@ int NCVar::ndim () const
6464
return ndims;
6565
}
6666

67+
/**
68+
* Error-checking function to get the name of each dimension from a NetCDF identity
69+
*/
70+
std::vector<std::string> NCVar::dimnames () const
71+
{
72+
int ndims = ndim();
73+
std::vector<int> dimids(ndims);
74+
std::vector<std::string> vnames(ndims);
75+
76+
for (int i = 0; i < ndims; ++i)
77+
check_nc_error(nc_inq_vardimid(ncid, varid, dimids.data()));
78+
79+
char buf[80];
80+
for (int i = 0; i < ndims; ++i) {
81+
check_nc_error(nc_inq_dimname(ncid, dimids[i], buf));
82+
vnames[i] = std::string(buf);
83+
}
84+
85+
return vnames;
86+
}
87+
6788
/**
6889
* Error-checking function to get the length of each dimension from a NetCDF identity
6990
*/

‎Source/IO/ERF_NCWpsFile.H

+64-7
Original file line numberDiff line numberDiff line change
@@ -128,6 +128,53 @@ int BuildFABsFromWRFBdyFile (const std::string &fname,
128128
amrex::Vector<amrex::Vector<amrex::FArrayBox>>& bdy_data_ylo,
129129
amrex::Vector<amrex::Vector<amrex::FArrayBox>>& bdy_data_yhi);
130130

131+
template<typename DType>
132+
void ReadTimeSliceFromNetCDFFile(const std::string& fname,
133+
const int tidx,
134+
amrex::Vector<std::string> names,
135+
amrex::Vector<NDArray<DType> >& arrays,
136+
amrex::Vector<int>& success)
137+
{
138+
amrex::Print() << "Reading time slice " << tidx << " from " << fname << std::endl;
139+
AMREX_ASSERT(arrays.size() == names.size());
140+
141+
if (amrex::ParallelDescriptor::IOProcessor())
142+
{
143+
auto ncf = ncutils::NCFile::open(fname, NC_CLOBBER | NC_NETCDF4);
144+
145+
int ntimes = ncf.dim("Time").len();
146+
AMREX_ALWAYS_ASSERT((tidx >= 0) && (tidx < ntimes));
147+
148+
for (auto n=0; n<arrays.size(); ++n)
149+
{
150+
std::string vname_to_write = names[n];
151+
std::string vname_to_read = names[n];
152+
if (vname_to_read.substr(0,2) == "R_") {
153+
vname_to_read = names[n+4]; // This allows us to read "T" instead -- we will over-write this later
154+
}
155+
156+
success[n] = ncf.has_var(vname_to_read);
157+
158+
if (success[n] == 1)
159+
{
160+
std::vector<std::string> dimnames = ncf.var(vname_to_read).dimnames();
161+
AMREX_ALWAYS_ASSERT(dimnames[0] == "Time");
162+
163+
std::vector<size_t> count = ncf.var(vname_to_read).shape();
164+
std::vector<size_t> start(count.size(), 0);
165+
start[0] = tidx;
166+
count[0] = 1;
167+
168+
arrays[n] = NDArray<DType>(vname_to_read, count);
169+
DType* dataPtr = arrays[n].get_data();
170+
171+
ncf.var(vname_to_read).get(dataPtr, start, count);
172+
} // has_var
173+
}
174+
ncf.close();
175+
}
176+
}
177+
131178
template<typename DType>
132179
void ReadNetCDFFile (const std::string& fname, amrex::Vector<std::string> names,
133180
amrex::Vector<NDArray<DType> >& arrays, amrex::Vector<int>& success)
@@ -172,19 +219,29 @@ void ReadNetCDFFile (const std::string& fname, amrex::Vector<std::string> names,
172219

173220
if (success[n] == 1) {
174221

222+
std::vector<std::string> dimnames = ncf.var(vname_to_read).dimnames();
223+
AMREX_ALWAYS_ASSERT(dimnames[0] == "Time");
224+
175225
std::vector<size_t> shape = ncf.var(vname_to_read).shape();
176226
arrays[n] = NDArray<DType>(vname_to_read,shape);
177227
DType* dataPtr = arrays[n].get_data();
178228

179229
std::vector<size_t> start(shape.size(), 0);
180230

181-
// auto numPts = arrays[n].ndim();
182-
// amrex::Print() << "NetCDF Variable name = " << vname_to_read << std::endl;
183-
// amrex::Print() << "numPts read from NetCDF file/var = " << numPts << std::endl;
184-
// amrex::Print() << "Dims of the variable = (";
185-
// for (auto &dim:shape)
186-
// amrex::Print() << dim << ", " ;
187-
// amrex::Print() << ")" << std::endl;
231+
#if 0
232+
auto numPts = arrays[n].ndim();
233+
amrex::Print() << "NetCDF Variable name = " << vname_to_read << std::endl;
234+
amrex::Print() << "numPts read from NetCDF file/var = " << numPts << std::endl;
235+
amrex::Print() << "Dims in var = " << ncf.var(vname_to_read).ndim() << std::endl;
236+
amrex::Print() << "Dim names = (";
237+
for (auto &dim:dimnames)
238+
amrex::Print() << dim << ", " ;
239+
amrex::Print() << ")" << std::endl;
240+
amrex::Print() << "Dims of the variable = (";
241+
for (auto &dim:shape)
242+
amrex::Print() << dim << ", " ;
243+
amrex::Print() << ")" << std::endl;
244+
#endif
188245

189246
ncf.var(vname_to_read).get(dataPtr, start, shape);
190247

‎Source/IO/ERF_ReadFromWRFBdy.cpp

+45
Original file line numberDiff line numberDiff line change
@@ -71,6 +71,7 @@ read_from_wrfbdy (const std::string& nc_bdy_file, const Box& domain,
7171
for (int nt(0); nt < ntimes; nt++) {
7272
std::string date(&timeStamps[nt][0], &timeStamps[nt][dateStrLen-1]+1);
7373
auto epochTime = getEpochTime(date, dateTimeFormat);
74+
Print() << " wrfbdy datetime " << nt << " : " << date << " " << epochTime << std::endl;
7475
epochTimes.push_back(epochTime);
7576

7677
if (nt == 1) {
@@ -120,11 +121,55 @@ read_from_wrfbdy (const std::string& nc_bdy_file, const Box& domain,
120121
if (ParallelDescriptor::IOProcessor())
121122
{
122123
Vector<int> success(nc_var_names.size());
124+
125+
#if 1
126+
for (int itime=0; itime < ntimes; ++itime) {
127+
Vector<RARRAY> tslice(nc_var_names.size());
128+
ReadTimeSliceFromNetCDFFile(nc_bdy_file, itime, nc_var_names, tslice, success);
129+
130+
/*
131+
* DEV TEST: Copy tslice into full array to verify result
132+
*/
133+
for (int n=0; n < nc_var_names.size(); ++n) {
134+
std::vector<size_t> vshape = tslice[n].get_vshape();
135+
size_t offset{1};
136+
for (auto &dim:vshape) offset *= dim;
137+
138+
Print() << "Time " << itime << " " << nc_var_names[n] << "(";
139+
for (auto &dim:vshape)
140+
amrex::Print() << dim << ",";
141+
Print() << ") offset=" << offset << std::endl;
142+
143+
if (itime==0) {
144+
// allocate full array
145+
vshape[0] = ntimes;
146+
arrays[n] = NDArray<float>(tslice[n].get_vname(), vshape);
147+
}
148+
149+
float* fullPtr = arrays[n].get_data();
150+
float* slicePtr = tslice[n].get_data();
151+
std::copy(slicePtr, slicePtr+offset, fullPtr+itime*offset);
152+
}
153+
}
154+
155+
for (auto &istat:success) {
156+
AMREX_ALWAYS_ASSERT(istat==1);
157+
}
158+
159+
//for (int n=0; n < nc_var_names.size(); ++n) {
160+
// std::vector<size_t> shape = arrays[n].get_vshape();
161+
// Print() << "Read " << nc_var_names[n] << "(";
162+
// for (auto &dim:shape)
163+
// amrex::Print() << dim << ",";
164+
// Print() << ")" << std::endl;
165+
//}
166+
#else
123167
ReadNetCDFFile(nc_bdy_file, nc_var_names, arrays, success);
124168
125169
// Assert that the data has the same number of time snapshots
126170
int itimes = static_cast<int>(arrays[0].get_vshape()[0]);
127171
AMREX_ALWAYS_ASSERT(itimes == ntimes);
172+
#endif
128173

129174
// Width of the boundary region
130175
width = arrays[0].get_vshape()[1];

‎Source/IO/ERF_SampleData.H

+1-1
Original file line numberDiff line numberDiff line change
@@ -300,7 +300,7 @@ struct LineSampler
300300
void
301301
write_line_ascii (amrex::Vector<amrex::Real>& time)
302302
{
303-
// same as write_1D_profiles
303+
// same as definitions in ERF.H
304304
constexpr int datwidth = 14;
305305
constexpr int datprecision = 6;
306306
constexpr int timeprecision = 13;

‎Source/IO/ERF_Write1DProfiles.cpp

-4
Original file line numberDiff line numberDiff line change
@@ -16,10 +16,6 @@ ERF::write_1D_profiles (Real time)
1616
{
1717
BL_PROFILE("ERF::write_1D_profiles()");
1818

19-
int datwidth = 14;
20-
int datprecision = 6;
21-
int timeprecision = 13; // e.g., 1-yr LES: 31,536,000 s with dt ~ 0.01 ==> min prec = 10
22-
2319
if (NumDataLogs() > 1)
2420
{
2521
// Define the 1d arrays we will need

‎Source/IO/ERF_Write1DProfiles_stag.cpp

-4
Original file line numberDiff line numberDiff line change
@@ -26,10 +26,6 @@ ERF::write_1D_profiles_stag (Real time)
2626
{
2727
BL_PROFILE("ERF::write_1D_profiles()");
2828

29-
int datwidth = 14;
30-
int datprecision = 6;
31-
int timeprecision = 13; // e.g., 1-yr LES: 31,536,000 s with dt ~ 0.01 ==> min prec = 10
32-
3329
if (NumDataLogs() > 1)
3430
{
3531
// Define the 1d arrays we will need

‎Source/IO/ERF_WriteScalarProfiles.cpp

+9-21
Original file line numberDiff line numberDiff line change
@@ -19,10 +19,6 @@ ERF::sum_integrated_quantities (Real time)
1919
if (verbose <= 0)
2020
return;
2121

22-
int datwidth = 14;
23-
int datprecision = 6;
24-
int timeprecision = 13; // e.g., 1-yr LES: 31,536,000 s with dt ~ 0.01 ==> min prec = 10
25-
2622
// Single level sum
2723
Real mass_sl;
2824

@@ -107,22 +103,23 @@ ERF::sum_integrated_quantities (Real time)
107103
scal_ml = foo[i++];
108104

109105
Print() << '\n';
106+
Print() << "TIME= " << std::setw(datwidth) << std::setprecision(timeprecision) << std::left << time << '\n';
110107
if (finest_level == 0) {
111108
#if 1
112-
Print() << "TIME= " << time << " MASS = " << mass_sl << '\n';
109+
Print() << " MASS = " << mass_sl << '\n';
113110
#else
114-
Print() << "TIME= " << time << " PERT MASS = " << mass_sl << '\n';
111+
Print() << " PERT MASS = " << mass_sl << '\n';
115112
#endif
116-
Print() << "TIME= " << time << " RHO THETA = " << rhth_sl << '\n';
117-
Print() << "TIME= " << time << " RHO SCALAR = " << scal_sl << '\n';
113+
Print() << " RHO THETA = " << rhth_sl << '\n';
114+
Print() << " RHO SCALAR = " << scal_sl << '\n';
118115
} else {
119116
#if 1
120-
Print() << "TIME= " << time << " MASS SL/ML = " << mass_sl << " " << mass_ml << '\n';
117+
Print() << " MASS SL/ML = " << mass_sl << " " << mass_ml << '\n';
121118
#else
122-
Print() << "TIME= " << time << " PERT MASS SL/ML = " << mass_sl << " " << mass_ml << '\n';
119+
Print() << " PERT MASS SL/ML = " << mass_sl << " " << mass_ml << '\n';
123120
#endif
124-
Print() << "TIME= " << time << " RHO THETA SL/ML = " << rhth_sl << " " << rhth_ml << '\n';
125-
Print() << "TIME= " << time << " RHO SCALAR SL/ML = " << scal_sl << " " << scal_ml << '\n';
121+
Print() << " RHO THETA SL/ML = " << rhth_sl << " " << rhth_ml << '\n';
122+
Print() << " RHO SCALAR SL/ML = " << scal_sl << " " << scal_ml << '\n';
126123
}
127124

128125
// The first data log only holds scalars
@@ -175,10 +172,6 @@ ERF::sum_derived_quantities (Real time)
175172

176173
int lev = 0;
177174

178-
int datwidth = 14;
179-
int datprecision = 6;
180-
int timeprecision = 13; // e.g., 1-yr LES: 31,536,000 s with dt ~ 0.01 ==> min prec = 10
181-
182175
AMREX_ALWAYS_ASSERT(lev == 0);
183176

184177
// ************************************************************************
@@ -332,8 +325,6 @@ ERF::cloud_fraction (Real /*time*/)
332325
void
333326
ERF::sample_points (int /*lev*/, Real time, IntVect cell, MultiFab& mf)
334327
{
335-
int datwidth = 14;
336-
337328
int ifile = 0;
338329

339330
//
@@ -370,9 +361,6 @@ ERF::sample_points (int /*lev*/, Real time, IntVect cell, MultiFab& mf)
370361
void
371362
ERF::sample_lines (int lev, Real time, IntVect cell, MultiFab& mf)
372363
{
373-
int datwidth = 14;
374-
int datprecision = 6;
375-
376364
int ifile = 0;
377365

378366
const int ncomp = mf.nComp(); // cell-centered state vars

‎Source/Initialization/ERF_InitFromWRFInput.cpp

+3
Original file line numberDiff line numberDiff line change
@@ -705,6 +705,9 @@ ERF::init_from_wrfinput (int lev)
705705
// Note we only have start_bdy_time if at level 0 and init_type == InitType:WRFInput
706706
//
707707
if (lev == 0) {
708+
Print() << "Setting start_time to "
709+
<< std::setprecision(timeprecision) << start_bdy_time
710+
<< " from wrfbdy" << std::endl;
708711
t_new[lev] = start_bdy_time;
709712
t_old[lev] = start_bdy_time - 1.e200;
710713
} else {

‎Source/TimeIntegration/ERF_TI_no_substep_fun.H

+1
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@
1717
const amrex::GpuArray<int, IntVars::NumTypes> ncomp_fast = {2,1,1,1};
1818

1919
if (verbose) amrex::Print() << " No-substepping time integration at level " << level
20+
<< std::setprecision(timeprecision)
2021
<< " to " << time_for_fp
2122
<< " with dt = " << slow_dt << std::endl;
2223

0 commit comments

Comments
 (0)