From 955b688eb463d2ac5aae395ffb2f19d2e6c1db4a Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Fri, 25 Mar 2022 22:42:44 +0000 Subject: [PATCH 1/7] Introducing MultiField --- src/atlas/CMakeLists.txt | 2 + src/atlas/field/MultiField.cc | 207 +++++++++++++++++++++++ src/atlas/field/MultiField.h | 149 ++++++++++++++++ src/tests/field/CMakeLists.txt | 5 + src/tests/field/test_multifield_ifs.cc | 224 +++++++++++++++++++++++++ 5 files changed, 587 insertions(+) create mode 100644 src/atlas/field/MultiField.cc create mode 100644 src/atlas/field/MultiField.h create mode 100644 src/tests/field/test_multifield_ifs.cc diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index 12f0db620..55aadeb6d 100644 --- a/src/atlas/CMakeLists.txt +++ b/src/atlas/CMakeLists.txt @@ -474,6 +474,8 @@ field/FieldSet.cc field/FieldSet.h field/MissingValue.cc field/MissingValue.h +field/MultiField.cc +field/MultiField.h field/State.cc field/State.h field/detail/FieldImpl.cc diff --git a/src/atlas/field/MultiField.cc b/src/atlas/field/MultiField.cc new file mode 100644 index 000000000..9220226af --- /dev/null +++ b/src/atlas/field/MultiField.cc @@ -0,0 +1,207 @@ +/* + * (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 "MultiField.h" + +#include +#include +#include +#include + +#include "eckit/thread/AutoLock.h" +#include "eckit/thread/Mutex.h" + +#include "atlas/field/Field.h" +#include "atlas/grid/Grid.h" +#include "atlas/mesh/Mesh.h" +#include "atlas/runtime/Exception.h" +#include "atlas/runtime/Log.h" + +namespace atlas { +namespace field { + +namespace { + +static eckit::Mutex* local_mutex = nullptr; +static std::map* m = nullptr; +static pthread_once_t once = PTHREAD_ONCE_INIT; + +static void init() { + local_mutex = new eckit::Mutex(); + m = new std::map(); +} + +template +void load_builder() { + MultiFieldCreatorBuilder("tmp"); +} + +void force_link() { + static struct Link { Link() = default; } link; + [](const Link&) {}(link); // disable unused warnings +} + +} // namespace + +void MultiField::initialize(const std::string& generator, const eckit::Parametrisation& params) { + std::unique_ptr FieldPool_generator(MultiFieldCreatorFactory::build(generator, params)); + FieldPool_generator->generate(*this, params); +} + +array::Array& MultiField::allocate(array::DataType datatype, array::ArraySpec&& spec) { + array_.reset(array::Array::create(datatype, std::move(spec))); + return *array_; +} + +//------------------------------------------------------------------------------------------------------ + +MultiField::MultiField() = default; + +MultiField::MultiField(const std::string& generator, const eckit::Parametrisation& params) { + initialize(generator, params); +} + +const util::Metadata& MultiField::metadata() const { + return metadata_; +} + +util::Metadata& MultiField::metadata() { + return metadata_; +} + +std::vector MultiField::field_names() const { + std::vector ret; + if (fields_.size()) { + ret.reserve(fields_.size()); + } + + for (auto it = field_index_.begin(); it != field_index_.end(); ++it) { + ret.push_back(it->first); + } + return ret; +} + +//----------------------------------------------------------------------------- + +MultiFieldCreator::MultiFieldCreator(const eckit::Parametrisation&) {} + +MultiFieldCreator::~MultiFieldCreator() = default; + +MultiFieldCreator* MultiFieldCreatorFactory::build(const std::string& name, const eckit::Parametrisation& param) { + pthread_once(&once, init); + + eckit::AutoLock lock(local_mutex); + + force_link(); + + std::map::const_iterator j = m->find(name); + + Log::debug() << "Looking for MultiFieldCreatorFactory [" << name << "]" << std::endl; + + if (j == m->end()) { + Log::error() << "No MultiFieldCreatorFactory for [" << name << "]" << std::endl; + Log::error() << "FieldPoolFactories are:" << std::endl; + for (j = m->begin(); j != m->end(); ++j) { + Log::error() << " " << (*j).first << std::endl; + } + throw_Exception(std::string("No MultiFieldCreatorFactory called ") + name, Here()); + } + + return (*j).second->make(param); +} + +void MultiFieldCreatorFactory::list(std::ostream& out) { + pthread_once(&once, init); + + eckit::AutoLock lock(local_mutex); + + force_link(); + + const char* sep = ""; + for (std::map::const_iterator j = m->begin(); j != m->end(); ++j) { + out << sep << (*j).first; + sep = ", "; + } +} + +bool MultiFieldCreatorFactory::has(const std::string& name) { + pthread_once(&once, init); + + eckit::AutoLock lock(local_mutex); + + force_link(); + + return (m->find(name) != m->end()); +} + +MultiFieldCreatorFactory::MultiFieldCreatorFactory(const std::string& name): name_(name) { + pthread_once(&once, init); + + eckit::AutoLock lock(local_mutex); + + ATLAS_ASSERT(m->find(name) == m->end()); + (*m)[name] = this; +} + +MultiFieldCreatorFactory::~MultiFieldCreatorFactory() { + eckit::AutoLock lock(local_mutex); + m->erase(name_); +} + +//----------------------------------------------------------------------------- +// C wrapper interfaces to C++ routines + +extern "C" { + +MultiField* atlas__FieldPool__new() { + return new MultiField; +} + +void atlas__FieldPool__initialize(MultiField* This, const char* generator, const eckit::Parametrisation* params) { + ATLAS_ASSERT(This != nullptr, "Reason: Use of uninitialised atlas_FieldPool"); + This->initialize(std::string(generator), *params); +} + +void atlas__FieldPool__delete(MultiField* This) { + ATLAS_ASSERT(This != nullptr, "Reason: Use of uninitialised atlas_FieldPool"); + delete This; +} + +int atlas__FieldPool__has(MultiField* This, const char* name) { + ATLAS_ASSERT(This != nullptr, "Reason: Use of uninitialised atlas_FieldPool"); + return This->has(name); +} + +FieldImpl* atlas__FieldPool__field_by_name(MultiField* This, const char* name) { + ATLAS_ASSERT(This != nullptr, "Reason: Use of uninitialised atlas_FieldPool"); + return This->field(std::string(name)).get(); +} + +FieldImpl* atlas__FieldPool__field_by_index(MultiField* This, idx_t index) { + ATLAS_ASSERT(This != nullptr, "Reason: Use of uninitialised atlas_FieldPool"); + return This->field(index).get(); +} + +idx_t atlas__FieldPool__size(const MultiField* This) { + ATLAS_ASSERT(This != nullptr, "Reason: Use of uninitialised atlas_FieldPool"); + return This->size(); +} + +util::Metadata* atlas__FieldPool__metadata(MultiField* This) { + ATLAS_ASSERT(This != nullptr, "Reason: Use of uninitialised atlas_FieldPool"); + return &This->metadata(); +} +} +//----------------------------------------------------------------------------- + +} // namespace field +} // namespace atlas diff --git a/src/atlas/field/MultiField.h b/src/atlas/field/MultiField.h new file mode 100644 index 000000000..963d16b13 --- /dev/null +++ b/src/atlas/field/MultiField.h @@ -0,0 +1,149 @@ +/* + * (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 "atlas/array/Array.h" +#include "atlas/field/Field.h" +#include "atlas/util/Config.h" +#include "atlas/util/Metadata.h" +#include "atlas/util/Object.h" + +namespace eckit { +class Parametrisation; +} + +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::Object { +public: // methods + //-- Constructors + MultiField(); + + MultiField(const std::string& generator, const eckit::Parametrisation& = util::Config()); + + //-- Accessors + + const Field& field(const std::string& name) const { return fields_[field_index_.at(name)]; } + Field& field(const std::string& name) { return fields_[field_index_.at(name)]; } + bool has(const std::string& name) const { return (field_index_.find(name) != field_index_.end()); } + std::vector field_names() const; + + const Field& field(const idx_t idx) const { return fields_[idx]; } + Field& field(const idx_t idx) { return fields_[idx]; } + idx_t size() const { return static_cast(fields_.size()); } + + const Field& operator[](const idx_t idx) const { return fields_[idx]; } + Field& operator[](const idx_t idx) { return fields_[idx]; } + + const Field& operator[](const std::string& name) const { return field(name); } + Field& operator[](const std::string& name) { return field(name); } + + const util::Metadata& metadata() const; + util::Metadata& metadata(); + + // -- Modifiers + + void initialize(const std::string& generator, const eckit::Parametrisation& = util::Config()); + + array::Array& allocate(array::DataType datatype, array::ArraySpec&&); + + /// @brief Implicit conversion to Array + operator const array::Array&() const { return *array_; } + operator array::Array&() { return *array_; } + + /// @brief Access contained Array + const array::Array& array() const { return *array_; } + array::Array& array() { return *array_; } + + +public: // temporary public for prototyping + std::map field_index_; + std::vector fields_; + std::unique_ptr array_; + util::Metadata metadata_; +}; + +//------------------------------------------------------------------------------------------------------ + +class MultiFieldCreator : public util::Object { +public: + MultiFieldCreator(const eckit::Parametrisation& = util::Config()); + + virtual ~MultiFieldCreator(); + + virtual void generate(MultiField&, const eckit::Parametrisation& = util::Config()) const = 0; +}; + +//------------------------------------------------------------------------------------------------------ + +class MultiFieldCreatorFactory { +public: + /*! + * \brief build FieldPoolCreator with options specified in parametrisation + * \return mesh generator + */ + static MultiFieldCreator* build(const std::string& FieldPool_generator, + const eckit::Parametrisation& = util::Config()); + + /*! + * \brief list all registered field creators + */ + static void list(std::ostream&); + static bool has(const std::string& name); + +private: + virtual MultiFieldCreator* make(const eckit::Parametrisation& = util::Config()) = 0; + + std::string name_; + +protected: + MultiFieldCreatorFactory(const std::string&); + virtual ~MultiFieldCreatorFactory(); +}; + +template +class MultiFieldCreatorBuilder : public MultiFieldCreatorFactory { + virtual MultiFieldCreator* make(const eckit::Parametrisation& param = util::Config()) { return new T(param); } + +public: + MultiFieldCreatorBuilder(const std::string& name): MultiFieldCreatorFactory(name) {} +}; + +// ------------------------------------------------------------------------------------ + +// C wrapper interfaces to C++ routines +extern "C" { +MultiField* atlas__FieldPool__new(); +void atlas__FieldPool__initialize(MultiField* This, const char* generator, const eckit::Parametrisation* params); +void atlas__FieldPool__delete(MultiField* This); +int atlas__FieldPool__has(MultiField* This, const char* name); +FieldImpl* atlas__FieldPool__field_by_name(MultiField* This, const char* name); +FieldImpl* atlas__FieldPool__field_by_index(MultiField* This, idx_t index); +idx_t atlas__FieldPool__size(const MultiField* This); +util::Metadata* atlas__FieldPool__metadata(MultiField* This); +} + +} // namespace field +} // namespace atlas diff --git a/src/tests/field/CMakeLists.txt b/src/tests/field/CMakeLists.txt index dc91d6f65..e103c2886 100644 --- a/src/tests/field/CMakeLists.txt +++ b/src/tests/field/CMakeLists.txt @@ -23,6 +23,11 @@ ecbuild_add_test( TARGET atlas_test_field_foreach LIBS atlas OMP 4 ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} + +ecbuild_add_test( TARGET atlas_test_multifield_ifs + SOURCES test_multifield_ifs.cc + LIBS atlas + ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) ecbuild_add_test( TARGET atlas_test_field_acc diff --git a/src/tests/field/test_multifield_ifs.cc b/src/tests/field/test_multifield_ifs.cc new file mode 100644 index 000000000..e458fe1a2 --- /dev/null +++ b/src/tests/field/test_multifield_ifs.cc @@ -0,0 +1,224 @@ +/* + * (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 "eckit/log/JSON.h" +#include "eckit/parser/JSONParser.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/mesh/Mesh.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::Parametrisation& p = util::Config()): MultiFieldCreator(p) {} + ~MultiFieldCreatorIFS() override = default; + void generate(MultiField& fieldpool, const eckit::Parametrisation& p = util::Config()) const override; +}; + +// --- Implementation (in .cc file) +void MultiFieldCreatorIFS::generate(MultiField& multifield, const eckit::Parametrisation& p) const { + const eckit::LocalConfiguration* params = dynamic_cast(&p); + if (!params) { + throw_Exception("Parametrisation has to be of atlas::util::Config type"); + } + + + long nproma = 0; + long ngptot = 0; + ; + long nfld = 0; + ; + long nlev = 0; + ; + long nblk = 0; + ; + + if (!p.get("ngptot", ngptot)) { + std::string grid_uid; + if (p.get("grid", grid_uid)) { + ngptot = Grid(grid_uid).size(); + } + else { + throw_Exception("Could not find 'ngptot' in Parametrisation"); + } + } + + if (!p.get("nproma", nproma)) { + throw_Exception("Could not find 'nproma' in Parametrisation"); + } + + if (!p.get("nlev", nlev)) { + throw_Exception("Could not find 'nlev' in Parametrisation"); + } + + array::DataType datatype = array::DataType::create(); + std::string datatype_str; + if (p.get("datatype", datatype_str)) { + datatype = array::DataType(datatype_str); + } + else { + array::DataType::kind_t kind(array::DataType::kind()); + p.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)); + + std::vector fields; + params->get("fields", fields); + nfld = fields.size(); + + auto multiarray_spec = array::make_shape(nblk, nfld, nlev, nproma); + auto& multiarray = multifield.allocate(datatype, multiarray_spec); + + 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); + + for (size_t i = 0; i < fields.size(); ++i) { + std::string name; + fields[i].get("name", name); + Field field; + constexpr auto all = array::Range::all(); + if (datatype.kind() == array::DataType::KIND_REAL64) { + auto slice = array::make_view(multiarray).slice(all, i, 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, i, all, all); + field = Field(name, slice.data(), field_array_spec); + } + else { + ATLAS_NOTIMPLEMENTED; + } + field.set_levels(nlev); + multifield.fields_.emplace_back(field); + multifield.field_index_[field.name()] = i; + } +} + +// Register in factory +MultiFieldCreatorBuilder __MultiFieldCreatorIFS("MultiFieldCreatorIFS"); + +// =================================================================== +// BEGIN TESTS +// =================================================================== + + +CASE("multifield_generator") { + EXPECT(MultiFieldCreatorFactory::has("MultiFieldCreatorIFS")); + std::unique_ptr MultiFieldCreator(MultiFieldCreatorFactory::build("MultiFieldCreatorIFS")); +} + + +CASE("multifield_create") { + int nproma = 16; + int nlev = 100; + int ngptot = 2000; + + auto json = [&]() -> std::string { + util::Config p; + p.set("ngptot", ngptot); + p.set("nproma", nproma); + p.set("nlev", nlev); + p.set("datatype", array::DataType::real64().str()); + + std::vector fields(3); + fields[0].set("name", "temperature"); + fields[1].set("name", "pressure"); + fields[2].set("name", "density"); + + p.set("fields", fields); + + // We can also translate parameters to a json: + std::stringstream json; + eckit::JSON js(json); + js << p; + std::string json_str = json.str(); + Log::info() << "json = " << json_str << std::endl; + return json_str; + }; + + + // And we can create back parameters from json: + std::stringstream json_stream; + json_stream << json(); + util::Config config(json_stream); + + MultiField multifield("MultiFieldCreatorIFS", config); + + EXPECT(multifield.has("temperature")); + EXPECT(multifield.has("pressure")); + EXPECT(multifield.has("density")); + + Log::info() << multifield.field("temperature") << std::endl; + Log::info() << multifield.field("pressure") << std::endl; + Log::info() << multifield.field("density") << 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")); + + temp(1, 2, 3) = 4; + pres(5, 6, 7) = 8; + dens(9, 10, 11) = 12; + + EXPECT_EQ(temp.stride(2), 1); + EXPECT_EQ(temp.stride(1), nproma); + EXPECT_EQ(temp.stride(0), nproma * nlev * multifield.size()); + + // Underlying array of MultiField + { + auto multiarray = array::make_view(multifield); + EXPECT_EQ(multiarray.stride(3), 1); + EXPECT_EQ(multiarray.stride(2), nproma); + EXPECT_EQ(multiarray.stride(1), nproma * nlev); + EXPECT_EQ(multiarray.stride(0), nproma * nlev * multifield.size()); + + EXPECT_EQ(multiarray(1, 0, 2, 3), 4.); + EXPECT_EQ(multiarray(5, 1, 6, 7), 8.); + EXPECT_EQ(multiarray(9, 2, 10, 11), 12.); + } +} + +//----------------------------------------------------------------------------- + +} // namespace test +} // namespace atlas + +int main(int argc, char** argv) { + return atlas::test::run(argc, argv); +} From 8a9245e31d51d472061afd6ba112987ae6c79542 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 28 Mar 2022 15:51:11 +0100 Subject: [PATCH 2/7] Support for vector fields in multifield --- src/tests/field/test_multifield_ifs.cc | 175 +++++++++++++++++-------- 1 file changed, 123 insertions(+), 52 deletions(-) diff --git a/src/tests/field/test_multifield_ifs.cc b/src/tests/field/test_multifield_ifs.cc index e458fe1a2..8ea78cd67 100644 --- a/src/tests/field/test_multifield_ifs.cc +++ b/src/tests/field/test_multifield_ifs.cc @@ -50,25 +50,14 @@ void MultiFieldCreatorIFS::generate(MultiField& multifield, const eckit::Paramet throw_Exception("Parametrisation has to be of atlas::util::Config type"); } - long nproma = 0; long ngptot = 0; - ; - long nfld = 0; - ; - long nlev = 0; - ; - long nblk = 0; - ; + long nfld = 0; + long nlev = 0; + long nblk = 0; if (!p.get("ngptot", ngptot)) { - std::string grid_uid; - if (p.get("grid", grid_uid)) { - ngptot = Grid(grid_uid).size(); - } - else { - throw_Exception("Could not find 'ngptot' in Parametrisation"); - } + throw_Exception("Could not find 'ngptot' in Parametrisation"); } if (!p.get("nproma", nproma)) { @@ -99,34 +88,69 @@ void MultiFieldCreatorIFS::generate(MultiField& multifield, const eckit::Paramet std::vector fields; params->get("fields", fields); - nfld = fields.size(); + nfld = fields.size(); + size_t nvar_tot = 0; + for (const auto& field_config : fields) { + size_t nvar = 1; + field_config.get("nvar", nvar); + nvar_tot += nvar; + } - auto multiarray_spec = array::make_shape(nblk, nfld, nlev, nproma); + auto multiarray_spec = array::make_shape(nblk, nvar_tot, nlev, nproma); auto& multiarray = multifield.allocate(datatype, multiarray_spec); - 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); - + 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; - constexpr auto all = array::Range::all(); - if (datatype.kind() == array::DataType::KIND_REAL64) { - auto slice = array::make_view(multiarray).slice(all, i, 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, i, all, all); - field = Field(name, slice.data(), field_array_spec); + 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 { - ATLAS_NOTIMPLEMENTED; + 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.fields_.emplace_back(field); + ATLAS_ASSERT(multifield.field_index_.find(field.name()) == multifield.field_index_.end(), + "Field with name already exists!"); multifield.field_index_[field.name()] = i; + multiarray_field_idx += field_vars; } } @@ -154,63 +178,110 @@ CASE("multifield_create") { p.set("ngptot", ngptot); p.set("nproma", nproma); p.set("nlev", nlev); - p.set("datatype", array::DataType::real64().str()); + p.set("datatype", array::make_datatype().str()); + p.set("fields", [] { + std::vector fields(5); + fields[0].set("name", "temperature"); + + fields[1].set("name", "pressure"); - std::vector fields(3); - fields[0].set("name", "temperature"); - fields[1].set("name", "pressure"); - fields[2].set("name", "density"); + fields[2].set("name", "density"); - p.set("fields", fields); + fields[3].set("name", "clv"); + fields[3].set("nvar", 5); + + fields[4].set("name", "wind_u"); + return fields; + }()); // We can also translate parameters to a json: std::stringstream json; eckit::JSON js(json); js << p; - std::string json_str = json.str(); - Log::info() << "json = " << json_str << std::endl; - return json_str; + return json.str(); }; // And we can create back parameters from json: + Log::info() << "json = " << json() << std::endl; std::stringstream json_stream; json_stream << json(); util::Config config(json_stream); MultiField multifield("MultiFieldCreatorIFS", config); + const auto nblk = multifield.array().shape(0); + const auto nfld = multifield.size(); + EXPECT_EQ(nfld, 5); + + 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 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")); - temp(1, 2, 3) = 4; - pres(5, 6, 7) = 8; - dens(9, 10, 11) = 12; + // or + { + auto temp = array::make_view(multifield[0]); + auto pres = array::make_view(multifield[1]); + auto dens = array::make_view(multifield[2]); + auto clv = array::make_view(multifield[3]); // note rank 4 + auto wind_u = array::make_view(multifield[4]); + } + + 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); - EXPECT_EQ(temp.stride(2), 1); - EXPECT_EQ(temp.stride(1), nproma); - EXPECT_EQ(temp.stride(0), nproma * nlev * multifield.size()); // Underlying array of MultiField { auto multiarray = array::make_view(multifield); - EXPECT_EQ(multiarray.stride(3), 1); - EXPECT_EQ(multiarray.stride(2), nproma); - EXPECT_EQ(multiarray.stride(1), nproma * nlev); - EXPECT_EQ(multiarray.stride(0), nproma * nlev * multifield.size()); + 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 * (nfld - 1 + 5) * nlev * nproma); } } From f2a2f5445fb6ca2a74d3fc4db73812e2a54df673 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Mon, 4 Apr 2022 15:56:40 +0100 Subject: [PATCH 3/7] Use util::Factory for MultiFieldCreatorFactory --- src/atlas/field/MultiField.cc | 142 ++++------------------------------ src/atlas/field/MultiField.h | 43 +++------- 2 files changed, 27 insertions(+), 158 deletions(-) diff --git a/src/atlas/field/MultiField.cc b/src/atlas/field/MultiField.cc index 9220226af..aeeac5fcb 100644 --- a/src/atlas/field/MultiField.cc +++ b/src/atlas/field/MultiField.cc @@ -8,8 +8,6 @@ * 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 "MultiField.h" #include @@ -30,31 +28,20 @@ namespace atlas { namespace field { namespace { - -static eckit::Mutex* local_mutex = nullptr; -static std::map* m = nullptr; -static pthread_once_t once = PTHREAD_ONCE_INIT; - -static void init() { - local_mutex = new eckit::Mutex(); - m = new std::map(); -} - -template -void load_builder() { - MultiFieldCreatorBuilder("tmp"); -} - void force_link() { - static struct Link { Link() = default; } link; - [](const Link&) {}(link); // disable unused warnings + static struct Link { + Link() { + ; + // For static linking add here something like + // MultiFieldCreatorBuilder(); + } + } link; } - } // namespace -void MultiField::initialize(const std::string& generator, const eckit::Parametrisation& params) { - std::unique_ptr FieldPool_generator(MultiFieldCreatorFactory::build(generator, params)); - FieldPool_generator->generate(*this, params); +void MultiField::initialize(const std::string& creator, const eckit::Parametrisation& params) { + std::unique_ptr MultiField_creator(MultiFieldCreatorFactory::build(creator, params)); + MultiField_creator->generate(*this, params); } array::Array& MultiField::allocate(array::DataType datatype, array::ArraySpec&& spec) { @@ -66,8 +53,8 @@ array::Array& MultiField::allocate(array::DataType datatype, array::ArraySpec&& MultiField::MultiField() = default; -MultiField::MultiField(const std::string& generator, const eckit::Parametrisation& params) { - initialize(generator, params); +MultiField::MultiField(const std::string& creator, const eckit::Parametrisation& params) { + initialize(creator, params); } const util::Metadata& MultiField::metadata() const { @@ -96,111 +83,12 @@ MultiFieldCreator::MultiFieldCreator(const eckit::Parametrisation&) {} MultiFieldCreator::~MultiFieldCreator() = default; -MultiFieldCreator* MultiFieldCreatorFactory::build(const std::string& name, const eckit::Parametrisation& param) { - pthread_once(&once, init); - - eckit::AutoLock lock(local_mutex); - - force_link(); - - std::map::const_iterator j = m->find(name); - - Log::debug() << "Looking for MultiFieldCreatorFactory [" << name << "]" << std::endl; - - if (j == m->end()) { - Log::error() << "No MultiFieldCreatorFactory for [" << name << "]" << std::endl; - Log::error() << "FieldPoolFactories are:" << std::endl; - for (j = m->begin(); j != m->end(); ++j) { - Log::error() << " " << (*j).first << std::endl; - } - throw_Exception(std::string("No MultiFieldCreatorFactory called ") + name, Here()); - } - - return (*j).second->make(param); -} - -void MultiFieldCreatorFactory::list(std::ostream& out) { - pthread_once(&once, init); - - eckit::AutoLock lock(local_mutex); - +MultiFieldCreator* MultiFieldCreatorFactory::build(const std::string& builder, const eckit::Parametrisation& config) { force_link(); - - const char* sep = ""; - for (std::map::const_iterator j = m->begin(); j != m->end(); ++j) { - out << sep << (*j).first; - sep = ", "; - } + auto factory = get(builder); + return factory->make(config); } -bool MultiFieldCreatorFactory::has(const std::string& name) { - pthread_once(&once, init); - - eckit::AutoLock lock(local_mutex); - - force_link(); - - return (m->find(name) != m->end()); -} - -MultiFieldCreatorFactory::MultiFieldCreatorFactory(const std::string& name): name_(name) { - pthread_once(&once, init); - - eckit::AutoLock lock(local_mutex); - - ATLAS_ASSERT(m->find(name) == m->end()); - (*m)[name] = this; -} - -MultiFieldCreatorFactory::~MultiFieldCreatorFactory() { - eckit::AutoLock lock(local_mutex); - m->erase(name_); -} - -//----------------------------------------------------------------------------- -// C wrapper interfaces to C++ routines - -extern "C" { - -MultiField* atlas__FieldPool__new() { - return new MultiField; -} - -void atlas__FieldPool__initialize(MultiField* This, const char* generator, const eckit::Parametrisation* params) { - ATLAS_ASSERT(This != nullptr, "Reason: Use of uninitialised atlas_FieldPool"); - This->initialize(std::string(generator), *params); -} - -void atlas__FieldPool__delete(MultiField* This) { - ATLAS_ASSERT(This != nullptr, "Reason: Use of uninitialised atlas_FieldPool"); - delete This; -} - -int atlas__FieldPool__has(MultiField* This, const char* name) { - ATLAS_ASSERT(This != nullptr, "Reason: Use of uninitialised atlas_FieldPool"); - return This->has(name); -} - -FieldImpl* atlas__FieldPool__field_by_name(MultiField* This, const char* name) { - ATLAS_ASSERT(This != nullptr, "Reason: Use of uninitialised atlas_FieldPool"); - return This->field(std::string(name)).get(); -} - -FieldImpl* atlas__FieldPool__field_by_index(MultiField* This, idx_t index) { - ATLAS_ASSERT(This != nullptr, "Reason: Use of uninitialised atlas_FieldPool"); - return This->field(index).get(); -} - -idx_t atlas__FieldPool__size(const MultiField* This) { - ATLAS_ASSERT(This != nullptr, "Reason: Use of uninitialised atlas_FieldPool"); - return This->size(); -} - -util::Metadata* atlas__FieldPool__metadata(MultiField* This) { - ATLAS_ASSERT(This != nullptr, "Reason: Use of uninitialised atlas_FieldPool"); - return &This->metadata(); -} -} //----------------------------------------------------------------------------- } // namespace field diff --git a/src/atlas/field/MultiField.h b/src/atlas/field/MultiField.h index 963d16b13..9b03c5b03 100644 --- a/src/atlas/field/MultiField.h +++ b/src/atlas/field/MultiField.h @@ -18,6 +18,7 @@ #include "atlas/array/Array.h" #include "atlas/field/Field.h" #include "atlas/util/Config.h" +#include "atlas/util/Factory.h" #include "atlas/util/Metadata.h" #include "atlas/util/Object.h" @@ -98,52 +99,32 @@ class MultiFieldCreator : public util::Object { //------------------------------------------------------------------------------------------------------ -class MultiFieldCreatorFactory { + +class MultiFieldCreatorFactory : public util::Factory { public: - /*! - * \brief build FieldPoolCreator with options specified in parametrisation - * \return mesh generator - */ - static MultiFieldCreator* build(const std::string& FieldPool_generator, - const eckit::Parametrisation& = util::Config()); + static std::string className() { return "MultiFieldCreatorFactory"; } /*! - * \brief list all registered field creators + * \brief build MultiFieldCreator with options specified in parametrisation + * \return MutliField creator */ - static void list(std::ostream&); - static bool has(const std::string& name); - -private: - virtual MultiFieldCreator* make(const eckit::Parametrisation& = util::Config()) = 0; + static MultiFieldCreator* build(const std::string&, const eckit::Parametrisation& = util::NoConfig()); - std::string name_; + using Factory::Factory; -protected: - MultiFieldCreatorFactory(const std::string&); - virtual ~MultiFieldCreatorFactory(); +private: + virtual MultiFieldCreator* make(const eckit::Parametrisation&) = 0; }; template class MultiFieldCreatorBuilder : public MultiFieldCreatorFactory { - virtual MultiFieldCreator* make(const eckit::Parametrisation& param = util::Config()) { return new T(param); } + virtual MultiFieldCreator* make(const eckit::Parametrisation& param) override { return new T(param); } public: - MultiFieldCreatorBuilder(const std::string& name): MultiFieldCreatorFactory(name) {} + using MultiFieldCreatorFactory::MultiFieldCreatorFactory; }; // ------------------------------------------------------------------------------------ -// C wrapper interfaces to C++ routines -extern "C" { -MultiField* atlas__FieldPool__new(); -void atlas__FieldPool__initialize(MultiField* This, const char* generator, const eckit::Parametrisation* params); -void atlas__FieldPool__delete(MultiField* This); -int atlas__FieldPool__has(MultiField* This, const char* name); -FieldImpl* atlas__FieldPool__field_by_name(MultiField* This, const char* name); -FieldImpl* atlas__FieldPool__field_by_index(MultiField* This, idx_t index); -idx_t atlas__FieldPool__size(const MultiField* This); -util::Metadata* atlas__FieldPool__metadata(MultiField* This); -} - } // namespace field } // namespace atlas From 0b1d96726e532512c798cf72cd647e960f15e4e7 Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 5 Apr 2022 22:30:50 +0100 Subject: [PATCH 4/7] MultiField with ArrayAllocator --- src/atlas/field/MultiField.cc | 51 ++------ src/atlas/field/MultiField.h | 160 ++++++++++++++++++++----- src/tests/field/test_multifield_ifs.cc | 142 +++++++++------------- 3 files changed, 194 insertions(+), 159 deletions(-) diff --git a/src/atlas/field/MultiField.cc b/src/atlas/field/MultiField.cc index aeeac5fcb..463b3bc65 100644 --- a/src/atlas/field/MultiField.cc +++ b/src/atlas/field/MultiField.cc @@ -39,56 +39,27 @@ void force_link() { } } // namespace -void MultiField::initialize(const std::string& creator, const eckit::Parametrisation& params) { - std::unique_ptr MultiField_creator(MultiFieldCreatorFactory::build(creator, params)); - MultiField_creator->generate(*this, params); -} - -array::Array& MultiField::allocate(array::DataType datatype, array::ArraySpec&& spec) { - array_.reset(array::Array::create(datatype, std::move(spec))); - return *array_; -} - -//------------------------------------------------------------------------------------------------------ - -MultiField::MultiField() = default; - -MultiField::MultiField(const std::string& creator, const eckit::Parametrisation& params) { - initialize(creator, params); -} - -const util::Metadata& MultiField::metadata() const { - return metadata_; -} - -util::Metadata& MultiField::metadata() { - return metadata_; -} - -std::vector MultiField::field_names() const { - std::vector ret; - if (fields_.size()) { - ret.reserve(fields_.size()); - } - - for (auto it = field_index_.begin(); it != field_index_.end(); ++it) { - ret.push_back(it->first); - } - return ret; -} - //----------------------------------------------------------------------------- -MultiFieldCreator::MultiFieldCreator(const eckit::Parametrisation&) {} +MultiFieldCreator::MultiFieldCreator(const eckit::Configuration&) {} MultiFieldCreator::~MultiFieldCreator() = default; -MultiFieldCreator* MultiFieldCreatorFactory::build(const std::string& builder, const eckit::Parametrisation& config) { +MultiFieldCreator* MultiFieldCreatorFactory::build(const std::string& builder, const eckit::Configuration& config) { force_link(); auto factory = get(builder); return factory->make(config); } +atlas::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)); +} + //----------------------------------------------------------------------------- } // namespace field diff --git a/src/atlas/field/MultiField.h b/src/atlas/field/MultiField.h index 9b03c5b03..773a098e3 100644 --- a/src/atlas/field/MultiField.h +++ b/src/atlas/field/MultiField.h @@ -14,13 +14,16 @@ #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; @@ -37,69 +40,162 @@ namespace field { * * Fields have to all be of same memory layout and data type */ -class MultiField : public util::Object { + +namespace detail { +class ArrayAllocator { +public: + virtual array::Array* allocate(const array::ArraySpec&) = 0; + virtual void deallocate() = 0; + virtual ~ArrayAllocator() = default; +}; + +class NoArrayAllocator : public ArrayAllocator { +public: + virtual array::Array* allocate(const array::ArraySpec& arrayspec) override { return nullptr; } + virtual void deallocate() override {} +}; + +class StandardArrayAllocator : public ArrayAllocator { +public: + virtual array::Array* allocate(const array::ArraySpec& arrayspec) override { + array::ArraySpec spec(arrayspec); + return array::Array::create(std::move(spec)); + } + virtual void deallocate() override { + // Nothing to deallocate. + // The Array owns the allocated data. + } +}; +} // namespace detail + +class MultiFieldImpl : public util::Object { public: // methods - //-- Constructors - MultiField(); + //-- Constructors + + MultiFieldImpl() { allocator_ = std::unique_ptr(new detail::NoArrayAllocator()); } + + MultiFieldImpl(const array::ArraySpec& spec, detail::ArrayAllocator* allocator = nullptr) { + if (allocator) { + allocator_ = std::unique_ptr(allocator); + } + else { + allocator_ = std::unique_ptr(new detail::StandardArrayAllocator()); + } + ATLAS_ASSERT(allocator_); + array_ = std::unique_ptr(allocator_->allocate(spec)); + } + + virtual ~MultiFieldImpl() { allocator_->deallocate(); } - MultiField(const std::string& generator, const eckit::Parametrisation& = util::Config()); //-- Accessors - const Field& field(const std::string& name) const { return fields_[field_index_.at(name)]; } - Field& field(const std::string& name) { return fields_[field_index_.at(name)]; } - bool has(const std::string& name) const { return (field_index_.find(name) != field_index_.end()); } - std::vector field_names() const; + 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 fields_[idx]; } - Field& field(const idx_t idx) { return fields_[idx]; } - idx_t size() const { return static_cast(fields_.size()); } + 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 fields_[idx]; } - Field& operator[](const idx_t idx) { return fields_[idx]; } + 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 field(name); } - Field& operator[](const std::string& name) { return field(name); } + 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; - util::Metadata& metadata(); + const util::Metadata& metadata() const { return metadata_; } + util::Metadata& metadata() { return metadata_; } // -- Modifiers - void initialize(const std::string& generator, const eckit::Parametrisation& = util::Config()); + /// @brief Implicit conversion to Array + operator const array::Array&() const { return array(); } + operator array::Array&() { return array(); } - array::Array& allocate(array::DataType datatype, array::ArraySpec&&); + operator const FieldSet&() const { return fieldset_; } - /// @brief Implicit conversion to Array - operator const array::Array&() const { return *array_; } - operator array::Array&() { return *array_; } + operator FieldSet&() { return fieldset_; } /// @brief Access contained Array - const array::Array& array() const { return *array_; } - array::Array& array() { return *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_; } public: // temporary public for prototyping - std::map field_index_; - std::vector fields_; + FieldSet fieldset_; std::unique_ptr array_; + std::unique_ptr allocator_; util::Metadata metadata_; }; + +class MultiField : public util::ObjectHandle { +public: // methods + //-- Constructors + using Handle::Handle; + + MultiField(const eckit::Configuration&); + + //-- 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 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& operator[](const idx_t idx) const { return get()->field(idx); } + Field& operator[](const idx_t idx) { return get()->field(idx); } + + const Field& operator[](const std::string& name) const { return get()->field(name); } + Field& operator[](const std::string& name) { return get()->field(name); } + + const util::Metadata& metadata() const { return get()->metadata(); } + util::Metadata& metadata() { return get()->metadata(); } + + // -- Modifiers + + /// @brief Implicit conversion to Array + operator const array::Array&() const { return get()->array(); } + operator array::Array&() { return get()->array(); } + + operator const FieldSet&() const { return get()->fieldset_; } + operator FieldSet&() { return get()->fieldset_; } + + /// @brief Access contained Array + const array::Array& array() const { return get()->array(); } + array::Array& array() { return get()->array(); } +}; + + //------------------------------------------------------------------------------------------------------ class MultiFieldCreator : public util::Object { public: - MultiFieldCreator(const eckit::Parametrisation& = util::Config()); + MultiFieldCreator(const eckit::Configuration& = util::Config()); virtual ~MultiFieldCreator(); - virtual void generate(MultiField&, const eckit::Parametrisation& = util::Config()) const = 0; + virtual MultiFieldImpl* create(const eckit::Configuration& = util::Config()) const = 0; }; //------------------------------------------------------------------------------------------------------ - class MultiFieldCreatorFactory : public util::Factory { public: static std::string className() { return "MultiFieldCreatorFactory"; } @@ -108,17 +204,17 @@ class MultiFieldCreatorFactory : public util::Factory * \brief build MultiFieldCreator with options specified in parametrisation * \return MutliField creator */ - static MultiFieldCreator* build(const std::string&, const eckit::Parametrisation& = util::NoConfig()); + static MultiFieldCreator* build(const std::string&, const eckit::Configuration& = util::NoConfig()); using Factory::Factory; private: - virtual MultiFieldCreator* make(const eckit::Parametrisation&) = 0; + virtual MultiFieldCreator* make(const eckit::Configuration&) = 0; }; template class MultiFieldCreatorBuilder : public MultiFieldCreatorFactory { - virtual MultiFieldCreator* make(const eckit::Parametrisation& param) override { return new T(param); } + virtual MultiFieldCreator* make(const eckit::Configuration& config) override { return new T(config); } public: using MultiFieldCreatorFactory::MultiFieldCreatorFactory; diff --git a/src/tests/field/test_multifield_ifs.cc b/src/tests/field/test_multifield_ifs.cc index 8ea78cd67..87db13e45 100644 --- a/src/tests/field/test_multifield_ifs.cc +++ b/src/tests/field/test_multifield_ifs.cc @@ -11,8 +11,7 @@ #include #include -#include "eckit/log/JSON.h" -#include "eckit/parser/JSONParser.h" +#include "eckit/config/YAMLConfiguration.h" #include "atlas/array/ArrayView.h" #include "atlas/array/DataType.h" @@ -20,7 +19,6 @@ #include "atlas/field/Field.h" #include "atlas/field/MultiField.h" #include "atlas/grid/Grid.h" -#include "atlas/mesh/Mesh.h" #include "atlas/runtime/Exception.h" #include "atlas/runtime/Log.h" @@ -38,44 +36,27 @@ namespace test { // --- Declaration (in .h file) class MultiFieldCreatorIFS : public MultiFieldCreator { public: - MultiFieldCreatorIFS(const eckit::Parametrisation& p = util::Config()): MultiFieldCreator(p) {} + MultiFieldCreatorIFS(const eckit::Configuration& config = util::Config()): MultiFieldCreator(config) {} ~MultiFieldCreatorIFS() override = default; - void generate(MultiField& fieldpool, const eckit::Parametrisation& p = util::Config()) const override; + MultiFieldImpl* create(const eckit::Configuration& = util::Config()) const override; }; // --- Implementation (in .cc file) -void MultiFieldCreatorIFS::generate(MultiField& multifield, const eckit::Parametrisation& p) const { - const eckit::LocalConfiguration* params = dynamic_cast(&p); - if (!params) { - throw_Exception("Parametrisation has to be of atlas::util::Config type"); - } - - long nproma = 0; - long ngptot = 0; - long nfld = 0; - long nlev = 0; +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; - if (!p.get("ngptot", ngptot)) { - throw_Exception("Could not find 'ngptot' in Parametrisation"); - } - - if (!p.get("nproma", nproma)) { - throw_Exception("Could not find 'nproma' in Parametrisation"); - } - - if (!p.get("nlev", nlev)) { - throw_Exception("Could not find 'nlev' in Parametrisation"); - } array::DataType datatype = array::DataType::create(); std::string datatype_str; - if (p.get("datatype", datatype_str)) { + if (config.get("datatype", datatype_str)) { datatype = array::DataType(datatype_str); } else { array::DataType::kind_t kind(array::DataType::kind()); - p.get("kind", kind); + config.get("kind", kind); if (!array::DataType::kind_valid(kind)) { std::stringstream msg; msg << "Could not create field. kind parameter unrecognized"; @@ -86,18 +67,20 @@ void MultiFieldCreatorIFS::generate(MultiField& multifield, const eckit::Paramet nblk = std::ceil(static_cast(ngptot) / static_cast(nproma)); - std::vector fields; - params->get("fields", fields); - nfld = fields.size(); - size_t nvar_tot = 0; + auto fields = config.getSubConfigurations("fields"); + long nfld = 0; for (const auto& field_config : fields) { - size_t nvar = 1; + long nvar = 1; field_config.get("nvar", nvar); - nvar_tot += nvar; + nfld += nvar; } - auto multiarray_spec = array::make_shape(nblk, nvar_tot, nlev, nproma); - auto& multiarray = multifield.allocate(datatype, multiarray_spec); + auto multiarray_shape = array::make_shape(nblk, nfld, nlev, nproma); + + MultiFieldImpl* multifield = new MultiFieldImpl{array::ArraySpec{datatype, multiarray_shape}}; + + auto& multiarray = multifield->array(); + auto& fieldset = multifield->fieldset(); size_t multiarray_field_idx = 0; for (size_t i = 0; i < fields.size(); ++i) { @@ -146,12 +129,13 @@ void MultiFieldCreatorIFS::generate(MultiField& multifield, const eckit::Paramet } } field.set_levels(nlev); - multifield.fields_.emplace_back(field); - ATLAS_ASSERT(multifield.field_index_.find(field.name()) == multifield.field_index_.end(), - "Field with name already exists!"); - multifield.field_index_[field.name()] = i; + ATLAS_ASSERT(not fieldset.has(field.name()), "Field with name \"" + field.name() + "\" already exists!"); + + fieldset.add(field); + multiarray_field_idx += field_vars; } + return multifield; } // Register in factory @@ -169,50 +153,37 @@ CASE("multifield_generator") { CASE("multifield_create") { - int nproma = 16; - int nlev = 100; - int ngptot = 2000; + 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 fields(5); - fields[0].set("name", "temperature"); - - fields[1].set("name", "pressure"); - - fields[2].set("name", "density"); - - fields[3].set("name", "clv"); - fields[3].set("nvar", 5); - - fields[4].set("name", "wind_u"); - return fields; - }()); - - // We can also translate parameters to a json: - std::stringstream json; - eckit::JSON js(json); - js << p; - return json.str(); + 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") // + }); + return p.json(); }; - - // And we can create back parameters from json: Log::info() << "json = " << json() << std::endl; - std::stringstream json_stream; - json_stream << json(); - util::Config config(json_stream); - MultiField multifield("MultiFieldCreatorIFS", 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, 5); + EXPECT_EQ(nvar, 9); EXPECT_EQ(multifield.size(), 5); EXPECT(multifield.has("temperature")); @@ -227,20 +198,17 @@ CASE("multifield_create") { 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")); + 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")); - // or - { - auto temp = array::make_view(multifield[0]); - auto pres = array::make_view(multifield[1]); - auto dens = array::make_view(multifield[2]); - auto clv = array::make_view(multifield[3]); // note rank 4 - auto wind_u = array::make_view(multifield[4]); - } + 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; @@ -266,9 +234,10 @@ CASE("multifield_create") { EXPECT_EQ(clv.size(), nblk * 5 * nlev * nproma); - // Underlying array of MultiField + // 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); + 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); @@ -280,8 +249,7 @@ CASE("multifield_create") { EXPECT_EQ(multiarray(13, 5, 14, 15), 16.); EXPECT_EQ(multiarray(17, 8, 18, 3), 19.); - - EXPECT_EQ(multiarray.size(), nblk * (nfld - 1 + 5) * nlev * nproma); + EXPECT_EQ(multiarray.size(), nblk * nvar * nlev * nproma); } } From 40e801ec6bbc892aacfb35ca0b3e14153a94af02 Mon Sep 17 00:00:00 2001 From: Slavko Brdar Date: Thu, 31 Aug 2023 17:52:59 +0200 Subject: [PATCH 5/7] add Fortran interface and a unit test for multifield --- src/atlas/CMakeLists.txt | 2 + src/atlas/field/detail/MultiFieldInterface.cc | 133 ++++++++++++++++++ src/atlas/field/detail/MultiFieldInterface.h | 38 +++++ src/atlas_f/CMakeLists.txt | 4 + .../field/atlas_MultiField_module.fypp | 102 ++++++++++++++ src/tests/field/CMakeLists.txt | 18 ++- src/tests/field/fctest_multifield_ifs.F90 | 77 ++++++++++ 7 files changed, 369 insertions(+), 5 deletions(-) create mode 100644 src/atlas/field/detail/MultiFieldInterface.cc create mode 100644 src/atlas/field/detail/MultiFieldInterface.h create mode 100644 src/atlas_f/field/atlas_MultiField_module.fypp create mode 100644 src/tests/field/fctest_multifield_ifs.F90 diff --git a/src/atlas/CMakeLists.txt b/src/atlas/CMakeLists.txt index 55aadeb6d..58fed83e8 100644 --- a/src/atlas/CMakeLists.txt +++ b/src/atlas/CMakeLists.txt @@ -482,6 +482,8 @@ field/detail/FieldImpl.cc field/detail/FieldImpl.h field/detail/FieldInterface.cc field/detail/FieldInterface.h +field/detail/MultiFieldInterface.cc +field/detail/MultiFieldInterface.h field/detail/MissingValue.cc field/detail/MissingValue.h ) diff --git a/src/atlas/field/detail/MultiFieldInterface.cc b/src/atlas/field/detail/MultiFieldInterface.cc new file mode 100644 index 000000000..685867685 --- /dev/null +++ b/src/atlas/field/detail/MultiFieldInterface.cc @@ -0,0 +1,133 @@ +/* + * (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 "atlas/library/config.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); + 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 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* field = new MultiFieldImpl{array::ArraySpec{datatype, multiarray_shape}}; + auto& multiarray = field->array(); + auto& fieldset = field->fieldset(); + + 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); + ATLAS_ASSERT(not fieldset.has(field.name()), "Field with name \"" + field.name() + "\" already exists!"); + + fieldset.add(field); + + multiarray_field_idx += field_vars; + } + ATLAS_ASSERT(field); + return field; +} + +void atlas__MultiField__delete(MultiFieldImpl* This) { + delete This; +} + +} + +// ------------------------------------------------------------------ + +} // 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..ac3b1ae22 --- /dev/null +++ b/src/atlas/field/detail/MultiFieldInterface.h @@ -0,0 +1,38 @@ +/* + * (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/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); +void atlas__MultiField__delete(MultiFieldImpl* This); +} + +//---------------------------------------------------------------------------------------------------------------------- + +} // namespace field +} // namespace atlas diff --git a/src/atlas_f/CMakeLists.txt b/src/atlas_f/CMakeLists.txt index 60f55318f..eefbbf647 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.fypp 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.fypp new file mode 100644 index 000000000..4d45836a5 --- /dev/null +++ b/src/atlas_f/field/atlas_MultiField_module.fypp @@ -0,0 +1,102 @@ +! (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" + +#: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 + +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 + +#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 +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 + +!------------------------------------------------------------------------------- + +#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 e103c2886..dbf51ee37 100644 --- a/src/tests/field/CMakeLists.txt +++ b/src/tests/field/CMakeLists.txt @@ -23,11 +23,6 @@ ecbuild_add_test( TARGET atlas_test_field_foreach LIBS atlas OMP 4 ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} - -ecbuild_add_test( TARGET atlas_test_multifield_ifs - SOURCES test_multifield_ifs.cc - LIBS atlas - ENVIRONMENT ${ATLAS_TEST_ENVIRONMENT} ) ecbuild_add_test( TARGET atlas_test_field_acc @@ -40,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..c6ca36f47 --- /dev/null +++ b/src/tests/field/fctest_multifield_ifs.F90 @@ -0,0 +1,77 @@ +! (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 + type(atlas_config) :: config + + integer, parameter :: nproma = 16; + integer, parameter :: nlev = 100; + integer, parameter :: ngptot = 2000; + type(atlas_Config), dimension(5) :: field_configs + + config = atlas_Config() + call config%set("type", "MultiFieldCreatorIFS"); + call config%set("ngptot", ngptot); + call config%set("nproma", nproma); + call config%set("nlev", nlev); + 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") + call config%set("fields", field_configs) + + mfield = atlas_MultiField(config) +END_TEST + +! ----------------------------------------------------------------------------- + +END_TESTSUITE From 39ee651c7770e87fcd509925a35990d9d071105a Mon Sep 17 00:00:00 2001 From: Willem Deconinck Date: Tue, 5 Sep 2023 18:43:10 +0200 Subject: [PATCH 6/7] Use MultiFieldArrayRegistry to avoid accidental deletion of multifield arrays --- src/atlas/field/MultiField.cc | 41 ++++- src/atlas/field/MultiField.h | 49 +----- src/atlas/field/detail/MultiFieldInterface.cc | 12 +- src/tests/field/test_multifield_ifs.cc | 165 +++++++++--------- 4 files changed, 137 insertions(+), 130 deletions(-) diff --git a/src/atlas/field/MultiField.cc b/src/atlas/field/MultiField.cc index 463b3bc65..f7a3ba88d 100644 --- a/src/atlas/field/MultiField.cc +++ b/src/atlas/field/MultiField.cc @@ -14,9 +14,7 @@ #include #include #include - -#include "eckit/thread/AutoLock.h" -#include "eckit/thread/Mutex.h" +#include #include "atlas/field/Field.h" #include "atlas/grid/Grid.h" @@ -41,6 +39,34 @@ void force_link() { //----------------------------------------------------------------------------- +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; @@ -51,7 +77,7 @@ MultiFieldCreator* MultiFieldCreatorFactory::build(const std::string& builder, c return factory->make(config); } -atlas::field::MultiField::MultiField(const eckit::Configuration& config) { +MultiField::MultiField(const eckit::Configuration& config) { std::string type; if (!config.get("type", type)) { ATLAS_THROW_EXCEPTION("Could not find \"type\" in configuration"); @@ -60,6 +86,13 @@ atlas::field::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_); + +} + //----------------------------------------------------------------------------- } // namespace field diff --git a/src/atlas/field/MultiField.h b/src/atlas/field/MultiField.h index 773a098e3..46d0adea0 100644 --- a/src/atlas/field/MultiField.h +++ b/src/atlas/field/MultiField.h @@ -41,51 +41,18 @@ namespace field { * Fields have to all be of same memory layout and data type */ -namespace detail { -class ArrayAllocator { -public: - virtual array::Array* allocate(const array::ArraySpec&) = 0; - virtual void deallocate() = 0; - virtual ~ArrayAllocator() = default; -}; - -class NoArrayAllocator : public ArrayAllocator { -public: - virtual array::Array* allocate(const array::ArraySpec& arrayspec) override { return nullptr; } - virtual void deallocate() override {} -}; - -class StandardArrayAllocator : public ArrayAllocator { -public: - virtual array::Array* allocate(const array::ArraySpec& arrayspec) override { - array::ArraySpec spec(arrayspec); - return array::Array::create(std::move(spec)); - } - virtual void deallocate() override { - // Nothing to deallocate. - // The Array owns the allocated data. - } -}; -} // namespace detail - class MultiFieldImpl : public util::Object { public: // methods //-- Constructors - MultiFieldImpl() { allocator_ = std::unique_ptr(new detail::NoArrayAllocator()); } - - MultiFieldImpl(const array::ArraySpec& spec, detail::ArrayAllocator* allocator = nullptr) { - if (allocator) { - allocator_ = std::unique_ptr(allocator); - } - else { - allocator_ = std::unique_ptr(new detail::StandardArrayAllocator()); - } - ATLAS_ASSERT(allocator_); - array_ = std::unique_ptr(allocator_->allocate(spec)); + MultiFieldImpl() { } + + MultiFieldImpl(const array::ArraySpec& spec) { + array::ArraySpec s(spec); + array_.reset(array::Array::create(std::move(s))); } - virtual ~MultiFieldImpl() { allocator_->deallocate(); } + virtual ~MultiFieldImpl() {} //-- Accessors @@ -132,11 +99,11 @@ class MultiFieldImpl : public util::Object { const FieldSet& fieldset() const { return fieldset_; } FieldSet& fieldset() { return fieldset_; } + void add(Field& field); public: // temporary public for prototyping FieldSet fieldset_; - std::unique_ptr array_; - std::unique_ptr allocator_; + std::shared_ptr array_; util::Metadata metadata_; }; diff --git a/src/atlas/field/detail/MultiFieldInterface.cc b/src/atlas/field/detail/MultiFieldInterface.cc index 685867685..a5e729b2a 100644 --- a/src/atlas/field/detail/MultiFieldInterface.cc +++ b/src/atlas/field/detail/MultiFieldInterface.cc @@ -59,9 +59,8 @@ MultiFieldImpl* atlas__MultiField__create(eckit::Configuration* config) { } auto multiarray_shape = array::make_shape(nblk, nfld, nlev, nproma); - MultiFieldImpl* field = new MultiFieldImpl{array::ArraySpec{datatype, multiarray_shape}}; - auto& multiarray = field->array(); - auto& fieldset = field->fieldset(); + 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) { @@ -111,14 +110,13 @@ MultiFieldImpl* atlas__MultiField__create(eckit::Configuration* config) { } field.set_levels(nlev); // field.set_blocks(nblk); - ATLAS_ASSERT(not fieldset.has(field.name()), "Field with name \"" + field.name() + "\" already exists!"); - fieldset.add(field); + multifield->add(field); multiarray_field_idx += field_vars; } - ATLAS_ASSERT(field); - return field; + ATLAS_ASSERT(multifield); + return multifield; } void atlas__MultiField__delete(MultiFieldImpl* This) { diff --git a/src/tests/field/test_multifield_ifs.cc b/src/tests/field/test_multifield_ifs.cc index 87db13e45..ce11afa35 100644 --- a/src/tests/field/test_multifield_ifs.cc +++ b/src/tests/field/test_multifield_ifs.cc @@ -80,7 +80,6 @@ MultiFieldImpl* MultiFieldCreatorIFS::create(const eckit::Configuration& config) MultiFieldImpl* multifield = new MultiFieldImpl{array::ArraySpec{datatype, multiarray_shape}}; auto& multiarray = multifield->array(); - auto& fieldset = multifield->fieldset(); size_t multiarray_field_idx = 0; for (size_t i = 0; i < fields.size(); ++i) { @@ -129,9 +128,8 @@ MultiFieldImpl* MultiFieldCreatorIFS::create(const eckit::Configuration& config) } } field.set_levels(nlev); - ATLAS_ASSERT(not fieldset.has(field.name()), "Field with name \"" + field.name() + "\" already exists!"); - fieldset.add(field); + multifield->add(field); multiarray_field_idx += field_vars; } @@ -175,81 +173,92 @@ CASE("multifield_create") { return p.json(); }; - Log::info() << "json = " << json() << std::endl; - - 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); + 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); + } + } + + SECTION("test registry") { + { + Field field = MultiField {eckit::YAMLConfiguration{json()}}.field("temperature"); + auto temp = array::make_view(field); + } } } From 83cd8cdd15fe54e19dae5a83b381b61d3fe6a526 Mon Sep 17 00:00:00 2001 From: Slavko Brdar Date: Wed, 6 Sep 2023 19:47:58 +0200 Subject: [PATCH 7/7] Improve and clean up MultiField --- src/atlas/CMakeLists.txt | 8 + src/atlas/field/MultiField.cc | 95 ++--- src/atlas/field/MultiField.h | 167 +++----- src/atlas/field/MultiFieldCreator.cc | 58 +++ src/atlas/field/MultiFieldCreator.h | 89 ++++ src/atlas/field/MultiFieldCreatorArray.cc | 150 +++++++ src/atlas/field/MultiFieldCreatorArray.h | 58 +++ src/atlas/field/MultiFieldCreatorIFS.cc | 230 +++++++++++ src/atlas/field/MultiFieldCreatorIFS.h | 63 +++ src/atlas/field/detail/MultiFieldImpl.cc | 32 ++ src/atlas/field/detail/MultiFieldImpl.h | 96 +++++ src/atlas/field/detail/MultiFieldInterface.cc | 125 ++---- src/atlas/field/detail/MultiFieldInterface.h | 5 + src/atlas_f/CMakeLists.txt | 2 +- ...odule.fypp => atlas_MultiField_module.F90} | 70 +++- src/tests/field/fctest_multifield_ifs.F90 | 187 ++++++++- src/tests/field/test_multifield_ifs.cc | 387 ++++++++++++------ 17 files changed, 1412 insertions(+), 410 deletions(-) create mode 100644 src/atlas/field/MultiFieldCreator.cc create mode 100644 src/atlas/field/MultiFieldCreator.h create mode 100644 src/atlas/field/MultiFieldCreatorArray.cc create mode 100644 src/atlas/field/MultiFieldCreatorArray.h create mode 100644 src/atlas/field/MultiFieldCreatorIFS.cc create mode 100644 src/atlas/field/MultiFieldCreatorIFS.h create mode 100644 src/atlas/field/detail/MultiFieldImpl.cc create mode 100644 src/atlas/field/detail/MultiFieldImpl.h rename src/atlas_f/field/{atlas_MultiField_module.fypp => atlas_MultiField_module.F90} (52%) 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