<commit_before>/* -*-  mode: c++; indent-tabs-mode: nil -*- */
/* -------------------------------------------------------------------------
   ATS

   License: see $ATS_DIR/COPYRIGHT
   Author: <NAME>

   Interface for a StrongMPC which uses a preconditioner in which the
   block-diagonal cell-local matrix is dense.  If the system looks something
   like:

   A( y1, y2, x, t ) = 0
   B( y1, y2, x, t ) = 0

   where y1,y2 are spatially varying unknowns that are discretized using the MFD
   method (and therefore have both cell and face unknowns), an approximation to
   the Jacobian is written as

   [  dA_c/dy1_c  dA_c/dy1_f   dA_c/dy2_c       0      ]
   [  dA_f/dy1_c  dA_f/dy1_f      0              0      ]
   [  dB_c/dy1_c     0          dB_c/dy2_c  dB_c/dy2_f ]
   [      0           0          dB_f/dy2_c  dB_f/dy2_f ]


   Note that the upper left block is the standard preconditioner for the A
   system, and the lower right block is the standard precon for the B system,
   and we have simply added cell-based couplings, dA_c/dy2_c and dB_c/dy1_c.

   In the temperature/pressure system, these correspond to d_water_content /
   d_temperature and d_energy / d_pressure.

   ------------------------------------------------------------------------- */

#include <fstream>
#include "EpetraExt_RowMatrixOut.h"

#include "LinearOperatorFactory.hh"
#include "FieldEvaluator.hh"
#include "Operator.hh"
#include "TreeOperator.hh"
#include "PDE_Accumulation.hh"

#include "mpc_coupled_cells.hh"

