diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index 12f0db620..ebecad3bb 100644 --- a/src/atlas/CMakeLists.txt +++ b/src/atlas/CMakeLists.txt @@ -474,12 +474,24 @@ field/FieldSet.cc field/FieldSet.h field/MissingValue.cc field/MissingValue.h +field/MultiField.cc +field/MultiField.h +field/MultiFieldCreator.cc +field/MultiFieldCreator.h +field/MultiFieldCreatorIFS.cc +field/MultiFieldCreatorIFS.h +field/MultiFieldCreatorArray.cc +field/MultiFieldCreatorArray.h field/State.cc field/State.h field/detail/FieldImpl.cc field/detail/FieldImpl.h field/detail/FieldInterface.cc field/detail/FieldInterface.h +field/detail/MultiFieldImpl.cc +field/detail/MultiFieldImpl.h +field/detail/MultiFieldInterface.cc +field/detail/MultiFieldInterface.h field/detail/MissingValue.cc field/detail/MissingValue.h ) diff --git a/src/atlas/field/MultiField.cc b/src/atlas/field/MultiField.cc new file mode 100644 index 000000000..e646f17af --- /dev/null +++ b/src/atlas/field/MultiField.cc @@ -0,0 +1,74 @@ +/* + * (C) Copyright 2013 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#include "atlas/field/MultiField.h" + +#include +#include +#include +#include +#include +#include + +#include "atlas/field/MultiFieldCreator.h" +#include "atlas/field/detail/MultiFieldImpl.h" +#include "atlas/runtime/Exception.h" + +namespace atlas { +namespace field { + +//----------------------------------------------------------------------------- + +MultiField::MultiField(const eckit::Configuration& config) { + std::string type; + if (!config.get("type", type)) { + ATLAS_THROW_EXCEPTION("Could not find \"type\" in configuration"); + } + std::unique_ptr creator(MultiFieldCreatorFactory::build(type, config)); + reset(creator->create(config)); +} + +MultiField::MultiField(const array::DataType datatype, const std::vector& shape, + const std::vector& var_names) { + std::unique_ptr creator(MultiFieldCreatorFactory::build("MultiFieldCreatorArray")); + reset(creator->create(datatype, shape, var_names)); +} + +const Field& MultiField::field(const std::string& name) const { return get()->field(name); } +Field& MultiField::field(const std::string& name) { return get()->field(name); } +bool MultiField::has(const std::string& name) const { return get()->has(name); } +std::vector MultiField::field_names() const { return get()->field_names(); } + +const Field& MultiField::field(const idx_t idx) const { return get()->field(idx); } +Field& MultiField::field(const idx_t idx) { return get()->field(idx); } +idx_t MultiField::size() const { return get()->size(); } + +const Field& MultiField::operator[](const idx_t idx) const { return get()->field(idx); } +Field& MultiField::operator[](const idx_t idx) { return get()->field(idx); } + +const Field& MultiField::operator[](const std::string& name) const { return get()->field(name); } +Field& MultiField::operator[](const std::string& name) { return get()->field(name); } + +const util::Metadata& MultiField::metadata() const { return get()->metadata(); } +util::Metadata& MultiField::metadata() { return get()->metadata(); } + +MultiField::operator const array::Array&() const { return get()->array(); } +MultiField::operator array::Array&() { return get()->array(); } + +MultiField::operator const FieldSet&() const { return get()->fieldset_; } +MultiField::operator FieldSet&() { return get()->fieldset_; } + +const array::Array& MultiField::array() const { return get()->array(); } +array::Array& MultiField::array() { return get()->array(); } + +//----------------------------------------------------------------------------- + +} // namespace field +} // namespace atlas diff --git a/src/atlas/field/MultiField.h b/src/atlas/field/MultiField.h new file mode 100644 index 000000000..2cc3c6903 --- /dev/null +++ b/src/atlas/field/MultiField.h @@ -0,0 +1,132 @@ +/* + * (C) Copyright 2013 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +/// @author Willem Deconinck +/// @date June 2015 + +#pragma once + +#include +#include + +#include "atlas/array/Array.h" +#include "atlas/field/Field.h" +#include "atlas/field/FieldSet.h" +#include "atlas/util/Config.h" +#include "atlas/util/Factory.h" +#include "atlas/util/Metadata.h" +#include "atlas/util/Object.h" +#include "atlas/util/ObjectHandle.h" + +namespace eckit { +class Parametrisation; +} + +namespace atlas { +namespace field { +class MultiFieldImpl; +} +} + +namespace atlas { +namespace field { + +/** + * \brief MultiField class that owns a collection of fields that are co-allocated + * + * Fields can only be described by parametrisation during the construction. + * Once setup, no additional fields can be added. + * + * Fields have to all be of same memory layout and data type + */ + +class MultiField : public util::ObjectHandle { +public: // methods + //-- Constructors + using Handle::Handle; + + MultiField(const eckit::Configuration&); + MultiField(const array::DataType datatype, const std::vector& shape, + const std::vector& var_names); + + //-- Accessors + + const Field& field(const std::string& name) const; + Field& field(const std::string& name); + bool has(const std::string& name) const; + std::vector field_names() const; + + const Field& field(const idx_t idx) const; + Field& field(const idx_t idx); + idx_t size() const; + + const Field& operator[](const idx_t idx) const; + Field& operator[](const idx_t idx); + + const Field& operator[](const std::string& name) const; + Field& operator[](const std::string& name); + + const util::Metadata& metadata() const; + util::Metadata& metadata(); + + // -- Modifiers + + /// @brief Implicit conversion to Array + operator const array::Array&() const; + operator array::Array&(); + + operator const FieldSet&() const; + operator FieldSet&(); + + /// @brief Access contained Array + const array::Array& array() const; + array::Array& array(); + +private: + template + void create(const std::vector shape, const std::vector var_names); +}; + +/** + * \brief MultiFieldArrayRegistry + */ + +class MultiFieldArrayRegistry : public field::FieldObserver { +private: + MultiFieldArrayRegistry() {} + +public: + static MultiFieldArrayRegistry& instance() { + static MultiFieldArrayRegistry inst; + return inst; + } + void onFieldDestruction(FieldImpl& field) override { + std::lock_guard guard(lock_); + map_.erase(&field); + } + + ~MultiFieldArrayRegistry() override = default; + + void add(Field& field, std::shared_ptr array) { + std::lock_guard guard(lock_); + map_.emplace(field.get(), array); + field->attachObserver(*this); + } + +public: + std::mutex lock_; + std::map> map_; + +}; + +// ------------------------------------------------------------------------------------ + +} // namespace field +} // namespace atlas diff --git a/src/atlas/field/MultiFieldCreator.cc b/src/atlas/field/MultiFieldCreator.cc new file mode 100644 index 000000000..72908f997 --- /dev/null +++ b/src/atlas/field/MultiFieldCreator.cc @@ -0,0 +1,58 @@ +/* + * (C) Copyright 2013 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +// file deepcode ignore CppMemoryLeak: static pointers for global registry are OK and will be cleaned up at end + +#include "atlas/field/MultiFieldCreator.h" + +#include +#include + +#include "eckit/thread/AutoLock.h" +#include "eckit/thread/Mutex.h" + +#include "atlas/field/MultiFieldCreatorIFS.h" +#include "atlas/field/MultiFieldCreatorArray.h" +#include "atlas/grid/Grid.h" +#include "atlas/runtime/Exception.h" +#include "atlas/runtime/Log.h" + + +namespace atlas { +namespace field { + + +namespace { + +void force_link() { + static struct Link { + Link() { + MultiFieldCreatorBuilder(); + MultiFieldCreatorBuilder(); + } + } link; +} + +} // namespace + +// ------------------------------------------------------------------ + +MultiFieldCreator::MultiFieldCreator() = default; + +MultiFieldCreator::~MultiFieldCreator() = default; + +MultiFieldCreator* MultiFieldCreatorFactory::build(const std::string& builder, const eckit::Configuration& config) { + force_link(); + auto factory = get(builder); + return factory->make(config); +} + +} // namespace field +} // namespace atlas diff --git a/src/atlas/field/MultiFieldCreator.h b/src/atlas/field/MultiFieldCreator.h new file mode 100644 index 000000000..84a0cf48e --- /dev/null +++ b/src/atlas/field/MultiFieldCreator.h @@ -0,0 +1,89 @@ +/* + * (C) Copyright 2013 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#ifndef atlas_field_MultiFieldCreator_h +#define atlas_field_MultiFieldCreator_h + +#include + +#include "atlas/util/Object.h" + +#include "atlas/field/MultiField.h" + +namespace eckit { +class Configuration; +} + +//------------------------------------------------------------------------------------------------------ + +namespace atlas { +namespace field { + +//------------------------------------------------------------------------------------------------------ + +/*! + * \brief Base class for creating new multifields based on Configuration + * + * \details + * Example to create field[100][3] of default type double: + * \code{.cpp} + * FieldImpl* field = Field::create( + * Config + * ("creator","ArraySpec") // ArraySpec MultiFieldCreator + * ("shape",array::make_shape(100,3)) // Rank 2 field with indexing [100][3] + * ); + * \endcode + */ +class MultiFieldCreator : public util::Object { +public: + MultiFieldCreator(); + MultiFieldCreator(const eckit::Configuration& config); + + virtual ~MultiFieldCreator(); + + virtual MultiFieldImpl* create(const eckit::Configuration& config = util::Config()) const = 0; + virtual MultiFieldImpl* create(const array::DataType datatype, const std::vector& shape, + const std::vector& var_names) const = 0; +}; + +//------------------------------------------------------------------------------------------------------ + +class MultiFieldCreatorFactory : public util::Factory { +public: + static std::string className() { return "MultiFieldCreatorFactory"; } + + /*! + * \brief build MultiFieldCreator with options specified in parametrisation + * \return mesh generator + */ + static MultiFieldCreator* build(const std::string&, const eckit::Configuration& = util::Config()); + + using Factory::Factory; + +private: + virtual MultiFieldCreator* make() = 0; + virtual MultiFieldCreator* make(const eckit::Configuration&) = 0; +}; + +template +class MultiFieldCreatorBuilder : public MultiFieldCreatorFactory { + virtual MultiFieldCreator* make() { return new T(); } + virtual MultiFieldCreator* make(const eckit::Configuration& config) { return new T(config); } + +public: + using MultiFieldCreatorFactory::MultiFieldCreatorFactory; +}; + +//------------------------------------------------------------------------------------------------------ + +} // namespace field +} // namespace atlas + +#endif diff --git a/src/atlas/field/MultiFieldCreatorArray.cc b/src/atlas/field/MultiFieldCreatorArray.cc new file mode 100644 index 000000000..5ecc5c2aa --- /dev/null +++ b/src/atlas/field/MultiFieldCreatorArray.cc @@ -0,0 +1,150 @@ +/* + * (C) Copyright 2013 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +/// @author Slavko Brdar +/// @author Willem Deconinck +/// @date September 2024 + +#include +#include + +#include "atlas/array/ArrayDataStore.h" +#include "atlas/array/DataType.h" +#include "atlas/array/Range.h" +#include "atlas/field/MultiFieldCreatorArray.h" +#include "atlas/field/detail/MultiFieldImpl.h" +#include "atlas/grid/Grid.h" +#include "atlas/runtime/Exception.h" +#include "atlas/runtime/Log.h" + +namespace atlas { +namespace field { + +MultiFieldCreatorArray::MultiFieldCreatorArray() {} + +MultiFieldCreatorArray::MultiFieldCreatorArray(const eckit::Configuration&) {} + +MultiFieldCreatorArray::~MultiFieldCreatorArray() = default; + +MultiFieldImpl* MultiFieldCreatorArray::create(const eckit::Configuration& config) const { + array::DataType datatype = array::DataType::create(); + std::string datatype_str; + if (config.get("datatype", datatype_str)) { + datatype = array::DataType(datatype_str); + } + else { + array::DataType::kind_t kind(array::DataType::kind()); + config.get("kind", kind); + if (!array::DataType::kind_valid(kind)) { + std::stringstream msg; + msg << "Could not create field. kind parameter unrecognized"; + throw_Exception(msg.str()); + } + datatype = array::DataType(kind); + } + std::vector shape; + config.get("shape", shape); + const auto fields = config.getSubConfigurations("fields"); + int nflds = 0; + for (size_t i = 0; i < fields.size(); ++i) { + long nvar = 1; + fields[i].get("nvar", nvar); + nflds += nvar; + } + std::vector var_names; + var_names.resize(nflds); + for (size_t i = 0, cnt = 0; i < fields.size(); ++i) { + std::string name; + fields[i].get("name", name); + long nvar = 1; + fields[i].get("nvar", nvar); + if (nvar > 1) { + for (int ivar = 0; ivar < nvar; ivar++) { + std::stringstream ss; + ss << name << "_" << ivar; + var_names[cnt++] = ss.str(); + } + + } + else { + var_names[cnt++] = name; + } + } + return create(datatype, shape, var_names); +} + +MultiFieldImpl* MultiFieldCreatorArray::create(const array::DataType datatype, const std::vector& shape, const std::vector& var_names) const { + const int dim = shape.size(); + const int nvar = var_names.size(); + ATLAS_ASSERT(nvar > 0 && dim > 2, "MultiField must have at least one field name."); + + int varidx = -1; + for (int i = 0; i < dim; i++) { + if (shape[i] == -1) { + varidx = i; + break; + } + } + + array::ArrayShape multiarray_shape = shape; + multiarray_shape[varidx] = nvar; + + MultiFieldImpl* multifield = new MultiFieldImpl{array::ArraySpec{datatype, multiarray_shape}}; + auto& multiarray = multifield->array(); + + array::ArrayShape field_shape; + field_shape.resize(dim - 1); + array::ArrayStrides field_strides; + field_strides.resize(dim - 1); + for (int i = 0, j = 0; i < dim; i++) { + if (i != varidx) { + field_shape[j] = multiarray.shape()[i]; + field_strides[j] = multiarray.strides()[i]; + ++j; + } + } + array::ArraySpec field_array_spec(field_shape, field_strides); + + for (size_t ifield = 0; ifield < nvar; ++ifield) { + idx_t start_index = multiarray.strides()[varidx] * ifield; + Field field; + if (datatype.kind() == array::DataType::KIND_REAL64) { + double* slice_begin = multiarray.data() + start_index; + field = Field(var_names[ifield], slice_begin, field_array_spec); + } + else if (datatype.kind() == array::DataType::KIND_REAL32) { + float* slice_begin = multiarray.data() + start_index; + field = Field(var_names[ifield], slice_begin, field_array_spec); + } + else if (datatype.kind() == array::DataType::KIND_INT32) { + int* slice_begin = multiarray.data() + start_index; + field = Field(var_names[ifield], slice_begin, field_array_spec); + } + else if (datatype.kind() == array::DataType::KIND_INT64) { + long* slice_begin = multiarray.data() + start_index; + field = Field(var_names[ifield], slice_begin, field_array_spec); + } + else { + ATLAS_NOTIMPLEMENTED; + } + multifield->add(field); + } + Log::debug() << "Creating multifield of " << datatype.str() << " type" << std::endl; + return multifield; +} + +// ------------------------------------------------------------------ + +namespace { +static MultiFieldCreatorBuilder __MultiFieldCreatorArray("MultiFieldCreatorArray"); +} + +} // namespace field +} // namespace atlas diff --git a/src/atlas/field/MultiFieldCreatorArray.h b/src/atlas/field/MultiFieldCreatorArray.h new file mode 100644 index 000000000..c66f2c144 --- /dev/null +++ b/src/atlas/field/MultiFieldCreatorArray.h @@ -0,0 +1,58 @@ +/* + * (C) Copyright 2013 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +/// @author Slavko Brdar +/// @date September 2024 + +#pragma once + +#include "atlas/field/MultiFieldCreator.h" + +namespace eckit { +class Configuration; +} +namespace atlas { +namespace field { +class MultiFieldImpl; +} +} // namespace atlas + +namespace atlas { +namespace field { + +// ------------------------------------------------------------------ + +/*! + * \brief MultiField creator using datatype, shape, variable names as arguments + * \details + * shape argument contains -1 at the position which gets filled with variable names + * Example use: + * \code{.cpp} + * MultiFieldImpl* multifield = MultiField::create( + * datatype, + * shape, + * var_names + * ); + * \endcode + */ +class MultiFieldCreatorArray : public MultiFieldCreator { +public: + MultiFieldCreatorArray(); + MultiFieldCreatorArray(const eckit::Configuration& config); + ~MultiFieldCreatorArray() override; + MultiFieldImpl* create(const eckit::Configuration& config = util::Config()) const override; + MultiFieldImpl* create(const array::DataType datatype, const std::vector& shape, + const std::vector& var_names) const override; +}; + +// ------------------------------------------------------------------ + +} // namespace field +} // namespace atlas diff --git a/src/atlas/field/MultiFieldCreatorIFS.cc b/src/atlas/field/MultiFieldCreatorIFS.cc new file mode 100644 index 000000000..c90681694 --- /dev/null +++ b/src/atlas/field/MultiFieldCreatorIFS.cc @@ -0,0 +1,230 @@ +/* + * (C) Copyright 2013 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +/// @author Willem Deconinck +/// @author Slavko Brdar +/// @date June 2024 + +#include +#include + +#include "eckit/config/Parametrisation.h" + +#include "atlas/array/ArrayDataStore.h" +#include "atlas/array/DataType.h" +#include "atlas/field/MultiFieldCreatorIFS.h" +#include "atlas/field/detail/MultiFieldImpl.h" +#include "atlas/grid/Grid.h" +#include "atlas/runtime/Exception.h" +#include "atlas/runtime/Log.h" + +namespace atlas { +namespace field { + +MultiFieldCreatorIFS::MultiFieldCreatorIFS() {} + +MultiFieldCreatorIFS::MultiFieldCreatorIFS(const eckit::Configuration& config) {} + +MultiFieldCreatorIFS::~MultiFieldCreatorIFS() = default; + +MultiFieldImpl* MultiFieldCreatorIFS::create(const array::DataType datatype, const std::vector& shape, + const std::vector& var_names) const { + ATLAS_NOTIMPLEMENTED; + return nullptr; +} + +MultiFieldImpl* MultiFieldCreatorIFS::create(const eckit::Configuration& config) const { + long nproma; + config.get("nproma", nproma); + long nlev; + config.get("nlev", nlev); + long nblk = 0; + if (config.has("nblk")) { + config.get("nblk", nblk); + } + else if (config.has("ngptot")) { + long ngptot; + config.get("ngptot", ngptot); + nblk = std::ceil(static_cast(ngptot) / static_cast(nproma)); + } + else { + ATLAS_THROW_EXCEPTION("Configuration not found: ngptot or nblk"); + } + array::DataType datatype = array::DataType::create(); + std::string datatype_str; + if (config.get("datatype", datatype_str)) { + datatype = array::DataType(datatype_str); + } + else { + array::DataType::kind_t kind(array::DataType::kind()); + config.get("kind", kind); + if (!array::DataType::kind_valid(kind)) { + std::stringstream msg; + msg << "Could not create field. kind parameter unrecognized"; + throw_Exception(msg.str()); + } + datatype = array::DataType(kind); + } + + auto fields = config.getSubConfigurations("fields"); + long nfld = 0; + for (const auto& field_params : fields) { + long nvar = 1; + field_params.get("nvar", nvar); + nfld += nvar; + } + array::ArrayShape multiarray_shape = ( (nlev > 0 and nfld > 0) ? array::make_shape(nblk, nfld, nlev, nproma) : + ( (nlev > 0) ? array::make_shape(nblk, nlev, nproma) : ( (nfld > 0) ? array::make_shape(nblk, nfld, nproma) : + array::make_shape(nblk, nproma) ) ) ); + + MultiFieldImpl* multifield = new MultiFieldImpl{array::ArraySpec{datatype, multiarray_shape}}; + auto& multiarray = multifield->array(); + + size_t multiarray_field_idx = 0; + for (size_t i = 0; i < fields.size(); ++i) { + std::string name; + fields[i].get("name", name); + Field field; + size_t field_vars = 1; + if (fields[i].get("nvar", field_vars)) { + array::ArrayShape field_shape = + ( (nlev > 0 and field_vars > 0) ? + array::make_shape(multiarray.shape(0), field_vars, multiarray.shape(2), multiarray.shape(3)) : + ( nlev > 0 ? + array::make_shape(multiarray.shape(0), multiarray.shape(1), multiarray.shape(2)) : + ( field_vars > 0 ? + array::make_shape(multiarray.shape(0), field_vars, multiarray.shape(2)) : + array::make_shape(multiarray.shape(0), multiarray.shape(1)) ) ) ); + array::ArrayShape multiarray_shape = + ( (nlev > 0 and field_vars > 0) ? array::make_shape(nblk, field_vars, nlev, nproma) : + ( (nlev > 0) ? array::make_shape(nblk, nlev, nproma) : ( (field_vars > 0) ? + array::make_shape(nblk, field_vars, nproma) : array::make_shape(nblk, nproma) ) ) ); + auto field_strides = multiarray.strides(); + auto field_array_spec = array::ArraySpec(field_shape, field_strides); + + constexpr auto all = array::Range::all(); + const auto range = array::Range(multiarray_field_idx, multiarray_field_idx + field_vars); + if (datatype.kind() == array::DataType::KIND_REAL64) { + if (nlev > 0 and field_vars > 0) { + auto slice = array::make_view(multiarray).slice(all, range, all, all); + field = Field(name, slice.data(), field_array_spec); + } + else if (nlev > 0) { + auto slice = array::make_view(multiarray).slice(all, all, all); + field = Field(name, slice.data(), field_array_spec); + } + else if (field_vars > 0) { + auto slice = array::make_view(multiarray).slice(all, range, all); + field = Field(name, slice.data(), field_array_spec); + } + else { + auto slice = array::make_view(multiarray).slice(all, all); + field = Field(name, slice.data(), field_array_spec); + } + } + else if (datatype.kind() == array::DataType::KIND_REAL32) { + if (nlev > 0 and field_vars > 0) { + auto slice = array::make_view(multiarray).slice(all, range, all, all); + field = Field(name, slice.data(), field_array_spec); + } + else if (nlev > 0) { + auto slice = array::make_view(multiarray).slice(all, all, all); + field = Field(name, slice.data(), field_array_spec); + } + else if (field_vars > 0) { + auto slice = array::make_view(multiarray).slice(all, range, all); + field = Field(name, slice.data(), field_array_spec); + } + else { + auto slice = array::make_view(multiarray).slice(all, all); + field = Field(name, slice.data(), field_array_spec); + } + } + else { + ATLAS_NOTIMPLEMENTED; + } + field.set_variables(field_vars); + } + else { + array::ArraySpec field_array_spec; + if (nlev > 0) { + auto field_shape = array::make_shape(multiarray.shape(0), multiarray.shape(2), multiarray.shape(3)); + auto field_strides = array::make_strides(multiarray.stride(0), multiarray.stride(2), multiarray.stride(3)); + field_array_spec = array::ArraySpec(field_shape, field_strides); + } + else if (field_vars > 0) { + auto field_shape = array::make_shape(multiarray.shape(0), multiarray.shape(2)); + auto field_strides = array::make_strides(multiarray.stride(0), multiarray.stride(2)); + field_array_spec = array::ArraySpec(field_shape, field_strides); + } + + constexpr auto all = array::Range::all(); + if (datatype.kind() == array::DataType::KIND_REAL64) { + if (nlev > 0 and field_vars > 0) { + auto slice = array::make_view(multiarray).slice(all, multiarray_field_idx, all, all); + field = Field(name, slice.data(), field_array_spec); + } + else if (nlev > 0) { + auto slice = array::make_view(multiarray).slice(all, all, all); + field = Field(name, slice.data(), field_array_spec); + } + else if (field_vars > 0) { + auto slice = array::make_view(multiarray).slice(all, multiarray_field_idx, all); + field = Field(name, slice.data(), field_array_spec); + } + else { + auto slice = array::make_view(multiarray).slice(all, all); + field = Field(name, slice.data(), field_array_spec); + } + } + else if (datatype.kind() == array::DataType::KIND_REAL32) { + if (nlev > 0 and field_vars > 0) { + auto slice = array::make_view(multiarray).slice(all, multiarray_field_idx, all, all); + field = Field(name, slice.data(), field_array_spec); + } + else if (nlev > 0) { + auto slice = array::make_view(multiarray).slice(all, all, all); + field = Field(name, slice.data(), field_array_spec); + } + else if (field_vars > 0) { + auto slice = array::make_view(multiarray).slice(all, multiarray_field_idx, all); + field = Field(name, slice.data(), field_array_spec); + } + else { + auto slice = array::make_view(multiarray).slice(all, all); + field = Field(name, slice.data(), field_array_spec); + } + } + else { + ATLAS_NOTIMPLEMENTED; + } + } + field.set_levels(nlev); + //field.set_blocks(nblk); + + multifield->add(field); + + multiarray_field_idx += field_vars; + } + std::string name; + config.get("name", name); + Log::debug() << "Creating IFS " << datatype.str() << " multifield: " << name << "[nblk=" << nblk << "][nvar=" << nfld + << "][nlev=" << nlev << "][nproma=" << nproma << "]\n"; + return multifield; +} + +// ------------------------------------------------------------------ + +namespace { +static MultiFieldCreatorBuilder __MultiFieldCreatorIFS("MultiFieldCreatorIFS"); +} + +} // namespace field +} // namespace atlas diff --git a/src/atlas/field/MultiFieldCreatorIFS.h b/src/atlas/field/MultiFieldCreatorIFS.h new file mode 100644 index 000000000..dddd7cd04 --- /dev/null +++ b/src/atlas/field/MultiFieldCreatorIFS.h @@ -0,0 +1,63 @@ +/* + * (C) Copyright 2013 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +/// @author Willem Deconinck +/// @date June 2015 + +#pragma once + +#include "atlas/field/MultiFieldCreator.h" + +namespace eckit { +class Configuration; +} +namespace atlas { +namespace field { +class MultiFieldImpl; +} +} // namespace atlas + +namespace atlas { +namespace field { + +// ------------------------------------------------------------------ + +/*! + * \brief MultiField creator using IFS parametrisation + * \details + * Ideally this class should belong to IFS. + * The only reference to IFS in Atlas::MultiField should be here. + * Example use: + * \code{.cpp} + * MultiFieldImpl* multifield = MultiField::create( + * Config + * ("creator","MultiFieldIFS") // MultiFieldIFS FieldCreator + * ("ngptot",ngptot) // Total number of grid points + * ("nproma",nproma) // Grouping of grid points for vectorlength + * ("nlev",nlev) // Number of levels + * ("nvar",nvar) // Number of variables + * ("kind",8) // Real kind in bytes + * ); + * \endcode + */ +class MultiFieldCreatorIFS : public MultiFieldCreator { +public: + MultiFieldCreatorIFS(); + MultiFieldCreatorIFS(const eckit::Configuration& config); + ~MultiFieldCreatorIFS() override; + MultiFieldImpl* create(const eckit::Configuration& config = util::Config()) const override; + MultiFieldImpl* create(const array::DataType datatype, const std::vector& shape, + const std::vector& var_names) const override; +}; + +// ------------------------------------------------------------------ + +} // namespace field +} // namespace atlas diff --git a/src/atlas/field/detail/MultiFieldImpl.cc b/src/atlas/field/detail/MultiFieldImpl.cc new file mode 100644 index 000000000..c3a6abb73 --- /dev/null +++ b/src/atlas/field/detail/MultiFieldImpl.cc @@ -0,0 +1,32 @@ +/* + * (C) Copyright 2013 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +// file deepcode ignore CppMemoryLeak: static pointers for global registry are OK and will be cleaned up at end + +#include "atlas/field/MultiField.h" +#include "atlas/field/MultiFieldCreator.h" + +#include +#include +#include + +#include "atlas/field/detail/MultiFieldImpl.h" + +namespace atlas { +namespace field { + +void MultiFieldImpl::add(Field& field) { + ATLAS_ASSERT(not fieldset_.has(field.name()), "Field with name \"" + field.name() + "\" already exists!"); + fieldset_.add(field); + MultiFieldArrayRegistry::instance().add(field, array_); +} + +} +} diff --git a/src/atlas/field/detail/MultiFieldImpl.h b/src/atlas/field/detail/MultiFieldImpl.h new file mode 100644 index 000000000..4a1445e1d --- /dev/null +++ b/src/atlas/field/detail/MultiFieldImpl.h @@ -0,0 +1,96 @@ +/* + * (C) Copyright 2013 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +/// @author Willem Deconinck +/// @date June 2015 + +#pragma once + +#include +#include + +#include "atlas/array/Array.h" +#include "atlas/field/Field.h" +#include "atlas/field/FieldSet.h" +#include "atlas/util/Config.h" +#include "atlas/util/Metadata.h" +#include "atlas/util/Object.h" + +namespace atlas { +namespace field { + +class MultiFieldImpl : public util::Object { +public: // methods + //-- Constructors + + MultiFieldImpl() {} + + MultiFieldImpl(const array::ArraySpec& spec) { + array::ArraySpec s(spec); + array_.reset(array::Array::create(std::move(s))); + } + + virtual ~MultiFieldImpl() {} + + + //-- Accessors + + const Field& field(const std::string& name) const { return fieldset_.field(name); } + Field& field(const std::string& name) { return fieldset_.field(name); } + bool has(const std::string& name) const { return fieldset_.has(name); } + std::vector field_names() const { return fieldset_.field_names(); } + + const Field& field(const idx_t idx) const { return fieldset_[idx]; } + Field& field(const idx_t idx) { return fieldset_[idx]; } + idx_t size() const { return fieldset_.size(); } + + const Field& operator[](const idx_t idx) const { return fieldset_[idx]; } + Field& operator[](const idx_t idx) { return fieldset_[idx]; } + + const Field& operator[](const std::string& name) const { return fieldset_.field(name); } + Field& operator[](const std::string& name) { return fieldset_.field(name); } + + const util::Metadata& metadata() const { return metadata_; } + util::Metadata& metadata() { return metadata_; } + + // -- Modifiers + + /// @brief Implicit conversion to Array + operator const array::Array&() const { return array(); } + operator array::Array&() { return array(); } + + operator const FieldSet&() const { return fieldset_; } + + operator FieldSet&() { return fieldset_; } + + /// @brief Access contained Array + const array::Array& array() const { + ATLAS_ASSERT(array_); + return *array_; + } + array::Array& array() { + ATLAS_ASSERT(array_); + return *array_; + } + + /// @brief Access contained FieldSet + const FieldSet& fieldset() const { return fieldset_; } + FieldSet& fieldset() { return fieldset_; } + + void add(Field& field); + +public: // temporary public for prototyping + FieldSet fieldset_; + std::shared_ptr array_; + util::Metadata metadata_; +}; + +} // namespace field +} // namespace atlas diff --git a/src/atlas/field/detail/MultiFieldInterface.cc b/src/atlas/field/detail/MultiFieldInterface.cc new file mode 100644 index 000000000..805f5e162 --- /dev/null +++ b/src/atlas/field/detail/MultiFieldInterface.cc @@ -0,0 +1,74 @@ +/* + * (C) Copyright 2013 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#include +#include +#include + +#include "atlas/library/config.h" +#include "atlas/field/MultiField.h" +#include "atlas/field/detail/MultiFieldImpl.h" +#include "atlas/field/detail/MultiFieldInterface.h" +#include "atlas/runtime/Exception.h" + +namespace atlas { +namespace field { + +extern "C" { +MultiFieldImpl* atlas__MultiField__create(eckit::Configuration* config) { + ATLAS_ASSERT(config != nullptr); + auto multifield = new MultiField(*config); + ATLAS_ASSERT(multifield); + return multifield->get(); +} + +MultiFieldImpl* atlas__MultiField__create_shape(int kind, int rank, int shapef[], const char* var_names, + size_t length, size_t size) { + array::ArrayShape shape; + shape.resize(rank); + array::ArrayStrides strides; + for (idx_t j = 0, jf = rank - 1; j < rank; ++j) { + shape[j] = shapef[jf--]; + } + + std::vector var_names_str; + for (size_t jj = 0; jj < size; ++jj) { + char str[length + 1]; + ATLAS_ASSERT(snprintf(str, sizeof(str), "%s", var_names + jj * length ) >= 0); + std::string sstr(str); + sstr.erase(std::find_if(sstr.rbegin(), sstr.rend(), [](unsigned char ch) { + return !std::isspace(ch); + }).base(), sstr.end()); + var_names_str.push_back(sstr); + } + + auto multifield = new MultiField(array::DataType{kind}, shape, var_names_str); + ATLAS_ASSERT(multifield); + return multifield->get(); +} + +void atlas__MultiField__delete(MultiFieldImpl* This) { + delete This; +} + +int atlas__MultiField__size(MultiFieldImpl* This) { + return This->size(); +} + +FieldSetImpl* atlas__MultiField__fieldset(MultiFieldImpl* This) { + return This->fieldset().get(); +} + +} + +// ------------------------------------------------------------------ + +} // namespace field +} // namespace atlas diff --git a/src/atlas/field/detail/MultiFieldInterface.h b/src/atlas/field/detail/MultiFieldInterface.h new file mode 100644 index 000000000..c63d4f755 --- /dev/null +++ b/src/atlas/field/detail/MultiFieldInterface.h @@ -0,0 +1,43 @@ +/* + * (C) Copyright 2013 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +/// @file MultiFieldInterface.h +/// @author Willem Deconinck +/// @date Sep 2014 + +#pragma once + +#include "atlas/field/FieldSet.h" +#include "atlas/field/MultiField.h" + +namespace atlas { +namespace functionspace { +} +} // namespace atlas + +namespace atlas { +namespace field { + +//---------------------------------------------------------------------------------------------------------------------- + +// C wrapper interfaces to C++ routines +extern "C" { +MultiFieldImpl* atlas__MultiField__create(eckit::Configuration* config); +MultiFieldImpl* atlas__MultiField__create_shape(int kind, int rank, int shapef[], const char* var_names, + size_t length, size_t size); +void atlas__MultiField__delete(MultiFieldImpl* This); +int atlas__MultiField__size(MultiFieldImpl* This); +FieldSetImpl* atlas__MultiField__fieldset(MultiFieldImpl* This); +} + +//---------------------------------------------------------------------------------------------------------------------- + +} // namespace field +} // namespace atlas diff --git a/src/atlas_f/CMakeLists.txt b/src/atlas_f/CMakeLists.txt index 60f55318f..2fd47e7d8 100644 --- a/src/atlas_f/CMakeLists.txt +++ b/src/atlas_f/CMakeLists.txt @@ -109,6 +109,9 @@ generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/field/State.h) generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/field/detail/FieldInterface.h MODULE atlas_field_c_binding OUTPUT field_c_binding.f90 ) +generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/field/detail/MultiFieldInterface.h + MODULE atlas_multifield_c_binding + OUTPUT multifield_c_binding.f90 ) generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/field/FieldSet.h) generate_fortran_bindings(FORTRAN_BINDINGS ../atlas/functionspace/detail/FunctionSpaceInterface.h MODULE atlas_functionspace_c_binding @@ -231,6 +234,7 @@ ecbuild_add_library( TARGET atlas_f field/atlas_FieldSet_module.fypp field/atlas_State_module.F90 field/atlas_Field_module.fypp + field/atlas_MultiField_module.F90 grid/atlas_Grid_module.F90 grid/atlas_GridDistribution_module.F90 grid/atlas_Vertical_module.F90 diff --git a/src/atlas_f/field/atlas_MultiField_module.F90 b/src/atlas_f/field/atlas_MultiField_module.F90 new file mode 100644 index 000000000..57168efaa --- /dev/null +++ b/src/atlas_f/field/atlas_MultiField_module.F90 @@ -0,0 +1,166 @@ +! (C) Copyright 2013 ECMWF. +! +! This software is licensed under the terms of the Apache Licence Version 2.0 +! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +! In applying this licence, ECMWF does not waive the privileges and immunities +! granted to it by virtue of its status as an intergovernmental organisation nor +! does it submit to any jurisdiction. + +#include "atlas/atlas_f.h" + +module atlas_multifield_module + +use fckit_owned_object_module, only : fckit_owned_object +use atlas_Config_module, only: atlas_Config +use atlas_field_module, only: atlas_field, array_c_to_f +use atlas_fieldset_module, only: atlas_fieldset + +implicit none + +public :: atlas_MultiField + +private + +!------------------------------------------------------------------------------ +TYPE, extends(fckit_owned_object) :: atlas_MultiField + +! Purpose : +! ------- +! *MultiField* : Object containing field data and Metadata + +! Methods : +! ------- +! name : The name or tag this field was created with +! data : Return the field as a fortran array of specified shape +! Metadata : Return object that can contain a variety of Metadata + +! Author : +! ------ +! 20-Nov-2013 Willem Deconinck *ECMWF* +! 29-Aug-2023 Slavko Brdar *ECMWF* + +!------------------------------------------------------------------------------ +contains + procedure, public :: MultiField__fieldset + procedure, public :: MultiField__size + generic :: fieldset => MultiField__fieldset + generic :: size => MultiField__size + +#if FCKIT_FINAL_NOT_INHERITING + final :: atlas_MultiField__final_auto +#endif + +END TYPE + +interface atlas_MultiField + module procedure atlas_MultiField__cptr + module procedure atlas_MultiField__create + module procedure atlas_MultiField__create_names +end interface + +private :: fckit_owned_object +private :: atlas_Config + +!======================================================== +contains +!======================================================== + +!------------------------------------------------------------------------------- + +function atlas_MultiField__cptr(cptr) result(field) + use, intrinsic :: iso_c_binding, only : c_ptr + type(atlas_MultiField) :: field + type(c_ptr), intent(in) :: cptr + call field%reset_c_ptr( cptr ) + call field%return() +end function + +!------------------------------------------------------------------------------- + +function atlas_MultiField__create(params) result(field) + use atlas_multifield_c_binding + type(atlas_MultiField) :: field + class(atlas_Config), intent(in) :: params + field = atlas_MultiField__cptr( atlas__MultiField__create(params%CPTR_PGIBUG_B) ) + call field%return() +end function + +!------------------------------------------------------------------------------- + +function atlas_MultiField__create_names(kind, shape, field_names) result(field) + use, intrinsic :: iso_c_binding, only : c_char, c_int, c_int32_t, c_size_t + use atlas_multifield_c_binding + type(atlas_MultiField) :: field + integer(c_int), intent(in) :: kind + integer, intent(in) :: shape(:) + character(*), intent(in) :: field_names(:) + character(kind=c_char,len=:), allocatable :: flat_field_names + integer(c_size_t) :: length + integer(c_int32_t) :: ii + integer(c_int32_t), allocatable :: field_name_sizes(:) + + if (size(field_names) == 0 .or. size(shape) < 3) then + print *, "atlas_MultiField__create_names: must have at least one field name, and the size of shape", & + & " is minimum 3, e.g. [nproma,-1,nblk]" + stop -1 + end if + + length = len(field_names(1)) + allocate(field_name_sizes(size(field_names))) + field_name_sizes = len(field_names(:)) + + if (any(field_name_sizes /= length)) then + print *, "atlas_MultiField__create_names: field_names have to have same length in characters" + stop -1 + end if + + allocate(character(len=length*size(field_names) ) :: flat_field_names) + do ii = 1, size(field_names) + flat_field_names((ii-1)*length+1:ii*length) = field_names(ii) + enddo + + field = atlas_MultiField__cptr( atlas__MultiField__create_shape(kind, size(shape), shape, & + & flat_field_names, length, size(field_names,kind=c_size_t)) ) + call field%return() +end function + +!------------------------------------------------------------------------------- + +function MultiField__size(this) result(size) + use atlas_multifield_c_binding + class(atlas_MultiField), intent(in) :: this + integer :: size + size = atlas__MultiField__size(this%CPTR_PGIBUG_B) +end function + +!------------------------------------------------------------------------------- +function MultiField__fieldset(this) result(fset) + use, intrinsic :: iso_c_binding, only : c_ptr + use atlas_multifield_c_binding + class(atlas_MultiField), intent(in) :: this + type(c_ptr) :: fset_cptr + type(atlas_FieldSet) :: fset + fset_cptr = atlas__MultiField__fieldset(this%CPTR_PGIBUG_B) + fset = atlas_FieldSet( fset_cptr ) + call fset%return() +end function + +!------------------------------------------------------------------------------- + +#if FCKIT_FINAL_NOT_INHERITING +ATLAS_FINAL subroutine atlas_MultiField__final_auto(this) + type(atlas_MultiField), intent(inout) :: this +#if FCKIT_FINAL_DEBUGGING + write(0,*) "atlas_MultiField__final_auto" +#endif +#if FCKIT_FINAL_NOT_PROPAGATING + call this%final() +#endif + FCKIT_SUPPRESS_UNUSED( this ) +end subroutine +#endif + +!------------------------------------------------------------------------------- + +end module atlas_multifield_module + diff --git a/src/tests/field/CMakeLists.txt b/src/tests/field/CMakeLists.txt index dc91d6f65..dbf51ee37 100644 --- a/src/tests/field/CMakeLists.txt +++ b/src/tests/field/CMakeLists.txt @@ -35,8 +35,21 @@ if( TEST atlas_test_field_acc ) set_tests_properties ( atlas_test_field_acc PROPERTIES LABELS "gpu;acc") endif() +ecbuild_add_test( TARGET atlas_test_multifield_ifs + SOURCES test_multifield_ifs.cc + LIBS atlas + ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} +) + if( HAVE_FCTEST ) + add_fctest( TARGET atlas_fctest_multifield_ifs + LINKER_LANGUAGE Fortran + SOURCES fctest_multifield_ifs.F90 + LIBS atlas_f + ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} + ) + add_fctest( TARGET atlas_fctest_field LINKER_LANGUAGE Fortran SOURCES fctest_field.F90 diff --git a/src/tests/field/fctest_multifield_ifs.F90 b/src/tests/field/fctest_multifield_ifs.F90 new file mode 100644 index 000000000..e71798c02 --- /dev/null +++ b/src/tests/field/fctest_multifield_ifs.F90 @@ -0,0 +1,228 @@ +! (C) Copyright 2013 ECMWF. +! This software is licensed under the terms of the Apache Licence Version 2.0 +! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +! In applying this licence, ECMWF does not waive the privileges and immunities +! granted to it by virtue of its status as an intergovernmental organisation nor +! does it submit to any jurisdiction. + +! This File contains Unit Tests for testing the +! C++ / Fortran Interfaces to the Mesh Datastructure +! @author Willem Deconinck +! @author Slavko Brdar + +#include "fckit/fctest.h" + +! ----------------------------------------------------------------------------- + +module fcta_MultiField_fixture +use atlas_module +use atlas_multifield_module +use, intrinsic :: iso_c_binding +implicit none +end module + +! ----------------------------------------------------------------------------- + +TESTSUITE_WITH_FIXTURE(fctest_atlas_MultiField, fcta_MultiField_fixture) + +! ----------------------------------------------------------------------------- + +TESTSUITE_INIT + call atlas_library%initialise() +END_TESTSUITE_INIT + +! ----------------------------------------------------------------------------- + +TESTSUITE_FINALIZE + call atlas_library%finalise() +END_TESTSUITE_FINALIZE + +! ----------------------------------------------------------------------------- + +TEST( test_multifield ) + implicit none + + type(atlas_MultiField) :: mfield(2) + type(atlas_FieldSet) :: fieldset(2) + type(atlas_Field) :: field + type(atlas_config) :: config + integer, pointer :: fdata_int_2d(:,:) + real(c_float), pointer :: fdata_f2d(:,:), fdata_f3d(:,:,:) + real(c_double), pointer :: fdata_d3d(:,:,:) + + integer, parameter :: nproma = 16; + integer, parameter :: nlev = 1; + integer, parameter :: ngptot = 2000; + integer, parameter :: nblk = (ngptot + nproma - 1) / nproma + type(atlas_Config), allocatable :: field_configs(:) + integer :: i + character(len=64), parameter, dimension(5) :: var_names = [ character(64) :: & + "temperature", "pressure", "density", "clv", "wind_u" ] + + config = atlas_Config() + call config%set("type", "MultiFieldCreatorIFS"); + call config%set("ngptot", ngptot); + call config%set("nproma", nproma); + allocate(field_configs(size(var_names))) + do i = 1, size(var_names) + field_configs(i) = atlas_Config() + call field_configs(i)%set("name", trim(var_names(i))) + end do + call field_configs(4)%set("nvar", 4) ! clv has four subvariables + call config%set("fields", field_configs) + + call config%set("nlev", 0); ! surface fields + call config%set("datatype", "real32"); + mfield(1) = atlas_MultiField(config) + + call config%set("nlev", 4); ! fields are 3d + call config%set("datatype", "real64"); + mfield(2) = atlas_MultiField(config) + + fieldset(1) = mfield(1)%fieldset() + FCTEST_CHECK_EQUAL(mfield(1)%size(), 5) + FCTEST_CHECK_EQUAL(fieldset(1)%size(), 5) + + fieldset(2) = atlas_FieldSet() + call fieldset(2)%add(mfield(1)%fieldset()) + field = fieldset(2)%field("density") + call field%data(fdata_f2d) + fdata_f2d(1,1) = 2. + call field%rename("dens") + + ! check data access directly though multifield + call fieldset(1)%data("dens", fdata_f2d) + fdata_f2d(1,1) = 3. + + ! check access to the renamed variable + field = fieldset(1)%field("dens") + call field%data(fdata_f2d) + FCTEST_CHECK_EQUAL(fdata_f2d(1,1), 3._c_float) + + ! check dimesionality + fieldset(2) = mfield(2)%fieldset() + call fieldset(2)%data("density", fdata_d3d) + fdata_d3d(1,1,1) = 4. + fieldset(2) = atlas_FieldSet() + call fieldset(2)%add(mfield(2)%fieldset()) + field = fieldset(2)%field("density") + call field%data(fdata_d3d) + FCTEST_CHECK_EQUAL(fdata_d3d(1,1,1), 4._c_double) +END_TEST + + +TEST( test_multifield_array_direct_constructor ) + implicit none + + type(atlas_MultiField) :: mfield(2) + type(atlas_FieldSet) :: fieldset(2), fset + type(atlas_Field) :: field + type(atlas_config) :: config + real(c_float), pointer :: fdata_f2d(:,:), fdata_f3d(:,:,:) + real(c_double), pointer :: fdata_d3d(:,:,:) + + integer, parameter :: nproma = 16; + integer, parameter :: nlev = 100; + integer, parameter :: ngptot = 2000; + integer, parameter :: nblk = (ngptot + nproma - 1) / nproma + integer :: i + character(len=64), parameter, dimension(5) :: var_names = [ character(64) :: & + "temperature ", "pressure", "density", "clv", "wind_u" ] + +return + + ! surface fields + mfield(1) = atlas_MultiField(atlas_real(c_float), [nproma, -1, nblk], var_names) + + ! 3d fields + mfield(2) = atlas_MultiField(atlas_real(c_double), [nproma, nlev, -1, nblk], var_names) + + FCTEST_CHECK_EQUAL(mfield(1)%size(), 5) + + fieldset(1) = mfield(1)%fieldset() + call fieldset(1)%data("density", fdata_f2d) + fdata_f2d(1,1) = 3. + fieldset(2) = mfield(2)%fieldset() + call fieldset(2)%data("density", fdata_d3d) + fdata_d3d(1,1,1) = 4. + + fset = atlas_FieldSet() + call fset%add(mfield(1)%fieldset()) + call fset%add(mfield(2)%fieldset()) + +END_TEST + + +TEST( test_multifield_array_config_constuctor ) + implicit none + + type(atlas_MultiField) :: mfield(2) + type(atlas_FieldSet) :: fieldset(2) + type(atlas_Field) :: field + type(atlas_config) :: config + integer, pointer :: fdata_int_2d(:,:) + real(c_float), pointer :: fdata_f2d(:,:), fdata_f3d(:,:,:) + real(c_double), pointer :: fdata_d3d(:,:,:) + + integer, parameter :: nproma = 16; + integer, parameter :: nlev = 1; + integer, parameter :: nblk = 200; + type(atlas_Config), allocatable :: field_configs(:) + integer :: i + character(len=64), parameter, dimension(5) :: var_names = [ character(64) :: & + "temperature", "pressure", "density", "clv", "wind_u" ] + + config = atlas_Config() + call config%set("type", "MultiFieldCreatorArray"); + allocate(field_configs(size(var_names))) + do i = 1, size(var_names) + field_configs(i) = atlas_Config() + call field_configs(i)%set("name", trim(var_names(i))) + end do + call field_configs(4)%set("nvar", 5) ! clv has four subvariables + call config%set("fields", field_configs) + + ! surface fields + call config%set("shape", [nproma, -1, nblk]); + call config%set("datatype", "real32"); + mfield(1) = atlas_MultiField(config) + + ! fields are 3d + call config%set("shape", [nproma, nlev, -1, nblk]); + call config%set("datatype", "real64"); + mfield(2) = atlas_MultiField(config) + + fieldset(1) = mfield(1)%fieldset() + FCTEST_CHECK_EQUAL(mfield(1)%size(), 9) + FCTEST_CHECK_EQUAL(fieldset(1)%size(), 9) + + fieldset(2) = atlas_FieldSet() + call fieldset(2)%add(mfield(1)%fieldset()) + field = fieldset(2)%field("density") + call field%data(fdata_f2d) + fdata_f2d(1,1) = 2. + call field%rename("dens") + + ! check data access directly though multifield + call fieldset(1)%data("dens", fdata_f2d) + fdata_f2d(1,1) = 3. + + ! check access to the renamed variable + field = fieldset(1)%field("dens") + call field%data(fdata_f2d) + FCTEST_CHECK_EQUAL(fdata_f2d(1,1), 3._c_float) + + ! check dimesionality + fieldset(2) = mfield(2)%fieldset() + call fieldset(2)%data("density", fdata_d3d) + fdata_d3d(1,1,1) = 4. + fieldset(2) = atlas_FieldSet() + call fieldset(2)%add(mfield(2)%fieldset()) + field = fieldset(2)%field("density") + call field%data(fdata_d3d) + FCTEST_CHECK_EQUAL(fdata_d3d(1,1,1), 4._c_double) +END_TEST + +! ----------------------------------------------------------------------------- + +END_TESTSUITE diff --git a/src/tests/field/test_multifield_ifs.cc b/src/tests/field/test_multifield_ifs.cc new file mode 100644 index 000000000..e4539bc28 --- /dev/null +++ b/src/tests/field/test_multifield_ifs.cc @@ -0,0 +1,413 @@ +/* + * (C) Copyright 2013 ECMWF. + * + * This software is licensed under the terms of the Apache Licence Version 2.0 + * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. + * In applying this licence, ECMWF does not waive the privileges and immunities + * granted to it by virtue of its status as an intergovernmental organisation + * nor does it submit to any jurisdiction. + */ + +#include +#include +#include + +#include "eckit/config/YAMLConfiguration.h" + +#include "atlas/array/ArrayView.h" +#include "atlas/array/MakeView.h" +#include "atlas/field/Field.h" +#include "atlas/field/MultiField.h" +#include "atlas/field/MultiFieldCreatorArray.h" +#include "atlas/field/detail/MultiFieldImpl.h" +#include "atlas/runtime/Exception.h" +#include "atlas/runtime/Log.h" + +#include "tests/AtlasTestEnvironment.h" + +using namespace atlas::field; + +namespace atlas { +namespace test { + +const std::vector make_shape(const std::initializer_list& list) { + return std::vector(list); +} + +CASE("multifield_generator") { + EXPECT(MultiFieldCreatorFactory::has("MultiFieldCreatorIFS")); + std::unique_ptr MultiFieldCreatorIFS(MultiFieldCreatorFactory::build("MultiFieldCreatorIFS")); + + EXPECT(MultiFieldCreatorFactory::has("MultiFieldCreatorArray")); + std::unique_ptr MultiFieldCreatorArray(MultiFieldCreatorFactory::build("MultiFieldCreatorArray")); +} + + +CASE("multifield_ifs_create") { + using Value = float; + int nproma = 16; + int nlev = 100; + int ngptot = 2000; + + auto json = [&]() -> std::string { + util::Config p; + p.set("type", "MultiFieldCreatorIFS"); + p.set("ngptot", ngptot); + p.set("nproma", nproma); + p.set("nlev", nlev); + p.set("datatype", array::make_datatype().str()); + p.set("fields", std::vector{ + util::Config("name", "temperature"), + util::Config("name", "pressure"), + util::Config("name", "density"), + util::Config("name", "clv")("nvar", 5), // 'clv' with 5 subvariables + util::Config("name", "wind_u") + }); + return p.json(); + }; + + SECTION("Print configuration") { + Log::info() << "json = " << json() << std::endl; + } + + SECTION("test") { + MultiField multifield{eckit::YAMLConfiguration{json()}}; + + const auto nblk = multifield.array().shape(0); + const auto nvar = multifield.array().shape(1); + const auto nfld = multifield.size(); + EXPECT_EQ(nfld, 5); + EXPECT_EQ(nvar, 9); + + EXPECT_EQ(multifield.size(), 5); + EXPECT(multifield.has("temperature")); + EXPECT(multifield.has("pressure")); + EXPECT(multifield.has("density")); + EXPECT(multifield.has("clv")); + EXPECT(multifield.has("wind_u")); + + Log::info() << multifield.field("temperature") << std::endl; + Log::info() << multifield.field("pressure") << std::endl; + Log::info() << multifield.field("density") << std::endl; + Log::info() << multifield.field("clv") << std::endl; + Log::info() << multifield.field("wind_u") << std::endl; + + auto temp = array::make_view(multifield.field("temperature")); + auto pres = array::make_view(multifield.field("pressure")); + auto dens = array::make_view(multifield.field("density")); + auto clv = array::make_view(multifield.field("clv")); // note rank 4 + auto wind_u = array::make_view(multifield.field("wind_u")); + + EXPECT_EQ(multifield[0].name(), "temperature"); + EXPECT_EQ(multifield[1].name(), "pressure"); + EXPECT_EQ(multifield[2].name(), "density"); + EXPECT_EQ(multifield[3].name(), "clv"); + EXPECT_EQ(multifield[4].name(), "wind_u"); + + auto block_stride = multifield.array().stride(0); + auto field_stride = nproma * nlev; + auto level_stride = nproma; + auto nproma_stride = 1; + + temp(1, 2, 3) = 4; + pres(5, 6, 7) = 8; + dens(9, 10, 11) = 12; + clv(13, 2, 14, 15) = 16; + wind_u(17, 18, 3) = 19; + + EXPECT_EQ(temp.stride(0), block_stride); + EXPECT_EQ(temp.stride(1), level_stride); + EXPECT_EQ(temp.stride(2), nproma_stride); + EXPECT_EQ(temp.size(), nblk * nlev * nproma); + + EXPECT_EQ(clv.stride(0), block_stride); + EXPECT_EQ(clv.stride(1), field_stride); + EXPECT_EQ(clv.stride(2), level_stride); + EXPECT_EQ(clv.stride(3), nproma_stride); + + EXPECT_EQ(clv.size(), nblk * 5 * nlev * nproma); + + + // Advanced usage, to access underlying array. This should only be used + // in a driver and not be exposed to algorithms. + { + auto multiarray = array::make_view(multifield); + EXPECT_EQ(multiarray.stride(0), block_stride); + EXPECT_EQ(multiarray.stride(1), field_stride); + EXPECT_EQ(multiarray.stride(2), level_stride); + EXPECT_EQ(multiarray.stride(3), nproma_stride); + + EXPECT_EQ(multiarray(1, 0, 2, 3), 4.); + EXPECT_EQ(multiarray(5, 1, 6, 7), 8.); + EXPECT_EQ(multiarray(9, 2, 10, 11), 12.); + EXPECT_EQ(multiarray(13, 5, 14, 15), 16.); + EXPECT_EQ(multiarray(17, 8, 18, 3), 19.); + + EXPECT_EQ(multiarray.size(), nblk * nvar * nlev * nproma); + } + + // access FieldSet through MultiField + auto fieldset = multifield->fieldset(); + auto field_v = array::make_view(fieldset.field("temperature")); + EXPECT_EQ(fieldset.size(), 5); + EXPECT(fieldset.has("temperature")); + EXPECT(fieldset.has("wind_u")); + EXPECT_EQ(field_v(1,2,3), 4); + } + + SECTION("test registry") { + { + Field field = MultiField {eckit::YAMLConfiguration{json()}}.field("temperature"); + auto temp = array::make_view(field); + } + } +} + +//----------------------------------------------------------------------------- + +CASE("multifield_array_create") { + using Value = float; + int nproma = 16; + int nlev = 100; + int ngptot = 2000; + + const int nblks = (ngptot + nproma - 1) / nproma; + const std::vector var_names = {"temperature", "pressure", "density", "clv", "wind_u"}; + + auto json = [&]() -> std::string { + util::Config p; + p.set("type", "MultiFieldCreatorArray"); + p.set("shape", {nblks, -1, nlev, nproma}); + p.set("datatype", array::make_datatype().str()); + p.set("fields", std::vector{ + util::Config("name", "temperature"), // + util::Config("name", "pressure"), // + util::Config("name", "density"), // + util::Config("name", "clv")("nvar", 5), // + util::Config("name", "wind_u") // + }); + return p.json(); + }; + + SECTION("test_MultiFieldArray_noconfig_3d") { + int nlev = 3; + const std::vector vshape = make_shape({nblks, -1, nlev, nproma}); + MultiField multifield(array::make_datatype(), vshape, var_names); + + const auto nblk = multifield.array().shape(0); + const auto nvar = multifield.array().shape(1); + nlev = multifield.array().shape(2); + const auto nfld = multifield.size(); + EXPECT_EQ(nfld, 5); + EXPECT_EQ(nvar, 5); + + EXPECT_EQ(multifield.size(), 5); + EXPECT(multifield.has("temperature")); + EXPECT(multifield.has("clv")); + + Log::info() << multifield.field("temperature") << std::endl; + Log::info() << multifield.field("clv") << std::endl; + + auto temp = array::make_view(multifield.field("temperature")); + auto clv = array::make_view(multifield.field("clv")); + EXPECT_EQ(multifield[0].name(), "temperature"); + EXPECT_EQ(multifield[3].name(), "clv"); + + auto block_stride = multifield.array().stride(0); + auto field_stride = nproma * nlev; + auto level_stride = nproma; + auto nproma_stride = 1; + + temp(1, 2, 3) = 4; + clv(13, 2, 14) = 16; + + EXPECT_EQ(temp.stride(0), block_stride); + EXPECT_EQ(temp.stride(1), level_stride); + EXPECT_EQ(temp.stride(2), nproma_stride); + EXPECT_EQ(temp.size(), nblk * nlev * nproma); + + // Advanced usage, to access underlying array. This should only be used + // in a driver and not be exposed to algorithms. + { + auto multiarray = array::make_view(multifield); + EXPECT_EQ(multiarray.stride(0), block_stride); + EXPECT_EQ(multiarray.stride(1), field_stride); + EXPECT_EQ(multiarray.stride(2), level_stride); + EXPECT_EQ(multiarray.stride(3), nproma_stride); + + EXPECT_EQ(multiarray(1, 0, 2, 3), 4.); + EXPECT_EQ(multiarray(13, 3, 2, 14), 16.); + + EXPECT_EQ(multiarray.size(), nblk * nvar * nlev * nproma); + } + + // access FieldSet through MultiField + auto fieldset = multifield->fieldset(); + auto field_v = array::make_view(fieldset.field("temperature")); + EXPECT_EQ(fieldset.size(), 5); + EXPECT(fieldset.has("temperature")); + EXPECT(fieldset.has("wind_u")); + EXPECT_EQ(field_v(1,2,3), 4); + } + + SECTION("test_MultiFieldArray_noconfig_2d") { + const std::vector vshape = make_shape({nblks, -1, nproma}); + MultiField multifield(array::make_datatype(), vshape, var_names); + + const auto nblk = multifield.array().shape(0); + const auto nvar = multifield.array().shape(1); + nlev = multifield.array().shape(2); + const auto nfld = multifield.size(); + EXPECT_EQ(nfld, 5); + EXPECT_EQ(nvar, 5); + + EXPECT_EQ(multifield.size(), 5); + EXPECT(multifield.has("temperature")); + EXPECT(multifield.has("clv")); + + Log::info() << multifield.field("temperature") << std::endl; + Log::info() << multifield.field("clv") << std::endl; + + auto temp = array::make_view(multifield.field("temperature")); + auto clv = array::make_view(multifield.field("clv")); + EXPECT_EQ(multifield[0].name(), "temperature"); + EXPECT_EQ(multifield[3].name(), "clv"); + + auto block_stride = multifield.array().stride(0); + auto field_stride = nproma; + auto nproma_stride = 1; + + temp(1, 3) = 4; + clv(13, 14) = 16; + + EXPECT_EQ(temp.stride(0), block_stride); + EXPECT_EQ(temp.stride(1), nproma_stride); + EXPECT_EQ(temp.size(), nblk * nproma); + + // Advanced usage, to access underlying array. This should only be used + // in a driver and not be exposed to algorithms. + { + auto multiarray = array::make_view(multifield); + EXPECT_EQ(multiarray.stride(0), block_stride); + EXPECT_EQ(multiarray.stride(1), field_stride); + EXPECT_EQ(multiarray.stride(2), nproma_stride); + + EXPECT_EQ(multiarray(1, 0, 3), 4.); + EXPECT_EQ(multiarray(13, 3, 14), 16.); + + EXPECT_EQ(multiarray.size(), nblk * nvar * nproma); + } + + // access FieldSet through MultiField + auto fieldset = multifield->fieldset(); + auto field_v = array::make_view(fieldset.field("temperature")); + EXPECT_EQ(fieldset.size(), 5); + EXPECT(fieldset.has("temperature")); + EXPECT(fieldset.has("wind_u")); + EXPECT_EQ(field_v(1,3), 4); + } + + SECTION("Print configuration") { + Log::info() << "json = " << json() << std::endl; + } + + SECTION("test_MultiFieldArray_config") { + MultiField multifield{eckit::YAMLConfiguration{json()}}; + + const auto nblk = multifield.array().shape(0); + const auto nvar = multifield.array().shape(1); + const auto nfld = multifield.size(); + EXPECT_EQ(nfld, 9); + EXPECT_EQ(nvar, 9); + + EXPECT_EQ(multifield.size(), 9); + EXPECT(multifield.has("temperature")); + EXPECT(multifield.has("pressure")); + EXPECT(multifield.has("density")); + EXPECT(multifield.has("clv_0")); + EXPECT(multifield.has("wind_u")); + + Log::info() << multifield.field("temperature") << std::endl; + Log::info() << multifield.field("pressure") << std::endl; + Log::info() << multifield.field("density") << std::endl; + Log::info() << multifield.field("clv_0") << std::endl; + Log::info() << multifield.field("wind_u") << std::endl; + + auto temp = array::make_view(multifield.field("temperature")); + auto pres = array::make_view(multifield.field("pressure")); + auto dens = array::make_view(multifield.field("density")); + auto clv = array::make_view(multifield.field("clv_0")); // note rank 4 + auto wind_u = array::make_view(multifield.field("wind_u")); + + EXPECT_EQ(multifield[0].name(), "temperature"); + EXPECT_EQ(multifield[1].name(), "pressure"); + EXPECT_EQ(multifield[2].name(), "density"); + EXPECT_EQ(multifield[3].name(), "clv_0"); + EXPECT_EQ(multifield[8].name(), "wind_u"); + + auto block_stride = multifield.array().stride(0); + auto field_stride = nproma * nlev; + auto level_stride = nproma; + auto nproma_stride = 1; + + temp(1, 2, 3) = 4; + pres(5, 6, 7) = 8; + dens(9, 10, 11) = 12; + clv(13, 14, 15) = 16; + wind_u(17, 18, 3) = 19; + + EXPECT_EQ(temp.stride(0), block_stride); + EXPECT_EQ(temp.stride(1), level_stride); + EXPECT_EQ(temp.stride(2), nproma_stride); + EXPECT_EQ(temp.size(), nblk * nlev * nproma); + + EXPECT_EQ(clv.stride(0), block_stride); + EXPECT_EQ(clv.stride(1), level_stride); + EXPECT_EQ(clv.stride(2), nproma_stride); + + EXPECT_EQ(clv.size(), nblk * nlev * nproma); + + + // Advanced usage, to access underlying array. This should only be used + // in a driver and not be exposed to algorithms. + { + auto multiarray = array::make_view(multifield); + EXPECT_EQ(multiarray.stride(0), block_stride); + EXPECT_EQ(multiarray.stride(1), field_stride); + EXPECT_EQ(multiarray.stride(2), level_stride); + EXPECT_EQ(multiarray.stride(3), nproma_stride); + + EXPECT_EQ(multiarray(1, 0, 2, 3), 4.); + EXPECT_EQ(multiarray(5, 1, 6, 7), 8.); + EXPECT_EQ(multiarray(9, 2, 10, 11), 12.); + EXPECT_EQ(multiarray(13, 3, 14, 15), 16.); + EXPECT_EQ(multiarray(17, 8, 18, 3), 19.); + + EXPECT_EQ(multiarray.size(), nblk * nvar * nlev * nproma); + } + + // access FieldSet through MultiField + auto fieldset = multifield->fieldset(); + auto field_v = array::make_view(fieldset.field("temperature")); + EXPECT_EQ(fieldset.size(), 9); + EXPECT(fieldset.has("temperature")); + EXPECT(fieldset.has("wind_u")); + EXPECT_EQ(field_v(1,2,3), 4); + } + + SECTION("test registry") { + { + Field field = MultiField {eckit::YAMLConfiguration{json()}}.field("temperature"); + auto temp = array::make_view(field); + } + } +} + +//----------------------------------------------------------------------------- + +} // namespace test +} // namespace atlas + +int main(int argc, char** argv) { + return atlas::test::run(argc, argv); +}