28#include "pism/util/array/Array.hh"
29#include "pism/util/array/Array_impl.hh"
31#include "pism/util/ConfigInterface.hh"
32#include "pism/util/Grid.hh"
33#include "pism/util/Time.hh"
35#include "pism/util/Context.hh"
36#include "pism/util/Logger.hh"
37#include "pism/util/Profiling.hh"
38#include "pism/util/VariableMetadata.hh"
39#include "pism/util/error_handling.hh"
40#include "pism/util/io/File.hh"
41#include "pism/util/io/IO_Flags.hh"
42#include "pism/util/io/io_helpers.hh"
43#include "pism/util/petscwrappers/VecScatter.hh"
44#include "pism/util/petscwrappers/Viewer.hh"
45#include "pism/util/pism_utilities.hh"
47#include "pism/util/InputInterpolation.hh"
56 ierr = DMGlobalToLocalBegin(dm, source, INSERT_VALUES, destination);
57 PISM_CHK(ierr,
"DMGlobalToLocalBegin");
59 ierr = DMGlobalToLocalEnd(dm, source, INSERT_VALUES, destination);
60 PISM_CHK(ierr,
"DMGlobalToLocalEnd");
63Array::Array(std::shared_ptr<const Grid> grid,
const std::string &name,
Kind ghostedp,
size_t dof,
64 size_t stencil_width,
const std::vector<double> &zlevels) {
74 const auto &config =
grid->ctx()->config();
76 auto max_stencil_width =
static_cast<size_t>(config->get_number(
"grid.max_stencil_width"));
85 auto system =
m_impl->
grid->ctx()->unit_system();
88 for (
unsigned int j = 0; j <
m_impl->
dof; ++j) {
96 if (zlevels.size() > 1) {
63Array::Array(std::shared_ptr<const Grid> grid,
const std::string &name,
Kind ghostedp,
size_t dof, {
…}
199 return NORM_INFINITY;
209 PetscErrorCode ierr = VecAXPY(
vec(), alpha, x.
vec());
217 PetscErrorCode ierr = VecShift(
vec(), alpha);
225 PetscErrorCode ierr = VecScale(
vec(), alpha);
245 this->
get_dof(destination_da, destination, 0, N);
249 unsigned int count)
const {
256 double ***result_a =
static_cast<double ***
>(tmp_res.get()),
257 ***source_a =
static_cast<double ***
>(tmp_v.
get());
261 for (
auto p =
m_impl->
grid->points(); p; p.next()) {
262 const int i = p.i(), j = p.j();
263 PetscErrorCode ierr =
264 PetscMemcpy(result_a[j][i], &source_a[j][i][start],
count *
sizeof(PetscScalar));
274 unsigned int count) {
281 double ***source_a =
static_cast<double ***
>(tmp_src.get()),
282 ***result_a =
static_cast<double ***
>(tmp_v.
get());
286 for (
auto p =
m_impl->
grid->points(); p; p.next()) {
287 const int i = p.i(), j = p.j();
288 PetscErrorCode ierr =
289 PetscMemcpy(&result_a[j][i][start], source_a[j][i],
count *
sizeof(PetscScalar));
312 PetscErrorCode ierr = 0;
315 PISM_CHK(ierr,
"DMCreateLocalVector");
318 PISM_CHK(ierr,
"DMCreateGlobalVector");
359 const Logger &log, Vec output) {
361 if (default_value.
exists()) {
363 log.
message(2,
" variable %s (%s) is not found: using the default constant\n",
366 PetscErrorCode ierr = VecSet(output, default_value);
380 assert(
ndims() == 2);
382 auto log =
grid()->ctx()->log();
386 auto dm_for_io =
ndof() == 1 ?
dm() :
grid()->get_dm(1, 0);
391 for (
unsigned int j = 0; j <
ndof(); ++j) {
394 auto V = file.
find_variable(variable.get_name(), variable[
"standard_name"]);
399 int last_record = -1;
400 interp->regrid(variable, file, last_record, *
grid(), tmp);
418 auto log =
grid()->ctx()->log();
419 log->message(4,
" Reading %s...\n",
m_impl->
name.c_str());
439 auto da2 =
grid()->get_dm(1, 0);
445 for (
unsigned int j = 0; j <
ndof(); ++j) {
464 for (
unsigned int j = 0; j <
ndof(); ++j) {
477 assert(N < m_impl->dof);
482 assert(N < m_impl->dof);
493 log->message(3,
"[%s] Writing %s...\n",
519 for (
unsigned int j = 0; j <
ndof(); ++j) {
523 log->message(3,
"[%s] Writing %s...\n",
548 PetscInt X_size, Y_size;
555 ierr = VecGetSize(
vec(), &X_size);
558 ierr = VecGetSize(other.
vec(), &Y_size);
561 if (X_size != Y_size) {
578 PISM_CHK(ierr,
"DMDAVecGetArrayDOF");
594 "Array::end_access(): a == NULL (looks like begin_access() was not called)");
605 PISM_CHK(ierr,
"DMDAVecRestoreArrayDOF");
608 PISM_CHK(ierr,
"DMDAVecRestoreArray");
621 ierr = DMLocalToLocalBegin(*
dm(),
vec(), INSERT_VALUES,
vec());
622 PISM_CHK(ierr,
"DMLocalToLocalBegin");
624 ierr = DMLocalToLocalEnd(*
dm(),
vec(), INSERT_VALUES,
vec());
625 PISM_CHK(ierr,
"DMLocalToLocalEnd");
630 PetscErrorCode ierr = VecSet(
vec(),c);
637 double ghost_width = 0;
647 bool out_of_range = (i <
m_impl->
grid->xs() - ghost_width) ||
649 (j < m_impl->
grid->ys() - ghost_width) ||
674 PetscErrorCode ierr = VecStrideNormAll(
vec(), type, result.data());
677 PetscErrorCode ierr = VecNorm(
vec(), type, result.data());
687 "Array::norm_all(...): NORM_1_AND_2"
688 " not implemented (called as %s.norm_all(...))",
704 case NORM_INFINITY: {
712 "Array::norm_all(...): unknown norm type"
713 " (called as %s.norm_all(...))",
733 this->
read(file, time);
738 this->
regrid(file, default_value);
743 const std::string &filename,
746 PetscErrorCode ierr = 0;
751 ierr = VecMin(v, NULL, &
min);
753 ierr = VecMax(v, NULL, &
max);
795 m_impl->
grid->ctx()->profiling().begin(
"io.regridding");
801 for (
unsigned k = 0;
k <
ndof(); ++
k) {
802 auto dm =
grid()->get_dm(1, 0);
820 m_impl->
grid->ctx()->profiling().end(
"io.regridding");
825 if (time_spent > 1.0) {
826 log->message(3,
" done in %f seconds.\n", time_spent);
828 log->message(3,
" done.\n");
848 megabyte = pow(2, 20),
849 N =
static_cast<double>(
size()),
850 mb_double =
static_cast<double>(
sizeof(
double)) * N / megabyte,
851 mb_float =
static_cast<double>(
sizeof(
float)) * N / megabyte;
854 std::string spacer(
timestamp.size(),
' ');
855 if (time_spent > 1) {
858 " [%s] Done writing %s (%f Mb double, %f Mb float)\n"
859 " %s in %f seconds (%f minutes).\n"
860 " %s Effective throughput: double: %f Mb/s, float: %f Mb/s.\n",
862 spacer.c_str(), time_spent, time_spent / minute,
864 mb_double / time_spent, mb_float / time_spent);
866 m_impl->
grid->ctx()->log()->message(3,
" done.\n");
875 while (not
m_vecs.empty()) {
877 m_vecs.back()->end_access();
886 for (
const auto *j : vecs) {
887 assert(j !=
nullptr);
902 for (
const auto *v : vecs) {
903 assert(v !=
nullptr);
921 return Mx * My * Mz * dof;
931 ierr = PetscObjectQuery((PetscObject)
dm()->get(),
"v_proc0", (PetscObject*)&v_proc0);
934 if (v_proc0 == NULL) {
941 ierr = DMDACreateNaturalVector(*
dm(), natural_work.
rawptr());
942 PISM_CHK(ierr,
"DMDACreateNaturalVector");
945 ierr = PetscObjectCompose((PetscObject)
dm()->get(),
"natural_work",
946 (PetscObject)((::Vec)natural_work));
947 PISM_CHK(ierr,
"PetscObjectCompose");
955 ierr = VecScatterCreateToZero(natural_work, scatter_to_zero.
rawptr(),
957 PISM_CHK(ierr,
"VecScatterCreateToZero");
960 ierr = PetscObjectCompose((PetscObject)
dm()->get(),
"scatter_to_zero",
961 (PetscObject)((::VecScatter)scatter_to_zero));
962 PISM_CHK(ierr,
"PetscObjectCompose");
965 ierr = PetscObjectCompose((PetscObject)
dm()->get(),
"v_proc0",
966 (PetscObject)v_proc0);
967 PISM_CHK(ierr,
"PetscObjectCompose");
975 ierr = VecDuplicate(v_proc0, &result);
978 return std::shared_ptr<petsc::Vec>(
new petsc::Vec(result));
982 PetscErrorCode ierr = 0;
983 VecScatter scatter_to_zero = NULL;
984 Vec natural_work = NULL;
986 ierr = PetscObjectQuery((PetscObject)
dm()->get(),
"scatter_to_zero",
987 (PetscObject*)&scatter_to_zero);
990 ierr = PetscObjectQuery((PetscObject)
dm()->get(),
"natural_work",
991 (PetscObject*)&natural_work);
994 if (natural_work == NULL || scatter_to_zero == NULL) {
998 ierr = DMDAGlobalToNaturalBegin(*
dm(), parallel, INSERT_VALUES, natural_work);
999 PISM_CHK(ierr,
"DMDAGlobalToNaturalBegin");
1001 ierr = DMDAGlobalToNaturalEnd(*
dm(), parallel, INSERT_VALUES, natural_work);
1002 PISM_CHK(ierr,
"DMDAGlobalToNaturalEnd");
1004 ierr = VecScatterBegin(scatter_to_zero, natural_work, onp0,
1005 INSERT_VALUES, SCATTER_FORWARD);
1008 ierr = VecScatterEnd(scatter_to_zero, natural_work, onp0,
1009 INSERT_VALUES, SCATTER_FORWARD);
1026 PetscErrorCode ierr;
1028 VecScatter scatter_to_zero = NULL;
1029 Vec natural_work = NULL;
1030 ierr = PetscObjectQuery((PetscObject)
dm()->get(),
"scatter_to_zero",
1031 (PetscObject*)&scatter_to_zero);
1032 PISM_CHK(ierr,
"PetscObjectQuery");
1034 ierr = PetscObjectQuery((PetscObject)
dm()->get(),
"natural_work",
1035 (PetscObject*)&natural_work);
1036 PISM_CHK(ierr,
"PetscObjectQuery");
1038 if (natural_work == NULL || scatter_to_zero == NULL) {
1042 ierr = VecScatterBegin(scatter_to_zero, onp0, natural_work,
1043 INSERT_VALUES, SCATTER_REVERSE);
1046 ierr = VecScatterEnd(scatter_to_zero, onp0, natural_work,
1047 INSERT_VALUES, SCATTER_REVERSE);
1050 ierr = DMDANaturalToGlobalBegin(*
dm(), natural_work, INSERT_VALUES, parallel);
1051 PISM_CHK(ierr,
"DMDANaturalToGlobalBegin");
1053 ierr = DMDANaturalToGlobalEnd(*
dm(), natural_work, INSERT_VALUES, parallel);
1054 PISM_CHK(ierr,
"DMDANaturalToGlobalEnd");
1083 MPI_Comm_rank(com, &rank);
1085 uint64_t result = 0;
1090 PetscErrorCode ierr = VecGetLocalSize(*v, &
size);
1095 MPI_Bcast(&result, 1, MPI_UINT64_T, 0, com);
1108 static_assert(
sizeof(
double) == 2 *
sizeof(uint32_t),
1109 "Cannot compile Array::fletcher64() (sizeof(double) != 2 * sizeof(uint32_t))");
1111 MPI_Status mpi_stat;
1112 const int checksum_tag = 42;
1117 MPI_Comm_rank(com, &rank);
1120 MPI_Comm_size(com, &comm_size);
1122 PetscInt local_size = 0;
1123 PetscErrorCode ierr = VecGetLocalSize(
vec(), &local_size);
PISM_CHK(ierr,
"VecGetLocalSize");
1132 std::vector<uint64_t> sums(comm_size);
1136 for (
int r = 1; r < comm_size; ++r) {
1137 MPI_Recv(&sums[r], 1, MPI_UINT64_T, r, checksum_tag, com, &mpi_stat);
1143 MPI_Send(&
sum, 1, MPI_UINT64_T, 0, checksum_tag, com);
1147 MPI_Bcast(&
sum, 1, MPI_UINT64_T, 0, com);
1164 log->message(1,
"%s%s: %s\n", prefix,
m_impl->
name.c_str(),
checksum(serial).c_str());
1168void Array::view(std::vector<std::shared_ptr<petsc::Viewer> > viewers)
const {
1169 PetscErrorCode ierr;
1173 "cannot 'view' a 3D field '%s'",
1183 for (
unsigned int i = 0; i <
ndof(); ++i) {
1190 output_units.c_str());
1194 PetscDraw draw = NULL;
1195 ierr = PetscViewerDrawGetDraw(v, 0, &draw);
1196 PISM_CHK(ierr,
"PetscViewerDrawGetDraw");
1198 ierr = PetscDrawSetTitle(draw, title.c_str());
1199 PISM_CHK(ierr,
"PetscDrawSetTitle");
1204 units, output_units);
1206 double bounds[2] = {0.0, 0.0};
1207 ierr = VecMin(tmp, NULL, &bounds[0]);
PISM_CHK(ierr,
"VecMin");
1208 ierr = VecMax(tmp, NULL, &bounds[1]);
PISM_CHK(ierr,
"VecMax");
1210 if (bounds[0] == bounds[1]) {
1215 ierr = PetscViewerDrawSetBounds(v, 1, bounds);
1216 PISM_CHK(ierr,
"PetscViewerDrawSetBounds");
1218 ierr = VecView(tmp, v);
1168void Array::view(std::vector<std::shared_ptr<petsc::Viewer> > viewers)
const {
…}
1226 const std::string &spec1,
const std::string &spec2) {
1230 PetscInt data_size = 0;
1231 PetscErrorCode ierr = VecGetLocalSize(v, &data_size);
VariableLookupData find_variable(const std::string &short_name, const std::string &std_name) const
Find a variable using its standard name and/or short name.
High-level PISM I/O class.
void message(int threshold, const char format[],...) const __attribute__((format(printf
Print a message to the log.
void failed()
Indicates a failure of a parallel section.
virtual void begin_access() const =0
void add_context(const std::string &message)
Add a message providing some context. This way we can (sort of) get a stack trace even though C++ exc...
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
void add(const PetscAccessible &v)
std::vector< const PetscAccessible * > m_vecs
void read(const std::string &filename, unsigned int time)
void write_impl(const File &file) const
Writes an Array to a NetCDF file.
virtual void end_access() const
Checks if an Array is allocated and calls DAVecRestoreArray.
virtual void regrid_impl(const File &file, io::Default default_value)
Gets an Array from a file file, interpolating onto the current grid.
Array(std::shared_ptr< const Grid > grid, const std::string &name, Kind ghostedp, size_t dof, size_t stencil_width, const std::vector< double > &zlevels)
unsigned int ndof() const
Returns the number of degrees of freedom per grid point.
void set_interpolation_type(InterpolationType type)
void set_name(const std::string &name)
Sets the variable name to name.
void dump(const char filename[]) const
Dumps a variable to a file, overwriting this file's contents (for debugging).
void define(const File &file, io::Type default_type) const
Define variables corresponding to an Array in a file opened using file.
virtual void begin_access() const
Checks if an Array is allocated and calls DAVecGetArray.
void shift(double alpha)
Result: v[j] <- v[j] + alpha for all j. Calls VecShift.
const std::vector< double > & levels() const
size_t size() const
Return the total number of elements in the owned part of an array.
void get_from_proc0(petsc::Vec &onp0)
Gets a local Array2 from processor 0.
void write(const std::string &filename) const
void copy_to_vec(std::shared_ptr< petsc::DM > destination_da, petsc::Vec &destination) const
Copies v to a global vector 'destination'. Ghost points are discarded.
const std::string & get_name() const
Get the name of an Array object.
uint64_t fletcher64_serial() const
void get_dof(std::shared_ptr< petsc::DM > da_result, petsc::Vec &result, unsigned int start, unsigned int count=1) const
void scale(double alpha)
Result: v <- v * alpha. Calls VecScale.
void print_checksum(const char *prefix="", bool serial=false) const
int state_counter() const
Get the object state counter.
std::vector< int > shape() const
std::string checksum(bool serial) const
void view(std::vector< std::shared_ptr< petsc::Viewer > > viewers) const
View a 2D vector field using existing PETSc viewers.
std::shared_ptr< const Grid > grid() const
void set(double c)
Result: v[j] <- c for all j.
void set_dof(std::shared_ptr< petsc::DM > da_source, petsc::Vec &source, unsigned int start, unsigned int count=1)
void inc_state_counter()
Increment the object state counter.
std::shared_ptr< petsc::DM > dm() const
unsigned int ndims() const
Returns the number of spatial dimensions.
void regrid(const std::string &filename, io::Default default_value)
std::shared_ptr< petsc::Vec > allocate_proc0_copy() const
uint64_t fletcher64() const
void put_on_proc0(petsc::Vec &onp0) const
Puts a local array::Scalar on processor 0.
void set_begin_access_use_dof(bool flag)
void add(double alpha, const Array &x)
Result: v <- v + alpha * x. Calls VecAXPY.
std::vector< double > norm(int n) const
Computes the norm of all the components of an Array.
void update_ghosts()
Updates ghost points.
unsigned int stencil_width() const
Get the stencil width of the current Array. Returns 0 if ghosts are not available.
void checkCompatibility(const char *function, const Array &other) const
Checks if two Arrays have compatible sizes, dimensions and numbers of degrees of freedom.
void read_impl(const File &file, unsigned int time)
Reads appropriate NetCDF variable(s) into an Array.
void check_array_indices(int i, int j, unsigned int k) const
Check array indices and warn if they are out of range.
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Abstract class for reading, writing, allocating, and accessing a DA-based PETSc Vec (2D and 3D fields...
Wrapper around VecGetArray and VecRestoreArray.
void convert_doubles(double *data, size_t length) const
std::shared_ptr< System > Ptr
#define PISM_CHK(errcode, name)
#define PISM_ERROR_LOCATION
static void check_range(petsc::Vec &v, const SpatialVariableMetadata &metadata, const std::string &filename, const Logger &log, bool report_range)
void global_to_local(petsc::DM &dm, Vec source, Vec destination)
double max(const array::Scalar &input)
Finds maximum over all the values in an array::Scalar object. Ignores ghosts.
void set_default_value_or_stop(const VariableMetadata &variable, io::Default default_value, const Logger &log, Vec output)
double sum(const array::Scalar &input)
Sums up all the values in an array::Scalar object. Ignores ghosts.
Kind
What "kind" of a vector to create: with or without ghosts.
static NormType int_to_normtype(int input)
double min(const array::Scalar &input)
Finds minimum over all the values in an array::Scalar object. Ignores ghosts.
void append_time(const File &file, const Config &config, double time_seconds)
Prepare a file for output.
void read_spatial_variable(const SpatialVariableMetadata &variable, const Grid &grid, const File &file, unsigned int time, double *output)
Read a variable from a file into an array output.
@ PISM_READWRITE_CLOBBER
create a file for writing, overwrite if present
@ PISM_READONLY
open an existing file for reading only
@ PISM_READWRITE
open an existing file for reading and writing
void write_spatial_variable(const SpatialVariableMetadata &metadata, const Grid &grid, const File &file, const double *input)
Write a double array to a file.
void define_spatial_variable(const SpatialVariableMetadata &metadata, const Grid &grid, const File &file, io::Type default_type)
Define a NetCDF variable corresponding to a VariableMetadata object.
void define_time(const File &file, const Context &ctx)
Prepare a file for output.
double get_time(MPI_Comm comm)
io::Backend string_to_backend(const std::string &backend)
void GlobalMax(MPI_Comm comm, double *local, double *result, int count)
std::string printf(const char *format,...)
uint64_t fletcher64(const uint32_t *data, size_t length)
static double end_time(const Config &config, double time_start, const std::string &calendar, const units::Unit &time_units)
static void report_range(const std::vector< double > &data, const VariableMetadata &metadata, const Logger &log)
Report the range of a time-series stored in data.
std::string timestamp(MPI_Comm com)
Creates a time-stamp used for the history NetCDF attribute.
void handle_fatal_errors(MPI_Comm com)
static double start_time(const Config &config, const File *file, const std::string &reference_date, const std::string &calendar, const units::Unit &time_units)
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)
void convert_vec(petsc::Vec &v, units::System::Ptr system, const std::string &spec1, const std::string &spec2)
bool begin_access_use_dof
If true, use DMDAVecGetArrayDOF() in begin_access()
int state_counter
Internal array::Array "revision number".
bool ghosted
true if this Array is ghosted
std::shared_ptr< const Grid > grid
The computational grid.
InterpolationType interpolation_type
bool report_range
If true, report range when regridding.
gsl_interp_accel * bsearch_accel
unsigned int da_stencil_width
stencil width supported by the DA
std::vector< SpatialVariableMetadata > metadata
Metadata (NetCDF variable attributes)
std::vector< double > zlevels
Vertical levels (for 3D fields)
std::shared_ptr< petsc::DM > da
unsigned int dof
number of "degrees of freedom" per grid point