From d1c9a069463a92760e522a9a499edcc25f98791b Mon Sep 17 00:00:00 2001 From: Anna Shlyaeva Date: Thu, 2 Apr 2026 15:04:33 -0600 Subject: [PATCH 1/3] Towards obsoleting VariableChange::inverse, updating yamls and references --- test/testinput/3dvar.yml | 1 + test/testinput/3dvar_lowres.yml | 2 +- test/testinput/3dvar_nicas.yml | 1 + test/testinput/3dvarfgat_pseudo.yml | 1 + test/testinput/4dvar_htlm.yml | 1 + test/testinput/4dvar_identity.yml | 1 + test/testref/3dvar_nicas.test | 4 ++-- test/testref/3dvarfgat_pseudo.test | 4 ++-- test/testref/4dvar_htlm.test | 4 ++-- test/testref/4dvar_identity.test | 4 ++-- 10 files changed, 14 insertions(+), 9 deletions(-) diff --git a/test/testinput/3dvar.yml b/test/testinput/3dvar.yml index 9fda7153a..c859c2dbc 100644 --- a/test/testinput/3dvar.yml +++ b/test/testinput/3dvar.yml @@ -52,6 +52,7 @@ cost function: - sea_water_depth background error: + full inverse: true covariance model: SABER saber central block: saber block name: diffusion diff --git a/test/testinput/3dvar_lowres.yml b/test/testinput/3dvar_lowres.yml index 398780299..24111d29a 100644 --- a/test/testinput/3dvar_lowres.yml +++ b/test/testinput/3dvar_lowres.yml @@ -54,7 +54,7 @@ cost function: background error: covariance model: SABER change background resolution: true - + full inverse: true saber central block: saber block name: diffusion read: diff --git a/test/testinput/3dvar_nicas.yml b/test/testinput/3dvar_nicas.yml index f73a41e16..48e3e0ebb 100644 --- a/test/testinput/3dvar_nicas.yml +++ b/test/testinput/3dvar_nicas.yml @@ -39,6 +39,7 @@ cost function: background error: covariance model: SABER + full inverse: true saber central block: saber block name: BUMP_NICAS read: diff --git a/test/testinput/3dvarfgat_pseudo.yml b/test/testinput/3dvarfgat_pseudo.yml index a31d4fcf5..90be9b8c9 100644 --- a/test/testinput/3dvarfgat_pseudo.yml +++ b/test/testinput/3dvarfgat_pseudo.yml @@ -66,6 +66,7 @@ cost function: background error: covariance model: SABER + full inverse: true saber central block: saber block name: diffusion read: diff --git a/test/testinput/4dvar_htlm.yml b/test/testinput/4dvar_htlm.yml index 4e964eec2..877669532 100644 --- a/test/testinput/4dvar_htlm.yml +++ b/test/testinput/4dvar_htlm.yml @@ -55,6 +55,7 @@ cost function: background error: covariance model: SABER + full inverse: true saber central block: saber block name: diffusion read: diff --git a/test/testinput/4dvar_identity.yml b/test/testinput/4dvar_identity.yml index 0fb4549ed..5d7bd67f7 100644 --- a/test/testinput/4dvar_identity.yml +++ b/test/testinput/4dvar_identity.yml @@ -55,6 +55,7 @@ cost function: background error: covariance model: SABER + full inverse: true saber central block: saber block name: diffusion read: diff --git a/test/testref/3dvar_nicas.test b/test/testref/3dvar_nicas.test index 77ddc1569..5e2fd752c 100644 --- a/test/testref/3dvar_nicas.test +++ b/test/testref/3dvar_nicas.test @@ -21,5 +21,5 @@ CostJo : Nonlinear Jo(SeaSurfaceSalinity) = 92.7854784143652296, nobs = 49, Jo 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 +CostJb : Nonlinear Jb = 29.9035929242882439 +CostFunction: Nonlinear J = 976.8080560330458866 diff --git a/test/testref/3dvarfgat_pseudo.test b/test/testref/3dvarfgat_pseudo.test index 31233c0e8..20876afd1 100644 --- a/test/testref/3dvarfgat_pseudo.test +++ b/test/testref/3dvarfgat_pseudo.test @@ -21,5 +21,5 @@ CostJo : Nonlinear Jo(SeaSurfaceSalinity) = 46.9330495529855583, nobs = 50, Jo CostJo : Nonlinear Jo(ADT) = 82.2554700491828896, nobs = 94, Jo/n = 0.8750581920125839, err = 0.1000000014901161 CostJo : Nonlinear Jo(InsituTemperature) = 364.1365038607749511, nobs = 206, Jo/n = 1.7676529313629852, err = 0.8991441572394548 CostJo : Nonlinear Jo(InsituSalinity) = 58.4283558070863904, nobs = 218, Jo/n = 0.2680199807664513, err = 0.6081969954577753 -CostJb : Nonlinear Jb = 28.9887911055684242 -CostFunction: Nonlinear J = 1070.9592406972394656 +CostJb : Nonlinear Jb = 29.0442063676153452 +CostFunction: Nonlinear J = 1071.0149672573315911 diff --git a/test/testref/4dvar_htlm.test b/test/testref/4dvar_htlm.test index 8be0170cd..572c0dc07 100644 --- a/test/testref/4dvar_htlm.test +++ b/test/testref/4dvar_htlm.test @@ -13,5 +13,5 @@ sea_water_potential_temperature min=-1.8883899372702533 max=31.7004645720658 ocean_mixed_layer_thickness min=2.2854716757984130 max=4593.1533423819937525 mean=192.4109073940401515 sea_water_depth min=2.2854716757984130 max=5658.3057467114012979 mean=1200.5229536158342398 CostJo : Nonlinear Jo(SeaSufaceTemp) = 594.7782934004992512, nobs = 90, Jo/n = 6.6086477044499921, err = 0.3727400724804861 -CostJb : Nonlinear Jb = 52.6241435560154684 -CostFunction: Nonlinear J = 647.4024369565147481 +CostJb : Nonlinear Jb = 52.6119624335335203 +CostFunction: Nonlinear J = 647.3902558340328142 diff --git a/test/testref/4dvar_identity.test b/test/testref/4dvar_identity.test index c0fd596d6..fc3cbeea7 100644 --- a/test/testref/4dvar_identity.test +++ b/test/testref/4dvar_identity.test @@ -13,5 +13,5 @@ sea_water_potential_temperature min=-1.8883899372702533 max=31.7004645720658 ocean_mixed_layer_thickness min=2.2854716757984130 max=4593.1533423819937525 mean=192.4109073940401515 sea_water_depth min=2.2854716757984130 max=5658.3057467114012979 mean=1200.5229536158342398 CostJo : Nonlinear Jo(SeaSufaceTemp) = 597.2848931579876535, nobs = 90, Jo/n = 6.6364988128665292, err = 0.3727400724804861 -CostJb : Nonlinear Jb = 48.8761841139292841 -CostFunction: Nonlinear J = 646.1610772719169518 +CostJb : Nonlinear Jb = 48.8661052114184997 +CostFunction: Nonlinear J = 646.1509983694011225 From c83abd76e0aa8be4fd8780ea5d2308ab57b93624 Mon Sep 17 00:00:00 2001 From: Anna Date: Fri, 10 Apr 2026 15:33:04 -0600 Subject: [PATCH 2/3] Remove obsoleted interfaces --- .../LinearVariableChange/Balance/Balance.h | 4 +- .../Base/LinearVariableChangeBase.h | 2 - .../LinearModel2GeoVaLs.cc | 16 ------- .../LinearModel2GeoVaLs/LinearModel2GeoVaLs.h | 2 - .../LinearVariableChange.cc | 42 ------------------- .../LinearVariableChange.h | 2 - .../VariableChange/Base/VariableChangeBase.h | 1 - .../VariableChange/Model2Ana/Model2Ana.cc | 2 - src/soca/VariableChange/Model2Ana/Model2Ana.h | 2 +- .../Model2GeoVaLs/Model2GeoVaLs.cc | 6 --- .../Model2GeoVaLs/Model2GeoVaLs.h | 1 - .../VariableChange/Soca2Cice/Soca2Cice.cc | 6 --- src/soca/VariableChange/Soca2Cice/Soca2Cice.h | 1 - .../Soca2Cice/soca_soca2cice_mod.F90 | 3 -- src/soca/VariableChange/VariableChange.cc | 8 ---- src/soca/VariableChange/VariableChange.h | 1 - 16 files changed, 3 insertions(+), 96 deletions(-) diff --git a/src/soca/LinearVariableChange/Balance/Balance.h b/src/soca/LinearVariableChange/Balance/Balance.h index daa785f34..c420d9373 100644 --- a/src/soca/LinearVariableChange/Balance/Balance.h +++ b/src/soca/LinearVariableChange/Balance/Balance.h @@ -40,9 +40,9 @@ class Balance: public LinearVariableChangeBase { /// Perform linear transforms void multiply(const Increment &, Increment &) const override; - void multiplyInverse(const Increment &, Increment &) const override; + void multiplyInverse(const Increment &, Increment &) const; void multiplyAD(const Increment &, Increment &) const override; - void multiplyInverseAD(const Increment &, Increment &) const override; + void multiplyInverseAD(const Increment &, Increment &) const; private: void print(std::ostream &) const override; diff --git a/src/soca/LinearVariableChange/Base/LinearVariableChangeBase.h b/src/soca/LinearVariableChange/Base/LinearVariableChangeBase.h index cda532cb9..935fc6917 100644 --- a/src/soca/LinearVariableChange/Base/LinearVariableChangeBase.h +++ b/src/soca/LinearVariableChange/Base/LinearVariableChangeBase.h @@ -58,9 +58,7 @@ class LinearVariableChangeBase : public util::Printable, LinearVariableChangeBase() {} virtual ~LinearVariableChangeBase() {} virtual void multiply(const Increment &, Increment &) const = 0; - virtual void multiplyInverse(const Increment &, Increment &) const = 0; virtual void multiplyAD(const Increment &, Increment &) const = 0; - virtual void multiplyInverseAD(const Increment &, Increment &) const = 0; private: virtual void print(std::ostream &) const = 0; diff --git a/src/soca/LinearVariableChange/LinearModel2GeoVaLs/LinearModel2GeoVaLs.cc b/src/soca/LinearVariableChange/LinearModel2GeoVaLs/LinearModel2GeoVaLs.cc index 15b7e3238..68f879ecf 100644 --- a/src/soca/LinearVariableChange/LinearModel2GeoVaLs/LinearModel2GeoVaLs.cc +++ b/src/soca/LinearVariableChange/LinearModel2GeoVaLs/LinearModel2GeoVaLs.cc @@ -69,14 +69,6 @@ void LinearModel2GeoVaLs::multiply(const Increment &dxin, // ----------------------------------------------------------------------------- -void LinearModel2GeoVaLs::multiplyInverse(const Increment &dxin, - Increment &dxout) const { - util::Timer timer("soca::LinearModel2GeoVaLs", "multiplyInverse"); - multiply(dxin, dxout); -} - -// ----------------------------------------------------------------------------- - void LinearModel2GeoVaLs::multiplyAD(const Increment &dxin, Increment &dxout) const { util::Timer timer("soca::LinearModel2GeoVaLs", "multiplyAD"); @@ -106,14 +98,6 @@ void LinearModel2GeoVaLs::multiplyAD(const Increment &dxin, // ----------------------------------------------------------------------------- -void LinearModel2GeoVaLs::multiplyInverseAD(const Increment &dxin, - Increment &dxout) const { - util::Timer timer("soca::LinearModel2GeoVaLs", "multiplyInverseAD"); - multiplyAD(dxin, dxout); -} - -// ----------------------------------------------------------------------------- - void LinearModel2GeoVaLs::print(std::ostream & os) const { os << "SOCA linear change variable: LinearModel2GeoVaLs"; } diff --git a/src/soca/LinearVariableChange/LinearModel2GeoVaLs/LinearModel2GeoVaLs.h b/src/soca/LinearVariableChange/LinearModel2GeoVaLs/LinearModel2GeoVaLs.h index ec301b495..2df59eed7 100644 --- a/src/soca/LinearVariableChange/LinearModel2GeoVaLs/LinearModel2GeoVaLs.h +++ b/src/soca/LinearVariableChange/LinearModel2GeoVaLs/LinearModel2GeoVaLs.h @@ -33,9 +33,7 @@ class LinearModel2GeoVaLs: public LinearVariableChangeBase { ~LinearModel2GeoVaLs(); void multiply(const Increment &, Increment &) const; - void multiplyInverse(const Increment &, Increment &) const; void multiplyAD(const Increment &, Increment &) const; - void multiplyInverseAD(const Increment &, Increment &) const; private: const Geometry & geom_; diff --git a/src/soca/LinearVariableChange/LinearVariableChange.cc b/src/soca/LinearVariableChange/LinearVariableChange.cc index 7d47e36c3..576556314 100644 --- a/src/soca/LinearVariableChange/LinearVariableChange.cc +++ b/src/soca/LinearVariableChange/LinearVariableChange.cc @@ -109,27 +109,6 @@ void LinearVariableChange::changeVarTL(Increment & dx, // ----------------------------------------------------------------------------- -void LinearVariableChange::changeVarInverseTL(Increment & dx, - const oops::Variables & vars) const { - Log::trace() << "LinearVariableChange::multiplyInverse starting" - << vars << std::endl; - - // Create output state - Increment dxout(dx.geometry(), vars, dx.validTime()); - - // Call variable change(s) - for (ircst_ it = linVarChas_.rbegin(); it != linVarChas_.rend(); ++it) { - dxout.zero(); - it->multiplyInverse(dx, dxout); - dx.updateFields(vars); - dx = dxout; - } - - Log::trace() << "LinearVariableChange::multiplyInverse done" << std::endl; -} - -// ----------------------------------------------------------------------------- - void LinearVariableChange::changeVarAD(Increment & dx, const oops::Variables & vars) const { Log::trace() << "LinearVariableChange::multiplyAD starting" << std::endl; @@ -148,27 +127,6 @@ void LinearVariableChange::changeVarAD(Increment & dx, // ----------------------------------------------------------------------------- -void LinearVariableChange::changeVarInverseAD(Increment & dx, - const oops::Variables & vars) const { - Log::trace() << "LinearVariableChange::multiplyInverseAD starting" - << std::endl; - - // Create output state - Increment dxout(dx.geometry(), vars, dx.validTime()); - - // Call variable change(s) - for (icst_ it = linVarChas_.begin(); it != linVarChas_.end(); ++it) { - dxout.zero(); - it->multiplyInverseAD(dx, dxout); - dx.updateFields(vars); - dx = dxout; - } - - Log::trace() << "LinearVariableChange::multiplyInverseAD done" << std::endl; -} - -// ----------------------------------------------------------------------------- - void LinearVariableChange::print(std::ostream & os) const { for (icst_ it = linVarChas_.begin(); it != linVarChas_.end(); ++it) { os << *it; diff --git a/src/soca/LinearVariableChange/LinearVariableChange.h b/src/soca/LinearVariableChange/LinearVariableChange.h index fa43d2601..d5005e3d1 100644 --- a/src/soca/LinearVariableChange/LinearVariableChange.h +++ b/src/soca/LinearVariableChange/LinearVariableChange.h @@ -59,9 +59,7 @@ class LinearVariableChange : public util::Printable { void changeVarTraj(const State &, const oops::Variables &); 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; private: void print(std::ostream &) const override; diff --git a/src/soca/VariableChange/Base/VariableChangeBase.h b/src/soca/VariableChange/Base/VariableChangeBase.h index 97896d2ab..9f57020a3 100644 --- a/src/soca/VariableChange/Base/VariableChangeBase.h +++ b/src/soca/VariableChange/Base/VariableChangeBase.h @@ -56,7 +56,6 @@ class VariableChangeBase : public util::Printable, private boost::noncopyable { virtual ~VariableChangeBase() {} virtual void changeVar(const State &, State &) const = 0; - virtual void changeVarInverse(const State &, State &) const = 0; virtual const std::string classname() = 0; diff --git a/src/soca/VariableChange/Model2Ana/Model2Ana.cc b/src/soca/VariableChange/Model2Ana/Model2Ana.cc index 46b8462c0..687283803 100644 --- a/src/soca/VariableChange/Model2Ana/Model2Ana.cc +++ b/src/soca/VariableChange/Model2Ana/Model2Ana.cc @@ -54,7 +54,6 @@ void Model2Ana::changeVar(const State & xm, std::endl; util::Timer timer("soca::Model2Ana", "changeVar"); - util::DateTime * vtime = &xa.validTime(); xa = xm; // Rotate from the logical grid to meridional/zonal @@ -81,7 +80,6 @@ void Model2Ana::changeVarInverse(const State & xa, std::endl; util::Timer timer("soca::Model2Ana", "changeVarInverse"); - util::DateTime * vtime = &xm.validTime(); xm = xa; // Rotate from meridional/zonal to the logical grid diff --git a/src/soca/VariableChange/Model2Ana/Model2Ana.h b/src/soca/VariableChange/Model2Ana/Model2Ana.h index d4c58944b..88ee2edc1 100644 --- a/src/soca/VariableChange/Model2Ana/Model2Ana.h +++ b/src/soca/VariableChange/Model2Ana/Model2Ana.h @@ -37,7 +37,7 @@ class Model2Ana: public VariableChangeBase { ~Model2Ana(); void changeVar(const State &, State &) const override; - void changeVarInverse(const State &, State &) const override; + void changeVarInverse(const State &, State &) const; std::vector initRotate(const eckit::Configuration & conf, const std::string & uv) const diff --git a/src/soca/VariableChange/Model2GeoVaLs/Model2GeoVaLs.cc b/src/soca/VariableChange/Model2GeoVaLs/Model2GeoVaLs.cc index 9622daff7..1c7c4cb5c 100644 --- a/src/soca/VariableChange/Model2GeoVaLs/Model2GeoVaLs.cc +++ b/src/soca/VariableChange/Model2GeoVaLs/Model2GeoVaLs.cc @@ -46,10 +46,4 @@ void Model2GeoVaLs::changeVar(const State & xin, State & xout) const { // ----------------------------------------------------------------------------- -void Model2GeoVaLs::changeVarInverse(const State &, State &) const { - util::abor1_cpp("Model2GeoVaLs::changeVarInverse not implemented"); -} - -// ----------------------------------------------------------------------------- - } // namespace soca diff --git a/src/soca/VariableChange/Model2GeoVaLs/Model2GeoVaLs.h b/src/soca/VariableChange/Model2GeoVaLs/Model2GeoVaLs.h index af5523321..fc96ad199 100644 --- a/src/soca/VariableChange/Model2GeoVaLs/Model2GeoVaLs.h +++ b/src/soca/VariableChange/Model2GeoVaLs/Model2GeoVaLs.h @@ -22,7 +22,6 @@ class Model2GeoVaLs: public VariableChangeBase { ~Model2GeoVaLs(); void changeVar(const State &, State &) const override; - void changeVarInverse(const State &, State &) const override; private: const Geometry & geom_; diff --git a/src/soca/VariableChange/Soca2Cice/Soca2Cice.cc b/src/soca/VariableChange/Soca2Cice/Soca2Cice.cc index 419f756e1..094d07248 100644 --- a/src/soca/VariableChange/Soca2Cice/Soca2Cice.cc +++ b/src/soca/VariableChange/Soca2Cice/Soca2Cice.cc @@ -70,10 +70,4 @@ void Soca2Cice::changeVar(const State & xin, State & xout) const // ----------------------------------------------------------------------------- -void Soca2Cice::changeVarInverse(const State &, State &) const { - util::Timer timer("soca::Soca2Cice", "changeVarInverse"); -} - -// ----------------------------------------------------------------------------- - } // namespace soca diff --git a/src/soca/VariableChange/Soca2Cice/Soca2Cice.h b/src/soca/VariableChange/Soca2Cice/Soca2Cice.h index 7849e9215..c89f1a47a 100644 --- a/src/soca/VariableChange/Soca2Cice/Soca2Cice.h +++ b/src/soca/VariableChange/Soca2Cice/Soca2Cice.h @@ -104,7 +104,6 @@ class Soca2Cice: public VariableChangeBase { ~Soca2Cice(); void changeVar(const State &, State &) const override; - void changeVarInverse(const State &, State &) const override; private: const Geometry & geom_; diff --git a/src/soca/VariableChange/Soca2Cice/soca_soca2cice_mod.F90 b/src/soca/VariableChange/Soca2Cice/soca_soca2cice_mod.F90 index 5759506e5..eb334bae0 100644 --- a/src/soca/VariableChange/Soca2Cice/soca_soca2cice_mod.F90 +++ b/src/soca/VariableChange/Soca2Cice/soca_soca2cice_mod.F90 @@ -33,9 +33,6 @@ module soca_soca2cice_mod !! !! - forward: deaggregates a 2D analysis of sea-ice and inserts !! analysis in CICE restarts -!! - inverse: TODO(G), aggregates seaice variables along CICE sea-ice -!! categories, save the aggregated variables in a file -!! readable by soca type, public :: soca_soca2cice_params real(kind=kind_real) :: seaice_edge diff --git a/src/soca/VariableChange/VariableChange.cc b/src/soca/VariableChange/VariableChange.cc index eddb7d26c..35505c3b1 100644 --- a/src/soca/VariableChange/VariableChange.cc +++ b/src/soca/VariableChange/VariableChange.cc @@ -125,14 +125,6 @@ void VariableChange::changeVar(State & x, const oops::Variables & vars) const { // ----------------------------------------------------------------------------- -void VariableChange::changeVarInverse(State & x, - const oops::Variables & vars) const { - util::Timer timer("soca::VariableChange", "changeVarInverse"); - changeVar(x, vars); -} - -// ----------------------------------------------------------------------------- - void VariableChange::print(std::ostream & os) const { os << *variableChange_; } diff --git a/src/soca/VariableChange/VariableChange.h b/src/soca/VariableChange/VariableChange.h index 8c7b865fd..6a37be052 100644 --- a/src/soca/VariableChange/VariableChange.h +++ b/src/soca/VariableChange/VariableChange.h @@ -48,7 +48,6 @@ class VariableChange : public util::Printable { ~VariableChange(); void changeVar(State &, const oops::Variables &) const; - void changeVarInverse(State &, const oops::Variables &) const; private: void print(std::ostream &) const override; From fc9d63a777543aeca4f5f847c1e5a7d9741c5c17 Mon Sep 17 00:00:00 2001 From: Anna Date: Thu, 2 Jul 2026 10:26:44 -0600 Subject: [PATCH 3/3] Add InverseBalance LVC and use it in some tests that compute Jb --- .../LinearVariableChange/Balance/Balance.cc | 15 ---- .../LinearVariableChange/Balance/Balance.h | 2 - .../Balance/CMakeLists.txt | 2 + .../Balance/InverseBalance.cc | 75 +++++++++++++++++++ .../Balance/InverseBalance.h | 51 +++++++++++++ test/testinput/3dvar.yml | 17 ++++- test/testinput/3dvarfgat_pseudo.yml | 18 ++++- test/testref/3dvarfgat_pseudo.test | 4 +- 8 files changed, 163 insertions(+), 21 deletions(-) create mode 100644 src/soca/LinearVariableChange/Balance/InverseBalance.cc create mode 100644 src/soca/LinearVariableChange/Balance/InverseBalance.h diff --git a/src/soca/LinearVariableChange/Balance/Balance.cc b/src/soca/LinearVariableChange/Balance/Balance.cc index a49e1ba4d..070cc8fa3 100644 --- a/src/soca/LinearVariableChange/Balance/Balance.cc +++ b/src/soca/LinearVariableChange/Balance/Balance.cc @@ -60,13 +60,6 @@ namespace soca { soca_balance_mult_f90(keyFtnConfig_, dxa.toFortran(), dxm.toFortran()); } // ----------------------------------------------------------------------------- - void Balance::multiplyInverse(const Increment & dxm, Increment & dxa) const { - // dxa = K^-1 dxm - oops::Log::trace() << "soca::Balance::multiplyInverse " << std::endl; - util::Timer timer("soca::Balance", "multiplyInverse"); - soca_balance_multinv_f90(keyFtnConfig_, dxm.toFortran(), dxa.toFortran()); - } - // ----------------------------------------------------------------------------- void Balance::multiplyAD(const Increment & dxm, Increment & dxa) const { // dxa = K^T dxm oops::Log::trace() << "soca::Balance::multiplyAD " << std::endl; @@ -74,14 +67,6 @@ namespace soca { soca_balance_multad_f90(keyFtnConfig_, dxm.toFortran(), dxa.toFortran()); } // ----------------------------------------------------------------------------- - void Balance::multiplyInverseAD(const Increment & dxa, - Increment & dxm) const { - // dxm = (K^-1)^T dxa - oops::Log::trace() << "soca::Balance::multiplyInverseAD " << std::endl; - util::Timer timer("soca::Balance", "multiplyInverseAD"); - soca_balance_multinvad_f90(keyFtnConfig_, dxa.toFortran(), dxm.toFortran()); - } - // ----------------------------------------------------------------------------- void Balance::print(std::ostream & os) const { os << "SOCA linear change variable: Balance"; } diff --git a/src/soca/LinearVariableChange/Balance/Balance.h b/src/soca/LinearVariableChange/Balance/Balance.h index c420d9373..a63689064 100644 --- a/src/soca/LinearVariableChange/Balance/Balance.h +++ b/src/soca/LinearVariableChange/Balance/Balance.h @@ -40,9 +40,7 @@ class Balance: public LinearVariableChangeBase { /// Perform linear transforms void multiply(const Increment &, Increment &) const override; - void multiplyInverse(const Increment &, Increment &) const; void multiplyAD(const Increment &, Increment &) const override; - void multiplyInverseAD(const Increment &, Increment &) const; private: void print(std::ostream &) const override; diff --git a/src/soca/LinearVariableChange/Balance/CMakeLists.txt b/src/soca/LinearVariableChange/Balance/CMakeLists.txt index c1f1756e9..d337abe13 100644 --- a/src/soca/LinearVariableChange/Balance/CMakeLists.txt +++ b/src/soca/LinearVariableChange/Balance/CMakeLists.txt @@ -2,6 +2,8 @@ soca_target_sources( Balance.cc Balance.h BalanceFortran.h + InverseBalance.cc + InverseBalance.h soca_balance_mod.F90 soca_balance.interface.F90 soca_ksshts_mod.F90 diff --git a/src/soca/LinearVariableChange/Balance/InverseBalance.cc b/src/soca/LinearVariableChange/Balance/InverseBalance.cc new file mode 100644 index 000000000..74810facd --- /dev/null +++ b/src/soca/LinearVariableChange/Balance/InverseBalance.cc @@ -0,0 +1,75 @@ +/* + * (C) Copyright 2017-2021 UCAR. + * + * 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 + +#include "eckit/config/Configuration.h" + +#include "oops/util/Logger.h" +#include "oops/util/Timer.h" + +#include "soca/Geometry/Geometry.h" +#include "soca/Increment/Increment.h" +#include "soca/State/State.h" +#include "soca/LinearVariableChange/Balance/InverseBalance.h" +#include "soca/LinearVariableChange/Balance/BalanceFortran.h" + +using oops::Log; + +namespace soca { + + // ----------------------------------------------------------------------------- + + static LinearVariableChangeMaker + makerLinearVariableChangeBalance_("Inverse BalanceSOCA"); + + // ----------------------------------------------------------------------------- + InverseBalance::InverseBalance(const State & bkg, + const State & traj, + const Geometry & geom, + const eckit::Configuration & conf) { + oops::Log::trace() << "soca::InverseBalance::setup " << std::endl; + util::Timer timer("soca::InverseBalance", "InverseBalance"); + + const eckit::Configuration * configc = &conf; + + // Interpolate trajectory to the geom resolution + State traj_at_geomres(geom, traj); + + // Compute Jacobians of the balance wrt traj + soca_balance_setup_f90(keyFtnConfig_, + &configc, + traj_at_geomres.toFortran(), + geom.toFortran()); + } + // ----------------------------------------------------------------------------- + InverseBalance::~InverseBalance() { + oops::Log::trace() << "soca::InverseBalance::delete " << std::endl; + soca_balance_delete_f90(keyFtnConfig_); + } + // ----------------------------------------------------------------------------- + void InverseBalance::multiply(const Increment & dxm, Increment & dxa) const { + // dxa = K^-1 dxm + oops::Log::trace() << "soca::InverseBalance::multiply " << std::endl; + util::Timer timer("soca::InverseBalance", "multiply"); + soca_balance_multinv_f90(keyFtnConfig_, dxm.toFortran(), dxa.toFortran()); + } + // ----------------------------------------------------------------------------- + void InverseBalance::multiplyAD(const Increment & dxa, + Increment & dxm) const { + // dxm = (K^-1)^T dxa + oops::Log::trace() << "soca::InverseBalance::multiplyAD " << std::endl; + util::Timer timer("soca::InverseBalance", "multiplyAD"); + soca_balance_multinvad_f90(keyFtnConfig_, dxa.toFortran(), dxm.toFortran()); + } + // ----------------------------------------------------------------------------- + void InverseBalance::print(std::ostream & os) const { + os << "SOCA linear change variable: inverse Balance"; + } + // ----------------------------------------------------------------------------- +} // namespace soca diff --git a/src/soca/LinearVariableChange/Balance/InverseBalance.h b/src/soca/LinearVariableChange/Balance/InverseBalance.h new file mode 100644 index 000000000..4a07ded70 --- /dev/null +++ b/src/soca/LinearVariableChange/Balance/InverseBalance.h @@ -0,0 +1,51 @@ +/* + * (C) Copyright 2017-2021 UCAR. + * + * 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 "oops/util/DateTime.h" +#include "oops/util/Printable.h" + +#include "soca/LinearVariableChange/Base/LinearVariableChangeBase.h" + +// Forward declarations +namespace eckit { + class Configuration; +} +namespace soca { + class State; + class Increment; + class Geometry; +} + +// ----------------------------------------------------------------------------- + +namespace soca { + +/// SOCA balance: inverse +class InverseBalance: public LinearVariableChangeBase { + public: + static const std::string classname() {return "soca::Balance";} + + explicit InverseBalance(const State &, const State &, + const Geometry &, const eckit::Configuration &); + ~InverseBalance(); + +/// Perform linear transforms + void multiply(const Increment &, Increment &) const override; + void multiplyAD(const Increment &, Increment &) const override; + + private: + void print(std::ostream &) const override; + int keyFtnConfig_; +}; +// ----------------------------------------------------------------------------- + +} // namespace soca diff --git a/test/testinput/3dvar.yml b/test/testinput/3dvar.yml index c859c2dbc..d61950ac0 100644 --- a/test/testinput/3dvar.yml +++ b/test/testinput/3dvar.yml @@ -52,7 +52,6 @@ cost function: - sea_water_depth background error: - full inverse: true covariance model: SABER saber central block: saber block name: diffusion @@ -121,6 +120,22 @@ cost function: dcdt: filename: data_static/72x35x25/dcdt.nc name: dcdt + inverse linear variable change: + input variables: *soca_an_vars + output variables: *soca_an_vars + linear variable changes: + + - linear variable change name: Inverse BalanceSOCA + kst: + dsdtmax: 0.1 + dsdzmin: 3.0e-6 + dtdzmin: 1.0e-6 + nlayers: 999 + ksshts: + nlayers: 10 + dcdt: + filename: data_static/72x35x25/dcdt.nc + name: dcdt observations: observers: diff --git a/test/testinput/3dvarfgat_pseudo.yml b/test/testinput/3dvarfgat_pseudo.yml index 90be9b8c9..cc0556756 100644 --- a/test/testinput/3dvarfgat_pseudo.yml +++ b/test/testinput/3dvarfgat_pseudo.yml @@ -66,7 +66,6 @@ cost function: background error: covariance model: SABER - full inverse: true saber central block: saber block name: diffusion read: @@ -104,6 +103,23 @@ cost function: filename: data_static/72x35x25/dcdt.nc name: dcdt + inverse linear variable change: + input variables: *soca_an_vars + output variables: *soca_an_vars + linear variable changes: + + - linear variable change name: Inverse BalanceSOCA + kst: + dsdtmax: 0.1 + dsdzmin: 3.0e-6 + dtdzmin: 1.0e-6 + nlayers: 2 + ksshts: + nlayers: 2 + dcdt: + filename: data_static/72x35x25/dcdt.nc + name: dcdt + observations: observers: - obs space: diff --git a/test/testref/3dvarfgat_pseudo.test b/test/testref/3dvarfgat_pseudo.test index 20876afd1..31233c0e8 100644 --- a/test/testref/3dvarfgat_pseudo.test +++ b/test/testref/3dvarfgat_pseudo.test @@ -21,5 +21,5 @@ CostJo : Nonlinear Jo(SeaSurfaceSalinity) = 46.9330495529855583, nobs = 50, Jo CostJo : Nonlinear Jo(ADT) = 82.2554700491828896, nobs = 94, Jo/n = 0.8750581920125839, err = 0.1000000014901161 CostJo : Nonlinear Jo(InsituTemperature) = 364.1365038607749511, nobs = 206, Jo/n = 1.7676529313629852, err = 0.8991441572394548 CostJo : Nonlinear Jo(InsituSalinity) = 58.4283558070863904, nobs = 218, Jo/n = 0.2680199807664513, err = 0.6081969954577753 -CostJb : Nonlinear Jb = 29.0442063676153452 -CostFunction: Nonlinear J = 1071.0149672573315911 +CostJb : Nonlinear Jb = 28.9887911055684242 +CostFunction: Nonlinear J = 1070.9592406972394656