<commit_before>/* -*-  mode: c++; c-default-style: "google"; indent-tabs-mode: nil -*- */
/* -------------------------------------------------------------------------
ATS

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

Derived MPC for flow and energy.  This couples using a block-diagonal coupler.
------------------------------------------------------------------------- */

#include "wrm_richards_evaluator.hh"
#include "mpc_frozen_prec_coupled_flow_energy.hh"

namespace Amanzi {

#define DEBUG_FLAG 1

RegisteredPKFactory<MPCFrozenCoupledFlowEnergy> MPCFrozenCoupledFlowEnergy::reg_("frozen energy-flow preconditioner coupled");

bool MPCFrozenCoupledFlowEnergy::modify_predictor(double h, Teuchos::RCP<TreeVector> u) {
  Teuchos::RCP<CompositeVector> temp_guess = u->SubVector("energy")->data();
  Teuchos::RCP<const CompositeVector> temp = S_next_->GetFieldData("temperature");

  bool update_faces(false);

  int ncells = temp->size("cell",false);
  std::vector<bool> changed(ncells, false);
  for (int c=0; c!=ncells; ++c) {
    if ((*temp)("cell",c) >= 273.15 && (*temp_guess)("cell",c) < 273.15) {
      // freezing
      (*temp_guess)("cell",c) = 273.15 - 1.e-3;
      changed[c] = true;
      update_faces = true;
    } else if ((*temp)("cell",c) <= 273.15 && (*temp_guess)("cell",c) > 273.15) {
      // thawing
      (*temp_guess)("cell",c) = 273.15 - 1.e-3;
      changed[c] = true;
      update_faces = true;
    } else if (273.15 > (*temp)("cell",c) &&
               (*temp)("cell",c) >= 273.1 &&
               ((*temp)("cell",c) - (*temp_guess)("cell",c)) > 1.e-2) {
      // catch the 2nd step in freezing -- after the 2nd step the
      // extrapolation should be ok?
      (*temp_guess)("cell",c) = (*temp)("cell",c);
      changed[c] = true;
      update_faces = true;
    }
  }

  if (update_faces) {
    AmanziMesh::Entity_ID_List cells;

    int f_owned = temp_guess->size("face");
    for (int f=0; f!=f_owned; ++f) {
      cells.clear();
      temp_guess->mesh()->face_get_cells(f, AmanziMesh::USED, &cells);
      if (cells.size() == 2) {
        if (changed[cells[0]] || changed[cells[1]]) {
          (*temp_guess)("face",f) = ((*temp_guess)("cell",cells[0])
                  + (*temp_guess)("cell",cells[1])) / 2.0;
        }
      }
    }
    return true;
  }
  return false;
}

bool MPCFrozenCoupledFlowEnergy::modify_predictor_temp(double h, Teuchos::RCP<TreeVector> u) {
  // modification of the initial guess occurs by updating T using this guess,
  // then calculating the p that would be required to keep the same water
  // mass.

  // Stuff temperature into state
  Teuchos::RCP<TreeVector> uT = u->SubVector("energy");
  sub_pks_[1]->changed_solution();

#if DEBUG_FLAG
  std::cout << std::endl;
  std::cout << "Modifying guess:" << std::endl;
  std::cout << "  T0: " << (*uT->data())("cell",0) << std::endl;
  std::cout << "  T1: " << (*uT->data())("cell",99) << std::endl;
#endif

  // update water content, which will get all the needed vals updated at the new temp.
  S_next_->GetFieldEvaluator("water_content")
      ->HasFieldChanged(S_next_.ptr(), "richards_pk");<commit_msg>updating the modification process to modify faces correctly as well<commit_after>/* -*-  mode: c++; c-default-style: "google"; indent-tabs-mode: nil -*- */
/* -------------------------------------------------------------------------
ATS

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

Derived MPC for flow and energy.  This couples using a block-diagonal coupler.
------------------------------------------------------------------------- */

#include "wrm_richards_evaluator.hh"
#include "mpc_frozen_prec_coupled_flow_energy.hh"

namespace Amanzi {

#define DEBUG_FLAG 0

RegisteredPKFactory<MPCFrozenCoupledFlowEnergy> MPCFrozenCoupledFlowEnergy::reg_("frozen energy-flow preconditioner coupled");

bool MPCFrozenCoupledFlowEnergy::modify_predictor(double h, Teuchos::RCP<TreeVector> u) {
  Teuchos::RCP<CompositeVector> temp_guess = u->SubVector("energy")->data();
  Epetra_MultiVector& guess_cells = *temp_guess->ViewComponent("cell",false);
  Epetra_MultiVector& guess_faces = *temp_guess->ViewComponent("face",false);

  Teuchos::RCP<const CompositeVector> temp = S_next_->GetFieldData("temperature");
  const Epetra_MultiVector& temp_cells = *temp->ViewComponent("cell",false);

  bool update_faces(false);

#if DEBUG_FLAG
  std::cout << "--- Modifying Guess: ---" << std::cout;
#endif

  int ncells = temp->size("cell",false);
  std::vector<bool> changed(ncells, false);
  for (int c=0; c!=ncells; ++c) {
    if (temp_cells[0][c] >= 273.15 && guess_cells[0][c] < 273.15) {
      // freezing
      guess_cells[0][c] = 273.15 - 1.e-4;
      changed[c] = true;
      update_faces = true;
    } else if (temp_cells[0][c] <= 273.15 && guess_cells[0][c] > 273.15) {
      // thawing
      guess_cells[0][c] = 273.15 - 1.e-3;
      changed[c] = true;
      update_faces = true;
    } else if (273.15 > temp_cells[0][c] &&
               temp_cells[0][c] >= 273.1 &&
               (temp_cells[0][c] - guess_cells[0][c]) > 1.e-2) {
      // catch the 2nd step in freezing -- after the 2nd step the
      // extrapolation should be ok?
      guess_cells[0][c] = temp_cells[0][c];
      changed[c] = true;
      update_faces = true;
    }
  }

  if (update_faces) {
    temp_guess->ScatterMasterToGhosted("cell");
    AmanziMesh::Entity_ID_List cells;

    int f_owned = temp_guess->size("face");
    for (int f=0; f!=f_owned; ++f) {
      cells.clear();
      temp_guess->mesh()->face_get_cells(f, AmanziMesh::USED, &cells);
      if (cells.size() == 2) {
        if (changed[cells[0]] || changed[cells[1]]) {
          guess_faces[0][f] = (guess_cells[0][cells[0]]
                  + guess_cells[0][cells[1]]) / 2.0;
        }
      }
    }
    return true;
  }
  return false;
}

bool MPCFrozenCoupledFlowEnergy::modify_predictor_temp(double h, Teuchos::RCP<TreeVector> u) {
  // modification of the initial guess occurs by updating T using this guess,
  // then calculating the p that would be required to keep the same water
  // mass.

  // Stuff temperature into state
  Teuchos::RCP<TreeVector> uT = u->SubVector("energy");
  sub_pks_[1]->changed_solution();

#if DEBUG_FLAG
  std::cout << std::endl;
  std::cout << "Modifying guess:" << std::endl;
  std::cout << "  T0: " << (*uT->data())("cell",0) << std::endl;
  std::cout << "  T1: " << (*uT->data())("cell",99) << std::endl;
#endif

  // update water content, which will get all the needed vals updated at the new temp.
  S_next_->GetFieldEvaluator("water_content")
      ->HasFieldChanged(S_next_.ptr(), "richards_pk");
{'fraction_non_alphanumeric': 0.20941140868512634, 'fraction_numerical': 0.021334543803903767, 'mean_word_length': 2.9368671828469326, 'pattern_counts': {'":': 0, '<': 63, '<?xml version=': 0, '>': 59, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4817708', 'n_tokens_mistral': 2429, 'n_tokens_neox': 2272, 'n_words': 476}