From 68d81508d9086d2561de415bfb96ae0d47112576 Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Mon, 13 Jan 2025 15:13:13 -0600 Subject: [PATCH] works, needs more testing so still wip --- utils/soca/gdas_ens_handler.h | 45 +++++++++++-------- utils/soca/gdas_postprocincr.h | 80 ++++++++++++++++++++++------------ utils/soca/gdas_soca_diagb.h | 4 +- 3 files changed, 83 insertions(+), 46 deletions(-) diff --git a/utils/soca/gdas_ens_handler.h b/utils/soca/gdas_ens_handler.h index 89c6e2a2c..741055ca1 100644 --- a/utils/soca/gdas_ens_handler.h +++ b/utils/soca/gdas_ens_handler.h @@ -75,14 +75,22 @@ namespace gdasapp { // ----------------------------------------------------------------------------- int execute(const eckit::Configuration & fullConfig) const { - // Setup the soca geometry + // Setup the native geometry of the ens. members, + // currently assumed to be the same as the deterministic const eckit::LocalConfiguration geomConfig(fullConfig, "geometry"); oops::Log::info() << "geometry: " << std::endl << geomConfig << std::endl; const soca::Geometry geom(geomConfig, this->getComm()); + // Setup the ensemble B (and output) soca geometry + oops::Log::info() << "====================== ens B geometry" << std::endl; + const std::string outputGeometryKey = fullConfig.has("output geometry") + ? "output geometry" // keep things backward compatible for now + : "geometry"; // and default to the input geometry + const eckit::LocalConfiguration geomOutConfig(fullConfig, outputGeometryKey); + const soca::Geometry geomOut(geomOutConfig, this->getComm()); // Initialize the post processing - PostProcIncr postProcIncr(fullConfig, geom, this->getComm()); + PostProcIncr postProcIncr(fullConfig, geom, this->getComm(), geomOut); oops::Log::info() << "soca increments: " << std::endl << postProcIncr.inputIncrConfig_ << std::endl; @@ -93,7 +101,9 @@ namespace gdasapp { oops::Log::info() << postProcIncr.inputIncrConfig_ << std::endl; // Assume z* output is the same for the trajectory and the state - ensMembers.push_back(postProcIncr.read(i)); + soca::Increment fullResIncr = postProcIncr.read(i); + soca::Increment lowResIncr(geomOut, fullResIncr); // interp to low resolution geometry + ensMembers.push_back(lowResIncr); } // Check if we need to recenter the increment around the deterministic @@ -103,9 +113,9 @@ namespace gdasapp { } // Compute ensemble moments - soca::Increment ensMean(geom, postProcIncr.socaIncrVar_, postProcIncr.dt_); - soca::Increment ensStd(geom, postProcIncr.socaIncrVar_, postProcIncr.dt_); - soca::Increment ensVariance(geom, postProcIncr.socaIncrVar_, postProcIncr.dt_); + soca::Increment ensMean(geomOut, postProcIncr.socaIncrVar_, postProcIncr.dt_); + soca::Increment ensStd(geomOut, postProcIncr.socaIncrVar_, postProcIncr.dt_); + soca::Increment ensVariance(geomOut, postProcIncr.socaIncrVar_, postProcIncr.dt_); gdasapp_ens_utils::ensMoments(ensMembers, ensMean, ensStd, ensVariance); oops::Log::info() << "mean: " << ensMean << std::endl; oops::Log::info() << "std: " << ensStd << std::endl; @@ -123,14 +133,15 @@ namespace gdasapp { // Initialize the trajectories used in the linear variable changes const eckit::LocalConfiguration trajConfig(fullConfig, "trajectory"); - soca::State determTraj(geom, trajConfig); // deterministic trajectory + soca::State determTrajNative(geom, trajConfig); // deterministic trajectory at full res + soca::State determTraj(geomOut, determTrajNative); // interpolated deterministic trajectory soca::State ensMeanTraj(determTraj); // trajectory defined as the ens. mean ensMeanTraj.zero(); ensMeanTraj += ensMean; // Compute the recentering increment as the difference between // the ensemble mean and the deterministic - soca::Increment recenteringIncr(geom, postProcIncr.socaIncrVar_, postProcIncr.dt_); + soca::Increment recenteringIncr(geomOut, postProcIncr.socaIncrVar_, postProcIncr.dt_); recenteringIncr.diff(determTraj, ensMeanTraj); postProcIncr.setToZero(recenteringIncr); @@ -192,7 +203,7 @@ namespace gdasapp { // Save total ssh oops::Log::info() << "ssh ensemble member " << i << std::endl; - soca::Increment ssh_tmp(geom, socaSshVar, postProcIncr.dt_); + soca::Increment ssh_tmp(geomOut, socaSshVar, postProcIncr.dt_); ssh_tmp = ensMembers[i]; sshTotal.push_back(ssh_tmp); @@ -257,9 +268,9 @@ namespace gdasapp { } // Compute ensemble moments for total ssh - soca::Increment sshMean(geom, socaSshVar, postProcIncr.dt_); - soca::Increment sshStd(geom, socaSshVar, postProcIncr.dt_); - soca::Increment sshTotalVariance(geom, socaSshVar, postProcIncr.dt_); + soca::Increment sshMean(geomOut, socaSshVar, postProcIncr.dt_); + soca::Increment sshStd(geomOut, socaSshVar, postProcIncr.dt_); + soca::Increment sshTotalVariance(geomOut, socaSshVar, postProcIncr.dt_); gdasapp_ens_utils::ensMoments(sshTotal, sshMean, sshStd, sshTotalVariance); oops::Log::info() << "mean ssh total: " << sshMean << std::endl; oops::Log::info() << "std ssh total: " << sshStd << std::endl; @@ -269,7 +280,7 @@ namespace gdasapp { // Compute ensemble moments for steric ssh sshMean.zero(); sshStd.zero(); - soca::Increment sshStericVariance(geom, socaSshVar, postProcIncr.dt_); + soca::Increment sshStericVariance(geomOut, socaSshVar, postProcIncr.dt_); gdasapp_ens_utils::ensMoments(sshSteric, sshMean, sshStd, sshStericVariance); oops::Log::info() << "mean steric ssh: " << sshMean << std::endl; oops::Log::info() << "std steric ssh: " << sshStd << std::endl; @@ -278,8 +289,8 @@ namespace gdasapp { // Compute ensemble moments for non-steric ssh sshMean.zero(); - soca::Increment sshNonStericVariance(geom, socaSshVar, postProcIncr.dt_); - soca::Increment sshNonStericStd(geom, socaSshVar, postProcIncr.dt_); + soca::Increment sshNonStericVariance(geomOut, socaSshVar, postProcIncr.dt_); + soca::Increment sshNonStericStd(geomOut, socaSshVar, postProcIncr.dt_); sshNonStericStd.zero(); gdasapp_ens_utils::ensMoments(sshNonSteric, sshMean, sshNonStericStd, sshNonStericVariance); oops::Log::info() << "mean non-steric ssh: " << sshMean << std::endl; @@ -308,7 +319,7 @@ namespace gdasapp { ensStd.write(bkgErrOutputConfig); // Explained variance by steric height R=1-SS(non-steric ssh)/SS(total ssh) - soca::Increment varianceRatio(geom, socaSshVar, postProcIncr.dt_); + soca::Increment varianceRatio(geomOut, socaSshVar, postProcIncr.dt_); varianceRatio = sshNonStericVariance; atlas::FieldSet varianceRatioFs; varianceRatio.toFieldSet(varianceRatioFs); @@ -319,7 +330,7 @@ namespace gdasapp { util::divideFieldSets(varianceRatioFs, sshTotalVarianceFs, sshTotalVarianceFs); varianceRatio.fromFieldSet(varianceRatioFs); - soca::Increment stericExplainedVariance(geom, socaSshVar, postProcIncr.dt_); + soca::Increment stericExplainedVariance(geomOut, socaSshVar, postProcIncr.dt_); stericExplainedVariance.ones(); stericExplainedVariance -= varianceRatio; diff --git a/utils/soca/gdas_postprocincr.h b/utils/soca/gdas_postprocincr.h index b8638a7b8..413e240ea 100644 --- a/utils/soca/gdas_postprocincr.h +++ b/utils/soca/gdas_postprocincr.h @@ -22,17 +22,38 @@ namespace gdasapp { +// ----------------------------------------------------------------------------- +/*! \class PostProcIncr + \brief This class handles the processing of increments in the GDAS application. + + The PostProcIncr class is responsible for managing the configuration and processing + of increments, including reading configuration parameters, initializing variables, + applying specified variable changes and handling input and output. +*/ +// ----------------------------------------------------------------------------- class PostProcIncr { public: // ----------------------------------------------------------------------------- - // Constructor - + // Constructors + + /** + * @brief Constructor for the PostProcIncr class. + * + * This constructor initializes the PostProcIncr object using the provided configuration, geometry, and MPI communicator. + * It sets up various parameters and configurations required for post-processing increments. + * + * @param fullConfig The full configuration object containing various settings. + * @param geom The native geometry of the increments output. + * @param comm The MPI communicator. + * @param geomProc The geometry to perform the processing on. + */ PostProcIncr(const eckit::Configuration & fullConfig, const soca::Geometry& geom, - const eckit::mpi::Comm & comm) + const eckit::mpi::Comm & comm, const soca::Geometry& geomProc) : dt_(getDate(fullConfig)), layerVar_(getLayerVar(fullConfig)), geom_(geom), - Layers_(getLayerThickness(fullConfig, geom)), + geomProc_(geomProc), + Layers_(getLayerThickness(fullConfig, geom, geomProc)), comm_(comm), ensSize_(1), pattern_() { @@ -66,6 +87,20 @@ class PostProcIncr { } } + /** + * @brief Constructor for the PostProcIncr class when the compute/processing geometry + * is the same as the native geometry. + * + * This constructor delegates the initialization to the main constructor + * + * @param fullConfig The full configuration object. + * @param geom The geometry object. + * @param comm The MPI communicator. + */ + PostProcIncr(const eckit::Configuration & fullConfig, const soca::Geometry & geom, + const eckit::mpi::Comm & comm) + : PostProcIncr(fullConfig, geom, comm, geom) {} + // ----------------------------------------------------------------------------- // Read ensemble member n @@ -75,7 +110,7 @@ class PostProcIncr { // initialize the soca increment soca::Increment socaIncr(geom_, socaIncrVar_, dt_); - eckit::LocalConfiguration memberConfig; // inputIncrConfig_); + eckit::LocalConfiguration memberConfig; memberConfig = inputIncrConfig_; // replace templated string if necessary @@ -88,8 +123,10 @@ class PostProcIncr { oops::Log::debug() << "-------------------- input increment: " << std::endl; oops::Log::debug() << socaIncr << std::endl; - return socaIncr; - } + + soca::Increment socaIncrOut(geomProc_, socaIncr); + return socaIncrOut; + } // ----------------------------------------------------------------------------- // Append variable to increment @@ -174,7 +211,7 @@ class PostProcIncr { const soca::State& xTraj) { oops::Log::info() << "==========================================" << std::endl; oops::Log::info() << "====== applying specified change of variables" << std::endl; - soca::LinearVariableChange lvc(this->geom_, lvcConfig); + soca::LinearVariableChange lvc(this->geomProc_, lvcConfig); lvc.changeVarTraj(xTraj, socaIncrVar_); lvc.changeVarTL(socaIncr, socaIncrVar_); oops::Log::info() << " in var change:" << socaIncr << std::endl; @@ -236,27 +273,15 @@ class PostProcIncr { } // Read the layer thickness from the relevant background soca::Increment getLayerThickness(const eckit::Configuration& fullConfig, - const soca::Geometry& geom) const { + const soca::Geometry& geom, + const soca::Geometry& geomProc) const { soca::Increment layerThick(geom, getLayerVar(fullConfig), getDate(fullConfig)); const eckit::LocalConfiguration vertGeomConfig(fullConfig, "vertical geometry"); layerThick.read(vertGeomConfig); - oops::Log::debug() << "layerThick: " << std::endl << layerThick << std::endl; - return layerThick; - } - // Initialize the trajectory - soca::State getTraj(const eckit::Configuration& fullConfig, - const soca::Geometry& geom) const { - if ( fullConfig.has("linear variable change") ) { - const eckit::LocalConfiguration trajConfig(fullConfig, "trajectory"); - soca::State traj(geom, trajConfig); - oops::Log::debug() << "traj:" << traj << std::endl; - return traj; - } else { - oops::Variables tmpVar(fullConfig, "layers variable"); - util::DateTime tmpDate(getDate(fullConfig)); - soca::State traj(geom, tmpVar, tmpDate); - return traj; - } + soca::Increment layerThickOut(geomProc, layerThick); + + oops::Log::debug() << "layerThickOut: " << std::endl << layerThickOut << std::endl; + return layerThickOut; } // ----------------------------------------------------------------------------- @@ -300,7 +325,8 @@ class PostProcIncr { util::DateTime dt_; // valid date of increment oops::Variables layerVar_; // layer variable const soca::Increment Layers_; // layer thicknesses - const soca::Geometry & geom_; + const soca::Geometry & geom_; // Native geometry + const soca::Geometry & geomProc_; // Geometry to perform processing on const eckit::mpi::Comm & comm_; // std::vector inputIncrConfig_; eckit::LocalConfiguration inputIncrConfig_; diff --git a/utils/soca/gdas_soca_diagb.h b/utils/soca/gdas_soca_diagb.h index 15cb7bce6..fc1ec08da 100644 --- a/utils/soca/gdas_soca_diagb.h +++ b/utils/soca/gdas_soca_diagb.h @@ -207,8 +207,8 @@ namespace gdasapp { // Setup the output soca geometry oops::Log::info() << "====================== output geometry" << std::endl; const std::string outputGeometryKey = fullConfig.has("output geometry") - ? "output geometry" // keep things backward compatible for now - : "geometry"; // and default to the input geometry + ? "output geometry" // keep things backward compatible for now + : "geometry"; // and default to the input geometry const eckit::LocalConfiguration geomOutConfig(fullConfig, outputGeometryKey); const soca::Geometry geomOut(geomOutConfig, this->getComm());