namespace Amanzi {

void MPCCoupledCells::Setup(const Teuchos::Ptr<State>& S) {
  StrongMPC<PK_PhysicalBDF_Default>::Setup(S);

  A_key_ = plist_->get<std::string>("conserved quantity A");
  B_key_ = plist_->get<std::string>("conserved quantity B");
  y1_key_ = plist_->get<std::string>("primary variable A");
  y2_key_ = plist_->get<std::string>("primary variable B");
  dA_dy2_key_ = std::string("d")+A_key_+std::string("_d")+y2_key_;
  dB_dy1_key_ = std::string("d")+B_key_+std::string("_d")+y1_key_;

  Key mesh_key = plist_->get<std::string>("domain name");
  mesh_ = S->GetMesh(mesh_key);

  // set up debugger
  db_ = sub_pks_[0]->debugger();

  // Get the sub-blocks from the sub-PK's preconditioners.
  Teuchos::RCP<Operators::Operator> pcA = sub_pks_[0]->preconditioner();
  Teuchos::RCP<Operators::Operator> pcB = sub_pks_[1]->preconditioner();

  // Create the combined operator
  Teuchos::RCP<TreeVectorSpace> tvs = Teuchos::rcp(new TreeVectorSpace());
  tvs->PushBack(Teuchos::rcp(new TreeVectorSpace(Teuchos::rcpFromRef(pcA->DomainMap()))));
  tvs->PushBack(Teuchos::rcp(new TreeVectorSpace(Teuchos::rcpFromRef(pcB->DomainMap()))));

  preconditioner_ = Teuchos::rcp(new Operators::TreeOperator(tvs));
  preconditioner_->SetOperatorBlock(0, 0, pcA);
  preconditioner_->SetOperatorBlock(1, 1, pcB);
  
  // create coupling blocks and push them into the preconditioner...
  S->RequireField(A_key_)->SetMesh(mesh_)->AddComponent("cell", AmanziMesh::CELL, 1);
  S->RequireField(y2_key_)->SetMesh(mesh_)->AddComponent("cell", AmanziMesh::CELL, 1);
  S->RequireFieldEvaluator(A_key_);
  S->RequireFieldEvaluator(y2_key_);

  if (!plist_->get<bool>("no dA/dy2 block", false)) {
    Teuchos::ParameterList& acc_pc_plist = plist_->sublist("dA_dy2 accumulation preconditioner");
    acc_pc_plist.set("entity kind", "cell");
    dA_dy2_ = Teuchos::rcp(new Operators::PDE_Accumulation(acc_pc_plist, mesh_));
    preconditioner_->SetOperatorBlock(0, 1, dA_dy2_->global_operator());
  }

  S->RequireField(B_key_)->SetMesh(mesh_)->AddComponent("cell", AmanziMesh::CELL, 1);
  S->RequireField(y1_key_)->SetMesh(mesh_)->AddComponent("cell", AmanziMesh::CELL, 1);
  S->RequireFieldEvaluator(B_key_);
  S->RequireFieldEvaluator(y1_key_);

  if (!plist_->get<bool>("no dB/dy1 block", false)) {
    Teuchos::ParameterList& acc_pc_plist = plist_->sublist("dB_dy1 accumulation preconditioner");
    acc_pc_plist.set("entity kind", "cell");
    dB_dy1_ = Teuchos::rcp(new Operators::PDE_Accumulation(acc_pc_plist, mesh_));
    preconditioner_->SetOperatorBlock(1, 0, dB_dy1_->global_operator());
  }

  // setup and initialize the preconditioner
  precon_used_ = plist_->isSublist("preconditioner");
  if (precon_used_) {
    preconditioner_->SymbolicAssembleMatrix();
    Teuchos::ParameterList& pc_sublist = plist_->sublist("preconditioner");
    preconditioner_->InitializePreconditioner(pc_sublist);
  }

  // setup and initialize the linear solver for the preconditioner
  if (plist_->isSublist("linear solver")) {
    Teuchos::ParameterList linsolve_sublist = plist_->sublist("linear solver");
    AmanziSolvers::LinearOperatorFactory<Operators::TreeOperator,TreeVector,TreeVectorSpace> fac;
    linsolve_preconditioner_ = fac.Create(linsolve_sublist, preconditioner_);
  } else {
    linsolve_preconditioner_ = preconditioner_;
  }
}


// updates the preconditioner
void MPCCoupledCells::UpdatePreconditioner(double t, Teuchos::RCP<const TreeVector> up,
        double h) {
  StrongMPC<PK_PhysicalBDF_Default>::UpdatePreconditioner(t,up,h);

  if (dA_dy2_ != Teuchos::null &&
      S_next_->GetFieldEvaluator(A_key_)->IsDependency(S_next_.ptr(), y2_key_)) {
    dA_dy2_->global_operator()->Init();
    S_next_->GetFieldEvaluator(A_key_)
        ->HasFieldDerivativeChanged(S_next_.ptr(), name_, y2_key_);
    Teuchos::RCP<const CompositeVector> dA_dy2_v = S_next_->GetFieldData(dA_dy2_key_);
    db_->WriteVector("  dwc_dT", dA_dy2_v.ptr());
    
    // -- update the cell-cell block
    dA_dy2_->AddAccumulationTerm(*dA_dy2_v, h, "cell", false);
  }

  if (dB_dy1_ != Teuchos::null &&
      S_next_->GetFieldEvaluator(B_key_)->IsDependency(S_next_.ptr(), y1_key_)) {
    dB_dy1_->global_operator()->Init();
    S_next_->GetFieldEvaluator(B_key_)
        ->HasFieldDerivativeChanged(S_next_.ptr(), name_, y1_key_);
    Teuchos::RCP<const CompositeVector> dB_dy1_v = S_next_->GetFieldData(dB_dy1_key_);
    db_->WriteVector("  dE_dp", dB_dy1_v.ptr());
    
    // -- update the cell-cell block
    dB_dy1_->AddAccumulationTerm(*dB_dy1_v, h, "cell", false);
  }

  if (precon_used_) {
    preconditioner_->AssembleMatrix();
    preconditioner_->UpdatePreconditioner();
  }
}


// applies preconditioner to u and returns the result in Pu
int MPCCoupledCells::ApplyPreconditioner(Teuchos::RCP<const TreeVector> u, Teuchos::RCP<TreeVector> Pu) {
  // write residuals
  if (vo_->os_OK(Teuchos::VERB_HIGH)) {
    *vo_->os() << "Residuals:" << std::endl;
    std::vector<std::string> vnames;
    vnames.push_back("  r_p"); vnames.push_back("  r_T"); 
    std::vector< Teuchos::Ptr<const CompositeVector> > vecs;
    vecs.push_back(u->SubVector(0)->Data().ptr()); 
    vecs.push_back(u->SubVector(1)->Data().ptr()); 
    db_->WriteVectors(vnames, vecs, true);
  }
  
  int ierr = linsolve_preconditioner_->ApplyInverse(*u, *Pu);

  if (vo_->os_OK(Teuchos::VERB_HIGH)) {
    *vo_->os() << "PC * residuals:" << std::endl;
    std::vector<std::string> vnames;
    vnames.push_back("  PC*r_p"); vnames.push_back("  PC*r_T"); 
    std::vector< Teuchos::Ptr<const CompositeVector> > vecs;
    vecs.push_back(Pu->SubVector(0)->Data().ptr()); 
    vecs.push_back(Pu->SubVector(1)->Data().ptr()); 
    db_->WriteVectors(vnames, vecs, true);
  }

  return (ierr > 0) ? 0 : 1;
}


} //  namespace
<commit_msg>pushes verbose object list to linear solver in coupled cells MPC for column runs<commit_after>/* -*-  mode: c++; indent-tabs-mode: nil -*- */
/* -------------------------------------------------------------------------
   ATS

   License: see $ATS_DIR/COPYRIGHT
   Author: <NAME>

   Interface for a StrongMPC which uses a preconditioner in which the
   block-diagonal cell-local matrix is dense.  If the system looks something
   like:

   A( y1, y2, x, t ) = 0
   B( y1, y2, x, t ) = 0

   where y1,y2 are spatially varying unknowns that are discretized using the MFD
   method (and therefore have both cell and face unknowns), an approximation to
   the Jacobian is written as

   [  dA_c/dy1_c  dA_c/dy1_f   dA_c/dy2_c       0      ]
   [  dA_f/dy1_c  dA_f/dy1_f      0              0      ]
   [  dB_c/dy1_c     0          dB_c/dy2_c  dB_c/dy2_f ]
   [      0           0          dB_f/dy2_c  dB_f/dy2_f ]


   Note that the upper left block is the standard preconditioner for the A
   system, and the lower right block is the standard precon for the B system,
   and we have simply added cell-based couplings, dA_c/dy2_c and dB_c/dy1_c.

   In the temperature/pressure system, these correspond to d_water_content /
   d_temperature and d_energy / d_pressure.

   ------------------------------------------------------------------------- */

