Skip to content

Commit

Permalink
works, needs more testing so still wip
Browse files Browse the repository at this point in the history
  • Loading branch information
@guillaumevernieres
guillaumevernieres committed Jan 13, 2025
1 parent 004aef2 commit 68d8150
Show file tree
Hide file tree
Showing 3 changed files with 83 additions and 46 deletions.
45 changes: 28 additions & 17 deletions utils/soca/gdas_ens_handler.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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
Expand All @@ -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;
Expand All @@ -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);

Expand Down Expand Up @@ -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);

Expand Down Expand Up @@ -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;
Expand All @@ -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;
Expand All @@ -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;
Expand Down Expand Up @@ -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);
Expand All @@ -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;

Expand Down
80 changes: 53 additions & 27 deletions utils/soca/gdas_postprocincr.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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_() {
Expand Down Expand Up @@ -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

Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
}

// -----------------------------------------------------------------------------
Expand Down Expand Up @@ -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<eckit::LocalConfiguration> inputIncrConfig_;
eckit::LocalConfiguration inputIncrConfig_;
Expand Down
4 changes: 2 additions & 2 deletions utils/soca/gdas_soca_diagb.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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());

Expand Down

0 comments on commit 68d8150

Please sign in to comment.