diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index 58fed83e8..ebecad3bb 100644 --- a/src/atlas/CMakeLists.txt +++ b/src/atlas/CMakeLists.txt @@ -476,12 +476,20 @@ 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 diff --git a/src/atlas/field/MultiField.cc b/src/atlas/field/MultiField.cc index f7a3ba88d..e646f17af 100644 --- a/src/atlas/field/MultiField.cc +++ b/src/atlas/field/MultiField.cc @@ -8,75 +8,24 @@ * nor does it submit to any jurisdiction. */ -#include "MultiField.h" +#include "atlas/field/MultiField.h" #include #include #include +#include #include #include -#include "atlas/field/Field.h" -#include "atlas/grid/Grid.h" -#include "atlas/mesh/Mesh.h" +#include "atlas/field/MultiFieldCreator.h" +#include "atlas/field/detail/MultiFieldImpl.h" #include "atlas/runtime/Exception.h" -#include "atlas/runtime/Log.h" namespace atlas { namespace field { -namespace { -void force_link() { - static struct Link { - Link() { - ; - // For static linking add here something like - // MultiFieldCreatorBuilder(); - } - } link; -} -} // namespace - //----------------------------------------------------------------------------- -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_; - -}; - -MultiFieldCreator::MultiFieldCreator(const eckit::Configuration&) {} - -MultiFieldCreator::~MultiFieldCreator() = default; - -MultiFieldCreator* MultiFieldCreatorFactory::build(const std::string& builder, const eckit::Configuration& config) { - force_link(); - auto factory = get(builder); - return factory->make(config); -} - MultiField::MultiField(const eckit::Configuration& config) { std::string type; if (!config.get("type", type)) { @@ -86,13 +35,39 @@ MultiField::MultiField(const eckit::Configuration& config) { reset(creator->create(config)); } -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_); - +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 diff --git a/src/atlas/field/MultiField.h b/src/atlas/field/MultiField.h index 46d0adea0..2cc3c6903 100644 --- a/src/atlas/field/MultiField.h +++ b/src/atlas/field/MultiField.h @@ -29,6 +29,12 @@ namespace eckit { class Parametrisation; } +namespace atlas { +namespace field { +class MultiFieldImpl; +} +} + namespace atlas { namespace field { @@ -41,150 +47,83 @@ namespace field { * Fields have to all be of same memory layout and data type */ -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_; -}; - - 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 { return get()->field(name); } - Field& field(const std::string& name) { return get()->field(name); } - bool has(const std::string& name) const { return get()->has(name); } - std::vector field_names() const { return get()->field_names(); } + 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 { return get()->field(idx); } - Field& field(const idx_t idx) { return get()->field(idx); } - idx_t size() const { return get()->size(); } + 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 { return get()->field(idx); } - Field& operator[](const idx_t idx) { return get()->field(idx); } + const Field& operator[](const idx_t idx) const; + Field& operator[](const idx_t idx); - const Field& operator[](const std::string& name) const { return get()->field(name); } - Field& operator[](const std::string& name) { return get()->field(name); } + const Field& operator[](const std::string& name) const; + Field& operator[](const std::string& name); - const util::Metadata& metadata() const { return get()->metadata(); } - util::Metadata& metadata() { return get()->metadata(); } + const util::Metadata& metadata() const; + util::Metadata& metadata(); // -- Modifiers /// @brief Implicit conversion to Array - operator const array::Array&() const { return get()->array(); } - operator array::Array&() { return get()->array(); } + operator const array::Array&() const; + operator array::Array&(); - operator const FieldSet&() const { return get()->fieldset_; } - operator FieldSet&() { return get()->fieldset_; } + operator const FieldSet&() const; + operator FieldSet&(); /// @brief Access contained Array - const array::Array& array() const { return get()->array(); } - array::Array& array() { return get()->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() {} -class MultiFieldCreator : public util::Object { public: - MultiFieldCreator(const eckit::Configuration& = util::Config()); - - virtual ~MultiFieldCreator(); + static MultiFieldArrayRegistry& instance() { + static MultiFieldArrayRegistry inst; + return inst; + } + void onFieldDestruction(FieldImpl& field) override { + std::lock_guard guard(lock_); + map_.erase(&field); + } - virtual MultiFieldImpl* create(const eckit::Configuration& = util::Config()) const = 0; -}; + ~MultiFieldArrayRegistry() override = default; -//------------------------------------------------------------------------------------------------------ + void add(Field& field, std::shared_ptr array) { + std::lock_guard guard(lock_); + map_.emplace(field.get(), array); + field->attachObserver(*this); + } -class MultiFieldCreatorFactory : public util::Factory { public: - static std::string className() { return "MultiFieldCreatorFactory"; } - - /*! - * \brief build MultiFieldCreator with options specified in parametrisation - * \return MutliField creator - */ - static MultiFieldCreator* build(const std::string&, const eckit::Configuration& = util::NoConfig()); + std::mutex lock_; + std::map> map_; - using Factory::Factory; - -private: - virtual MultiFieldCreator* make(const eckit::Configuration&) = 0; -}; - -template -class MultiFieldCreatorBuilder : public MultiFieldCreatorFactory { - virtual MultiFieldCreator* make(const eckit::Configuration& config) override { return new T(config); } - -public: - using MultiFieldCreatorFactory::MultiFieldCreatorFactory; }; // ------------------------------------------------------------------------------------ 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 index a5e729b2a..805f5e162 100644 --- a/src/atlas/field/detail/MultiFieldInterface.cc +++ b/src/atlas/field/detail/MultiFieldInterface.cc @@ -8,10 +8,13 @@ * 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" @@ -21,108 +24,48 @@ namespace field { extern "C" { MultiFieldImpl* atlas__MultiField__create(eckit::Configuration* config) { ATLAS_ASSERT(config != nullptr); - long nproma = config->getLong("nproma"); - long nlev = config->getLong("nlev"); - long nblk = 0; - if (config->has("nblk")) { - nblk = config->getLong("nblk"); - } - else if (config->has("ngptot")) { - long ngptot = config->getLong("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 multifield = new MultiField(*config); + ATLAS_ASSERT(multifield); + return multifield->get(); +} - auto fields = config->getSubConfigurations("fields"); - long nfld = 0; - for (const auto& field_config : fields) { - long nvar = 1; - field_config.get("nvar", nvar); - nfld += nvar; +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--]; } - auto multiarray_shape = array::make_shape(nblk, nfld, nlev, 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)) { - auto field_shape = - array::make_shape(multiarray.shape(0), field_vars, multiarray.shape(2), multiarray.shape(3)); - 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) { - auto slice = array::make_view(multiarray).slice(all, range, all, all); - field = Field(name, slice.data(), field_array_spec); - } - else if (datatype.kind() == array::DataType::KIND_REAL32) { - auto slice = array::make_view(multiarray).slice(all, range, all, all); - field = Field(name, slice.data(), field_array_spec); - } - else { - ATLAS_NOTIMPLEMENTED; - } - field.set_variables(field_vars); - } - else { - 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)); - auto field_array_spec = array::ArraySpec(field_shape, field_strides); - - constexpr auto all = array::Range::all(); - if (datatype.kind() == array::DataType::KIND_REAL64) { - auto slice = array::make_view(multiarray).slice(all, multiarray_field_idx, all, all); - field = Field(name, slice.data(), field_array_spec); - } - else if (datatype.kind() == array::DataType::KIND_REAL32) { - auto slice = array::make_view(multiarray).slice(all, multiarray_field_idx, 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::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; + 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(); +} + } // ------------------------------------------------------------------ diff --git a/src/atlas/field/detail/MultiFieldInterface.h b/src/atlas/field/detail/MultiFieldInterface.h index ac3b1ae22..c63d4f755 100644 --- a/src/atlas/field/detail/MultiFieldInterface.h +++ b/src/atlas/field/detail/MultiFieldInterface.h @@ -14,6 +14,7 @@ #pragma once +#include "atlas/field/FieldSet.h" #include "atlas/field/MultiField.h" namespace atlas { @@ -29,7 +30,11 @@ 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); } //---------------------------------------------------------------------------------------------------------------------- diff --git a/src/atlas_f/CMakeLists.txt b/src/atlas_f/CMakeLists.txt index eefbbf647..2fd47e7d8 100644 --- a/src/atlas_f/CMakeLists.txt +++ b/src/atlas_f/CMakeLists.txt @@ -234,7 +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.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.fypp b/src/atlas_f/field/atlas_MultiField_module.F90 similarity index 52% rename from src/atlas_f/field/atlas_MultiField_module.fypp rename to src/atlas_f/field/atlas_MultiField_module.F90 index 4d45836a5..57168efaa 100644 --- a/src/atlas_f/field/atlas_MultiField_module.fypp +++ b/src/atlas_f/field/atlas_MultiField_module.F90 @@ -8,13 +8,12 @@ #include "atlas/atlas_f.h" -#:include "atlas/atlas_f.fypp" -#:include "internals/atlas_generics.fypp" - 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 @@ -42,6 +41,10 @@ module atlas_multifield_module !------------------------------------------------------------------------------ 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 @@ -52,6 +55,7 @@ module atlas_multifield_module 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 @@ -83,6 +87,66 @@ function atlas_MultiField__create(params) result(field) !------------------------------------------------------------------------------- +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 diff --git a/src/tests/field/fctest_multifield_ifs.F90 b/src/tests/field/fctest_multifield_ifs.F90 index c6ca36f47..e71798c02 100644 --- a/src/tests/field/fctest_multifield_ifs.F90 +++ b/src/tests/field/fctest_multifield_ifs.F90 @@ -23,7 +23,7 @@ module fcta_MultiField_fixture ! ----------------------------------------------------------------------------- -TESTSUITE_WITH_FIXTURE(fctest_atlas_MultiField,fcta_MultiField_fixture) +TESTSUITE_WITH_FIXTURE(fctest_atlas_MultiField, fcta_MultiField_fixture) ! ----------------------------------------------------------------------------- @@ -42,34 +42,185 @@ module fcta_MultiField_fixture TEST( test_multifield ) implicit none - type(atlas_MultiField) :: mfield - type(atlas_config) :: config + 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 = 100; + integer, parameter :: nlev = 1; integer, parameter :: ngptot = 2000; - type(atlas_Config), dimension(5) :: field_configs + 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); - call config%set("nlev", nlev); + 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"); - field_configs(1) = atlas_Config() - field_configs(2) = atlas_Config() - field_configs(3) = atlas_Config() - field_configs(4) = atlas_Config() - field_configs(5) = atlas_Config() - call field_configs(1)%set("name", "temperature") - call field_configs(2)%set("name", "pressure") - call field_configs(3)%set("name", "density") - call field_configs(4)%set("name", "clv") - call field_configs(4)%set("nvar", 5) - call field_configs(5)%set("name", "wind_u") + 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) - mfield = atlas_MultiField(config) + ! 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 ! ----------------------------------------------------------------------------- diff --git a/src/tests/field/test_multifield_ifs.cc b/src/tests/field/test_multifield_ifs.cc index ce11afa35..e4539bc28 100644 --- a/src/tests/field/test_multifield_ifs.cc +++ b/src/tests/field/test_multifield_ifs.cc @@ -10,147 +10,40 @@ #include #include +#include #include "eckit/config/YAMLConfiguration.h" #include "atlas/array/ArrayView.h" -#include "atlas/array/DataType.h" #include "atlas/array/MakeView.h" #include "atlas/field/Field.h" #include "atlas/field/MultiField.h" -#include "atlas/grid/Grid.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; using namespace atlas::field; namespace atlas { namespace test { -// ------------------------------------------------------------------- -// Example IFS MultiField creato - -// --- Declaration (in .h file) -class MultiFieldCreatorIFS : public MultiFieldCreator { -public: - MultiFieldCreatorIFS(const eckit::Configuration& config = util::Config()): MultiFieldCreator(config) {} - ~MultiFieldCreatorIFS() override = default; - MultiFieldImpl* create(const eckit::Configuration& = util::Config()) const override; -}; - -// --- Implementation (in .cc file) -MultiFieldImpl* MultiFieldCreatorIFS::create(const eckit::Configuration& config) const { - long ngptot = config.getLong("ngptot"); - long nproma = config.getLong("nproma"); - long nlev = config.getLong("nlev"); - long nblk = 0; - - - 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); - } - - nblk = std::ceil(static_cast(ngptot) / static_cast(nproma)); - - auto fields = config.getSubConfigurations("fields"); - long nfld = 0; - for (const auto& field_config : fields) { - long nvar = 1; - field_config.get("nvar", nvar); - nfld += nvar; - } - - auto multiarray_shape = array::make_shape(nblk, nfld, nlev, 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)) { - auto field_shape = - array::make_shape(multiarray.shape(0), field_vars, multiarray.shape(2), multiarray.shape(3)); - 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) { - auto slice = array::make_view(multiarray).slice(all, range, all, all); - field = Field(name, slice.data(), field_array_spec); - } - else if (datatype.kind() == array::DataType::KIND_REAL32) { - auto slice = array::make_view(multiarray).slice(all, range, all, all); - field = Field(name, slice.data(), field_array_spec); - } - else { - ATLAS_NOTIMPLEMENTED; - } - field.set_variables(field_vars); - } - else { - 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)); - auto field_array_spec = array::ArraySpec(field_shape, field_strides); - - constexpr auto all = array::Range::all(); - if (datatype.kind() == array::DataType::KIND_REAL64) { - auto slice = array::make_view(multiarray).slice(all, multiarray_field_idx, all, all); - field = Field(name, slice.data(), field_array_spec); - } - else if (datatype.kind() == array::DataType::KIND_REAL32) { - auto slice = array::make_view(multiarray).slice(all, multiarray_field_idx, all, all); - field = Field(name, slice.data(), field_array_spec); - } - else { - ATLAS_NOTIMPLEMENTED; - } - } - field.set_levels(nlev); - - multifield->add(field); - - multiarray_field_idx += field_vars; - } - return multifield; +const std::vector make_shape(const std::initializer_list& list) { + return std::vector(list); } -// Register in factory -MultiFieldCreatorBuilder __MultiFieldCreatorIFS("MultiFieldCreatorIFS"); - -// =================================================================== -// BEGIN TESTS -// =================================================================== - - CASE("multifield_generator") { EXPECT(MultiFieldCreatorFactory::has("MultiFieldCreatorIFS")); - std::unique_ptr MultiFieldCreator(MultiFieldCreatorFactory::build("MultiFieldCreatorIFS")); + std::unique_ptr MultiFieldCreatorIFS(MultiFieldCreatorFactory::build("MultiFieldCreatorIFS")); + + EXPECT(MultiFieldCreatorFactory::has("MultiFieldCreatorArray")); + std::unique_ptr MultiFieldCreatorArray(MultiFieldCreatorFactory::build("MultiFieldCreatorArray")); } -CASE("multifield_create") { +CASE("multifield_ifs_create") { using Value = float; int nproma = 16; int nlev = 100; @@ -163,12 +56,12 @@ CASE("multifield_create") { p.set("nproma", nproma); p.set("nlev", nlev); p.set("datatype", array::make_datatype().str()); - p.set("fields", { - util::Config("name", "temperature"), // - util::Config("name", "pressure"), // - util::Config("name", "density"), // - util::Config("name", "clv")("nvar", 5), // - util::Config("name", "wind_u") // + 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(); }; @@ -252,6 +145,14 @@ CASE("multifield_create") { 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") { @@ -264,6 +165,246 @@ CASE("multifield_create") { //----------------------------------------------------------------------------- +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