From 005578854f182e1e5ef6bda2e5719ba8956865fc Mon Sep 17 00:00:00 2001 From: frld <37186816+frld@users.noreply.github.com> Date: Thu, 10 Oct 2024 11:20:52 +0100 Subject: [PATCH] Add 3DVar and error covariance toolbox (dirac test) main programs. (#106) * Add 3DVar and error covariance toolbox (dirac test) main programs. * Add interpolator increment apply and adjoint methods and add increment to state method. * Update ci container (#117) * Add reference data for 3dvar and dirac application tests. Turn off creation of json schema again for 3dvar and errorcovariancetoolbox applications. Co-authored-by: Toby Searle <14909402+twsearle@users.noreply.github.com> --- .github/workflows/ci.yml | 17 ++- CMakeLists.txt | 3 + VERSION | 2 +- ci/CMakeLists.txt | 2 + ci/hpccm_recipe_almalinux9.py | 8 ++ src/mains/CMakeLists.txt | 21 +++ src/mains/orcamodel3DVar.cc | 36 +++++ src/mains/orcamodelErrorCovarianceToolbox.cc | 25 ++++ src/orca-jedi/CMakeLists.txt | 7 +- src/orca-jedi/increment/Increment.cc | 25 +++- src/orca-jedi/interpolator/Interpolator.cc | 123 ++++++++++++++++++ src/orca-jedi/interpolator/Interpolator.h | 12 +- .../interpolator/InterpolatorParameters.h | 1 + src/orca-jedi/state/State.cc | 73 ++++++++++- src/orca-jedi/state/State.h | 3 + src/orca-jedi/utilities/ModelData.cc | 35 +++++ src/orca-jedi/utilities/ModelData.h | 40 ++++++ src/orca-jedi/utilities/OrcaModelTraits.h | 4 + .../variablechanges/LinearVariableChange.cc | 42 ++++++ .../variablechanges/LinearVariableChange.h | 55 ++++++++ .../LinearVariableChangeParameters.h | 28 ++++ src/tests/CMakeLists.txt | 10 ++ src/tests/Data/CMakeLists.txt | 1 + src/tests/Data/sic_obs_ideal.nc | 3 + src/tests/orca-jedi/test_interpolator.cc | 43 ++++++ src/tests/orca-jedi/test_state.cc | 5 + src/tests/testinput/3dvar_sic.yaml | 106 +++++++++++++++ src/tests/testinput/CMakeLists.txt | 2 + src/tests/testinput/dirac.yaml | 41 ++++++ src/tests/testoutput/CMakeLists.txt | 2 + src/tests/testoutput/test_3dvar_sic.ref | 16 +++ src/tests/testoutput/test_dirac.ref | 10 ++ 32 files changed, 777 insertions(+), 24 deletions(-) create mode 100644 src/mains/orcamodel3DVar.cc create mode 100644 src/mains/orcamodelErrorCovarianceToolbox.cc create mode 100644 src/orca-jedi/utilities/ModelData.cc create mode 100644 src/orca-jedi/utilities/ModelData.h create mode 100644 src/orca-jedi/variablechanges/LinearVariableChange.cc create mode 100644 src/orca-jedi/variablechanges/LinearVariableChange.h create mode 100644 src/orca-jedi/variablechanges/LinearVariableChangeParameters.h create mode 100644 src/tests/Data/sic_obs_ideal.nc create mode 100644 src/tests/testinput/3dvar_sic.yaml create mode 100644 src/tests/testinput/dirac.yaml create mode 100644 src/tests/testoutput/test_3dvar_sic.ref create mode 100644 src/tests/testoutput/test_dirac.ref diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1571097..e261dfc 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -5,9 +5,10 @@ on: branches: [develop] pull_request: branches: [develop] + workflow_dispatch: env: REGISTRY: ghcr.io - IMAGE_NAME: twsearle/orca-jedi/ci-almalinux9:v1.3.0 + IMAGE_NAME: twsearle/orca-jedi/ci-almalinux9:v1.4.0 jobs: build: runs-on: ubuntu-latest @@ -32,6 +33,20 @@ jobs: repository: JCSDA-internal/oops token: ${{ secrets.GHCR_PAT }} + - name: checkout vader + uses: actions/checkout@v3 + with: + path: ci/vader + repository: JCSDA-internal/vader + token: ${{ secrets.GHCR_PAT }} + + - name: checkout saber + uses: actions/checkout@v3 + with: + path: ci/saber + repository: JCSDA-internal/saber + token: ${{ secrets.GHCR_PAT }} + - name: checkout ioda uses: actions/checkout@v3 with: diff --git a/CMakeLists.txt b/CMakeLists.txt index 381d177..335d308 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -36,6 +36,9 @@ ecbuild_debug( " atlas_orca_FEATURES: [${atlas_orca_FEATURES}]" ) find_package( oops REQUIRED ) ecbuild_debug( " oops_FEATURES : [${oops_FEATURES}]" ) +find_package( saber REQUIRED) +ecbuild_debug( " saber_FEATURES : [${saber_FEATURES}]" ) + find_package( ufo REQUIRED) ecbuild_debug( " ufo_FEATURES : [${ufo_FEATURES}]" ) diff --git a/VERSION b/VERSION index f0bb29e..88c5fb8 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -1.3.0 +1.4.0 diff --git a/ci/CMakeLists.txt b/ci/CMakeLists.txt index 1a76f25..9118472 100644 --- a/ci/CMakeLists.txt +++ b/ci/CMakeLists.txt @@ -18,6 +18,8 @@ if(NOT DEFINED jedicmake_DIR) endif() add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/oops" EXCLUDE_FROM_ALL) +add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/vader" EXCLUDE_FROM_ALL) +add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/saber" EXCLUDE_FROM_ALL) add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/ioda" EXCLUDE_FROM_ALL) add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/ufo" EXCLUDE_FROM_ALL) add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/orca-jedi") diff --git a/ci/hpccm_recipe_almalinux9.py b/ci/hpccm_recipe_almalinux9.py index c9f0e29..18925ca 100644 --- a/ci/hpccm_recipe_almalinux9.py +++ b/ci/hpccm_recipe_almalinux9.py @@ -50,6 +50,7 @@ def gitlab_url(repo, vn): odc_vn = USERARG.get('odc_vn', '1.5.2') openmpi_vn = USERARG.get('openmpi_vn', '4.1.5') pycodestyle_vn = USERARG.get('pycodestyle_vn', '2.10') +qhull_vn = USERARG.get('qhull_vn', '8.0.2') # no qhull-devel rpm udunits_vn = USERARG.get('udunits_vn', '2.2.28') yaxt_vn = USERARG.get('yaxt_vn', '528-0.10.0') # URL has a number and a version @@ -203,6 +204,13 @@ def gitlab_url(repo, vn): cmake_opts=['-DCMAKE_BUILD_TYPE=Release', '-DJSON_BuildTests=OFF'], ) +Stage0 += generic_cmake( + prefix='/usr/local', + url=github_url('qhull/qhull', f'v{qhull_vn}'), + directory=f'qhull-{qhull_vn}', + cmake_opts=['-DCMAKE_BUILD_TYPE=Release'], +) + Stage0 += generic_cmake( prefix='/usr/local', url=github_url('NOAA-EMC/NCEPLIBS-bufr', f'v{nceplibs_bufr_vn}'), diff --git a/src/mains/CMakeLists.txt b/src/mains/CMakeLists.txt index a20224e..a233ad2 100644 --- a/src/mains/CMakeLists.txt +++ b/src/mains/CMakeLists.txt @@ -11,6 +11,12 @@ if( ${nemo-feedback_FOUND} ) LIBS orcamodel nemo_feedback DEFINITIONS NEMO_FEEDBACK_EXISTS ) + + ecbuild_add_executable( TARGET orcamodel_3DVar.x + SOURCES orcamodel3DVar.cc + LIBS orcamodel nemo_feedback + DEFINITIONS NEMO_FEEDBACK_EXISTS + ) else() message(" compiling without feedback file support") ecbuild_add_executable( TARGET orcamodel_hofx.x @@ -22,6 +28,21 @@ else() SOURCES orcamodelHofX3D.cc LIBS orcamodel ) + + ecbuild_add_executable( TARGET orcamodel_3DVar.x + SOURCES orcamodel3DVar.cc + LIBS orcamodel + ) endif() + +ecbuild_add_executable( TARGET orcamodel_ErrorCovarianceToolbox.x + SOURCES orcamodelErrorCovarianceToolbox.cc + LIBS orcamodel + ) + oops_output_json_schema( "orcamodel_hofx.x" ) oops_output_json_schema( "orcamodel_hofx3D.x" ) +# The applications below throw an error generating json schema. +# This appears to be a generic issue see also oops/qg/mains/CMakeLists.txt +#oops_output_json_schema( "orcamodel_3DVar.x" ) +#oops_output_json_schema( "orcamodel_ErrorCovarianceToolbox.x" ) diff --git a/src/mains/orcamodel3DVar.cc b/src/mains/orcamodel3DVar.cc new file mode 100644 index 0000000..9249cc6 --- /dev/null +++ b/src/mains/orcamodel3DVar.cc @@ -0,0 +1,36 @@ +/* + * (C) British Crown Copyright 2024 MetOffice + * + * 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. + */ + +#include "atlas/library/Library.h" +#include "oops/generic/instantiateModelFactory.h" + +#include "oops/runs/Run.h" +#include "oops/runs/Variational.h" +#include "saber/oops/instantiateCovarFactory.h" + +#include "ufo/ObsTraits.h" +#include "ufo/instantiateObsFilterFactory.h" +#if defined(NEMO_FEEDBACK_EXISTS) + #include "nemo-feedback/instantiateObsFilterFactory.h" +#endif + +#include "orca-jedi/utilities/OrcaModelTraits.h" + +int main(int argc, char ** argv) { + oops::Run run(argc, argv); + oops::instantiateModelFactory(); + atlas::Library::instance().initialise(); + saber::instantiateCovarFactory(); + ufo::instantiateObsFilterFactory(); +#if defined(NEMO_FEEDBACK_EXISTS) + nemo_feedback::instantiateObsFilterFactory(); +#endif + oops::Variational var; + int i = run.execute(var); + atlas::Library::instance().finalise(); + return i; +} diff --git a/src/mains/orcamodelErrorCovarianceToolbox.cc b/src/mains/orcamodelErrorCovarianceToolbox.cc new file mode 100644 index 0000000..f46e972 --- /dev/null +++ b/src/mains/orcamodelErrorCovarianceToolbox.cc @@ -0,0 +1,25 @@ +/* + * (C) British Crown Copyright 2024 MetOffice + * + * 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. + */ + +#include "atlas/library/Library.h" + +#include "saber/oops/ErrorCovarianceToolbox.h" +#include "saber/oops/instantiateCovarFactory.h" + +#include "oops/runs/Run.h" +#include "orca-jedi/utilities/OrcaModelTraits.h" + +int main(int argc, char ** argv) { + oops::Run run(argc, argv); + atlas::Library::instance().initialise(); + saber::instantiateCovarFactory(); + + saber::ErrorCovarianceToolbox var; + int i = run.execute(var); + atlas::Library::instance().finalise(); + return i; +} diff --git a/src/orca-jedi/CMakeLists.txt b/src/orca-jedi/CMakeLists.txt index 83362ad..02b36f7 100644 --- a/src/orca-jedi/CMakeLists.txt +++ b/src/orca-jedi/CMakeLists.txt @@ -35,15 +35,20 @@ utilities/Types.h utilities/Types.cc utilities/IOUtils.h utilities/IOUtils.cc +utilities/ModelData.h +utilities/ModelData.cc variablechanges/VariableChangeParameters.h variablechanges/VariableChange.h +variablechanges/LinearVariableChange.h +variablechanges/LinearVariableChange.cc +variablechanges/LinearVariableChangeParameters.h ) ecbuild_add_library( TARGET orcamodel SOURCES ${oops_orcamodel_src_files} - PUBLIC_LIBS oops ufo atlas atlas-orca eckit netcdf + PUBLIC_LIBS oops saber ufo atlas atlas-orca eckit netcdf INSTALL_HEADERS LISTED LINKER_LANGUAGE ${OOPS_LINKER_LANGUAGE} ) diff --git a/src/orca-jedi/increment/Increment.cc b/src/orca-jedi/increment/Increment.cc index d4a1327..a42b65d 100644 --- a/src/orca-jedi/increment/Increment.cc +++ b/src/orca-jedi/increment/Increment.cc @@ -207,7 +207,7 @@ Increment & Increment::operator*=(const double & zz) { return *this; } -/// \brief Create increment from the different of two state objects. +/// \brief Create increment from the difference of two state objects. /// \param x1 State object. /// \param x2 State object subtracted. void Increment::diff(const State & x1, const State & x2) { @@ -222,6 +222,12 @@ void Increment::diff(const State & x1, const State & x2) { atlas::Field field2 = x2.getField(i); atlas::Field fieldi = incrementFields_[i]; + atlas::field::MissingValue mv1(field1); + bool has_mv1 = static_cast(mv1); + atlas::field::MissingValue mv2(field2); + bool has_mv2 = static_cast(mv2); + bool has_mv = has_mv1 || has_mv2; + std::string fieldName1 = field1.name(); std::string fieldName2 = field2.name(); std::string fieldNamei = fieldi.name(); @@ -234,8 +240,15 @@ void Increment::diff(const State & x1, const State & x2) { auto field_viewi = atlas::array::make_view(fieldi); for (atlas::idx_t j = 0; j < field_viewi.shape(0); ++j) { for (atlas::idx_t k = 0; k < field_viewi.shape(1); ++k) { - if (!ghost(j)) { field_viewi(j, k) = field_view1(j, k) - field_view2(j, k); - } else { field_viewi(j, k) = 0; } + if (!ghost(j)) { + if (!has_mv1 || (has_mv1 && !mv1(field_view1(j, k)))) { + if (!has_mv2 || (has_mv2 && !mv2(field_view2(j, k)))) { + field_viewi(j, k) = field_view1(j, k) - field_view2(j, k); + } + } + } else { + field_viewi(j, k) = 0; + } } } } @@ -532,7 +545,7 @@ void Increment::read(const eckit::Configuration & conf) { throw eckit::NotImplemented(err_message, Here()); } -/// \brief write out increments fields to a file using params specified filename. +/// \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; @@ -545,10 +558,12 @@ void Increment::write(const OrcaIncrementParameters & params) const { oops::Log::debug() << "Increment::write to filename " << nemo_field_path << std::endl; +incrementFields_.haloExchange(); + writeFieldsToFile(nemo_field_path, *geom_, time_, incrementFields_); } -/// \brief write out increments fields to a file using config specified filename. +/// \brief Write out increments fields to a file using config specified filename. void Increment::write(const eckit::Configuration & config) const { write(oops::validateAndDeserialize(config)); } diff --git a/src/orca-jedi/interpolator/Interpolator.cc b/src/orca-jedi/interpolator/Interpolator.cc index 796e06c..4bb8251 100644 --- a/src/orca-jedi/interpolator/Interpolator.cc +++ b/src/orca-jedi/interpolator/Interpolator.cc @@ -27,6 +27,7 @@ #include "orca-jedi/geometry/Geometry.h" #include "orca-jedi/state/State.h" +#include "orca-jedi/increment/Increment.h" #include "orca-jedi/interpolator/Interpolator.h" @@ -177,6 +178,128 @@ template void Interpolator::executeInterpolation( const std::vector & mask, std::vector::iterator& iter) const; +/// \brief Interpolate from model space to observation space +/// \param vars Oops variables +/// \param inc Increment object (input) +/// \param mask Mask vector in observation space +/// \param result Result (output) vector in observation space +void Interpolator::apply(const oops::Variables& vars, const Increment& inc, + const std::vector & mask, + std::vector& result) const { + // input is inc output is result + const size_t nvars = vars.size(); + + for (size_t j=0; j < nvars; ++j) { + if (!inc.variables().has(vars[j])) { + std::stringstream err_stream; + err_stream << "orcamodel::Interpolator::apply varname \" " + << "\" " << vars[j] + << " not found in the model increment." << std::endl; + err_stream << " add the variable to the increment variables and " + << "add a mapping from the geometry to that variable." + << std::endl; + throw eckit::BadParameter(err_stream.str(), Here()); + } + } + + const std::vector varSizes = + inc.geometry()->variableSizes(vars); + size_t nvals = 0; + for (size_t jvar=0; jvar < nvars; ++jvar) nvals += nlocs_ * varSizes[jvar]; + result.resize(nvals); + + std::size_t out_idx = 0; + for (size_t jvar=0; jvar < nvars; ++jvar) { + auto gv_varname = vars[jvar].name(); + atlas::Field tgt_field = atlasObsFuncSpace_.createField( + atlas::option::name(gv_varname) | + atlas::option::levels(varSizes[jvar])); + interpolator_.execute(inc.incrementFields()[gv_varname], tgt_field); + auto field_view = atlas::array::make_view(tgt_field); + atlas::field::MissingValue mv(inc.incrementFields()[gv_varname]); + bool has_mv = static_cast(mv); + for (std::size_t klev=0; klev < varSizes[jvar]; ++klev) { + for (std::size_t iloc=0; iloc < nlocs_; iloc++) { + if (has_mv && mv(field_view(iloc, klev))) { + result[out_idx] = util::missingValue(); + } else { + result[out_idx] = field_view(iloc, klev); + } + ++out_idx; + } + } + } +} + +/// \brief Interpolate from observation space to model space +/// \param vars Oops variables +/// \param inc Increment object (output) +/// \param mask Mask (observation space) vector +/// \param resultin Values (observation space) vector (input) +void Interpolator::applyAD(const oops::Variables& vars, Increment& inc, + const std::vector & mask, + const std::vector & resultin) const +{ + // input is resultin output is inc + + oops::Log::trace() << "orcamodel::Interpolator::applyAD start " + << std::endl; + + const size_t nvars = vars.size(); + + for (size_t j=0; j < nvars; ++j) { + if (!inc.variables().has(vars[j])) { + std::stringstream err_stream; + err_stream << "orcamodel::Interpolator::apply varname \" " + << "\" " << vars[j] + << " not found in the model increment." << std::endl; + err_stream << " add the variable to the increment variables and " + << "add a mapping from the geometry to that variable." + << std::endl; + throw eckit::BadParameter(err_stream.str(), Here()); + } + } + + const std::vector varSizes = + inc.geometry()->variableSizes(vars); + size_t nvals = 0; + + for (size_t jvar=0; jvar < nvars; ++jvar) nvals += nlocs_ * varSizes[jvar]; + + std::size_t out_idx = 0; + for (size_t jvar=0; jvar < nvars; ++jvar) { + auto gv_varname = vars[jvar].name(); + auto tgt_field = atlasObsFuncSpace_.createField( + atlas::option::name(gv_varname) | + atlas::option::levels(varSizes[jvar])); + +// Copying observation array vector to an atlas observation field (tgt_field) + auto field_view = atlas::array::make_view(tgt_field); + atlas::field::MissingValue mv(inc.incrementFields()[gv_varname]); + bool has_mv = static_cast(mv); + + for (std::size_t klev=0; klev < varSizes[jvar]; ++klev) { + for (std::size_t iloc=0; iloc < nlocs_; iloc++) { + if (has_mv && mv(field_view(iloc, klev))) { + field_view(iloc, klev) = util::missingValue(); + } else { + field_view(iloc, klev) = resultin[out_idx]; + } + ++out_idx; + } + } + +// halo exchange update ghost points + std::shared_ptr geom = inc.geometry(); + geom->functionSpace().haloExchange(inc.incrementFields()[gv_varname]); + + interpolator_.execute_adjoint(inc.incrementFields()[gv_varname], tgt_field); + } // jvar + + oops::Log::trace() << "orcamodel::Interpolator::applyAD done " + << std::endl; +} + void Interpolator::print(std::ostream & os) const { os << "orcamodel::Interpolator: " << std::endl; os << " Obs function space " << atlasObsFuncSpace_ << std::endl; diff --git a/src/orca-jedi/interpolator/Interpolator.h b/src/orca-jedi/interpolator/Interpolator.h index c1b92c9..4695d28 100644 --- a/src/orca-jedi/interpolator/Interpolator.h +++ b/src/orca-jedi/interpolator/Interpolator.h @@ -58,16 +58,10 @@ class Interpolator : public util::Printable, std::vector& result) const; void apply(const oops::Variables& vars, const Increment& inc, const std::vector & mask, - std::vector& result) const { - throw eckit::NotImplemented("Increment interpolation not implemented", - Here()); - } - void applyAD(const oops::Variables& vars, const Increment& inc, + std::vector& result) const; + void applyAD(const oops::Variables& vars, Increment& inc, const std::vector & mask, - const std::vector &) const { - throw eckit::NotImplemented("Adjoint interpolation not implemented", - Here()); - } + const std::vector &) const; private: template void executeInterpolation( diff --git a/src/orca-jedi/interpolator/InterpolatorParameters.h b/src/orca-jedi/interpolator/InterpolatorParameters.h index 16d752b..833fa35 100644 --- a/src/orca-jedi/interpolator/InterpolatorParameters.h +++ b/src/orca-jedi/interpolator/InterpolatorParameters.h @@ -22,6 +22,7 @@ class OrcaAtlasInterpolatorParameters : public oops::Parameters { oops::OptionalParameter non_linear{"non_linear", this}; oops::OptionalParameter max_fraction_elems_to_try{"max_fraction_elems_to_try", this}; + oops::OptionalParameter adjoint{"adjoint", this}; }; class OrcaInterpolatorParameters : public oops::Parameters { diff --git a/src/orca-jedi/state/State.cc b/src/orca-jedi/state/State.cc index bd359a8..349ac3a 100644 --- a/src/orca-jedi/state/State.cc +++ b/src/orca-jedi/state/State.cc @@ -83,10 +83,12 @@ State::State(const Geometry & geom, readFieldsFromFile(nemo_file_name, *geom_, validTime(), "background", stateFields_); nemo_file_name = params.errorFieldFile.value().value_or(""); - readFieldsFromFile(nemo_file_name, *geom_, validTime(), "background error standard deviation", - stateFields_); - readFieldsFromFile(nemo_file_name, *geom_, validTime(), "background error variance", - stateFields_); + if (params.errorFieldFile.value()) { + readFieldsFromFile(nemo_file_name, *geom_, validTime(), "background error standard deviation", + stateFields_); + readFieldsFromFile(nemo_file_name, *geom_, validTime(), "background error variance", + stateFields_); + } } geom_->log_status(); oops::Log::trace() << "State(ORCA)::State created." << std::endl; @@ -165,11 +167,44 @@ State & State::operator=(const State & rhs) { // Interactions with Increments +/// \brief Add increment to state. +/// \brief Requires increment and state to have the same field names (currently). +/// \param Increment. State & State::operator+=(const Increment & dx) { - std::string err_message = - "orcamodel::State::State::operator+= not implemented"; - throw eckit::NotImplemented(err_message, Here()); oops::Log::trace() << "State(ORCA)::add increment starting" << std::endl; + + ASSERT(this->validTime() == dx.validTime()); + + auto ghost = atlas::array::make_view( + geom_->mesh().nodes().ghost()); + for (int i = 0; i< stateFields_.size(); i++) + { + atlas::Field field = stateFields_[i]; + atlas::field::MissingValue mv(field); + bool has_mv = static_cast(mv); + + atlas::Field fieldi = dx.incrementFields()[i]; + + std::string fieldName = field.name(); + std::string fieldNamei = fieldi.name(); + oops::Log::debug() << "orcamodel::Increment::add:: state field name = " << fieldName + << " increment field name = " << fieldNamei + << std::endl; + ASSERT(fieldName == fieldNamei); + + auto field_view = atlas::array::make_view(field); + auto field_viewi = atlas::array::make_view(fieldi); + for (atlas::idx_t j = 0; j < field_view.shape(0); ++j) { + for (atlas::idx_t k = 0; k < field_view.shape(1); ++k) { + if (!ghost(j)) { + if (!has_mv || (has_mv && !mv(field_view(j, k)))) { + field_view(j, k) += field_viewi(j, k); + } + } + } + } + } + oops::Log::trace() << "State(ORCA)::add increment done" << std::endl; return *this; } @@ -361,4 +396,28 @@ atlas::Field State::getField(int i) const { return stateFields_[i]; } +/// \brief Output state fieldset as an atlas fieldset. +/// \param fset Atlas fieldset to output to. +void State::toFieldSet(atlas::FieldSet & fset) const { + oops::Log::debug() << "State toFieldSet starting" << std::endl; + + fset = atlas::FieldSet(); + + for (size_t i=0; i < vars_.size(); ++i) { + // copy variable from increments to new field set + atlas::Field field = stateFields_[i]; + std::string fieldName = field.name(); + oops::Log::debug() << "Copy state toFieldSet " << fieldName << std::endl; + + fset->add(field); + } + oops::Log::debug() << "State toFieldSet done" << std::endl; +} + +void State::accumul(const double & zz, const State & xx) { + std::string err_message = + "orcamodel::State::accumul not implemented"; + throw eckit::NotImplemented(err_message, Here()); +} + } // namespace orcamodel diff --git a/src/orca-jedi/state/State.h b/src/orca-jedi/state/State.h index 8d9d8ef..f070b49 100644 --- a/src/orca-jedi/state/State.h +++ b/src/orca-jedi/state/State.h @@ -106,6 +106,9 @@ class State : public util::Printable, oops::Variables & variables() {return vars_;} atlas::Field getField(int) const; + void toFieldSet(atlas::FieldSet &) const; + + void accumul(const double &, const State &); private: void setupStateFields(); diff --git a/src/orca-jedi/utilities/ModelData.cc b/src/orca-jedi/utilities/ModelData.cc new file mode 100644 index 0000000..8988580 --- /dev/null +++ b/src/orca-jedi/utilities/ModelData.cc @@ -0,0 +1,35 @@ +/* + * (C) Crown Copyright 2024 Met Office. + * + * 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. + */ + + +#include + +#include "orca-jedi/utilities/ModelData.h" + + +namespace orcamodel { + + +const eckit::LocalConfiguration ModelData::modelData() const { + eckit::LocalConfiguration model_data; + + // + // retrieve data that need to be shared with other system components + // and stored them into the data structure 'model_data' + // + + return model_data; +} + + +void ModelData::print(std::ostream & os) const { + os << "orcamodel::ModelData::ModelData(): " << modelData(); +} + + +} // namespace orcamodel + diff --git a/src/orca-jedi/utilities/ModelData.h b/src/orca-jedi/utilities/ModelData.h new file mode 100644 index 0000000..4c3055d --- /dev/null +++ b/src/orca-jedi/utilities/ModelData.h @@ -0,0 +1,40 @@ +/* + * (C) Crown Copyright 2024 Met Office. + * + * 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. + */ + + +#pragma once + +#include +#include + +#include "eckit/config/LocalConfiguration.h" + +#include "oops/util/Printable.h" + +#include "orca-jedi/geometry/Geometry.h" + + +namespace orcamodel { + +class Geometry; + + +class ModelData : public util::Printable { + public: + explicit ModelData(const Geometry &) {} + ~ModelData() {} + + static const std::string classname() {return "orcamodel::ModelData";} + + const eckit::LocalConfiguration modelData() const; + + private: + void print(std::ostream &) const override; +}; + + +} // namespace orcamodel diff --git a/src/orca-jedi/utilities/OrcaModelTraits.h b/src/orca-jedi/utilities/OrcaModelTraits.h index 182c70d..9a58cc7 100644 --- a/src/orca-jedi/utilities/OrcaModelTraits.h +++ b/src/orca-jedi/utilities/OrcaModelTraits.h @@ -15,7 +15,9 @@ #include "orca-jedi/model/ModelBiasIncrement.h" #include "orca-jedi/model/ModelBiasCovariance.h" #include "orca-jedi/state/State.h" +#include "orca-jedi/utilities/ModelData.h" #include "orca-jedi/variablechanges/VariableChange.h" +#include "orca-jedi/variablechanges/LinearVariableChange.h" namespace orcamodel { @@ -35,6 +37,8 @@ struct OrcaModelTraits { typedef orcamodel::ModelBiasCovariance ModelAuxCovariance; typedef orcamodel::State State; typedef orcamodel::VariableChange VariableChange; + typedef orcamodel::LinearVariableChange LinearVariableChange; + typedef orcamodel::ModelData ModelData; }; } // namespace orcamodel diff --git a/src/orca-jedi/variablechanges/LinearVariableChange.cc b/src/orca-jedi/variablechanges/LinearVariableChange.cc new file mode 100644 index 0000000..3c754c7 --- /dev/null +++ b/src/orca-jedi/variablechanges/LinearVariableChange.cc @@ -0,0 +1,42 @@ +/* + * (C) British Crown Copyright 2024 Met Office + * + * 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. +*/ + +#include + +#include "eckit/exception/Exceptions.h" + +#include "orca-jedi/variablechanges/LinearVariableChange.h" + +#include "oops/util/Printable.h" +#include "orca-jedi/geometry/Geometry.h" +#include "orca-jedi/state/State.h" +#include "orca-jedi/variablechanges/LinearVariableChangeParameters.h" + +namespace orcamodel { + +// ----------------------------------------------------------------------------- +LinearVariableChange::LinearVariableChange(const Geometry &, const Parameters_ &) {} +LinearVariableChange::LinearVariableChange(const Geometry &, const eckit::Configuration &) {} +// ----------------------------------------------------------------------------- +// LinearVariableChange::~LinearVariableChange() {} + +void LinearVariableChange::changeVarTL(Increment & dx, const oops::Variables & vars) +const {} + +void LinearVariableChange::changeVarInverseTL(Increment & dx, const oops::Variables & vars) +const {} + +void LinearVariableChange::changeVarAD(Increment & dx, const oops::Variables & vars) +const {} + +void LinearVariableChange::changeVarInverseAD(Increment & dx, const oops::Variables & vars) +const {} + +void LinearVariableChange::changeVarTraj(const State & xbg, const oops::Variables & xfg) +{} + +} // namespace orcamodel diff --git a/src/orca-jedi/variablechanges/LinearVariableChange.h b/src/orca-jedi/variablechanges/LinearVariableChange.h new file mode 100644 index 0000000..73c4bcd --- /dev/null +++ b/src/orca-jedi/variablechanges/LinearVariableChange.h @@ -0,0 +1,55 @@ +/* + * (C) British Crown Copyright 2024 Met Office + * + * 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. +*/ + +#pragma once + +#include +#include +#include + +#include "eckit/config/Configuration.h" + +#include "oops/util/Printable.h" +#include "orca-jedi/geometry/Geometry.h" +#include "orca-jedi/state/State.h" +#include "orca-jedi/variablechanges/LinearVariableChangeParameters.h" + +// Forward declarations +namespace oops { +class Variables; +} + +namespace orcamodel { + +// ----------------------------------------------------------------------------- + +class LinearVariableChange: public util::Printable { + public: + typedef LinearVariableChangeParameters Parameters_; + static const std::string classname() { + return "orcamodel::LinearVariableChange"; + } + + LinearVariableChange(const Geometry &, const Parameters_ &); + LinearVariableChange(const Geometry &, const eckit::Configuration &); + + void changeVarTL(Increment &, const oops::Variables &) const; + void changeVarInverseTL(Increment &, const oops::Variables &) const; + void changeVarAD(Increment &, const oops::Variables &) const; + void changeVarInverseAD(Increment &, const oops::Variables &) const; + + void changeVarTraj(const State &, const oops::Variables &); + + + private: + void print(std::ostream & out) const override { + out << "orcamodel::LinearVariableChange::print Not Implemented"; }; +}; + +} // namespace orcamodel + + diff --git a/src/orca-jedi/variablechanges/LinearVariableChangeParameters.h b/src/orca-jedi/variablechanges/LinearVariableChangeParameters.h new file mode 100644 index 0000000..8d853ec --- /dev/null +++ b/src/orca-jedi/variablechanges/LinearVariableChangeParameters.h @@ -0,0 +1,28 @@ +/* + * (C) Crown Copyright 2021, the Met Office. All rights reserved. + * + * Refer to COPYRIGHT.txt of this distribution for details. + */ + +#pragma once + +#include +#include + +#include "eckit/exception/Exceptions.h" +#include "oops/util/DateTime.h" +#include "oops/base/Variables.h" +#include "oops/base/LinearVariableChangeParametersBase.h" +#include "oops/util/parameters/Parameter.h" +#include "oops/util/parameters/RequiredParameter.h" +#include "oops/util/parameters/OptionalParameter.h" + + +class LinearVariableChangeParameters : + public oops::LinearVariableChangeParametersBase { + OOPS_CONCRETE_PARAMETERS(LinearVariableChangeParameters, + oops::LinearVariableChangeParametersBase) + public: + // No linear variable change. No additional parameters +}; + diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt index 337d068..fb930b2 100644 --- a/src/tests/CMakeLists.txt +++ b/src/tests/CMakeLists.txt @@ -90,3 +90,13 @@ ecbuild_add_test( TARGET test_orcamodel_hofx3D_prof_2var OMP 1 ARGS testinput/hofx3d_nc_prof_2vars.yaml COMMAND orcamodel_hofx3D.x ) + +ecbuild_add_test( TARGET test_orcamodel_3DVar_sic + OMP 1 + ARGS testinput/3dvar_sic.yaml + COMMAND orcamodel_3DVar.x ) + +ecbuild_add_test( TARGET test_orcamodel_Dirac + OMP 1 + ARGS testinput/dirac.yaml + COMMAND orcamodel_ErrorCovarianceToolbox.x ) diff --git a/src/tests/Data/CMakeLists.txt b/src/tests/Data/CMakeLists.txt index e141869..e1c09ad 100644 --- a/src/tests/Data/CMakeLists.txt +++ b/src/tests/Data/CMakeLists.txt @@ -26,6 +26,7 @@ amm1r_nemo.nc amm1r_atlas_grid_spec.yaml hofx_potm_amm1r_obs.nc hofx_ssh_amm_obs.nc +sic_obs_ideal.nc ) foreach(FILENAME ${orcajedi_test_data}) diff --git a/src/tests/Data/sic_obs_ideal.nc b/src/tests/Data/sic_obs_ideal.nc new file mode 100644 index 0000000..fd58d8c --- /dev/null +++ b/src/tests/Data/sic_obs_ideal.nc @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1f04618a2934a648e4289450471e75f999168b4f6dae663ca45c72baa870a0ba +size 16702 diff --git a/src/tests/orca-jedi/test_interpolator.cc b/src/tests/orca-jedi/test_interpolator.cc index d99c593..d64187b 100644 --- a/src/tests/orca-jedi/test_interpolator.cc +++ b/src/tests/orca-jedi/test_interpolator.cc @@ -34,6 +34,7 @@ struct InterpTestSettingsFixture { public: eckit::LocalConfiguration geometry_config; eckit::LocalConfiguration state_config; + eckit::LocalConfiguration increment_config; eckit::LocalConfiguration interpolator_config; size_t nlocs, nlevs; std::vector lons; @@ -76,6 +77,7 @@ CASE("test interpolator") { eckit::LocalConfiguration interp_conf; interp_conf.set("type", "unstructured-bilinear-lonlat"); interp_conf.set("non_linear", "missing-if-all-missing-real32"); + interp_conf.set("adjoint", true); settings_map["ORCA2_T"].interpolator_config.set("atlas-interpolator", interp_conf); std::vector state_variables { @@ -87,6 +89,10 @@ CASE("test interpolator") { settings_map["ORCA2_T"].state_config.set("nemo field file", "../Data/orca2_t_nemo.nc"); settings_map["ORCA2_T"].state_config.set("nemo error field file", "../Data/orca2_t_bkg_var.nc"); + settings_map["ORCA2_T"].increment_config.set("date", "2021-06-30T00:00:00Z"); + settings_map["ORCA2_T"].increment_config.set("output path", + "../testoutput/orca2_t_increment_interptest.nc"); + settings_map["ORCA2_T"].lons = std::vector{0, 120, 270}; settings_map["ORCA2_T"].lats = std::vector{88, 0, 30}; @@ -126,6 +132,7 @@ CASE("test interpolator") { eckit::LocalConfiguration interp_conf; interp_conf.set("type", "unstructured-bilinear-lonlat"); interp_conf.set("non_linear", "missing-if-all-missing-real32"); + interp_conf.set("adjoint", true); settings_map["AMM1"].interpolator_config.set("atlas-interpolator", interp_conf); std::vector state_variables { @@ -137,6 +144,10 @@ CASE("test interpolator") { settings_map["AMM1"].state_config.set("nemo field file", "../Data/amm1_nemo.nc"); settings_map["AMM1"].state_config.set("nemo error field file", "../Data/amm1_nemo.nc"); + settings_map["AMM1"].increment_config.set("date", "2021-06-30T00:00:00Z"); + settings_map["AMM1"].increment_config.set("output path", + "../testoutput/amm1_t_increment_interptest.nc"); + settings_map["AMM1"].lons = std::vector{-17.5, -6.78, -16.1}; settings_map["AMM1"].lats = std::vector{58.16, 58.91, 63.55}; @@ -195,6 +206,7 @@ CASE("test interpolator") { EXPECT(std::abs(vals[i] - settings.surf_values[i]) < ATOL); } } + SECTION("test " + key + " interpolator.apply multiple levels") { std::vector vals(settings.nlevs*settings.nlocs); std::vector mask(settings.nlocs, true); @@ -209,6 +221,37 @@ CASE("test interpolator") { EXPECT(std::abs(vals[i] - settings.vol_values[i]) < ATOL); } } + + SECTION("test " + key + " interpolator.apply/applyAD with increment") { + OrcaIncrementParameters incrementParams; + incrementParams.validateAndDeserialize(settings.increment_config); + + Increment increment(geometry, settings.surf_vars, incrementParams.date); + increment.ones(); + + // two variables at n locations + std::vector vals(2*settings.nlocs); + std::vector mask(settings.nlocs, true); + + // increment -> observation space + interpolator.apply(settings.surf_vars, increment, mask, vals); + + double kgo_value = 1; + + for (size_t i=0; i < settings.surf_values.size(); ++i) { + std::cout << "vals[" << i << "] " << std::setprecision(12) << vals[i] + << " kgo_values[" << i << "] " << kgo_value << std::endl; + } + for (size_t i=0; i < settings.surf_values.size(); ++i) { + EXPECT(std::abs(vals[i] - kgo_value) < ATOL); + } + + Increment incrementout(geometry, settings.surf_vars, incrementParams.date); + // observation space -> increment + interpolator.applyAD(settings.surf_vars, incrementout, mask, vals); + + incrementout.write(incrementParams); + } } } diff --git a/src/tests/orca-jedi/test_state.cc b/src/tests/orca-jedi/test_state.cc index da65300..3d4b609 100644 --- a/src/tests/orca-jedi/test_state.cc +++ b/src/tests/orca-jedi/test_state.cc @@ -135,6 +135,11 @@ CASE("test basic state") { std::cout << field.name() << std::endl; EXPECT(field.name() == "sea_ice_area_fraction"); } + SECTION("test state to fieldset") { + atlas::FieldSet statefset = atlas::FieldSet(); + state.State::toFieldSet(statefset); + EXPECT(statefset[0].name() == "sea_ice_area_fraction"); + } } } // namespace test diff --git a/src/tests/testinput/3dvar_sic.yaml b/src/tests/testinput/3dvar_sic.yaml new file mode 100644 index 0000000..c84f108 --- /dev/null +++ b/src/tests/testinput/3dvar_sic.yaml @@ -0,0 +1,106 @@ +cost function: + cost type: 3D-Var + time window: + begin: '2021-06-29T00:00:00Z' + length: P1D + analysis variables: [ice_area_fraction] + background error: + covariance model: SABER + saber central block: + saber block name: ID + geometry : + nemo variables: + - name: ice_area_fraction + nemo field name: iiceconc + model space: surface + field precision: double + - name: sea_surface_temperature + nemo field name: votemper + model space: surface + field precision: double + - name: sea_water_potential_temperature + nemo field name: votemper + model space: volume + field precision: double + - name: ice_area_fraction_background_error + nemo field name: sic_tot_var + model space: surface + variable type: background variance + field precision: double + grid name: ORCA2_T + number levels: 10 + background : + date: 2021-06-29T12:00:00Z + state variables: [ ice_area_fraction ] + nemo field file: Data/orca2_t_nemo.nc + observations: + observers: + - obs space: + name: Sea Ice + obsdatain: + engine: + type: H5File + obsfile: Data/sic_obs_ideal.nc + obsdataout: + engine: + type: H5File + obsfile: testoutput/obsdataout_3dvar_sic.nc + simulated variables: [ice_area_fraction] + get values: + time interpolation: linear + atlas-interpolator: + type: unstructured-bilinear-lonlat + non_linear: missing-if-all-missing + max_fraction_elems_to_try: 0.0 + adjoint: true + obs operator: + name: Composite + components: + - name: Identity + obs error: + covariance model: diagonal +variational: + minimizer: + algorithm: DRPCG + iterations: + - ninner: 1 + gradient norm reduction: 1e-30 + geometry : + nemo variables: + - name: ice_area_fraction + nemo field name: iiceconc + model space: surface + field precision: double + - name: sea_surface_temperature + nemo field name: votemper + model space: surface + field precision: double + - name: sea_water_potential_temperature + nemo field name: votemper + model space: volume + field precision: double + - name: ice_area_fraction_background_error + nemo field name: sic_tot_var + model space: surface + variable type: background variance + field precision: double + grid name: ORCA2_T + number levels: 10 + test: on + online diagnostics: + write increment: true + increment: + state component: + output path: testoutput/increments_3dvar_sic.nc + date: 2021-06-30T00:00:00Z +final: + diagnostics: + departures: oman +output: + date: 2021-06-30T00:00:00Z + state variables: [ ice_area_fraction ] + nemo field file: testoutput/3dvar_sic2.nc + output nemo field file: testoutput/analysis_3dvar_sic.nc +test: + reference filename: testoutput/test_3dvar_sic.ref + test output filename: testoutput/test_3dvar_sic.out diff --git a/src/tests/testinput/CMakeLists.txt b/src/tests/testinput/CMakeLists.txt index 4470af6..9109a97 100644 --- a/src/tests/testinput/CMakeLists.txt +++ b/src/tests/testinput/CMakeLists.txt @@ -18,6 +18,8 @@ list( APPEND orcajedi_test_input odb_sic_query.yaml odb_sic_mapping.yaml test_name_map.yaml + dirac.yaml + 3dvar_sic.yaml ) foreach(FILENAME ${orcajedi_test_input}) diff --git a/src/tests/testinput/dirac.yaml b/src/tests/testinput/dirac.yaml new file mode 100644 index 0000000..fe124a2 --- /dev/null +++ b/src/tests/testinput/dirac.yaml @@ -0,0 +1,41 @@ +background error: + covariance model: SABER + saber central block: + saber block name: ID + +dirac: + x indices: [20] + y indices: [10] + z indices: [0] + +geometry: + nemo variables: + - name: ice_area_fraction + nemo field name: iiceconc + model space: surface + field precision: double + - name: sea_surface_temperature + nemo field name: votemper + model space: surface + field precision: double + - name: sea_water_potential_temperature + nemo field name: votemper + model space: volume + field precision: double + - name: ice_area_fraction_background_error + nemo field name: sic_tot_var + model space: surface + variable type: background variance + field precision: double + grid name: ORCA2_T + number levels: 10 +background: + date: 2021-06-30T12:00:00Z + state variables: [ ice_area_fraction ] + nemo field file: Data/orca2_t_nemo.nc +output dirac: + date: 2021-06-30T12:00:00Z + output path: testoutput/dirac_cov_%id%.nc +test: + reference filename: testoutput/test_dirac.ref + test output filename: testoutput/test_dirac.out diff --git a/src/tests/testoutput/CMakeLists.txt b/src/tests/testoutput/CMakeLists.txt index 2ecd402..340085a 100644 --- a/src/tests/testoutput/CMakeLists.txt +++ b/src/tests/testoutput/CMakeLists.txt @@ -12,6 +12,8 @@ list( APPEND orcajedi_test_output test_hofx3d_nc_prof_2vars.ref test_hofx3d_nc_ssh_amm1.ref test_hofx3d_nc_ssh_amm1r.ref + test_3dvar_sic.ref + test_dirac.ref ) foreach(FILENAME ${orcajedi_test_output}) diff --git a/src/tests/testoutput/test_3dvar_sic.ref b/src/tests/testoutput/test_3dvar_sic.ref new file mode 100644 index 0000000..32c8f87 --- /dev/null +++ b/src/tests/testoutput/test_3dvar_sic.ref @@ -0,0 +1,16 @@ +CostJb : Nonlinear Jb = 0.0000000000000000e+00 +CostJo : Nonlinear Jo(Sea Ice) = 1.5000000000000000e+00, nobs = 48, Jo/n = 3.1250000000000000e-02, err = 2.0000000000000000e+00 +CostFunction: Nonlinear J = 1.5000000000000000e+00 +DRIPCGMinimizer: reduction in residual norm = 7.2455060481720379e-01 +CostFunction::addIncrement: Analysis: + Model state valid at time: 2021-06-29T12:00:00Z + 1 variables: ice_area_fraction + atlas field norms: + ice_area_fraction: 3.20155e-03 + + + + +CostJb : Nonlinear Jb = 0.00000e+00 +CostJo : Nonlinear Jo(Sea Ice) = 1.13468e+00, nobs = 48, Jo/n = 2.36392e-02, err = 2.00000e+00 +CostFunction: Nonlinear J = 1.13468e+00 diff --git a/src/tests/testoutput/test_dirac.ref b/src/tests/testoutput/test_dirac.ref new file mode 100644 index 0000000..0c2c41f --- /dev/null +++ b/src/tests/testoutput/test_dirac.ref @@ -0,0 +1,10 @@ +Input Dirac increment:Increment valid at time: 2021-06-30T12:00:00Z + 1 variables: ice_area_fraction + atlas field: + ice_area_fraction num: 26553 mean: 3.76605e-05 rms: 6.13682e-03 min: 0.00000e+00 max: 1.00000e+00 + +Covariance(SABER) * Increment:Increment valid at time: 2021-06-30T12:00:00Z + 1 variables: ice_area_fraction + atlas field: + ice_area_fraction num: 26553 mean: 3.76605e-05 rms: 6.13682e-03 min: 0.00000e+00 max: 1.00000e+00 +