#include <fstream>
#include "EpetraExt_RowMatrixOut.h"

#include "LinearOperatorFactory.hh"
#include "FieldEvaluator.hh"
#include "Operator.hh"
#include "TreeOperator.hh"
#include "PDE_Accumulation.hh"

#include "mpc_coupled_cells.hh"

namespace Amanzi {

void MPCCoupledCells::Setup(const Teuchos::Ptr<State>& S) {
  StrongMPC<PK_PhysicalBDF_Default>::Setup(S);

  A_key_ = plist_->get<std::string>("conserved quantity A");
  B_key_ = plist_->get<std::string>("conserved quantity B");
  y1_key_ = plist_->get<std::string>("primary variable A");
  y2_key_ = plist_->get<std::string>("primary variable B");
  dA_dy2_key_ = std::string("d")+A_key_+std::string("_d")+y2_key_;
  dB_dy1_key_ = std::string("d")+B_key_+std::string("_d")+y1_key_;

  Key mesh_key = plist_->get<std::string>("domain name");
  mesh_ = S->GetMesh(mesh_key);

  // set up debugger
  db_ = sub_pks_[0]->debugger();

  // Get the sub-blocks from the sub-PK's preconditioners.
  Teuchos::RCP<Operators::Operator> pcA = sub_pks_[0]->preconditioner();
  Teuchos::RCP<Operators::Operator> pcB = sub_pks_[1]->preconditioner();

  // Create the combined operator
  Teuchos::RCP<TreeVectorSpace> tvs = Teuchos::rcp(new TreeVectorSpace());
  tvs->PushBack(Teuchos::rcp(new TreeVectorSpace(Teuchos::rcpFromRef(pcA->DomainMap()))));
  tvs->PushBack(Teuchos::rcp(new TreeVectorSpace(Teuchos::rcpFromRef(pcB->DomainMap()))));

  preconditioner_ = Teuchos::rcp(new Operators::TreeOperator(tvs));
  preconditioner_->SetOperatorBlock(0, 0, pcA);
  preconditioner_->SetOperatorBlock(1, 1, pcB);
  
  // create coupling blocks and push them into the preconditioner...
  S->RequireField(A_key_)->SetMesh(mesh_)->AddComponent("cell", AmanziMesh::CELL, 1);
  S->RequireField(y2_key_)->SetMesh(mesh_)->AddComponent("cell", AmanziMesh::CELL, 1);
  S->RequireFieldEvaluator(A_key_);
  S->RequireFieldEvaluator(y2_key_);

  if (!plist_->get<bool>("no dA/dy2 block", false)) {
    Teuchos::ParameterList& acc_pc_plist = plist_->sublist("dA_dy2 accumulation preconditioner");
    acc_pc_plist.set("entity kind", "cell");
    dA_dy2_ = Teuchos::rcp(new Operators::PDE_Accumulation(acc_pc_plist, mesh_));
    preconditioner_->SetOperatorBlock(0, 1, dA_dy2_->global_operator());
  }

  S->RequireField(B_key_)->SetMesh(mesh_)->AddComponent("cell", AmanziMesh::CELL, 1);
  S->RequireField(y1_key_)->SetMesh(mesh_)->AddComponent("cell", AmanziMesh::CELL, 1);
  S->RequireFieldEvaluator(B_key_);
  S->RequireFieldEvaluator(y1_key_);

  if (!plist_->get<bool>("no dB/dy1 block", false)) {
    Teuchos::ParameterList& acc_pc_plist = plist_->sublist("dB_dy1 accumulation preconditioner");
    acc_pc_plist.set("entity kind", "cell");
    dB_dy1_ = Teuchos::rcp(new Operators::PDE_Accumulation(acc_pc_plist, mesh_));
    preconditioner_->SetOperatorBlock(1, 0, dB_dy1_->global_operator());
  }

  // setup and initialize the preconditioner
  precon_used_ = plist_->isSublist("preconditioner");
  if (precon_used_) {
    preconditioner_->SymbolicAssembleMatrix();
    Teuchos::ParameterList& pc_sublist = plist_->sublist("preconditioner");
    preconditioner_->InitializePreconditioner(pc_sublist);
  }

  // setup and initialize the linear solver for the preconditioner
  if (plist_->isSublist("linear solver")) {
    Teuchos::ParameterList linsolve_sublist = plist_->sublist("linear solver");
    if (!linsolve_sublist.isSublist("verbose object"))
      linsolve_sublist.set("verbose object", plist_->sublist("verbose object"));
    AmanziSolvers::LinearOperatorFactory<Operators::TreeOperator,TreeVector,TreeVectorSpace> fac;
    linsolve_preconditioner_ = fac.Create(linsolve_sublist, preconditioner_);
  } else {
    linsolve_preconditioner_ = preconditioner_;
  }
}


// updates the preconditioner
void MPCCoupledCells::UpdatePreconditioner(double t, Teuchos::RCP<const TreeVector> up,
        double h) {
  StrongMPC<PK_PhysicalBDF_Default>::UpdatePreconditioner(t,up,h);

  if (dA_dy2_ != Teuchos::null &&
      S_next_->GetFieldEvaluator(A_key_)->IsDependency(S_next_.ptr(), y2_key_)) {
    dA_dy2_->global_operator()->Init();
    S_next_->GetFieldEvaluator(A_key_)
        ->HasFieldDerivativeChanged(S_next_.ptr(), name_, y2_key_);
    Teuchos::RCP<const CompositeVector> dA_dy2_v = S_next_->GetFieldData(dA_dy2_key_);
    db_->WriteVector("  dwc_dT", dA_dy2_v.ptr());
    
    // -- update the cell-cell block
    dA_dy2_->AddAccumulationTerm(*dA_dy2_v, h, "cell", false);
  }

  if (dB_dy1_ != Teuchos::null &&
      S_next_->GetFieldEvaluator(B_key_)->IsDependency(S_next_.ptr(), y1_key_)) {
    dB_dy1_->global_operator()->Init();
    S_next_->GetFieldEvaluator(B_key_)
        ->HasFieldDerivativeChanged(S_next_.ptr(), name_, y1_key_);
    Teuchos::RCP<const CompositeVector> dB_dy1_v = S_next_->GetFieldData(dB_dy1_key_);
    db_->WriteVector("  dE_dp", dB_dy1_v.ptr());
    
    // -- update the cell-cell block
    dB_dy1_->AddAccumulationTerm(*dB_dy1_v, h, "cell", false);
  }

  if (precon_used_) {
    preconditioner_->AssembleMatrix();
    preconditioner_->UpdatePreconditioner();
  }
}


// applies preconditioner to u and returns the result in Pu
int MPCCoupledCells::ApplyPreconditioner(Teuchos::RCP<const TreeVector> u, Teuchos::RCP<TreeVector> Pu) {
  // write residuals
  if (vo_->os_OK(Teuchos::VERB_HIGH)) {
    *vo_->os() << "Residuals:" << std::endl;
    std::vector<std::string> vnames;
    vnames.push_back("  r_p"); vnames.push_back("  r_T"); 
    std::vector< Teuchos::Ptr<const CompositeVector> > vecs;
    vecs.push_back(u->SubVector(0)->Data().ptr()); 
    vecs.push_back(u->SubVector(1)->Data().ptr()); 
    db_->WriteVectors(vnames, vecs, true);
  }
  
  int ierr = linsolve_preconditioner_->ApplyInverse(*u, *Pu);

  if (vo_->os_OK(Teuchos::VERB_HIGH)) {
    *vo_->os() << "PC * residuals:" << std::endl;
    std::vector<std::string> vnames;
    vnames.push_back("  PC*r_p"); vnames.push_back("  PC*r_T"); 
    std::vector< Teuchos::Ptr<const CompositeVector> > vecs;
    vecs.push_back(Pu->SubVector(0)->Data().ptr()); 
    vecs.push_back(Pu->SubVector(1)->Data().ptr()); 
    db_->WriteVectors(vnames, vecs, true);
  }

  return (ierr > 0) ? 0 : 1;
}


} //  namespace

{'fraction_non_alphanumeric': 0.16480692500169067, 'fraction_numerical': 0.011361330898762427, 'mean_word_length': 4.141863699582754, 'pattern_counts': {'":': 0, '<': 73, '<?xml version=': 0, '>': 220, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2342024', 'n_tokens_mistral': 5491, 'n_tokens_neox': 5324, 'n_words': 950}