diff --git a/src/orca-jedi/CMakeLists.txt b/src/orca-jedi/CMakeLists.txt index d4aad4e..83362ad 100644 --- a/src/orca-jedi/CMakeLists.txt +++ b/src/orca-jedi/CMakeLists.txt @@ -21,8 +21,6 @@ model/ModelBiasIncrement.h model/Model.h state/State.h state/State.cc -state/StateIOUtils.h -state/StateIOUtils.cc nemo_io/NemoFieldReader.h nemo_io/NemoFieldReader.cc nemo_io/NemoFieldWriter.h @@ -35,6 +33,8 @@ nemo_io/WriteServer.cc utilities/OrcaModelTraits.h utilities/Types.h utilities/Types.cc +utilities/IOUtils.h +utilities/IOUtils.cc variablechanges/VariableChangeParameters.h variablechanges/VariableChange.h ) diff --git a/src/orca-jedi/increment/Increment.cc b/src/orca-jedi/increment/Increment.cc index 5020103..d4a1327 100644 --- a/src/orca-jedi/increment/Increment.cc +++ b/src/orca-jedi/increment/Increment.cc @@ -28,10 +28,11 @@ #include "ufo/GeoVaLs.h" +#include "orca-jedi/utilities/IOUtils.h" #include "orca-jedi/geometry/Geometry.h" #include "orca-jedi/state/State.h" -#include "orca-jedi/state/StateIOUtils.h" #include "orca-jedi/increment/Increment.h" +#include "orca-jedi/increment/IncrementParameters.h" #include "atlas/mesh.h" #include "atlas-orca/grid/OrcaGrid.h" @@ -397,11 +398,16 @@ void Increment::random() { } /// \brief Apply Dirac delta functions to configuration specified points. -void Increment::dirac(const eckit::Configuration & conf) { +void Increment::dirac(const eckit::Configuration & config) { + dirac(oops::validateAndDeserialize(config)); +} + +/// \brief Apply Dirac delta functions to params specified points. +void Increment::dirac(const OrcaDiracParameters & params) { // Adding a delta function at points specified by ixdir, iydir, izdir - const std::vector & ixdir = conf.getIntVector("ixdir"); - const std::vector & iydir = conf.getIntVector("iydir"); - const std::vector & izdir = conf.getIntVector("izdir"); + const std::vector & ixdir = params.ixdir; + const std::vector & iydir = params.iydir; + const std::vector & izdir = params.izdir; ASSERT(ixdir.size() == iydir.size() && ixdir.size() == izdir.size()); int ndir = ixdir.size(); @@ -526,10 +532,25 @@ void Increment::read(const eckit::Configuration & conf) { throw eckit::NotImplemented(err_message, Here()); } -void Increment::write(const eckit::Configuration & conf) const { - std::string err_message = - "orcamodel::Increment::write not implemented"; - throw eckit::NotImplemented(err_message, Here()); +/// \brief write out increments fields to a file using params specified filename. +void Increment::write(const OrcaIncrementParameters & params) const { + oops::Log::debug() << "orcamodel::increment::write" << std::endl; + + std::string output_filename = params.nemoFieldFile.value(); + if (output_filename == "") + throw eckit::BadValue(std::string("orcamodel::writeIncrementFieldsToFile:: ") + + "file name not specified", Here()); + + auto nemo_field_path = eckit::PathName(output_filename); + oops::Log::debug() << "Increment::write to filename " + << nemo_field_path << std::endl; + + writeFieldsToFile(nemo_field_path, *geom_, time_, incrementFields_); +} + +/// \brief write out increments fields to a file using config specified filename. +void Increment::write(const eckit::Configuration & config) const { + write(oops::validateAndDeserialize(config)); } /// \brief Print some basic information about the self increment object. diff --git a/src/orca-jedi/increment/Increment.h b/src/orca-jedi/increment/Increment.h index f9143f3..29bb98e 100644 --- a/src/orca-jedi/increment/Increment.h +++ b/src/orca-jedi/increment/Increment.h @@ -21,6 +21,7 @@ #include "orca-jedi/geometry/Geometry.h" #include "orca-jedi/state/State.h" +#include "orca-jedi/increment/IncrementParameters.h" #include "eckit/exception/Exceptions.h" @@ -78,6 +79,7 @@ class Increment : public util::Printable, double dot_product_with(const Increment &) const; void schur_product_with(const Increment &); void random(); + void dirac(const OrcaDiracParameters &); void dirac(const eckit::Configuration &); /// ATLAS @@ -96,7 +98,9 @@ class Increment : public util::Printable, }; void read(const eckit::Configuration &); + void write(const OrcaIncrementParameters &) const; void write(const eckit::Configuration &) const; + void print(std::ostream & os) const override; struct stats stats(const std::string & field_name) const; double norm() const; diff --git a/src/orca-jedi/increment/IncrementParameters.h b/src/orca-jedi/increment/IncrementParameters.h new file mode 100644 index 0000000..3084e05 --- /dev/null +++ b/src/orca-jedi/increment/IncrementParameters.h @@ -0,0 +1,35 @@ +/* + * (C) British Crown Copyright 2024 Met Office + */ + +#pragma once + +#include +#include + +#include "eckit/exception/Exceptions.h" +#include "oops/util/DateTime.h" +#include "oops/base/Variables.h" +#include "oops/base/ParameterTraitsVariables.h" +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/RequiredParameter.h" + +namespace orcamodel { + +class OrcaIncrementParameters : public oops::Parameters { + OOPS_CONCRETE_PARAMETERS(OrcaIncrementParameters, oops::Parameters) + + public: + oops::RequiredParameter nemoFieldFile{"output path", this}; + oops::RequiredParameter date{"date", this}; +}; + +class OrcaDiracParameters : public oops::Parameters { + OOPS_CONCRETE_PARAMETERS(OrcaDiracParameters, oops::Parameters) + + public: + oops::RequiredParameter> ixdir{"x indices", this}; + oops::RequiredParameter> iydir{"y indices", this}; + oops::RequiredParameter> izdir{"z indices", this}; +}; +} // namespace orcamodel diff --git a/src/orca-jedi/state/State.cc b/src/orca-jedi/state/State.cc index b7d53b8..9334079 100644 --- a/src/orca-jedi/state/State.cc +++ b/src/orca-jedi/state/State.cc @@ -39,7 +39,7 @@ #include "orca-jedi/variablechanges/VariableChange.h" #include "orca-jedi/increment/Increment.h" #include "orca-jedi/state/State.h" -#include "orca-jedi/state/StateIOUtils.h" +#include "orca-jedi/utilities/IOUtils.h" #include "orca-jedi/utilities/Types.h" @@ -78,9 +78,12 @@ State::State(const Geometry & geom, if (params_.analyticInit.value().value_or(false)) { this->analytic_init(*geom_); } else { - readFieldsFromFile(params_, *geom_, validTime(), "background", - stateFields_); - readFieldsFromFile(params_, *geom_, validTime(), "background variance", + std::string nemo_file_name; + nemo_file_name = params.nemoFieldFile.value(); + readFieldsFromFile(nemo_file_name, *geom_, validTime(), "background", + stateFields_); + nemo_file_name = params.varianceFieldFile.value().value_or(""); + readFieldsFromFile(nemo_file_name, *geom_, validTime(), "background variance", stateFields_); } geom_->log_status(); @@ -185,7 +188,8 @@ void State::read(const OrcaStateParameters & params) { throw eckit::UserError(msg.str(), Here()); } - readFieldsFromFile(params, *geom_, validTime(), "background", + std::string nemo_file_name = params.nemoFieldFile.value(); + readFieldsFromFile(nemo_file_name, *geom_, validTime(), "background", stateFields_); oops::Log::trace() << "State(ORCA)::read done" << std::endl; } @@ -225,7 +229,8 @@ void State::setupStateFields() { void State::write(const OrcaStateParameters & params) const { oops::Log::trace() << "State(ORCA)::write starting" << std::endl; - writeFieldsToFile(params, *geom_, validTime(), stateFields_); + writeFieldsToFile(params.outputNemoFieldFile.value().value_or(""), *geom_, validTime(), + stateFields_); } void State::write(const eckit::Configuration & config) const { diff --git a/src/orca-jedi/state/StateIOUtils.cc b/src/orca-jedi/utilities/IOUtils.cc similarity index 78% rename from src/orca-jedi/state/StateIOUtils.cc rename to src/orca-jedi/utilities/IOUtils.cc index 6498e01..995c614 100644 --- a/src/orca-jedi/state/StateIOUtils.cc +++ b/src/orca-jedi/utilities/IOUtils.cc @@ -2,7 +2,7 @@ * (C) British Crown Copyright 2024 Met Office */ -#include "orca-jedi/state/StateIOUtils.h" +#include "orca-jedi/utilities/IOUtils.h" #include #include @@ -15,15 +15,21 @@ #include "atlas/field/MissingValue.h" #include "orca-jedi/geometry/Geometry.h" -#include "orca-jedi/state/State.h" #include "orca-jedi/nemo_io/WriteServer.h" +#include "orca-jedi/utilities/Types.h" #define NEMO_FILL_TOL 1e-6 namespace orcamodel { +/// \brief Read atlas fields to a file. +/// \param nemo_file_name The name of the netCDF file to read. +/// \param geom The orcamodel::Geometry of the model data. +/// \param valid_date The datetime of data to write to the file. +/// \param variable_type The type of variable (e.g "background" for background data). +/// \param fs The atlas FieldSet collection of fields to write. void readFieldsFromFile( - const OrcaStateParameters & params, + const std::string & nemo_file_name, const Geometry & geom, const util::DateTime & valid_date, const std::string & variable_type, @@ -31,15 +37,6 @@ void readFieldsFromFile( oops::Log::trace() << "orcamodel::readFieldsFromFile:: start for valid_date" << " " << valid_date << std::endl; - // Open Nemo Feedback file - std::string nemo_file_name; - if (variable_type == "background") { - nemo_file_name = params.nemoFieldFile.value(); - } - if (variable_type == "background variance") { - nemo_file_name = params.varianceFieldFile.value().value_or(""); - } - auto nemo_field_path = eckit::PathName(nemo_file_name); oops::Log::debug() << "orcamodel::readFieldsFromFile:: nemo_field_path " << nemo_field_path << std::endl; @@ -77,7 +74,7 @@ void readFieldsFromFile( }; ApplyForFieldType(populate, geom.fieldPrecision(fieldName), - std::string("State(ORCA)::readFieldsFromFile '") + std::string("orcamodel::readFieldsFromFile '") + nemoName + "' field type not recognised"); // Add a halo exchange following read to fill out halo points geom.functionSpace().haloExchange(field); @@ -125,21 +122,29 @@ template void populateField( ReadServer & nemo_reader, atlas::Field & field); +/// \brief write out atlas fields to a file. +/// \param output_file_name The name of the netCDF file. +/// \param geom The orcamodel::Geometry of the model data. +/// \param valid_date The datetime of data to write to the file. +/// \param fs The atlas FieldSet collection of fields to write. void writeFieldsToFile( - const OrcaStateParameters & params, + const std::string & output_file_name, const Geometry & geom, const util::DateTime & valid_date, const atlas::FieldSet & fs) { - oops::Log::trace() << "orcamodel::writeFieldsToFile:: start for valid_date" + oops::Log::trace() << "orcamodel::writeStateFieldsToFile:: start for valid_date" << " " << valid_date << std::endl; - std::string output_filename = - params.outputNemoFieldFile.value().value_or(""); - if (output_filename == "") - throw eckit::BadValue(std::string("orcamodel::writeFieldsToFile:: ") + if (output_file_name == "") + throw eckit::BadValue(std::string("orcamodel::writeStateFieldsToFile:: ") + "file name not specified", Here()); + auto nemo_field_path = eckit::PathName(output_file_name); + oops::Log::debug() << "orcamodel::writeStateFieldsToFile:: " + << nemo_field_path << std::endl; + std::map varCoordTypeMap; + { const oops::Variables vars = geom.variables(); const std::vector coordSpaces = @@ -148,7 +153,6 @@ void writeFieldsToFile( varCoordTypeMap[vars[i].name()] = coordSpaces[i]; } - auto nemo_field_path = eckit::PathName(output_filename); oops::Log::debug() << "orcamodel::writeFieldsToFile:: nemo_field_path " << nemo_field_path << std::endl; std::vector datetimes = {valid_date}; @@ -170,10 +174,13 @@ void writeFieldsToFile( writer.write_vol_var(nemoName, 0, field_mv, field_view); } }; + + oops::Log::debug() << "field " << fieldName + << " data type " << field.datatype().str() << std::endl; ApplyForFieldType(write, - geom.fieldPrecision(fieldName), - std::string("State(ORCA)::writeFieldsToFile '") - + nemoName + "' field type not recognised."); + field.datatype(), + std::string("orcamodel::writeFieldsToFile '") + + nemoName + "' field data type not recognised."); } } diff --git a/src/orca-jedi/state/StateIOUtils.h b/src/orca-jedi/utilities/IOUtils.h similarity index 78% rename from src/orca-jedi/state/StateIOUtils.h rename to src/orca-jedi/utilities/IOUtils.h index ab439d7..aeaacb0 100644 --- a/src/orca-jedi/state/StateIOUtils.h +++ b/src/orca-jedi/utilities/IOUtils.h @@ -10,21 +10,18 @@ #include "atlas/field.h" #include "orca-jedi/geometry/Geometry.h" -#include "orca-jedi/state/State.h" -#include "orca-jedi/state/StateParameters.h" -#include "orca-jedi/utilities/Types.h" #include "orca-jedi/nemo_io/ReadServer.h" namespace orcamodel { void readFieldsFromFile( - const OrcaStateParameters & params, + const std::string & nemo_file_name, const Geometry & geom, const util::DateTime & valid_date, const std::string & variable_type, atlas::FieldSet & fs); void writeFieldsToFile( - const OrcaStateParameters & params, + const std::string & output_file_name, const Geometry & geom, const util::DateTime & valid_date, const atlas::FieldSet & fs); diff --git a/src/orca-jedi/utilities/Types.h b/src/orca-jedi/utilities/Types.h index 7ac80f8..15ecb2d 100644 --- a/src/orca-jedi/utilities/Types.h +++ b/src/orca-jedi/utilities/Types.h @@ -10,6 +10,8 @@ #include "eckit/exception/Exceptions.h" +#include "atlas/array/DataType.h" + namespace orcamodel { //// \brief Enum type for obs variable data types @@ -31,6 +33,19 @@ void ApplyForFieldType(const Functor& functor, FieldDType field_type, } } +/// \brief Apply a function for a given Atlas DataType +template +void ApplyForFieldType(const Functor& functor, atlas::DataType datatype, + const std::string& error_message) { + if (datatype.str() == "real32") { + functor(float{}); + } else if (datatype.str() == "real64") { + functor(double{}); + } else { + throw eckit::BadParameter(error_message); + } +} + // create a mapping between C++ types and NetCDF type objects template struct NetCDFTypeMap { diff --git a/src/tests/orca-jedi/test_increment.cc b/src/tests/orca-jedi/test_increment.cc index f556f0b..e44f474 100644 --- a/src/tests/orca-jedi/test_increment.cc +++ b/src/tests/orca-jedi/test_increment.cc @@ -16,6 +16,7 @@ #include "atlas/library/Library.h" #include "orca-jedi/increment/Increment.h" +#include "orca-jedi/utilities/IOUtils.h" #include "tests/orca-jedi/OrcaModelTestEnvironment.h" @@ -49,6 +50,11 @@ CASE("test increment") { config.set("number levels", 10); Geometry geometry(config, eckit::mpi::comm()); + eckit::LocalConfiguration inc_config; + inc_config.set("date", "2021-06-30T00:00:00Z"); + + OrcaIncrementParameters params; + oops::Variables oops_vars2{{oops::Variable{"sea_ice_area_fraction"}, oops::Variable{"sea_water_potential_temperature"}}}; @@ -56,6 +62,15 @@ CASE("test increment") { util::DateTime datetime("2021-06-30T00:00:00Z"); + SECTION("test increment parameters") { + inc_config.set("output path", "../Data/orca2_t_nemo.nc"); + params.validateAndDeserialize(inc_config); + EXPECT(params.nemoFieldFile.value() == + inc_config.getString("output path")); + auto datetime = static_cast(inc_config.getString("date")); + EXPECT(params.date.value() == datetime); + } + SECTION("test constructor") { Increment increment(geometry, oops_vars2, datetime); // copy @@ -82,14 +97,14 @@ CASE("test increment") { std::vector ix = {20, 30}; std::vector iy = {10, 40}; std::vector iz = {1, 3}; - dirac_config.set("ixdir", ix); - dirac_config.set("iydir", iy); - dirac_config.set("izdir", iz); + dirac_config.set("x indices", ix); + dirac_config.set("y indices", iy); + dirac_config.set("z indices", iz); Increment increment(geometry, oops_vars, datetime); EXPECT_THROWS_AS(increment.dirac(dirac_config), eckit::BadValue); - dirac_config.set("izdir", std::vector{0, 0}); + dirac_config.set("z indices", std::vector{0, 0}); increment.dirac(dirac_config); increment.print(std::cout); EXPECT(std::abs(increment.norm() - 0.0086788) < 1e-6); @@ -155,6 +170,23 @@ CASE("test increment") { std::cout << "increment (diff state1 state2):" << std::endl; increment.print(std::cout); } + + SECTION("test increment write") { + Increment increment1(geometry, oops_vars, datetime); + increment1.ones(); + inc_config.set("output path", "../testoutput/orca2_t_increment_output.nc"); + params.validateAndDeserialize(inc_config); + increment1.write(params); + } + + SECTION("test orca atlas fieldset write") { + Increment increment1(geometry, oops_vars, datetime); + increment1.ones(); + increment1 *=2; + atlas::FieldSet incfset = atlas::FieldSet(); + increment1.Increment::toFieldSet(incfset); + writeFieldsToFile("../testoutput/orca2_t_inc_fs_output.nc", geometry, datetime, incfset); + } } } // namespace test