Skip to content

add linear variable change in GetValuesParameters #26

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 2 commits into
base: develop
Choose a base branch
from
Draft

add linear variable change in GetValuesParameters #26

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
3f9f101
add linearvariablechange group in GetValuesParameters and use it in…
TingLei-daprediction Jul 16, 2025
a7c4057
WIP : a mpasjedi Ens3DVar using tv as an analysis variable could run…
TingLei-daprediction Aug 4, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions src/oops/assimilation/BMatrix.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
#include "oops/assimilation/ControlIncrement.h"
#include "oops/assimilation/CostFunction.h"
#include "oops/base/ObsAuxCovariances.h"
#include "oops/util/Logger.h"

namespace oops {

Expand All @@ -33,7 +34,13 @@ template<typename MODEL, typename OBS> class BMatrix : private boost::noncopyabl
public:
explicit BMatrix(const CostFct_ & j): j_(j) {}
void multiply(const CtrlInc_ & x, CtrlInc_ & bx) const {
//tothink
Log::trace() << "BMatrix:: multiply start" <<std::endl;
Log::trace() << "BMatrix:: multiply start bx " <<bx<<std::endl;
Log::trace() << "BMatrix:: multiply start const x " <<x<<std::endl;
j_.jb().multiplyB(x, bx);
Log::trace() << "BMatrix:: multiply start after bx " <<bx<<std::endl;
Log::trace() << "BMatrix:: multiply start after const x " <<x<<std::endl;
}

/// Return ObsBias block of control variable error covariances
Expand Down
5 changes: 4 additions & 1 deletion src/oops/assimilation/CostFunction.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -280,7 +280,8 @@ double CostFunction<MODEL, OBS>::evaluate(CtrlVar_ & fguess,

template<typename MODEL, typename OBS>
void CostFunction<MODEL, OBS>::computeGradientFG(CtrlInc_ & grad) const {
Log::trace() << "CostFunction::computeGradientFG start" << std::endl;
Log::trace() << "CostFunction::computeGradientFG start clt" << std::endl;
Log::trace() << "CostFunction::computeGradientFG start clt grad.var " << grad.state().variables()<<std::endl;
PostProcessor<Increment_> pp;
PostProcessorTLAD<MODEL> costad;
this->zeroAD(grad);
Expand All @@ -290,7 +291,9 @@ void CostFunction<MODEL, OBS>::computeGradientFG(CtrlInc_ & grad) const {
jterms_[jj]->computeCostAD(tmp, grad, costad);
}

Log::trace() << "CostFunction::computeGradientFG start clt beforerunADJ grad " << grad.state().variables()<<std::endl;
this->runADJ(grad, costad, pp);
Log::trace() << "CostFunction::computeGradientFG start clt afterunADJ grad " << grad.state().variables()<<std::endl;

for (size_t jj = 0; jj < jterms_.size(); ++jj) {
jterms_[jj]->setPostProcAD();
Expand Down
22 changes: 22 additions & 0 deletions src/oops/assimilation/DRPCGMinimizer.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -172,14 +172,17 @@ double DRPCGMinimizer<MODEL, OBS>::solve(CtrlInc_ & dx, CtrlInc_ & dxh, CtrlInc_
scals.push_back(1.0/dotRr0);

Log::info() << std::endl;
Log::trace() << " DRPCG Starting Iteration0 0.0dxh" <<dxh<< std::endl;
for (int jiter = 0; jiter < maxiter; ++jiter) {
Log::info() << " DRPCG Starting Iteration " << jiter+1 << std::endl;
Log::trace() << " DRPCG Starting Iteration0 " << jiter << std::endl;

if (jiter < 5 || (jiter + 1) % 5 == 0 || jiter + 1 == maxiter) {
util::update_workflow_meter("iteration", jiter+1);
util::printRunStats("DRPCG iteration " + std::to_string(jiter+1));
}

Log::trace() << " DRPCG Starting Iteration0 1" << jiter << std::endl;
if (jiter > 0) {
// beta_{i} = r_{i+1}^T z_{i+1} / r_{i}^T z_{i}
double beta = rdots/rdots_old;
Expand All @@ -193,6 +196,8 @@ double DRPCGMinimizer<MODEL, OBS>::solve(CtrlInc_ & dx, CtrlInc_ & dxh, CtrlInc_
hh += pr;
}

Log::trace() << " DRPCG Starting Iteration0 0.7pp" <<pp<< jiter << std::endl;
Log::trace() << " DRPCG Starting Iteration0 0.7qq" <<qq<< jiter << std::endl;
// q_{i} = h_{i} + H^T R^{-1} H p_{i}
HtRinvH.multiply(pp, qq);
qq += hh;
Expand All @@ -202,10 +207,16 @@ double DRPCGMinimizer<MODEL, OBS>::solve(CtrlInc_ & dx, CtrlInc_ & dxh, CtrlInc_
double alpha = rdots/rho;

// dx_{i+1} = dx_{i} + alpha * p_{i}
Log::trace() << " DRPCG Starting Iteration0 1.9dx" <<dx<< jiter << std::endl;
Log::trace() << " DRPCG Starting Iteration0 1.9pp" <<pp<< jiter << std::endl;
dx.axpy(alpha, pp);
// dxh_{i+1} = dxh_{i} + alpha * h_{i} ! for diagnosing Jb
Log::trace() << " DRPCG Starting Iteration0 1.9hh" <<hh<< jiter << std::endl;
dxh.axpy(alpha, hh);
Log::trace() << " DRPCG Starting Iteration0 1.9dxh" <<dxh<< jiter << std::endl;
// r_{i+1} = r_{i} - alpha * q_{i}
Log::trace() << " DRPCG Starting Iteration0 1.9qq" <<qq<< jiter << std::endl;
Log::trace() << " DRPCG Starting Iteration0 1.9rr" <<rr<< jiter << std::endl;
rr.axpy(-alpha, qq);

// Compute the quadratic cost function
Expand All @@ -215,16 +226,27 @@ double DRPCGMinimizer<MODEL, OBS>::solve(CtrlInc_ & dx, CtrlInc_ & dxh, CtrlInc_
double costJb = costJ0Jb + dot_product(dx, gradJb) + 0.5 * dot_product(dx, dxh);
// Jo[dx_{i}] + Jc[dx_{i}] = J[dx_{i}] - Jb[dx_{i}]
double costJoJc = costJ - costJb;
Log::trace() << " DRPCG Starting Iteration0 2.0rr" <<rr<< jiter << std::endl;

// Re-orthogonalization
for (int jj = 0; jj < jiter; ++jj) {
double proj = scals[jj] * dot_product(rr, zvecs[jj]);
rr.axpy(-proj, vvecs[jj]);
}

Log::trace() << " DRPCG Starting Iteration0 2.5" << jiter << std::endl;
Log::trace() << " DRPCG Starting Iteration0 2.5pr" <<pr<< jiter << std::endl;
Log::trace() << " DRPCG Starting Iteration0 2.5rr" <<rr<< jiter << std::endl;
// z_{i+1} = B LMP r_{i+1}
lmp_.multiply(rr, pr);
Log::trace() << " DRPCG Starting Iteration0 3pr" << jiter << std::endl;
Log::trace() << pr << jiter << std::endl;
Log::trace() << " DRPCG Starting Iteration0 3zz" << jiter << std::endl;
Log::trace() << zz << jiter << std::endl;

Log::trace() << " DRPCG Starting Iteration0 3Before tothink B.multiply" << jiter << std::endl;
B.multiply(pr, zz);
Log::trace() << " DRPCG Starting Iteration0 3after B.multiply" << jiter << std::endl;

// r_{i}^T z_{i}
rdots_old = rdots;
Expand Down
6 changes: 6 additions & 0 deletions src/oops/base/EnsembleCovariance.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,11 @@ EnsembleCovariance<MODEL>::EnsembleCovariance(const Geometry_ & resol, const Var
ens_.reset(new IncrementSet<MODEL>(resol, vars, tmp, true));
*ens_ -= ens_->ens_mean();

for (size_t i = 0; i < ens_->size(); ++i) { //cltthink
auto member = (*ens_)[i];
Log::trace()<<"EnsembleCovariance::EnsembleCovariance start ith dx "<<i <<" "<<member<<std::endl;
}

// Setup ensemble transform
if (conf.has("ensemble transform")) {
const eckit::LocalConfiguration confEns(conf, "ensemble transform");
Expand Down Expand Up @@ -159,6 +164,7 @@ EnsembleCovariance<MODEL>::EnsembleCovariance(const Geometry_ & resol, const Var
}
loc_.reset(new Localization_(resol, locVars, locconf));
}
Log::trace() << "EnsembleCovariance::EnsembleCovariance vars " <<vars<< std::endl;

size_t current = eckit::system::ResourceUsage().maxResidentSetSize();
this->setObjectSize(current - init);
Expand Down
2 changes: 2 additions & 0 deletions src/oops/base/GetValuePosts.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,8 @@ class GetValuesParameters : public oops::Parameters {
public:
Parameter<eckit::LocalConfiguration> variableChange{"variable change",
eckit::LocalConfiguration(), this};
Parameter<eckit::LocalConfiguration> linearvariableChange{"linear variable change",
eckit::LocalConfiguration(), this};
};

/// \brief Fills GeoVaLs with requested variables at requested locations during model run
Expand Down
9 changes: 7 additions & 2 deletions src/oops/base/GetValueTLADs.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,15 +64,17 @@ class GetValueTLADs : public PostBaseTLAD<MODEL> {
std::vector<GetValPtr_> getvals_;
std::map<util::DateTime, CVarPtr_> chvartlad_;
const eckit::LocalConfiguration cvConf_;
const eckit::LocalConfiguration linearcvConf_;
};

// -----------------------------------------------------------------------------
template <typename MODEL, typename OBS>
GetValueTLADs<MODEL, OBS>::GetValueTLADs(const GetValuesParameters<MODEL> & params,
const util::DateTime & bgn, const util::DateTime & end)
: PostBaseTLAD<MODEL>(bgn, end), geovars_(), linvars_(), getvals_(), chvartlad_(),
cvConf_(params.variableChange.value())
linearcvConf_(params.linearvariableChange.value()),cvConf_(params.variableChange.value())
{
Log::trace() << "GetValueTLADs::GetValueTLADs clt linarcvConf_ " <<linearcvConf_ <<std::endl;
Log::trace() << "GetValueTLADs::GetValueTLADs" << std::endl;
}
// -----------------------------------------------------------------------------
Expand All @@ -82,6 +84,8 @@ void GetValueTLADs<MODEL, OBS>::append(GetValPtr_ getval) {
getvals_.push_back(getval);
geovars_ += getval->requiredVariables();
linvars_ += getval->linearVariables();
Log::trace() << "GetValuePosts::append start geovars_ " <<geovars_<< std::endl;
Log::trace() << "GetValuePosts::append start linvars_ " <<linvars_<< std::endl;
Log::trace() << "GetValuePosts::append done" << std::endl;
}

Expand Down Expand Up @@ -110,7 +114,7 @@ void GetValueTLADs<MODEL, OBS>::doProcessingTraj(const State_ & xx) {

for (GetValPtr_ getval : getvals_) getval->process(zz);

CVarPtr_ cvtlad(new LinearVariableChange<MODEL>(xx.geometry(), cvConf_));
CVarPtr_ cvtlad(new LinearVariableChange<MODEL>(xx.geometry(), linearcvConf_));
cvtlad->changeVarTraj(xx, linvars_);
chvartlad_[xx.validTime()] = std::move(cvtlad);

Expand Down Expand Up @@ -170,6 +174,7 @@ void GetValueTLADs<MODEL, OBS>::doProcessingAD(Increment_ & dx) {
Log::trace() << "GetValueTLADs::doProcessingAD start" << std::endl;
const util::DateTime now = dx.validTime();
ASSERT(chvartlad_.find(now) != chvartlad_.end());
Log::trace() << "GetValueTLADs::doProcessingAD linearvars" <<linvars_<<std::endl;
Increment_ dz(dx.geometry(), linvars_, dx.validTime());
dz.zero();

Expand Down
1 change: 1 addition & 0 deletions src/oops/base/GetValues.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -690,6 +690,7 @@ template <typename MODEL, typename OBS>
void GetValues<MODEL, OBS>::processAD(Increment_ & dx) {
Log::trace() << "GetValues::processAD start" << std::endl;
util::Timer timer("oops::GetValues", "processAD");
Log::trace() << "GetValues::processAD startclt variables " <<dx.variables()<< std::endl;

for (size_t jtask = 0; jtask < ntasks_; ++jtask) {
// Mask obs outside time slot
Expand Down
1 change: 1 addition & 0 deletions src/oops/base/ObserversTLAD.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -112,6 +112,7 @@ template <typename MODEL, typename OBS>
void ObserversTLAD<MODEL, OBS>::initializeTraj(const Geometry_ & geom, const ObsAuxCtrls_ & ybias,
PostProcTLAD_ & pp) {
Log::trace() << "ObserversTLAD<MODEL, OBS>::initializeTraj start" << std::endl;
Log::trace() << "ObserversTLAD<MODEL, OBS>::initializeTraj clt params " << getValuesParams_<<std::endl;
posts_.reset(new GetValueTLADs_(getValuesParams_, winbgn_, winend_));
for (size_t jj = 0; jj < observers_.size(); ++jj) {
if (observers_[jj]) {
Expand Down
1 change: 1 addition & 0 deletions src/oops/interface/Increment.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -155,6 +155,7 @@ Increment<MODEL>::Increment(const Geometry_ & resol, const Variables & vars,
{
Log::trace() << "Increment<MODEL>::Increment starting" << std::endl;
util::Timer timer(classname(), "Increment");
Log::trace() << "Increment<MODEL>::Increment vars " <<vars<< std::endl;
increment_.reset(new Increment_(resol.geometry(), vars, time));
this->setObjectSize(increment_->serialSize()*sizeof(double));
Log::trace() << "Increment<MODEL>::Increment done" << std::endl;
Expand Down