From 94faae82e8396bc9edb55bedea10ce45d1475397 Mon Sep 17 00:00:00 2001 From: Travis Sluka Date: Mon, 28 Aug 2023 13:01:06 -0600 Subject: [PATCH 1/4] change one 3dvar test, it will break --- test/testinput/3dvar_godas.yml | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/test/testinput/3dvar_godas.yml b/test/testinput/3dvar_godas.yml index 1852d0d24..ebb8baf7b 100644 --- a/test/testinput/3dvar_godas.yml +++ b/test/testinput/3dvar_godas.yml @@ -192,7 +192,13 @@ cost function: obsfile: data_static/obs/prof.nc simulated variables: [waterTemperature] obs operator: - name: InsituTemperature + name: VertInterp + observation alias file: testinput/obsop_name_map.yml + variables: + - waterTemperature + vertical coordinate: sea_water_depth + observation vertical coordinate: depth + interpolation method: linear obs error: covariance model: diagonal obs filters: From 01bf1fb3defb36c16438883d9796e582a3ab4c34 Mon Sep 17 00:00:00 2001 From: Travis Sluka Date: Mon, 2 Oct 2023 15:27:30 -0600 Subject: [PATCH 2/4] add vader to linearvariable change --- .../LinearModel2GeoVaLs.F90 | 35 +++-- .../LinearVariableChange.cc | 120 ++++++++++++++++-- .../LinearVariableChange.h | 5 + 3 files changed, 140 insertions(+), 20 deletions(-) diff --git a/src/soca/LinearVariableChange/LinearModel2GeoVaLs/LinearModel2GeoVaLs.F90 b/src/soca/LinearVariableChange/LinearModel2GeoVaLs/LinearModel2GeoVaLs.F90 index 390db804d..c8ddfb8fb 100644 --- a/src/soca/LinearVariableChange/LinearModel2GeoVaLs/LinearModel2GeoVaLs.F90 +++ b/src/soca/LinearVariableChange/LinearModel2GeoVaLs/LinearModel2GeoVaLs.F90 @@ -46,18 +46,29 @@ subroutine soca_model2geovals_linear_changevar_f90(c_key_geom, c_key_dxin, c_key ! identity operators do i=1, size(dxout%fields) - call dxin%get(dxout%fields(i)%metadata%name, field) - if (dxout%fields(i)%name == field%metadata%name .or. & - dxout%fields(i)%name == field%metadata%getval_name) then - dxout%fields(i)%val(:,:,:) = field%val(:,:,:) !< full field - elseif (field%metadata%getval_name_surface == dxout%fields(i)%name) then - dxout%fields(i)%val(:,:,1) = field%val(:,:,1) !< surface only of a 3D field - - else - call abor1_ftn( 'error in soca_model2geovals_linear_changevar_f90 processing ' & - // dxout%fields(i)%name ) - endif - + select case (dxout%fields(i)%name) + case ( 'latitude' ) + dxout%fields(i)%val(:,:,1) = real(geom%lat, kind=kind_real) + + case ( 'longitude' ) + dxout%fields(i)%val(:,:,1) = real(geom%lon, kind=kind_real) + + case ( 'sea_water_depth') + call geom%thickness2depth(geom%h, dxout%fields(i)%val) + + case default + call dxin%get(dxout%fields(i)%metadata%name, field) + if (dxout%fields(i)%name == field%metadata%name .or. & + dxout%fields(i)%name == field%metadata%getval_name) then + dxout%fields(i)%val(:,:,:) = field%val(:,:,:) !< full field + elseif (field%metadata%getval_name_surface == dxout%fields(i)%name) then + dxout%fields(i)%val(:,:,1) = field%val(:,:,1) !< surface only of a 3D field + + else + call abor1_ftn( 'error in soca_model2geovals_linear_changevar_f90 processing ' & + // dxout%fields(i)%name ) + endif + end select end do end subroutine diff --git a/src/soca/LinearVariableChange/LinearVariableChange.cc b/src/soca/LinearVariableChange/LinearVariableChange.cc index f4ef49436..03b3e8e59 100644 --- a/src/soca/LinearVariableChange/LinearVariableChange.cc +++ b/src/soca/LinearVariableChange/LinearVariableChange.cc @@ -5,10 +5,12 @@ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. */ +#include #include #include #include "soca/LinearVariableChange/LinearVariableChange.h" +#include "soca/VariableChange/Model2GeoVaLs/Model2GeoVaLs.h" #include "soca/Geometry/Geometry.h" #include "soca/Increment/Increment.h" @@ -22,10 +24,26 @@ namespace soca { // ----------------------------------------------------------------------------- +std::map> SocaLinVaderCookbook { + {"sea_water_temperature", {"SeaWaterTemperature_A"}}, +}; + +// ----------------------------------------------------------------------------- + LinearVariableChange::LinearVariableChange(const Geometry & geom, const eckit::Configuration & config) - : geom_(geom), params_(), linVarChas_() { + : geom_(geom), params_(), linVarChas_(), vader_(), default_(false) { params_.deserialize(config); + + // setup vader + eckit::LocalConfiguration vaderConfig, vaderCookbookConfig; + for (auto kv : SocaLinVaderCookbook) vaderCookbookConfig.set(kv.first, kv.second); + vaderConfig.set(vader::configCookbookKey, vaderCookbookConfig); + vader_.reset(new vader::Vader(params_.vader, vaderConfig)); + + const boost::optional> & + linVarChgs = params_.linearVariableChangesWrapper; + default_ = linVarChgs == boost::none; } // ----------------------------------------------------------------------------- @@ -36,14 +54,53 @@ LinearVariableChange::~LinearVariableChange() {} void LinearVariableChange::changeVarTraj(const State & xfg, const oops::Variables & vars) { - Log::trace() << "LinearVariableChange::setTrajectory starting" << std::endl; + Log::trace() << "LinearVariableChange::changeVarTraj starting, default: "<{ + "latitude", + "longitude", + "sea_water_potential_temperature", + "sea_water_salinity", + "sea_water_depth"}); + preVaderVars += xfg.variables(); + xfg2.updateFields(preVaderVars); + Model2GeoVaLs m2gv(xfg.geometry(), params_.toConfiguration()); + m2gv.changeVar(xfg, xfg2); + Log::debug() << "LinearVariableChange::changeVarTraj variables after var change: " + << xfg2.variables() << std::endl; + } - // TODO(travis): do something with vars? + // call Vader + // ---------------------------------------------------------------------------- + Log::debug() << "LinearVariableChange::changeVarTraj VADER variable changes. " << std::endl; + oops::Variables varsFilled = xfg2.variables(); + oops::Variables varsVader = vars; + varsVader -= varsFilled; // pass only the needed variables + atlas::FieldSet xfs; + xfg2.toFieldSet(xfs); + varsFilled += vader_->changeVarTraj(xfs, varsVader); + xfg2.updateFields(varsFilled); + xfg2.fromFieldSet(xfs); + Log::debug() << "LinearVariableChange::changeVarTraj variables after var change: " + << xfg2.variables() << std::endl; // TODO(travis): this is not ideal. We are saving the first trajectory and // assuming it is the background. This should all get ripped out when the // variable changes that rely on the background are dealt with properly. - if (!bkg_) { bkg_.reset(new State(xfg)); } + if (!bkg_) { bkg_.reset(new State(xfg2)); } const boost::optional> & linVarChgs = params_.linearVariableChangesWrapper; @@ -59,13 +116,13 @@ void LinearVariableChange::changeVarTraj(const State & xfg, linVarChaParWra.linearVariableChangeParameters; // Add linear variable change to vector linVarChas_.push_back( - LinearVariableChangeFactory::create(*bkg_, xfg, geom_, linVarChaPar)); + LinearVariableChangeFactory::create(*bkg_, xfg2, geom_, linVarChaPar)); } - } else { + } else { // No variable changes were specified, use the default (LinearModel2GeoVaLs) eckit::LocalConfiguration conf; conf.set("linear variable change name", "default"); - linVarChas_.push_back(LinearVariableChangeFactory::create(*bkg_, xfg, geom_, + linVarChas_.push_back(LinearVariableChangeFactory::create(*bkg_, xfg2, geom_, oops::validateAndDeserialize( conf))); } @@ -76,7 +133,54 @@ void LinearVariableChange::changeVarTraj(const State & xfg, void LinearVariableChange::changeVarTL(Increment & dx, const oops::Variables & vars) const { - Log::trace() << "LinearVariableChange::multiply starting" << std::endl; + Log::trace() << "LinearVariableChange::changeVarTL starting, default: " << default_ << std::endl; + + Log::debug() << "LinearVariableChange::changeVarTL vars in: " + << dx.variables() << std::endl; + Log::debug() << "LinearVariableChange::changeVarTL vars out: " + << vars << std::endl; + + + // The following is TEMPORARY. + // ---------------------------------------------------------------------------- + // We need to do some variable renaming BEFORE we run VADER. + // Eventually, we will internally rename these variables when they are + // first loaded in so that we won't have to worry about it here. + if (default_ && vars.has("sea_water_temperature")) { + Log::debug() << "LinearVariableChange::changeVarTL Pre-VADER variable changes. " << std::endl; + oops::Variables preVaderVars(std::vector{ + "latitude", + "longitude", + "sea_water_potential_temperature", + "sea_water_salinity", + "sea_water_depth"}); + preVaderVars += dx.variables(); + Increment preVader(dx.geometry(), preVaderVars, dx.time()); + + for (icst_ it = linVarChas_.begin(); it != linVarChas_.end(); ++it) { + it->multiply(dx, preVader); + dx.updateFields(preVaderVars); + dx = preVader; + } + Log::debug() << "LinearVariableChange::changeVarTL variables after var change: " + << dx.variables() << std::endl; + } + + + // call Vader + // ---------------------------------------------------------------------------- + Log::debug() << "LinearVariableChange::changeVarTL VADER variable changes. " << std::endl; + atlas::FieldSet xfs; + dx.toFieldSet(xfs); + oops::Variables varsFilled = dx.variables(); + oops::Variables varsVader = vars; + varsVader -= varsFilled; + varsFilled += vader_->changeVarTL(xfs, varsVader); + Log::debug() << "DBG " << varsFilled <> linearVariableChangesWrapper{"linear variable changes", this}; + oops::Parameter vader{"vader", {}, this}; }; // ----------------------------------------------------------------------------- @@ -68,6 +71,8 @@ class LinearVariableChange : public util::Printable { const Geometry & geom_; std::unique_ptr bkg_; LinVarChaVec_ linVarChas_; + std::unique_ptr vader_; + bool default_; }; // ----------------------------------------------------------------------------- From e30ac92d30e791fde7f35860ca6e979106a45c98 Mon Sep 17 00:00:00 2001 From: Travis Sluka Date: Tue, 3 Oct 2023 10:38:56 -0600 Subject: [PATCH 3/4] linearvariable change working for adjoint --- .../LinearModel2GeoVaLs.F90 | 13 ++-- .../LinearVariableChange.cc | 63 +++++++++++++------ .../LinearVariableChange.h | 2 + src/soca/VariableChange/VariableChange.cc | 3 +- 4 files changed, 55 insertions(+), 26 deletions(-) diff --git a/src/soca/LinearVariableChange/LinearModel2GeoVaLs/LinearModel2GeoVaLs.F90 b/src/soca/LinearVariableChange/LinearModel2GeoVaLs/LinearModel2GeoVaLs.F90 index c8ddfb8fb..7ba0e20bb 100644 --- a/src/soca/LinearVariableChange/LinearModel2GeoVaLs/LinearModel2GeoVaLs.F90 +++ b/src/soca/LinearVariableChange/LinearModel2GeoVaLs/LinearModel2GeoVaLs.F90 @@ -47,14 +47,17 @@ subroutine soca_model2geovals_linear_changevar_f90(c_key_geom, c_key_dxin, c_key ! identity operators do i=1, size(dxout%fields) select case (dxout%fields(i)%name) - case ( 'latitude' ) - dxout%fields(i)%val(:,:,1) = real(geom%lat, kind=kind_real) + ! TODO: These should NOT be needed in an increment. Need to make some changes to vader... + case ( 'sea_area_fraction' ) + dxout%fields(i)%val(:,:,1) = 0.0 + case ( 'latitude' ) + dxout%fields(i)%val(:,:,1) = 0.0 case ( 'longitude' ) - dxout%fields(i)%val(:,:,1) = real(geom%lon, kind=kind_real) - + dxout%fields(i)%val(:,:,1) = 0.0 case ( 'sea_water_depth') - call geom%thickness2depth(geom%h, dxout%fields(i)%val) + dxout%fields(i)%val(:,:,1) = 0.0 + ! call geom%thickness2depth(geom%h, dxout%fields(i)%val) case default call dxin%get(dxout%fields(i)%metadata%name, field) diff --git a/src/soca/LinearVariableChange/LinearVariableChange.cc b/src/soca/LinearVariableChange/LinearVariableChange.cc index 03b3e8e59..cbf49e8c1 100644 --- a/src/soca/LinearVariableChange/LinearVariableChange.cc +++ b/src/soca/LinearVariableChange/LinearVariableChange.cc @@ -74,7 +74,9 @@ void LinearVariableChange::changeVarTraj(const State & xfg, "longitude", "sea_water_potential_temperature", "sea_water_salinity", - "sea_water_depth"}); + "sea_water_depth", + "sea_area_fraction", + }); preVaderVars += xfg.variables(); xfg2.updateFields(preVaderVars); Model2GeoVaLs m2gv(xfg.geometry(), params_.toConfiguration()); @@ -91,7 +93,8 @@ void LinearVariableChange::changeVarTraj(const State & xfg, varsVader -= varsFilled; // pass only the needed variables atlas::FieldSet xfs; xfg2.toFieldSet(xfs); - varsFilled += vader_->changeVarTraj(xfs, varsVader); + varsVaderPopulates_ = vader_->changeVarTraj(xfs, varsVader); + varsFilled += varsVaderPopulates_; xfg2.updateFields(varsFilled); xfg2.fromFieldSet(xfs); Log::debug() << "LinearVariableChange::changeVarTraj variables after var change: " @@ -153,7 +156,8 @@ void LinearVariableChange::changeVarTL(Increment & dx, "longitude", "sea_water_potential_temperature", "sea_water_salinity", - "sea_water_depth"}); + "sea_water_depth", + "sea_area_fraction"}); preVaderVars += dx.variables(); Increment preVader(dx.geometry(), preVaderVars, dx.time()); @@ -176,21 +180,11 @@ void LinearVariableChange::changeVarTL(Increment & dx, oops::Variables varsVader = vars; varsVader -= varsFilled; varsFilled += vader_->changeVarTL(xfs, varsVader); - Log::debug() << "DBG " << varsFilled <{ + "sea_water_potential_temperature"}); + oops::Variables varsToDrop(std::vector{ + "sea_water_temperature"}); + + // run Vader + oops::Variables vaderVars = dx.variables(); + vaderVars += varsToAdj; + dx.updateFields(vaderVars); + + atlas::FieldSet xfs; + dx.toFieldSet(xfs); + + oops::Variables varsVaderWillAdjoint = varsVaderPopulates_; + vader_->changeVarAD(xfs, varsVaderWillAdjoint); + ASSERT(varsVaderWillAdjoint.size() == 0); + + vaderVars -= varsToDrop; + dx.updateFields(vaderVars); + dx.fromFieldSet(xfs); + Log::debug() << "LinearVariableChange::changeVarTL variables after var change: " + << dx.variables() << std::endl; + } + + Increment dxout(dx.geometry(), vars, dx.time()); // Call variable change(s) diff --git a/src/soca/LinearVariableChange/LinearVariableChange.h b/src/soca/LinearVariableChange/LinearVariableChange.h index bad638158..a54170f50 100644 --- a/src/soca/LinearVariableChange/LinearVariableChange.h +++ b/src/soca/LinearVariableChange/LinearVariableChange.h @@ -71,7 +71,9 @@ class LinearVariableChange : public util::Printable { const Geometry & geom_; std::unique_ptr bkg_; LinVarChaVec_ linVarChas_; + std::unique_ptr vader_; + oops::Variables varsVaderPopulates_; bool default_; }; diff --git a/src/soca/VariableChange/VariableChange.cc b/src/soca/VariableChange/VariableChange.cc index 34f899438..1d84139f7 100644 --- a/src/soca/VariableChange/VariableChange.cc +++ b/src/soca/VariableChange/VariableChange.cc @@ -76,7 +76,8 @@ void VariableChange::changeVar(State & x, const oops::Variables & vars) const { "longitude", "sea_water_potential_temperature", "sea_water_salinity", - "sea_water_depth"}); + "sea_water_depth", + "sea_area_fraction"}); preVaderVars += x.variables(); State preVader(x.geometry(), preVaderVars, x.time()); variableChange_->changeVar(x, preVader); From b9f37ccb3675b708034f650ff5972b8a2c35fd7d Mon Sep 17 00:00:00 2001 From: Anna Shlyaeva Date: Tue, 30 Dec 2025 10:41:51 -0700 Subject: [PATCH 4/4] Updates to work with updated vader, some cleanup --- .../LinearVariableChange.cc | 213 ++++++++---------- .../LinearVariableChange.h | 7 +- src/soca/VariableChange/VariableChange.cc | 5 +- test/testinput/3dvar.yml | 8 +- test/testinput/3dvar_nicas.yml | 8 +- test/testref/3dvar_nicas.test | 32 +-- 6 files changed, 128 insertions(+), 145 deletions(-) diff --git a/src/soca/LinearVariableChange/LinearVariableChange.cc b/src/soca/LinearVariableChange/LinearVariableChange.cc index a4994a947..5fef9ae86 100644 --- a/src/soca/LinearVariableChange/LinearVariableChange.cc +++ b/src/soca/LinearVariableChange/LinearVariableChange.cc @@ -16,6 +16,7 @@ #include "soca/Increment/Increment.h" #include "soca/State/State.h" +#include "oops/util/FieldSetHelpers.h" #include "oops/util/Logger.h" using oops::Log; @@ -25,7 +26,7 @@ namespace soca { // ----------------------------------------------------------------------------- std::map> SocaLinVaderCookbook { - {"sea_water_temperature", {"SeaWaterTemperature_A"}}, + {"sea_water_temperature", {"SeaWaterTemperature_A", "SeaWaterTemperature_B"}} }; // ----------------------------------------------------------------------------- @@ -41,9 +42,7 @@ LinearVariableChange::LinearVariableChange(const Geometry & geom, vaderConfig.set(vader::configCookbookKey, vaderCookbookConfig); vader_.reset(new vader::Vader(params_.vader, vaderConfig)); - const boost::optional> & - linVarChgs = params_.linearVariableChangesWrapper; - default_ = linVarChgs == boost::none; + default_ = (params_.linearVariableChangesWrapper.value() == boost::none); } // ----------------------------------------------------------------------------- @@ -54,56 +53,42 @@ LinearVariableChange::~LinearVariableChange() {} void LinearVariableChange::changeVarTraj(const State & xfg, const oops::Variables & vars) { - Log::trace() << "LinearVariableChange::changeVarTraj starting, default: "<{ "latitude", "longitude", - "sea_water_potential_temperature", - "sea_water_salinity", "sea_water_depth", "sea_area_fraction", }); preVaderVars += xfg.variables(); - xfg2.updateFields(preVaderVars); - Model2GeoVaLs m2gv(xfg.geometry(), params_.toConfiguration()); - m2gv.changeVar(xfg, xfg2); + xfg_vadertraj.updateFields(preVaderVars); + Model2GeoVaLs m2gv(geom_, params_.toConfiguration()); + m2gv.changeVar(xfg, xfg_vadertraj); Log::debug() << "LinearVariableChange::changeVarTraj variables after var change: " - << xfg2.variables() << std::endl; + << xfg_vadertraj.variables() << std::endl; } // call Vader // ---------------------------------------------------------------------------- - Log::debug() << "LinearVariableChange::changeVarTraj VADER variable changes. " << std::endl; - oops::Variables varsFilled = xfg2.variables(); + oops::Variables varsFilled = xfg_vadertraj.variables(); oops::Variables varsVader = vars; - varsVader -= varsFilled; // pass only the needed variables + varsVader -= varsFilled; // pass only the needed variables atlas::FieldSet xfs; - xfg2.toFieldSet(xfs); - varsVaderPopulates_ = vader_->changeVarTraj(xfs, varsVader); - varsFilled += varsVaderPopulates_; - xfg2.updateFields(varsFilled); - xfg2.fromFieldSet(xfs); - Log::debug() << "LinearVariableChange::changeVarTraj variables after var change: " - << xfg2.variables() << std::endl; + xfg_vadertraj.toFieldSet(xfs); + xfs.haloExchange(); + vader_->changeVarTraj(xfs, varsVader); // TODO(travis): this is not ideal. We are saving the first trajectory and // assuming it is the background. This should all get ripped out when the // variable changes that rely on the background are dealt with properly. - if (!bkg_) { bkg_.reset(new State(xfg2)); } + if (!bkg_) { bkg_.reset(new State(xfg)); } const boost::optional> & linVarChgs = params_.linearVariableChangesWrapper; @@ -119,71 +104,63 @@ void LinearVariableChange::changeVarTraj(const State & xfg, linVarChaParWra.linearVariableChangeParameters; // Add linear variable change to vector linVarChas_.push_back( - LinearVariableChangeFactory::create(*bkg_, xfg2, geom_, linVarChaPar)); + LinearVariableChangeFactory::create(*bkg_, xfg, geom_, linVarChaPar)); } - } else { + } else { // No variable changes were specified, use the default (LinearModel2GeoVaLs) eckit::LocalConfiguration conf; conf.set("linear variable change name", "default"); - linVarChas_.push_back(LinearVariableChangeFactory::create(*bkg_, xfg2, geom_, + linVarChas_.push_back(LinearVariableChangeFactory::create(*bkg_, xfg, geom_, oops::validateAndDeserialize( conf))); } Log::trace() << "LinearVariableChange::setTrajectory done" << std::endl; } +// ------------------------------------------------------------------------------------------------- + +void LinearVariableChange::initVaderTLAD(oops::Variables & ingredientVars) const { + oops::Log::trace() << "LinearVariableChange::initVaderTLAD starting" << std::endl; + oops::Variables originalIngredientVars = ingredientVars; + varsVaderPopulates_ = vader_->initTLAD(ingredientVars); + varsVaderPopulates_ -= originalIngredientVars; + oops::Log::trace() << "LinearVariableChange::initVaderTLAD done" << std::endl; +} + // ----------------------------------------------------------------------------- void LinearVariableChange::changeVarTL(Increment & dx, const oops::Variables & vars) const { - Log::trace() << "LinearVariableChange::changeVarTL starting, default: " << default_ << std::endl; - - Log::debug() << "LinearVariableChange::changeVarTL vars in: " - << dx.variables() << std::endl; - Log::debug() << "LinearVariableChange::changeVarTL vars out: " - << vars << std::endl; - - - // The following is TEMPORARY. - // ---------------------------------------------------------------------------- - // We need to do some variable renaming BEFORE we run VADER. - // Eventually, we will internally rename these variables when they are - // first loaded in so that we won't have to worry about it here. - if (default_ && vars.has("sea_water_temperature")) { - Log::debug() << "LinearVariableChange::changeVarTL Pre-VADER variable changes. " << std::endl; - oops::Variables preVaderVars(std::vector{ - "latitude", - "longitude", - "sea_water_potential_temperature", - "sea_water_salinity", - "sea_water_depth", - "sea_area_fraction"}); - preVaderVars += dx.variables(); - Increment preVader(dx.geometry(), preVaderVars, dx.time()); - - for (icst_ it = linVarChas_.begin(); it != linVarChas_.end(); ++it) { - it->multiply(dx, preVader); - dx.updateFields(preVaderVars); - dx = preVader; - } - Log::debug() << "LinearVariableChange::changeVarTL variables after var change: " - << dx.variables() << std::endl; + Log::trace() << "LinearVariableChange::changeVarTL starting" << std::endl; + + // If all variables already in incoming state just remove the no longer + // needed fields + // if (hasAllFields) { + // dx.updateFields(vars); + // oops::Log::trace() << "VariableChange::changeVar done (identity)" + // << std::endl; + // return; + // } + + // Make sure vader is fully initialized + if (vader_->needsTLADInit()) { + oops::Variables ingredientVars = dx.variables(); + initVaderTLAD(ingredientVars); + oops::Log::info() << " changevarTL: vader will populate variables: " + << varsVaderPopulates_ << std::endl; + } + // If Vader is doing anything, call Vader + if (varsVaderPopulates_.size() > 0) { + atlas::FieldSet dxfs; + dx.toFieldSet(dxfs); + vader_->changeVarTL(dxfs); + // Set intermediate state for the Increment containing original fields plus the ones + // Vader has done + oops::Variables varsVader = dx.variables(); + varsVader += varsVaderPopulates_; + dx.updateFields(varsVader); + dx.fromFieldSet(dxfs); } - - - // call Vader - // ---------------------------------------------------------------------------- - Log::debug() << "LinearVariableChange::changeVarTL VADER variable changes. " << std::endl; - atlas::FieldSet xfs; - dx.toFieldSet(xfs); - oops::Variables varsFilled = dx.variables(); - oops::Variables varsVader = vars; - varsVader -= varsFilled; - varsFilled += vader_->changeVarTL(xfs, varsVader); - dx.updateFields(varsFilled); - dx.fromFieldSet(xfs); - Log::debug() << "LinearVariableChange::changeVarTL variables after var change: " - << dx.variables() << std::endl; // Create output state Increment dxout(dx.geometry(), vars, dx.validTime()); @@ -226,43 +203,27 @@ void LinearVariableChange::changeVarAD(Increment & dx, const oops::Variables & vars) const { Log::trace() << "LinearVariableChange::changeVarAD starting" << std::endl; - Log::debug() << "LinearVariableChange::changeVarAD vars in: " - << dx.variables() << std::endl; - Log::debug() << "LinearVariableChange::changeVarAD vars out: " - << vars << std::endl; - - // NOTE: the IF is temporary, we need to do some variable renaming afterward - if (default_ && dx.variables().has("sea_water_temperature")) { - // call vader - Log::debug() << "LinearVariableChange::changeVarAD VADER variable changes. " << std::endl; - oops::Variables varsToAdj(std::vector{ - "sea_water_potential_temperature"}); - oops::Variables varsToDrop(std::vector{ - "sea_water_temperature"}); - - // run Vader - oops::Variables vaderVars = dx.variables(); - vaderVars += varsToAdj; - dx.updateFields(vaderVars); - - atlas::FieldSet xfs; - dx.toFieldSet(xfs); - - oops::Variables varsVaderWillAdjoint = varsVaderPopulates_; - vader_->changeVarAD(xfs, varsVaderWillAdjoint); - ASSERT(varsVaderWillAdjoint.size() == 0); - - vaderVars -= varsToDrop; - dx.updateFields(vaderVars); - dx.fromFieldSet(xfs); - Log::debug() << "LinearVariableChange::changeVarTL variables after var change: " - << dx.variables() << std::endl; + // Make sure vader is fully initialized + if (vader_->needsTLADInit()) { + oops::Variables ingredientVars(vars); + initVaderTLAD(ingredientVars); + oops::Log::info() << " changevarAD: vader will populate variables: " + << varsVaderPopulates_ << std::endl; } + // Create dx_for_vader as a copy of dx, with just the variables created by Vader + // (in the forward direction). Remove the variables created by vader from dx. + // This way we ensure the model code will not be able to do the adjoint for these vars + Increment dx_for_vader(dx, true); // true => full copy + oops::Variables varsVaderDidntPopulate = dx.variables(); + varsVaderDidntPopulate -= varsVaderPopulates_; + dx_for_vader.updateFields(varsVaderPopulates_); + dx.updateFields(varsVaderDidntPopulate); - Increment dxout(dx.geometry(), vars, dx.time()); + // Create empty output state + Increment dxout(dx.geometry(), vars, dx.validTime()); - // Call variable change(s) + // Call model's adjoint variable change. for (ircst_ it = linVarChas_.rbegin(); it != linVarChas_.rend(); ++it) { dxout.zero(); it->multiplyAD(dx, dxout); @@ -270,6 +231,30 @@ void LinearVariableChange::changeVarAD(Increment & dx, dx = dxout; } + // dxout needs to temporarily have the variables that Vader populated put into it before + // being passed into vader_.changeVarAD, so Vader can do its adjoints. + if (varsVaderPopulates_.size() > 0) { + atlas::FieldSet dxfs; + dx.toFieldSet(dxfs); + oops::Variables varsVaderWillAdjoint = varsVaderPopulates_; + atlas::FieldSet dxvaderfs; + dx_for_vader.toFieldSet(dxvaderfs); + for (const auto field : dxvaderfs) { + dxfs.add(field); + } + + oops::Variables varsAdjointed = vader_->changeVarAD(dxfs); + varsVaderWillAdjoint -= varsAdjointed; + // After changeVarAD, vader should have removed everything from varsVaderWillAdjoint, + // indicating it did all the adjoints we expected it to. + ASSERT(varsVaderWillAdjoint.size() == 0); + // remove the fields that were adjointed by vader from dxout_fs + util::removeFieldsFromFieldSet(dxfs, varsAdjointed.variables()); + // Copy dxout into dx for return + dx.updateFields(vars); + dx.fromFieldSet(dxfs); + } + Log::trace() << "LinearVariableChange::multiplyAD done" << std::endl; } diff --git a/src/soca/LinearVariableChange/LinearVariableChange.h b/src/soca/LinearVariableChange/LinearVariableChange.h index dc6cf9149..b671a7253 100644 --- a/src/soca/LinearVariableChange/LinearVariableChange.h +++ b/src/soca/LinearVariableChange/LinearVariableChange.h @@ -40,7 +40,7 @@ class LinearVariableChangeParameters : public oops::Parameters { oops::OptionalParameter outputVariables{"output variables", this}; oops::OptionalParameter> linearVariableChangesWrapper{"linear variable changes", this}; - oops::Parameter vader{"vader", {}, this}; + oops::Parameter vader{"vader", {}, this}; }; // ----------------------------------------------------------------------------- @@ -68,13 +68,14 @@ class LinearVariableChange : public util::Printable { private: void print(std::ostream &) const override; + void initVaderTLAD(oops::Variables & ingredientVars) const; Parameters_ params_; const Geometry & geom_; std::unique_ptr bkg_; LinVarChaVec_ linVarChas_; - + std::unique_ptr vader_; - oops::Variables varsVaderPopulates_; + mutable oops::Variables varsVaderPopulates_; bool default_; }; diff --git a/src/soca/VariableChange/VariableChange.cc b/src/soca/VariableChange/VariableChange.cc index 73d7d4e58..eddb7d26c 100644 --- a/src/soca/VariableChange/VariableChange.cc +++ b/src/soca/VariableChange/VariableChange.cc @@ -78,10 +78,7 @@ void VariableChange::changeVar(State & x, const oops::Variables & vars) const { oops::Variables preVaderVars(std::vector{ "latitude", "longitude", - "sea_water_potential_temperature", - "sea_water_salinity", - "sea_water_depth", - "sea_area_fraction"}); + "sea_water_depth"}); preVaderVars += x.variables(); State preVader(x.geometry(), preVaderVars, x.validTime()); variableChange_->changeVar(x, preVader); diff --git a/test/testinput/3dvar.yml b/test/testinput/3dvar.yml index ea500b6a9..9fda7153a 100644 --- a/test/testinput/3dvar.yml +++ b/test/testinput/3dvar.yml @@ -230,13 +230,7 @@ cost function: obsfile: data_static/obs/prof.nc simulated variables: [waterTemperature] obs operator: - name: VertInterp - observation alias file: testinput/obsop_name_map.yml - variables: - - waterTemperature - vertical coordinate: sea_water_depth - observation vertical coordinate: depth - interpolation method: linear + name: InsituTemperature obs error: covariance model: diagonal obs filters: diff --git a/test/testinput/3dvar_nicas.yml b/test/testinput/3dvar_nicas.yml index f73a41e16..3c87d057d 100644 --- a/test/testinput/3dvar_nicas.yml +++ b/test/testinput/3dvar_nicas.yml @@ -176,7 +176,13 @@ cost function: obsfile: data_static/obs/prof.nc simulated variables: [waterTemperature] obs operator: - name: InsituTemperature + name: VertInterp + observation alias file: testinput/obsop_name_map.yml + variables: + - waterTemperature + vertical coordinate: sea_water_depth + observation vertical coordinate: depth + interpolation method: linear obs error: covariance model: diagonal obs filters: diff --git a/test/testref/3dvar_nicas.test b/test/testref/3dvar_nicas.test index 77ddc1569..ee62512e8 100644 --- a/test/testref/3dvar_nicas.test +++ b/test/testref/3dvar_nicas.test @@ -1,25 +1,25 @@ -CostJo : Nonlinear Jo(SeaSurfaceTemp) = 6.8893641637637177e+02, nobs = 113, Jo/n = 6.0967824458086000e+00, err = 3.6655656239312961e-01 +CostJo : Nonlinear Jo(SeaSurfaceTemp) = 6.8893641637637154e+02, nobs = 113, Jo/n = 6.0967824458085973e+00, err = 3.6655656239312961e-01 CostJo : Nonlinear Jo(SeaSurfaceSalinity) = 9.2074193097584697e+01, nobs = 49, Jo/n = 1.8790651652568306e+00, err = 1.0000000000000000e+00 -CostJo : Nonlinear Jo(ADT) = 2.2281591911422575e+02, nobs = 89, Jo/n = 2.5035496529688288e+00, err = 1.0000000149011612e-01 -CostJo : Nonlinear Jo(InsituTemperature) = 3.4206375767190184e+02, nobs = 203, Jo/n = 1.6850431412408957e+00, err = 9.0034813852418905e-01 -CostJo : Nonlinear Jo(InsituSalinity) = 3.0723498032711791e+02, nobs = 218, Jo/n = 1.4093347721427427e+00, err = 6.0819699545777528e-01 +CostJo : Nonlinear Jo(ADT) = 2.2281591911422569e+02, nobs = 89, Jo/n = 2.5035496529688279e+00, err = 1.0000000149011612e-01 +CostJo : Nonlinear Jo(InsituTemperature) = 3.4207282026717877e+02, nobs = 203, Jo/n = 1.6850877845673831e+00, err = 9.0034813852418905e-01 +CostJo : Nonlinear Jo(InsituSalinity) = 3.0723498032711768e+02, nobs = 218, Jo/n = 1.4093347721427416e+00, err = 6.0819699545777528e-01 CostJb : Nonlinear Jb = 0.0000000000000000e+00 -CostFunction: Nonlinear J = 1.6531252665872018e+03 -RPCGMinimizer: reduction in residual norm = 1.2006776473098589e-01 +CostFunction: Nonlinear J = 1.6531343291824785e+03 +RPCGMinimizer: reduction in residual norm = 1.2414458192527290e-01 CostFunction::addIncrement: Analysis: Valid time: 2018-04-15T00:00:00Z sea_water_cell_thickness min=0.0009999999999977 max=1345.6400000000003274 mean=112.4337347359288941 - sea_water_salinity min=0.0000000000000000 max=38.6959149218921183 mean=29.0069721693647047 -sea_water_potential_temperature min=-1.8883899372702533 max=31.0125853407566474 mean=5.2973741711447007 + sea_water_salinity min=0.0000000000000000 max=38.6959149218921183 mean=29.0074916345918368 +sea_water_potential_temperature min=-1.8883899372702533 max=31.7246904174557400 mean=5.3128599069475166 eastward_sea_water_velocity min=-0.6975788196638589 max=0.4840625031242498 mean=0.0010009071341291 northward_sea_water_velocity min=-0.7661101215480253 max=0.9402890903466523 mean=0.0021201094019194 - sea_surface_height_above_geoid min=-2.0936283531745894 max=0.8190071369678245 mean=-0.2049739891582961 + sea_surface_height_above_geoid min=-2.1508977733523640 max=0.8678867820476882 mean=-0.2004871633702467 ocean_mixed_layer_thickness min=0.0005000000000000 max=4593.1523423819935488 mean=165.9225478094618893 sea_water_depth min=0.0005000000000000 max=5587.8978052886013757 mean=1053.5232951172040430 -CostJo : Nonlinear Jo(SeaSurfaceTemp) = 275.5343698338789409, nobs = 113, Jo/n = 2.4383572551670705, err = 0.3665565623931296 -CostJo : Nonlinear Jo(SeaSurfaceSalinity) = 92.7854784143652296, nobs = 49, Jo/n = 1.8935811921299026, err = 1.0000000000000000 -CostJo : Nonlinear Jo(ADT) = 178.4631772821379343, nobs = 89, Jo/n = 2.0052042391251454, err = 0.1000000014901161 -CostJo : Nonlinear Jo(InsituTemperature) = 280.3382013630258029, nobs = 203, Jo/n = 1.3809763613942159, err = 0.9003481385241890 -CostJo : Nonlinear Jo(InsituSalinity) = 119.7825517174736092, nobs = 218, Jo/n = 0.5494612464104294, err = 0.6081969954577753 -CostJb : Nonlinear Jb = 29.8302471366160731 -CostFunction: Nonlinear J = 976.7340257474976397 +CostJo : Nonlinear Jo(SeaSurfaceTemp) = 218.3295992979803088, nobs = 113, Jo/n = 1.9321203477697373, err = 0.3665565623931296 +CostJo : Nonlinear Jo(SeaSurfaceSalinity) = 91.0601724565712800, nobs = 49, Jo/n = 1.8583708664606384, err = 1.0000000000000000 +CostJo : Nonlinear Jo(ADT) = 169.9471585510565319, nobs = 89, Jo/n = 1.9095186354051297, err = 0.1000000014901161 +CostJo : Nonlinear Jo(InsituTemperature) = 420.0581423685579239, nobs = 203, Jo/n = 2.0692519328500389, err = 0.9003481385241890 +CostJo : Nonlinear Jo(InsituSalinity) = 121.1026570349107629, nobs = 218, Jo/n = 0.5555167753894990, err = 0.6081969954577753 +CostJb : Nonlinear Jb = 39.6288485214363888 +CostFunction: Nonlinear J = 1060.1265782305131324