diff --git a/src/orca-jedi/increment/Increment.cc b/src/orca-jedi/increment/Increment.cc index b51e167..9db3a28 100644 --- a/src/orca-jedi/increment/Increment.cc +++ b/src/orca-jedi/increment/Increment.cc @@ -37,11 +37,6 @@ #include "atlas/mesh.h" #include "atlas-orca/grid/OrcaGrid.h" -#include // uuid class DJL -#include // generators DJL -#include // streaming operators etc. DJL - - namespace orcamodel { // ----------------------------------------------------------------------------- @@ -232,8 +227,6 @@ void Increment::diff(const State & x1, const State & x2) { atlas::field::MissingValue mv2(field2); bool has_mv2 = static_cast(mv2); bool has_mv = has_mv1 || has_mv2; - oops::Log::debug() << "DJL Increment::diff mv1 " << mv1 << " mv2 " << mv2 << " has_mv " - << has_mv << std::endl; std::string fieldName1 = field1.name(); std::string fieldName2 = field2.name(); @@ -284,18 +277,11 @@ void Increment::setval(const double & val) { if (!has_mv || (has_mv && !mv(field_view(j, k)))) { field_view(j, k) = val; } - } else { - field_view(j, k) = -99; /// DJL } } } } - // DJL write debug fields to file - boost::uuids::uuid uuid = boost::uuids::random_generator()(); - writeFieldsToFile("incsetval"+ boost::uuids::to_string(uuid) +".nc", *geom_, - validTime(), incrementFields_); - oops::Log::trace() << "Increment(ORCA)::setval done" << std::endl; } @@ -559,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; @@ -572,10 +558,16 @@ void Increment::write(const OrcaIncrementParameters & params) const { oops::Log::debug() << "Increment::write to filename " << nemo_field_path << std::endl; + for (size_t i=0; i < vars_.size(); ++i) { + auto gv_varname = vars_[i].name(); +// halo exchange update ghost points + geom_->functionSpace().haloExchange(incrementFields_[gv_varname]); + } + 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 a335b20..4bb8251 100644 --- a/src/orca-jedi/interpolator/Interpolator.cc +++ b/src/orca-jedi/interpolator/Interpolator.cc @@ -31,13 +31,6 @@ #include "orca-jedi/interpolator/Interpolator.h" -#include // uuid class DJL -#include // generators DJL -#include // streaming operators etc. DJL - -#include "orca-jedi/utilities/IOUtils.h" // DJL - - namespace eckit { class Configuration; } @@ -193,7 +186,6 @@ template void Interpolator::executeInterpolation( void Interpolator::apply(const oops::Variables& vars, const Increment& inc, const std::vector & mask, std::vector& result) const { - // DJL question can the atlas templates help? // input is inc output is result const size_t nvars = vars.size(); @@ -226,7 +218,6 @@ void Interpolator::apply(const oops::Variables& vars, const Increment& inc, auto field_view = atlas::array::make_view(tgt_field); atlas::field::MissingValue mv(inc.incrementFields()[gv_varname]); bool has_mv = static_cast(mv); - oops::Log::debug() << "DJL Interpolator::apply mv " << mv << " has_mv " << has_mv << std::endl; 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))) { @@ -254,10 +245,6 @@ void Interpolator::applyAD(const oops::Variables& vars, Increment& inc, oops::Log::trace() << "orcamodel::Interpolator::applyAD start " << std::endl; - oops::Log::debug() << "DJL ** Interpolator::applyAD this needs checking **" << std::endl; - -// ** Trying to do the opposite of apply DJL - const size_t nvars = vars.size(); for (size_t j=0; j < nvars; ++j) { @@ -277,65 +264,38 @@ void Interpolator::applyAD(const oops::Variables& vars, Increment& inc, inc.geometry()->variableSizes(vars); size_t nvals = 0; -// DJL write debug fields to file - boost::uuids::uuid uuid = boost::uuids::random_generator()(); - std::shared_ptr geom = inc.geometry(); - writeFieldsToFile("applyADpre"+ boost::uuids::to_string(uuid) +".nc", *geom, - inc.validTime(), inc.incrementFields()); - 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) { - oops::Log::debug() << "DJL ** jvar " << jvar << " " << nvars - << "varSizes " << varSizes[jvar] - << std::endl; auto gv_varname = vars[jvar].name(); -// atlas::Field tgt_field = atlasObsFuncSpace_.createField( -// atlas::option::name(gv_varname) | -// atlas::option::levels(varSizes[jvar])); -// atlas::Field incField = inc.incrementFields()[jvar]; - 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) -// DJL not sure if the missing value aspect does anything here auto field_view = atlas::array::make_view(tgt_field); -// field_view.assign(0.0); atlas::field::MissingValue mv(inc.incrementFields()[gv_varname]); bool has_mv = static_cast(mv); - oops::Log::debug() << "DJL Interpolator::applyAD mv " - << mv << " has_mv " << has_mv << std::endl; 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 { - oops::Log::debug() << "DJL iloc " << iloc << " klev " << klev << " out_idx " - << out_idx << " resultin[out_idx] " << resultin[out_idx] << std::endl; field_view(iloc, klev) = resultin[out_idx]; } ++out_idx; } } -// halo exchange update ghost points DJL +// 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 -// DJL write debug fields to file -// boost::uuids::uuid uuid = boost::uuids::random_generator()(); -// std::shared_ptr geom = inc.geometry(); - writeFieldsToFile("applyAD"+ boost::uuids::to_string(uuid) +".nc", *geom, inc.validTime(), - inc.incrementFields()); - oops::Log::trace() << "orcamodel::Interpolator::applyAD done " << std::endl; } diff --git a/src/orca-jedi/state/State.cc b/src/orca-jedi/state/State.cc index ade4b26..d60a8f3 100644 --- a/src/orca-jedi/state/State.cc +++ b/src/orca-jedi/state/State.cc @@ -164,14 +164,14 @@ 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) { oops::Log::trace() << "State(ORCA)::add increment starting" << std::endl; ASSERT(this->validTime() == dx.validTime()); -// DJL assumes state and increments have the same fields in the field set -// DJL add something to enforce/check this - auto ghost = atlas::array::make_view( geom_->mesh().nodes().ghost()); for (int i = 0; i< stateFields_.size(); i++) @@ -187,6 +187,8 @@ State & State::operator+=(const Increment & dx) { 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) { diff --git a/src/tests/testinput/3dvar_ice_obs.yaml b/src/tests/testinput/3dvar_ice_obs.yaml index 1280805..f3d90f1 100644 --- a/src/tests/testinput/3dvar_ice_obs.yaml +++ b/src/tests/testinput/3dvar_ice_obs.yaml @@ -44,7 +44,7 @@ cost function: obsdataout: engine: type: H5File - obsfile: testoutput/test_3dvar_obsdataout.nc + obsfile: testoutput/obsdataout_3dvar_ice.nc simulated variables: [ice_area_fraction] get values: time interpolation: linear @@ -91,9 +91,8 @@ variational: write increment: true increment: state component: - output path: testoutput/increments_ice.nc + output path: testoutput/increments_3dvar_ice.nc date: 2021-06-30T00:00:00Z -# type: in final: diagnostics: departures: oman @@ -101,7 +100,4 @@ output: date: 2021-06-30T00:00:00Z state variables: [ ice_area_fraction ] nemo field file: testoutput/3dvar_ice2.nc - output nemo field file: testoutput/3dvar_ice_analysis.nc -# steps: ["2021-06-29T12:00:00Z"] -#test: -# reference filename: testoutput/test_3dvar_ice_obs.ref + output nemo field file: testoutput/analysis_3dvar_ice.nc diff --git a/src/tests/testinput/dirac.yaml b/src/tests/testinput/dirac.yaml index 6fc2a45..582aa02 100644 --- a/src/tests/testinput/dirac.yaml +++ b/src/tests/testinput/dirac.yaml @@ -36,6 +36,3 @@ background: output dirac: date: 2021-06-30T12:00:00Z output path: testoutput/dirac_cov_%id%.nc - -#test: -# reference filename: testoutput/test_dirac.ref