text:
/* -*-  mode: c++; indent-tabs-mode: nil -*- */
/*
  ATS is released under the three-clause BSD License.
  The terms of use and "as is" disclaimer for this license are
  provided in the top-level COPYRIGHT file.

  Authors: Ethan Coon (coonet@ornl.gov)
*/

#include "Epetra_Import.h"

#include "flow_bc_factory.hh"

#include "Point.hh"

#include "upwind_cell_centered.hh"
#include "upwind_arithmetic_mean.hh"
#include "upwind_total_flux.hh"
#include "upwind_gravity_flux.hh"

#include "CompositeVectorFunction.hh"
#include "CompositeVectorFunctionFactory.hh"

#include "predictor_delegate_bc_flux.hh"
#include "wrm_evaluator.hh"
#include "rel_perm_evaluator.hh"
#include "richards_water_content_evaluator.hh"
#include "OperatorDefs.hh"
#include "BoundaryFlux.hh"
#include "pk_helpers.hh"

#include "richards.hh"

namespace Amanzi {
namespace Flow {

// -------------------------------------------------------------
// Constructor
// -------------------------------------------------------------

Richards::Richards(Teuchos::ParameterList& pk_tree,
                   const Teuchos::RCP<Teuchos::ParameterList>& glist,
                   const Teuchos::RCP<State>& S,
                   const Teuchos::RCP<TreeVector>& solution) :
    PK(pk_tree, glist,  S, solution),
    PK_PhysicalBDF_Default(pk_tree, glist,  S, solution),
    coupled_to_surface_via_head_(false),
    coupled_to_surface_via_flux_(false),
    infiltrate_only_if_unfrozen_(false),
    modify_predictor_with_consistent_faces_(false),
    modify_predictor_wc_(false),
    modify_predictor_bc_flux_(false),
    modify_predictor_first_bc_flux_(false),
    upwind_from_prev_flux_(false),
    dynamic_mesh_(false),
    clobber_boundary_flux_dir_(false),
    vapor_diffusion_(false),
    perm_scale_(1.),
    jacobian_(false),
    jacobian_lag_(0),
    iter_(0),
    iter_counter_time_(0.)
{
  // set a default absolute tolerance
  if (!plist_->isParameter("absolute error tolerance"))
    plist_->set("absolute error tolerance", .5 * .1 * 55000.); // phi * s * nl

  // get field names
  conserved_key_ = Keys::readKey(*plist_, domain_, "conserved", "water_content");
  mass_dens_key_ = Keys::readKey(*plist_, domain_, "mass density", "mass_density_liquid");
  molar_dens_key_ = Keys::readKey(*plist_, domain_, "molar density", "molar_density_liquid");
  perm_key_ = Keys::readKey(*plist_, domain_, "permeability", "permeability");
  coef_key_ = Keys::readKey(*plist_, domain_, "conductivity", "relative_permeability");
  uw_coef_key_ = Keys::readKey(*plist_, domain_, "upwinded conductivity", "upwind_relative_permeability");
  flux_key_ = Keys::readKey(*plist_, domain_, "darcy flux", "water_flux");
  flux_dir_key_ = Keys::readKey(*plist_, domain_, "darcy flux direction", "water_flux_direction");
  velocity_key_ = Keys::readKey(*plist_, domain_, "darcy velocity", "darcy_velocity");
  sat_key_ = Keys::readKey(*plist_, domain_, "saturation", "saturation_liquid");
  sat_gas_key_ = Keys::readKey(*plist_, domain_, "saturation gas", "saturation_gas");
  sat_ice_key_ = Keys::readKey(*plist_, domain_, "saturation ice", "saturation_ice");
  capillary_pressure_gas_liq_key_ = Keys::readKey(*plist_, domain_,
      "capillary_pressure_gas_liq", "capillary_pressure_gas_liq");
  capillary_pressure_liq_ice_key_ = Keys::readKey(*plist_, domain_,
      "capillary_pressure_liq_ice", "capillary_pressure_liq_ice");
}

// -------------------------------------------------------------
// Setup data
// -------------------------------------------------------------
void Richards::Setup()
{
  PK_PhysicalBDF_Default::Setup();
  SetupRichardsFlow_();
  SetupPhysicalEvaluators_();
};

// -------------------------------------------------------------
// Pieces of the construction process that are common to all
// Richards-like PKs.
// -------------------------------------------------------------
void Richards::SetupRichardsFlow_()
{
  // Get data for special-case entities.
  S_->Require<CompositeVector,CompositeVectorSpace>(cell_vol_key_, tag_next_)
    .SetMesh(mesh_)
    ->AddComponent("cell", AmanziMesh::CELL, 1);
  S_->RequireEvaluator(cell_vol_key_, tag_next_);

  // Set up Operators
  // -- boundary conditions
  Teuchos::ParameterList bc_plist = plist_->sublist("boundary conditions", true);
  FlowBCFactory bc_factory(mesh_, bc_plist);
  bc_pressure_ = bc_factory.CreatePressure();
  bc_head_ = bc_factory.CreateHead();
  bc_level_ = bc_factory.CreateFixedLevel();
  bc_flux_ = bc_factory.CreateMassFlux();
  bc_seepage_ = bc_factory.CreateSeepageFacePressure();
  bc_seepage_->Compute(0.); // compute at t=0 to set up
  std::tie(bc_seepage_infilt_explicit_, bc_seepage_infilt_) =
    bc_factory.CreateSeepageFacePressureWithInfiltration();
  bc_seepage_infilt_->Compute(0.); // compute at t=0 to set up
  bc_rho_water_ = bc_plist.get<double>("hydrostatic water density [kg m^-3]",1000.);

  // -- linear tensor coefficients
  // scaling for permeability
  perm_scale_ = plist_->get<double>("permeability rescaling", 1.e7);

  // permeability type - scalar or tensor?
  Teuchos::ParameterList& perm_list = S_->GetEvaluatorList(perm_key_);
  std::string perm_type = perm_list.get<std::string>("permeability type", "scalar");
  if (perm_type == "scalar") {
    perm_tensor_rank_ = 1;
    num_perm_vals_ = 1;
  } else if (perm_type == "horizontal and vertical") {
    perm_tensor_rank_ = 2;
    num_perm_vals_ = 2;
  } else if (perm_type == "diagonal tensor") {
    perm_tensor_rank_ = 2;
    num_perm_vals_ = mesh_->space_dimension();
  } else if (perm_type == "full tensor") {
    perm_tensor_rank_ = 2;
    num_perm_vals_ = (mesh_->space_dimension() == 3) ? 6 : 3;
  } else {
    Errors::Message message("`permeability type` must be one of the following: \"scalar\", \"diagonal tensor\", \"full tensor\", or \"horizontal and vertical\".");
    Exceptions::amanzi_throw(message);
  }

  // data allocation -- move to State!
  unsigned int c_owned = mesh_->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED);
  K_ = Teuchos::rcp(new std::vector<WhetStone::Tensor>(c_owned));
  for (unsigned int c=0; c!=c_owned; ++c) {
    (*K_)[c].Init(mesh_->space_dimension(),perm_tensor_rank_);
  }

  // -- nonlinear coefficients/upwinding
  // if coupled to the surface, how do we deal with the surface face
  if (plist_->isParameter("surface rel perm strategy")) {
    clobber_policy_ = plist_->get<std::string>("surface rel perm strategy");
  } else if (plist_->get<bool>("clobber surface rel perm", false)) {
    clobber_policy_ = "clobber";
  } else if (plist_->get<bool>("max surface rel perm", false)) {
    clobber_policy_ = "max";
  } else if (plist_->get<bool>("unsaturated clobber surface rel perm", false)) {
    clobber_policy_ = "unsaturated";
  } else {
    clobber_policy_ = "none";
  }
  clobber_boundary_flux_dir_ = plist_->get<bool>("clobber boundary flux direction for upwinding", false);

  // what upwinding method to use
  std::string method_name = plist_->get<std::string>("relative permeability method", "upwind with Darcy flux");
  if (method_name == "upwind with gravity") {
    upwinding_ = Teuchos::rcp(new Operators::UpwindGravityFlux(name_,
            tag_next_, K_));
    Krel_method_ = Operators::UPWIND_METHOD_GRAVITY;
  } else if (method_name == "cell centered") {
    upwinding_ = Teuchos::rcp(new Operators::UpwindCellCentered(name_, tag_next_));
    Krel_method_ = Operators::UPWIND_METHOD_CENTERED;
  } else if (method_name == "upwind with Darcy flux") {
    upwinding_ = Teuchos::rcp(new Operators::UpwindTotalFlux(name_,
            tag_next_, flux_dir_key_, 1.e-5));
    Krel_method_ = Operators::UPWIND_METHOD_TOTAL_FLUX;
  } else if (method_name == "arithmetic mean") {
    upwinding_ = Teuchos::rcp(new Operators::UpwindArithmeticMean(name_,
            tag_next_));
    Krel_method_ = Operators::UPWIND_METHOD_ARITHMETIC_MEAN;
  } else {
    std::stringstream messagestream;
    messagestream << "Richards Flow PK has no upwinding method named: " << method_name;
    Errors::Message message(messagestream.str());
    Exceptions::amanzi_throw(message);
  }

  // require the data on appropriate locations
  std::string coef_location = upwinding_->CoefficientLocation();
  if (coef_location == "upwind: face") {
    S_->Require<CompositeVector,CompositeVectorSpace>(uw_coef_key_, tag_next_,  name_).SetMesh(mesh_)
        ->SetGhosted()->SetComponent("face", AmanziMesh::FACE, 1);
  } else if (coef_location == "standard: cell") {
    S_->Require<CompositeVector,CompositeVectorSpace>(uw_coef_key_, tag_next_,  name_).SetMesh(mesh_)
        ->SetGhosted()->SetComponent("cell", AmanziMesh::CELL, 1);
  } else {
    Errors::Message message("Unknown upwind coefficient location in Richards flow.");
    Exceptions::amanzi_throw(message);
  }
  S_->GetRecordW(uw_coef_key_, tag_next_, name_).set_io_vis(false);

  // -- create the forward operator for the diffusion term
  Teuchos::ParameterList& mfd_plist = plist_->sublist("diffusion");
  mfd_plist.set("nonlinear coefficient", coef_location);
  mfd_plist.set("gravity", true);

  Operators::PDE_DiffusionFactory opfactory;
  matrix_diff_ = opfactory.CreateWithGravity(mfd_plist, mesh_, bc_);
  matrix_ = matrix_diff_->global_operator();

  // -- create the operator, data for flux directions
  Teuchos::ParameterList face_diff_list(mfd_plist);
  face_diff_list.set("nonlinear coefficient", "none");
  face_matrix_diff_ = opfactory.CreateWithGravity(face_diff_list, mesh_, bc_);

  S_->Require<CompositeVector,CompositeVectorSpace>(flux_dir_key_, tag_next_,  name_)
    .SetMesh(mesh_)->SetGhosted()
    ->SetComponent("face", AmanziMesh::FACE, 1);

  // -- create the operators for the preconditioner
  //    diffusion
  // NOTE: Can this be a clone of the primary operator? --etc
  //      get the discretiation type
  Teuchos::ParameterList& mfd_pc_plist = plist_->sublist("diffusion preconditioner");
  mfd_pc_plist.set("nonlinear coefficient", coef_location);
  mfd_pc_plist.set("gravity", true);
  if (!mfd_pc_plist.isParameter("discretization primary"))
    mfd_pc_plist.set("discretization primary", mfd_plist.get<std::string>("discretization primary"));
  if (!mfd_pc_plist.isParameter("discretization secondary") && mfd_plist.isParameter("discretization secondary"))
    mfd_pc_plist.set("discretization secondary", mfd_plist.get<std::string>("discretization secondary"));
  if (!mfd_pc_plist.isParameter("schema") && mfd_plist.isParameter("schema"))
    mfd_pc_plist.set("schema", mfd_plist.get<Teuchos::Array<std::string> >("schema"));
  if (mfd_pc_plist.get<bool>("include Newton correction", false)) {
    if (mfd_pc_plist.get<std::string>("discretization primary") == "fv: default") {
      mfd_pc_plist.set("Newton correction", "true Jacobian");
    } else {
      mfd_pc_plist.set("Newton correction", "approximate Jacobian");
    }
  }

  //    get the inverse method
  precon_used_ = plist_->isSublist("preconditioner") ||
    plist_->isSublist("inverse") ||
    plist_->isSublist("linear solver");
  if (precon_used_) {
    mfd_pc_plist.set("inverse", plist_->sublist("inverse"));
    // old style... deprecate me!
    mfd_pc_plist.sublist("inverse").setParameters(plist_->sublist("preconditioner"));
    mfd_pc_plist.sublist("inverse").setParameters(plist_->sublist("linear solver"));
  }

  //    create the operator
  preconditioner_diff_ = opfactory.CreateWithGravity(mfd_pc_plist, mesh_, bc_);
  preconditioner_ = preconditioner_diff_->global_operator();

  //    If using approximate Jacobian for the preconditioner, we also need
  //    derivative information.  For now this means upwinding the derivative.
  jacobian_ = mfd_pc_plist.get<std::string>("Newton correction", "none") != "none";
  if (jacobian_) {
    jacobian_lag_ = mfd_pc_plist.get<int>("Newton correction lag", 0);

    // require the derivative drel_perm/dp
    S_->RequireDerivative<CompositeVector,CompositeVectorSpace>(coef_key_,
            tag_next_, key_, tag_next_);
    if (mfd_pc_plist.get<std::string>("discretization primary") != "fv: default"){
      // MFD -- upwind required, require data
      duw_coef_key_ = Keys::getDerivKey(uw_coef_key_, key_);
      S_->Require<CompositeVector,CompositeVectorSpace>(duw_coef_key_, tag_next_,  name_)
        .SetMesh(mesh_)->SetGhosted()
        ->SetComponent("face", AmanziMesh::FACE, 1);
      upwinding_deriv_ = Teuchos::rcp(new Operators::UpwindTotalFlux(name_,
              tag_next_, flux_dir_key_, 1.e-8));
    } else {
      // FV -- no upwinding of derivative
      duw_coef_key_ = std::string();
    }
  }

  // -- accumulation terms
  Teuchos::ParameterList& acc_pc_plist = plist_->sublist("accumulation preconditioner");
  acc_pc_plist.set<std::string>("entity kind", "cell");
  preconditioner_acc_ = Teuchos::rcp(new Operators::PDE_Accumulation(acc_pc_plist, preconditioner_));

  // // -- vapor diffusion terms
  // vapor_diffusion_ = plist_->get<bool>("include vapor diffusion", false);
  // if (vapor_diffusion_){
  //   AMANZI_ASSERT(0); // untested!

  //   // Create the vapor diffusion vectors
  //   S_->Require<CompositeVector,CompositeVectorSpace>("vapor_diffusion_pressure", tag_next_,  name_).SetMesh(mesh_)->SetGhosted()
  //       ->SetComponent("cell", AmanziMesh::CELL, 1);
  //   S_->GetRecordW("vapor_diffusion_pressure",name_)->set_io_vis(true);

  //   S_->Require<CompositeVector,CompositeVectorSpace>("vapor_diffusion_temperature", tag_next_,  name_).SetMesh(mesh_)->SetGhosted()
  //     ->SetComponent("cell", AmanziMesh::CELL, 1);
  //   S_->GetRecordW("vapor_diffusion_temperature",name_)->set_io_vis(true);

  //   // operator for the vapor diffusion terms
  //   matrix_vapor_ = Operators::CreateMatrixMFD(mfd_plist, mesh_);
  // }

  // -- source terms
  is_source_term_ = plist_->get<bool>("source term", false);
  if (is_source_term_) {
    if (source_key_.empty())
      source_key_ = Keys::readKey(*plist_, domain_, "source", "water_source");
    source_term_is_differentiable_ =
        plist_->get<bool>("source term is differentiable", true);
    explicit_source_ = plist_->get<bool>("explicit source term", false);

    requireAtNext(source_key_, tag_next_, *S_)
      .SetMesh(mesh_)
      ->AddComponent("cell", AmanziMesh::CELL, 1);
    if (source_term_is_differentiable_) {
      // require derivative of source
      S_->RequireDerivative<CompositeVector,CompositeVectorSpace>(source_key_,
              tag_next_, key_, tag_next_);
    }
  }

  // coupling to the surface
  // -- coupling done by a Neumann condition
  coupled_to_surface_via_flux_ = plist_->get<bool>("coupled to surface via flux", false);
  if (coupled_to_surface_via_flux_) {
    Key domain_surf = Keys::readDomainHint(*plist_, domain_, "subsurface", "surface");
    ss_flux_key_ = Keys::readKey(*plist_, domain_surf, "surface-subsurface flux", "surface_subsurface_flux");
    S_->Require<CompositeVector,CompositeVectorSpace>(ss_flux_key_, tag_next_)
      .SetMesh(S_->GetMesh(domain_surf))
      ->AddComponent("cell", AmanziMesh::CELL, 1);
  }

  // -- coupling done by a Dirichlet condition
  coupled_to_surface_via_head_ = plist_->get<bool>("coupled to surface via head", false);
  if (coupled_to_surface_via_head_) {
    Key domain_surf = Keys::readDomainHint(*plist_, domain_, "subsurface", "surface");
    ss_primary_key_ = Keys::readKey(*plist_, domain_surf, "pressure", "pressure");
    S_->Require<CompositeVector,CompositeVectorSpace>(ss_primary_key_, tag_next_)
      .SetMesh(S_->GetMesh(domain_surf))
      ->AddComponent("cell", AmanziMesh::Entity_kind::CELL, 1);
  }

  // -- Make sure coupling isn't flagged multiple ways.
  if (coupled_to_surface_via_flux_ && coupled_to_surface_via_head_) {
    Errors::Message message("Richards PK requested both flux and head coupling -- choose one.");
    Exceptions::amanzi_throw(message);
  }

  // Require fields and evaluators
  // -- pressure, the primary variable
  //  NOTE: no need to require evaluator for p here, either at the old or new
  //  times, as this was done in pk_physical.  All we have to do is set the
  //  structure.
  CompositeVectorSpace matrix_cvs = matrix_->RangeMap();
  compute_boundary_values_ = plist_->get<bool>("compute boundary values", false);
  if (compute_boundary_values_)
    matrix_cvs.AddComponent("boundary_face", AmanziMesh::BOUNDARY_FACE, 1);
  S_->Require<CompositeVector,CompositeVectorSpace>(key_, tag_next_, name_)
    .Update(matrix_cvs)->SetGhosted();

  // -- flux is managed here as a primary variable
  requireAtNext(flux_key_, tag_next_, *S_, name_)
    .SetMesh(mesh_)->SetGhosted()
    ->SetComponent("face", AmanziMesh::FACE, 1);

  // -- also need a velocity, but only for vis/diagnostics, so might as well
  // -- only keep at NEXT
  requireAtNext(velocity_key_, Tags::NEXT, *S_, name_)
    .SetMesh(mesh_)->SetGhosted()
    ->SetComponent("cell", AmanziMesh::CELL, 3);

  // Globalization and other timestep control flags
  // -- predictors
  modify_predictor_with_consistent_faces_ =
    plist_->get<bool>("modify predictor with consistent faces", false);
  modify_predictor_bc_flux_ =
    plist_->get<bool>("modify predictor for flux BCs", false);
  modify_predictor_first_bc_flux_ =
    plist_->get<bool>("modify predictor for initial flux BCs", false);
  modify_predictor_wc_ =
    plist_->get<bool>("modify predictor via water content", false);

  // -- correctors
  p_limit_ = plist_->get<double>("limit correction to pressure change [Pa]", -1.);
  patm_limit_ = plist_->get<double>("limit correction to pressure change when crossing atmospheric [Pa]", -1.);

  // -- valid step controls
  sat_change_limit_ = plist_->get<double>("max valid change in saturation in a time step [-]", -1.);
  sat_ice_change_limit_ = plist_->get<double>("max valid change in ice saturation in a time step [-]", -1.);
}

// -------------------------------------------------------------
// Create the physical evaluators for water content, water
// retention, rel perm, etc, that are specific to Richards.
// -------------------------------------------------------------
void Richards::SetupPhysicalEvaluators_()
{
  // -- Absolute permeability.
  //       For now, we assume scalar permeability.  This will change.
  requireAtNext(perm_key_, tag_next_, *S_)
    .SetMesh(mesh_)->SetGhosted()
    ->AddComponent("cell", AmanziMesh::CELL, num_perm_vals_);

  // -- water content, and evaluator, and derivative for PC
  requireAtNext(conserved_key_, tag_next_, *S_)
    .SetMesh(mesh_)->SetGhosted()
    ->AddComponent("cell", AmanziMesh::CELL, 1);
  S_->RequireDerivative<CompositeVector,CompositeVectorSpace>(conserved_key_,
          tag_next_, key_, tag_next_);

  //    and at the current time, where it is a copy evaluator
  requireAtCurrent(conserved_key_, tag_current_, *S_, name_);

  // -- Water retention evaluators
  // This deals with deprecated location for the WRM list (in the PK).  Move it
  // to state.
  if (plist_->isSublist("water retention evaluator")) {
    auto& wrm_plist = S_->GetEvaluatorList(sat_key_);
    wrm_plist.setParameters(plist_->sublist("water retention evaluator"));
    wrm_plist.set("evaluator type", "WRM");
  }
  if (!S_->HasEvaluator(coef_key_, tag_next_) &&
      (S_->GetEvaluatorList(coef_key_).numParams() == 0)) {
    Teuchos::ParameterList& kr_plist = S_->GetEvaluatorList(coef_key_);
    kr_plist.setParameters(S_->GetEvaluatorList(sat_key_));
    kr_plist.set<std::string>("evaluator type", "WRM rel perm");
  }

  // -- saturation
  requireAtNext(sat_key_, tag_next_, *S_)
    .SetMesh(mesh_)->SetGhosted()
    ->AddComponent("cell", AmanziMesh::CELL, 1);
  requireAtNext(sat_gas_key_, tag_next_, *S_)
    .SetMesh(mesh_)->SetGhosted()
    ->AddComponent("cell", AmanziMesh::CELL, 1);
  auto& wrm = S_->RequireEvaluator(sat_key_, tag_next_);

  //    and at the current time, where it is a copy evaluator
  requireAtCurrent(sat_key_, tag_current_, *S_, name_);

  // -- rel perm
  S_->GetEvaluatorList(coef_key_).set<double>("permeability rescaling", perm_scale_);
  requireAtNext(coef_key_, tag_next_, *S_)
    .SetMesh(mesh_)->SetGhosted()
    ->AddComponent("cell", AmanziMesh::CELL, 1)
    ->AddComponent("boundary_face", AmanziMesh::BOUNDARY_FACE, 1);

  // -- get the WRM models
  auto wrm_eval = dynamic_cast<Flow::WRMEvaluator*>(&wrm);
  AMANZI_ASSERT(wrm_eval != nullptr);
  wrms_ = wrm_eval->get_WRMs();

  // -- molar density used to infer liquid Darcy velocity from flux
  requireAtNext(molar_dens_key_, tag_next_, *S_)
    .SetMesh(mesh_)->SetGhosted()
    ->AddComponent("cell", AmanziMesh::CELL, 1);

  // -- liquid mass density for the gravity fluxes
  requireAtNext(mass_dens_key_, tag_next_, *S_)
    .SetMesh(mesh_)->SetGhosted()
    ->AddComponent("cell", AmanziMesh::CELL, 1);
}

// -------------------------------------------------------------
// Initialize PK
// -------------------------------------------------------------
void Richards::Initialize()
{
  // Initialize via hydrostatic balance
  if (!S_->GetRecordW(key_, tag_next_, name_).initialized()) InitializeHydrostatic_(tag_next_);

  // Initialize in the standard ways
  PK_PhysicalBDF_Default::Initialize();

  // check whether this is a dynamic mesh problem
  if (S_->IsDeformableMesh(domain_)) dynamic_mesh_ = true;

  // Set extra fields as initialized -- these don't currently have evaluators,
  // and will be initialized in the call to commit_state()
  S_->GetW<CompositeVector>(uw_coef_key_, tag_next_, name_).PutScalar(1.0);
  S_->GetRecordW(uw_coef_key_, tag_next_, name_).set_initialized();

  if (!duw_coef_key_.empty()) {
    S_->GetW<CompositeVector>(duw_coef_key_, tag_next_, name_).PutScalar(1.0);
    S_->GetRecordW(duw_coef_key_, tag_next_, name_).set_initialized();
  }

  S_->GetW<CompositeVector>(flux_key_, tag_next_, name()).PutScalar(0.0);
  S_->GetRecordW(flux_key_, tag_next_, name()).set_initialized();
  changedEvaluatorPrimary(flux_key_, tag_next_, *S_);

  S_->GetW<CompositeVector>(flux_dir_key_, tag_next_, name()).PutScalar(0.0);
  S_->GetRecordW(flux_dir_key_, tag_next_, name()).set_initialized();
  S_->GetW<CompositeVector>(velocity_key_, Tags::NEXT, name()).PutScalar(0.0);
  S_->GetRecordW(velocity_key_, Tags::NEXT, name()).set_initialized();

  // absolute perm
  SetAbsolutePermeabilityTensor_(tag_next_);

  // operators
  const AmanziGeometry::Point& g = S_->Get<AmanziGeometry::Point>("gravity", Tags::DEFAULT);
  matrix_diff_->SetGravity(g);
  matrix_diff_->SetBCs(bc_, bc_);
  matrix_diff_->SetTensorCoefficient(K_);

  preconditioner_diff_->SetGravity(g);
  preconditioner_diff_->SetBCs(bc_, bc_);
  preconditioner_diff_->SetTensorCoefficient(K_);

  face_matrix_diff_->SetGravity(g);
  face_matrix_diff_->SetBCs(bc_, bc_);
  face_matrix_diff_->SetTensorCoefficient(K_);
  face_matrix_diff_->SetScalarCoefficient(Teuchos::null, Teuchos::null);

  // if (vapor_diffusion_){
  //   //vapor diffusion
  //   matrix_vapor_->CreateMFDmassMatrices(Teuchos::null);
  //   // residual vector for vapor diffusion
  //   res_vapor = Teuchos::rcp(new CompositeVector(*S_->GetPtr<CompositeVector>(key_)));
  // }
};

void Richards::InitializeHydrostatic_(const Tag& tag)
{
  // constant head over the surface
  if (plist_->sublist("initial condition").isParameter("hydrostatic head [m]")) {
    double head_wt = plist_->sublist("initial condition").get<double>("hydrostatic head [m]");
    double rho = plist_->sublist("initial condition").get<double>("hydrostatic water density [kg m^-3]");
    int ncols = mesh_->num_columns(false);

    const auto& gvec = S_->Get<AmanziGeometry::Point>("gravity", Tags::DEFAULT);
    int z_index = mesh_->space_dimension() - 1;
    double g = -gvec[z_index];
    double p_atm = S_->Get<double>("atmospheric_pressure", Tags::DEFAULT);

    // set pressure on the column of faces and cells
    Teuchos::RCP<CompositeVector> pres = S_->GetPtrW<CompositeVector>(key_, tag, name());
    Teuchos::RCP<Epetra_IntVector> flags = Teuchos::null;
    if (pres->HasComponent("face"))
      flags = Teuchos::rcp(new Epetra_IntVector(*pres->Map().Map("face",false)));

    { // context for viewcomponent -- do cells
      Epetra_MultiVector& pres_c = *pres->ViewComponent("cell", false);
      Teuchos::RCP<Epetra_MultiVector> pres_f = Teuchos::null;
      if (pres->HasComponent("face")) {
        pres_f = pres->ViewComponent("face", false);
      }

      int ncells_per = -1;
      for (int col=0; col!=ncols; ++col) {
        auto& col_cells = mesh_->cells_of_column(col);
        if (ncells_per < 0) ncells_per = col_cells.size();
        AMANZI_ASSERT(col_cells.size() == ncells_per);
        auto& col_faces = mesh_->faces_of_column(col);
        AMANZI_ASSERT(col_faces.size() == col_cells.size()+1);
        double z_wt = mesh_->face_centroid(col_faces[0])[z_index] + head_wt;

        if (pres_f.get()) {
          (*pres_f)[0][col_faces[0]] = p_atm + rho * g * head_wt;
          (*flags)[col_faces[0]] = 1;
        }
        for (int lcv_c=0; lcv_c!=col_cells.size(); ++lcv_c) {
          AmanziMesh::Entity_ID c = col_cells[lcv_c];
          AmanziMesh::Entity_ID f = col_faces[lcv_c+1];
          pres_c[0][c] = p_atm + rho * g * (z_wt - mesh_->cell_centroid(c)[z_index]);
          if (pres_f.get()) {
            (*pres_f)[0][f] = p_atm + rho * g * (z_wt - mesh_->face_centroid(f)[z_index]);
            (*flags)[f] = 1;
          }
        }
      }
    }

    if (pres->HasComponent("face")) {
      // communicate, then deal with horizontal-normal faces
      pres->ScatterMasterToGhosted("cell");
      {
        const Epetra_MultiVector& pres_c = *pres->ViewComponent("cell", false);
        Epetra_MultiVector& pres_f = *pres->ViewComponent("face", false);
        for (AmanziMesh::Entity_ID f=0; f!=pres_f.MyLength(); ++f) {
          if (!(*flags)[f]) {
            AmanziMesh::Entity_ID_List f_cells;
            mesh_->face_get_cells(f, AmanziMesh::Parallel_type::ALL, &f_cells);
            if (f_cells.size() == 1) {
              // boundary face, use the cell value as the water table is
              // assumed to parallel the cell structure
              pres_f[0][f] = pres_c[0][f_cells[0]];
              (*flags)[f] = 1;
            } else {
              AMANZI_ASSERT(f_cells.size() == 2);
              // interpolate between cells
              pres_f[0][f] = (pres_c[0][f_cells[0]] + pres_c[0][f_cells[1]]) / 2.;
              (*flags)[f] = 1;
            }
          }
        }
      }
    }
    S_->GetRecordW(key_, tag, name()).set_initialized();
  }

  // constant head datum
  if (plist_->sublist("initial condition").isParameter("hydrostatic water level [m]")) {
    double z_wt = plist_->sublist("initial condition").get<double>("hydrostatic water level [m]");
    double rho = plist_->sublist("initial condition").get<double>("hydrostatic water density [kg m^-3]");

    int z_index = mesh_->space_dimension() - 1;
    const auto& gravity = S_->Get<AmanziGeometry::Point>("gravity", Tags::DEFAULT);
    double g = -gravity[z_index];

    double p_atm = S_->Get<double>("atmospheric_pressure", Tags::DEFAULT);

    Teuchos::RCP<CompositeVector> pres = S_->GetPtrW<CompositeVector>(key_, tag, name());
    Epetra_MultiVector& pres_c = *pres->ViewComponent("cell", false);
    for (int c=0; c!=pres_c.MyLength(); ++c) {
      pres_c[0][c] = p_atm + rho * g * (z_wt - mesh_->cell_centroid(c)[z_index]);
    }

    if (pres->HasComponent("face")) {
      Epetra_MultiVector& pres_f = *pres->ViewComponent("face", false);
      for (int f=0; f!=pres_f.MyLength(); ++f) {
        pres_f[0][f] = p_atm + rho * g * (z_wt - mesh_->face_centroid(f)[z_index]);
      }
    }
    S_->GetRecordW(key_, tag, name()).set_initialized();
  }
}

// -----------------------------------------------------------------------------
// Update any secondary (dependent) variables given a solution.
//
//   After a timestep is evaluated (or at ICs), there is no way of knowing if
//   secondary variables have been updated to be consistent with the new
//   solution.
// -----------------------------------------------------------------------------
void Richards::CommitStep(double t_old, double t_new, const Tag& tag_next)
{
  // saves primary variable
  PK_PhysicalBDF_Default::CommitStep(t_old, t_new, tag_next);

  AMANZI_ASSERT(tag_next == tag_next_ || tag_next == Tags::NEXT);
  Tag tag_current = tag_next == tag_next_ ? tag_current_ : Tags::CURRENT;

  // also save saturation
  assign(sat_key_, tag_current, tag_next, *S_);
  if (S_->HasRecordSet(sat_ice_key_)) {
    assign(sat_ice_key_, tag_current, tag_next, *S_);
  }
};

// -----------------------------------------------------------------------------
// Check for controls on saturation
// -----------------------------------------------------------------------------
bool
Richards::ValidStep()
{
  Teuchos::OSTab tab = vo_->getOSTab();
  if (vo_->os_OK(Teuchos::VERB_EXTREME))
    *vo_->os() << "Validating time step." << std::endl;

  if (sat_change_limit_ > 0.0) {
    const Epetra_MultiVector& sl_new = *S_->GetPtr<CompositeVector>(sat_key_, tag_next_)
        ->ViewComponent("cell",false);
    const Epetra_MultiVector& sl_old = *S_->GetPtr<CompositeVector>(sat_key_, tag_current_)
        ->ViewComponent("cell",false);
    Epetra_MultiVector dsl(sl_new);
    dsl.Update(-1., sl_old, 1.);
    double change = 0.;
    dsl.NormInf(&change);

    if (change > sat_change_limit_) {
      if (vo_->os_OK(Teuchos::VERB_LOW))
        *vo_->os() << "Invalid time step, max sl change="
                   << change << " > limit=" << sat_change_limit_ << std::endl;
      return false;
    }
  }

  if (S_->HasRecordSet(sat_ice_key_) && (sat_ice_change_limit_ > 0.0)) {
    const Epetra_MultiVector& si_new = *S_->GetPtr<CompositeVector>(sat_ice_key_, tag_next_)
        ->ViewComponent("cell",false);
    const Epetra_MultiVector& si_old = *S_->GetPtr<CompositeVector>(sat_ice_key_, tag_current_)
        ->ViewComponent("cell",false);
    Epetra_MultiVector dsi(si_new);
    dsi.Update(-1., si_old, 1.);
    double change = 0.;
    dsi.NormInf(&change);

    if (change > sat_ice_change_limit_) {
      if (vo_->os_OK(Teuchos::VERB_LOW))
        *vo_->os() << "Invalid time step, max si change="
                   << change << " > limit=" << sat_ice_change_limit_ << std::endl;
      return false;
    }
  }
  return PK_PhysicalBDF_Default::ValidStep();
}

// -----------------------------------------------------------------------------
// Update any diagnostic variables prior to vis (in this case velocity field).
// -----------------------------------------------------------------------------
void Richards::CalculateDiagnostics(const Tag& tag)
{
  AMANZI_ASSERT(tag == Tags::NEXT); // what else would this be?
  Teuchos::OSTab tab = vo_->getOSTab();
  if (vo_->os_OK(Teuchos::VERB_EXTREME))
    *vo_->os() << "Calculating diagnostic variables." << std::endl;

  // update the cell velocities
  UpdateBoundaryConditions_(tag_next_);

  Teuchos::RCP<const CompositeVector> pres = S_->GetPtr<CompositeVector>(key_, tag_next_);
  Teuchos::RCP<const CompositeVector> rel_perm = S_->GetPtr<CompositeVector>(uw_coef_key_, tag_next_);
  Teuchos::RCP<const CompositeVector> rho = S_->GetPtr<CompositeVector>(mass_dens_key_, tag_next_);
  // update the stiffness matrix
  matrix_diff_->SetDensity(rho);
  matrix_diff_->SetScalarCoefficient(rel_perm, Teuchos::null);
  matrix_diff_->UpdateMatrices(Teuchos::null, pres.ptr());
  matrix_diff_->ApplyBCs(true, true, true);

  // derive fluxes
  Teuchos::RCP<CompositeVector> flux = S_->GetPtrW<CompositeVector>(flux_key_, tag_next_, name_);
  matrix_diff_->UpdateFlux(pres.ptr(), flux.ptr());
  UpdateVelocity_(tag);
};

// -----------------------------------------------------------------------------
// Use the physical rel perm (on cells) to update a work vector for rel perm.
//
//   This deals with upwinding, etc.
// -----------------------------------------------------------------------------
bool Richards::UpdatePermeabilityData_(const Tag& tag)
{
  Teuchos::OSTab tab = vo_->getOSTab();
  if (vo_->os_OK(Teuchos::VERB_EXTREME))
    *vo_->os() << "  Updating permeability?";

  Teuchos::RCP<const CompositeVector> rel_perm = S_->GetPtr<CompositeVector>(coef_key_, tag);
  bool update_perm = S_->GetEvaluator(coef_key_, tag)
      .Update(*S_, name_);

  // requirements due to the upwinding method
  if (Krel_method_ == Operators::UPWIND_METHOD_TOTAL_FLUX) {
    bool update_dir = S_->GetEvaluator(mass_dens_key_, tag)
        .Update(*S_, name_);
    update_dir |= S_->GetEvaluator(key_, tag).Update(*S_, name_);

    if (update_dir) {
      // update the direction of the flux -- note this is NOT the flux
      Teuchos::RCP<const CompositeVector> rho = S_->GetPtr<CompositeVector>(mass_dens_key_, tag);
      Teuchos::RCP<CompositeVector> flux_dir = S_->GetPtrW<CompositeVector>(flux_dir_key_, tag, name_);
      Teuchos::RCP<const CompositeVector> pres = S_->GetPtr<CompositeVector>(key_, tag);

      face_matrix_diff_->SetDensity(rho);
      face_matrix_diff_->UpdateMatrices(Teuchos::null, pres.ptr());
      face_matrix_diff_->ApplyBCs(true, true, true);
      face_matrix_diff_->UpdateFlux(pres.ptr(), flux_dir.ptr());

      if (clobber_boundary_flux_dir_) {
        Epetra_MultiVector& flux_dir_f = *flux_dir->ViewComponent("face",false);

        auto& markers = bc_markers();
        auto& values = bc_values();

        for (int f=0; f!=markers.size(); ++f) {
          if (markers[f] == Operators::OPERATOR_BC_NEUMANN) {
            AmanziMesh::Entity_ID_List cells;
            mesh_->face_get_cells(f, AmanziMesh::Parallel_type::ALL, &cells);
            AMANZI_ASSERT(cells.size() == 1);
            int c = cells[0];
            AmanziMesh::Entity_ID_List faces;
            std::vector<int> dirs;
            mesh_->cell_get_faces_and_dirs(c, &faces, &dirs);
            int i = std::find(faces.begin(), faces.end(), f) - faces.begin();

            flux_dir_f[0][f] = values[f]*dirs[i];
          }
        }
      }
    }

    update_perm |= update_dir;
  }

  if (update_perm) {
    Teuchos::RCP<CompositeVector> uw_rel_perm = S_->GetPtrW<CompositeVector>(uw_coef_key_, tag, name_);

    // Move rel perm on boundary_faces into uw_rel_perm on faces
    const Epetra_Import& vandelay = mesh_->exterior_face_importer();
    const Epetra_MultiVector& rel_perm_bf =
        *rel_perm->ViewComponent("boundary_face",false);
    {
      Epetra_MultiVector& uw_rel_perm_f = *uw_rel_perm->ViewComponent("face",false);
      uw_rel_perm_f.Export(rel_perm_bf, vandelay, Insert);
    }

    // Upwind, only overwriting boundary faces if the wind says to do so.
    upwinding_->Update(*rel_perm, *uw_rel_perm, *S_);

    if (clobber_policy_ == "clobber") {
      Epetra_MultiVector& uw_rel_perm_f = *uw_rel_perm->ViewComponent("face",false);
      uw_rel_perm_f.Export(rel_perm_bf, vandelay, Insert);
    } else if (clobber_policy_ == "max") {
      Epetra_MultiVector& uw_rel_perm_f = *uw_rel_perm->ViewComponent("face",false);
      const auto& fmap = mesh_->face_map(true);
      const auto& bfmap = mesh_->exterior_face_map(true);
      for (int bf=0; bf!=rel_perm_bf.MyLength(); ++bf) {
        auto f = fmap.LID(bfmap.GID(bf));
        if (rel_perm_bf[0][bf] > uw_rel_perm_f[0][f]) {
          uw_rel_perm_f[0][f] = rel_perm_bf[0][bf];
        }
      }
    } else if (clobber_policy_ == "unsaturated") {
      // clobber only when the interior cell is unsaturated
      Epetra_MultiVector& uw_rel_perm_f = *uw_rel_perm->ViewComponent("face",false);
      const Epetra_MultiVector& pres = *S_->Get<CompositeVector>(key_, tag).ViewComponent("cell",false);
      const auto& fmap = mesh_->face_map(true);
      const auto& bfmap = mesh_->exterior_face_map(true);
      for (int bf=0; bf!=rel_perm_bf.MyLength(); ++bf) {
        auto f = fmap.LID(bfmap.GID(bf));
        AmanziMesh::Entity_ID_List fcells;
        mesh_->face_get_cells(f, AmanziMesh::Parallel_type::ALL, &fcells);
        AMANZI_ASSERT(fcells.size() == 1);
        if (pres[0][fcells[0]] < 101225.) {
          uw_rel_perm_f[0][f] = rel_perm_bf[0][bf];
        } else if (pres[0][fcells[0]] < 101325.) {
          double frac = (101325. - pres[0][fcells[0]])/100.;
          uw_rel_perm_f[0][f] = rel_perm_bf[0][bf] * frac + uw_rel_perm_f[0][f] * (1-frac);
        }
      }
    }

    if (uw_rel_perm->HasComponent("face"))
      uw_rel_perm->ScatterMasterToGhosted("face");
  }

  // debugging
  if (vo_->os_OK(Teuchos::VERB_EXTREME)) {
    *vo_->os() << " " << update_perm << std::endl;
  }
  return update_perm;
};

bool Richards::UpdatePermeabilityDerivativeData_(const Tag& tag)
{
  Teuchos::OSTab tab = vo_->getOSTab();
  if (vo_->os_OK(Teuchos::VERB_EXTREME))
    *vo_->os() << "  Updating permeability derivatives?";

  bool update_perm = S_->GetEvaluator(coef_key_, tag).UpdateDerivative(*S_, name_, key_, tag);
  if (update_perm) {
    const CompositeVector& drel_perm = S_->GetDerivative<CompositeVector>(coef_key_, tag, key_, tag);

    if (!duw_coef_key_.empty()) {
      // must also upwind
      CompositeVector& duw_rel_perm = S_->GetW<CompositeVector>(duw_coef_key_, tag, name_);
      duw_rel_perm.PutScalar(0.);

      // Upwind, only overwriting boundary faces if the wind says to do so.
      upwinding_deriv_->Update(drel_perm, duw_rel_perm, *S_);
      duw_rel_perm.ScatterMasterToGhosted("face");
    } else {
      drel_perm.ScatterMasterToGhosted("cell");
    }
  }

  // debugging
  if (vo_->os_OK(Teuchos::VERB_EXTREME)) {
    *vo_->os() << " " << update_perm << std::endl;
  }
  return update_perm;
};

// -----------------------------------------------------------------------------
// Compute boundary condition functions at the current time.
// -----------------------------------------------------------------------------
void Richards::ComputeBoundaryConditions_(const Tag& tag)
{
  bc_pressure_->Compute(S_->get_time(tag));
  bc_head_->Compute(S_->get_time(tag));
  bc_level_->Compute(S_->get_time(tag));
  bc_flux_->Compute(S_->get_time(tag));
  bc_seepage_->Compute(S_->get_time(tag));
  bc_seepage_infilt_->Compute(S_->get_time(tag));
}

// -----------------------------------------------------------------------------
// Push boundary conditions into the global array.
// -----------------------------------------------------------------------------
void Richards::UpdateBoundaryConditions_(const Tag& tag, bool kr)
{
  Teuchos::OSTab tab = vo_->getOSTab();
  if (vo_->os_OK(Teuchos::VERB_EXTREME))
    *vo_->os() << "  Updating BCs." << std::endl;

  auto& markers = bc_markers();
  auto& values = bc_values();

  // initialize all to 0
  for (unsigned int n=0; n!=markers.size(); ++n) {
    markers[n] = Operators::OPERATOR_BC_NONE;
    values[n] = 0.0;
  }

  // count for debugging
  std::vector<int> bc_counts;
  std::vector<std::string> bc_names;

  // Dirichlet-type boundary conditions
  // -------------------------------------
  // pressure boundary conditions -- the primary
  bc_counts.push_back(bc_pressure_->size());
  bc_names.push_back(key_);
  for (const auto& bc : *bc_pressure_) {
    int f = bc.first;
#ifdef ENABLE_DBC
    AmanziMesh::Entity_ID_List cells;
    mesh_->face_get_cells(f, AmanziMesh::Parallel_type::ALL, &cells);
    AMANZI_ASSERT(cells.size() == 1);
#endif
    markers[f] = Operators::OPERATOR_BC_DIRICHLET;
    values[f] = bc.second;
  }

  // head boundary conditions
  bc_counts.push_back(bc_head_->size());
  bc_names.push_back("head");
  if (bc_head_->size() > 0) {
    double p_atm = S_->Get<double>("atmospheric_pressure", Tags::DEFAULT);
    int z_index = mesh_->space_dimension() - 1;
    const auto& gravity = S_->Get<AmanziGeometry::Point>("gravity", Tags::DEFAULT);
    double g = -gravity[z_index];

    for (const auto& bc : *bc_head_) {
      int f = bc.first;

      AmanziMesh::Entity_ID_List cells;
      mesh_->face_get_cells(f, AmanziMesh::Parallel_type::ALL, &cells);
      AMANZI_ASSERT(cells.size() == 1);

      // we need to find the elevation of the surface, but finding the top edge
      // of this stack of faces is not possible currently.  The best approach
      // is instead to work with the cell.
      int col = mesh_->column_ID(cells[0]);
      double z_surf = mesh_->face_centroid(mesh_->faces_of_column(col)[0])[z_index];
      double z_wt = bc.second + z_surf;

      markers[f] = Operators::OPERATOR_BC_DIRICHLET;
      // note, here the cell centroid's z is used to relate to the column's top
      // face centroid, specifically NOT the boundary face's centroid.
      values[f] = p_atm + bc_rho_water_ * g *
        (z_wt - mesh_->cell_centroid(cells[0])[z_index]);
    }
  }

  // fixed level head boundary conditions
  bc_counts.push_back(bc_head_->size());
  bc_names.push_back("head");
  if (bc_level_->size() > 0) {
    double p_atm = S_->Get<double>("atmospheric_pressure", Tags::DEFAULT);
    int z_index = mesh_->space_dimension() - 1;
    const auto& gravity = S_->Get<AmanziGeometry::Point>("gravity", Tags::DEFAULT);
    double g = -gravity[z_index];

    for (const auto& bc : *bc_level_) {
      int f = bc.first;
      markers[f] = Operators::OPERATOR_BC_DIRICHLET;
      values[f] = p_atm + bc_rho_water_ * g *
        (bc.second - mesh_->face_centroid(f)[z_index]);
    }
  }

  // Neumann type boundary conditions
  // -------------------------------------
  const Epetra_MultiVector& rel_perm =
    *S_->Get<CompositeVector>(uw_coef_key_, tag).ViewComponent("face",false);

  // standard Neumann flux BCs
  bc_counts.push_back(bc_flux_->size());
  bc_names.push_back(flux_key_);

  if (!infiltrate_only_if_unfrozen_) {
    for (const auto& bc : *bc_flux_) {
      int f = bc.first;
#ifdef ENABLE_DBC
      AmanziMesh::Entity_ID_List cells;
      mesh_->face_get_cells(f, AmanziMesh::Parallel_type::ALL, &cells);
      AMANZI_ASSERT(cells.size() == 1);
#endif
      markers[f] = Operators::OPERATOR_BC_NEUMANN;
      values[f] = bc.second;
      if (!kr && rel_perm[0][f] > 0.) values[f] /= rel_perm[0][f];
    }

  } else {
    // Neumann boundary conditions that turn off if temp < freezing
    const Epetra_MultiVector& temp =
      *S_->GetPtr<CompositeVector>(Keys::getKey(domain_,"temperature"), tag)
         ->ViewComponent("face");
    for (const auto& bc : *bc_flux_) {
      int f = bc.first;
#ifdef ENABLE_DBC
      AmanziMesh::Entity_ID_List cells;
      mesh_->face_get_cells(f, AmanziMesh::Parallel_type::ALL, &cells);
      AMANZI_ASSERT(cells.size() == 1);
#endif
      markers[f] = Operators::OPERATOR_BC_NEUMANN;
      if (temp[0][f] > 273.15) {
        values[f] = bc.second;
        if (!kr && rel_perm[0][f] > 0.) values[f] /= rel_perm[0][f];
      } else {
        values[f] = 0.;
      }
    }
  }

  // seepage face -- pressure <= specified value (usually 101325), outward water flux >= 0
  S_->Get<CompositeVector>(flux_key_, tag).ScatterMasterToGhosted("face");
  const Epetra_MultiVector& flux = *S_->Get<CompositeVector>(flux_key_, tag)
      .ViewComponent("face", true);

  const double& p_atm = S_->Get<double>("atmospheric_pressure", Tags::DEFAULT);
  Teuchos::RCP<const CompositeVector> u = S_->GetPtr<CompositeVector>(key_, tag);
  double seepage_tol = 10.;

  bc_counts.push_back(bc_seepage_->size());
  bc_names.push_back("standard seepage");
  for (const auto& bc : *bc_seepage_) {
    int f = bc.first;
#ifdef ENABLE_DBC
    AmanziMesh::Entity_ID_List cells;
    mesh_->face_get_cells(f, AmanziMesh::Parallel_type::ALL, &cells);
    AMANZI_ASSERT(cells.size() == 1);
#endif

    //double boundary_pressure = std::max(getFaceOnBoundaryValue(f, *u, *bc_), 101325.); // does not make sense to seep from nonsaturated cells
    double boundary_pressure = getFaceOnBoundaryValue(f, *u, *bc_); // does not make sense to seep from nonsaturated cells
    double boundary_flux = flux[0][f]*getBoundaryDirection(*mesh_, f);
    if (boundary_pressure > bc.second) {
      markers[f] = Operators::OPERATOR_BC_DIRICHLET;
      values[f] = bc.second;
    } else if (boundary_pressure < bc.second - seepage_tol) {
      markers[f] = Operators::OPERATOR_BC_NEUMANN;
      values[f] = 0.;
    } else if (boundary_flux >= 0.) {
      markers[f] = Operators::OPERATOR_BC_DIRICHLET;
      values[f] = bc.second;
    } else {
      markers[f] = Operators::OPERATOR_BC_NEUMANN;
      values[f] = 0.;
    }
  }

  // seepage face -- pressure <= p_atm, outward water flux is specified
  bc_counts.push_back(bc_seepage_infilt_->size());
  bc_names.push_back("seepage with infiltration");
  {
    Tag seepage_tag;
    if (bc_seepage_infilt_explicit_) {
      seepage_tag = tag_current_;
    } else {
      seepage_tag = tag;
    }
    const Epetra_MultiVector& flux = *S_->Get<CompositeVector>(flux_key_, seepage_tag).ViewComponent("face", true);
    Teuchos::RCP<const CompositeVector> u = S_->GetPtr<CompositeVector>(key_, seepage_tag);

    int i = 0;
    for (const auto& bc : *bc_seepage_infilt_) {
      int f = bc.first;
#ifdef ENABLE_DBC
      AmanziMesh::Entity_ID_List cells;
      mesh_->face_get_cells(f, AmanziMesh::Parallel_type::ALL, &cells);
      AMANZI_ASSERT(cells.size() == 1);
#endif

      double flux_seepage_tol = std::abs(bc.second) * .001;
      double boundary_pressure = getFaceOnBoundaryValue(f, *u, *bc_);
      double boundary_flux = flux[0][f]*getBoundaryDirection(*mesh_, f);
      if (i == 0)
        std::cout << "BFlux = " << boundary_flux << " with constraint = " << bc.second - flux_seepage_tol << std::endl;

      if (boundary_flux < bc.second - flux_seepage_tol &&
          boundary_pressure > p_atm + seepage_tol) {
        // both constraints are violated, either option should push things in the right direction
        markers[f] = Operators::OPERATOR_BC_DIRICHLET;
        values[f] = p_atm;
        if (i == 0)
          std::cout << "BC PRESSURE ON SEEPAGE = " << boundary_pressure << " with flux " << flux[0][f]*getBoundaryDirection(*mesh_, f) << " resulted in DIRICHLET pressure " << p_atm << std::endl;

      } else if (boundary_flux >= bc.second - flux_seepage_tol &&
                 boundary_pressure > p_atm - seepage_tol) {
        // max pressure condition violated
        markers[f] = Operators::OPERATOR_BC_DIRICHLET;
        values[f] = p_atm;
        if (i == 0)
          std::cout << "BC PRESSURE ON SEEPAGE = " << boundary_pressure << " with flux " << flux[0][f]*getBoundaryDirection(*mesh_, f) << " resulted in DIRICHLET pressure " << p_atm << std::endl;

      } else if (boundary_flux < bc.second - flux_seepage_tol &&
                 boundary_pressure <= p_atm + seepage_tol) {
        // max infiltration violated
        markers[f] = Operators::OPERATOR_BC_NEUMANN;
        values[f] = bc.second;
        if (i == 0)
          std::cout << "BC PRESSURE ON SEEPAGE = " << boundary_pressure << " with flux " << flux[0][f]*getBoundaryDirection(*mesh_, f) << " resulted in NEUMANN flux " << bc.second << std::endl;

      } else if (boundary_flux >= bc.second - flux_seepage_tol &&
                 boundary_pressure <= p_atm - seepage_tol) {
        // both conditions are valid
        markers[f] = Operators::OPERATOR_BC_NEUMANN;
        values[f] = bc.second;
        if (i == 0)
          std::cout << "BC PRESSURE ON SEEPAGE = " << boundary_pressure << " with flux " << flux[0][f]*getBoundaryDirection(*mesh_, f) << " resulted in NEUMANN flux " << bc.second << std::endl;

      } else {
        AMANZI_ASSERT(0);
      }
      i++;
    }
  }

  // surface coupling
  bc_counts.push_back(0);
  bc_names.push_back("surface coupling (head)");

  if (coupled_to_surface_via_head_) {
    // Face is Dirichlet with value of surface head
    Teuchos::RCP<const AmanziMesh::Mesh> surface = S_->GetMesh("surface");
    const Epetra_MultiVector& head = *S_->GetPtr<CompositeVector>("surface_pressure", tag)
        ->ViewComponent("cell",false);

    unsigned int ncells_surface = head.MyLength();
    bc_counts[bc_counts.size()-1] = ncells_surface;

    for (unsigned int c=0; c!=ncells_surface; ++c) {
      // -- get the surface cell's equivalent subsurface face
      AmanziMesh::Entity_ID f = surface->entity_get_parent(AmanziMesh::CELL, c);
#ifdef ENABLE_DBC
      AmanziMesh::Entity_ID_List cells;
      mesh_->face_get_cells(f, AmanziMesh::Parallel_type::ALL, &cells);
      AMANZI_ASSERT(cells.size() == 1);
#endif
      // -- set that value to dirichlet
      markers[f] = Operators::OPERATOR_BC_DIRICHLET;
      values[f] = head[0][c];
    }
  }

  // surface coupling
  bc_counts.push_back(0);
  bc_names.push_back("surface coupling (flux)");
  if (coupled_to_surface_via_flux_) {
    // Face is Neumann with value of surface residual
    Teuchos::RCP<const AmanziMesh::Mesh> surface = S_->GetMesh(Keys::getDomain(ss_flux_key_));
    const Epetra_MultiVector& ss_flux = *S_->Get<CompositeVector>(ss_flux_key_, tag).ViewComponent("cell",false);
    unsigned int ncells_surface = ss_flux.MyLength();
    bc_counts[bc_counts.size()-1] = ncells_surface;
    for (unsigned int c=0; c!=ncells_surface; ++c) {
      // -- get the surface cell's equivalent subsurface face
      AmanziMesh::Entity_ID f = surface->entity_get_parent(AmanziMesh::CELL, c);
#ifdef ENABLE_DBC
      AmanziMesh::Entity_ID_List cells;
      mesh_->face_get_cells(f, AmanziMesh::Parallel_type::ALL, &cells);
      AMANZI_ASSERT(cells.size() == 1);
#endif
      // -- set that value to Neumann
      markers[f] = Operators::OPERATOR_BC_NEUMANN;

      // NOTE: the flux provided by the coupler is in units of mols / s, where
      //       as Neumann BCs are in units of mols / s / A.  The right A must
      //       be chosen, as it is the subsurface mesh's face area, not the
      //       surface mesh's cell area.
      values[f] = ss_flux[0][c] / mesh_->face_area(f);

      if (!kr && rel_perm[0][f] > 0.) values[f] /= rel_perm[0][f];
    }
  }

  // mark all remaining boundary conditions as zero flux conditions
  AmanziMesh::Entity_ID_List cells;
  int n_default = 0;
  int nfaces_owned = mesh_->num_entities(AmanziMesh::FACE, AmanziMesh::Parallel_type::OWNED);
  for (int f = 0; f < nfaces_owned; f++) {
    if (markers[f] == Operators::OPERATOR_BC_NONE) {
      mesh_->face_get_cells(f, AmanziMesh::Parallel_type::ALL, &cells);
      int ncells = cells.size();

      if (ncells == 1) {
        n_default++;
        markers[f] = Operators::OPERATOR_BC_NEUMANN;
        values[f] = 0.0;
      }
    }
  }
  bc_names.push_back("default (zero flux)");
  bc_counts.push_back(n_default);

  // report on counts
  if (vo_->os_OK(Teuchos::VERB_HIGH)) {
    std::vector<int> bc_counts_global(bc_counts.size(), 0);
    mesh_->get_comm()->SumAll(&bc_counts[0], &bc_counts_global[0], bc_counts.size());

    *vo_->os() << "  BCs applied:" << std::endl;

    for (int i=0; i!=bc_counts_global.size(); ++i) {
      *vo_->os() << "    " << bc_names[i] << ": " << bc_counts_global[i] << std::endl;
    }
  }
};

bool Richards::ModifyPredictor(double h, Teuchos::RCP<const TreeVector> u0,
        Teuchos::RCP<TreeVector> u)
{
  Teuchos::OSTab tab = vo_->getOSTab();
  if (vo_->os_OK(Teuchos::VERB_EXTREME))
    *vo_->os() << "Modifying predictor:" << std::endl;

  // update boundary conditions
  ComputeBoundaryConditions_(tag_next_);
  UpdateBoundaryConditions_(tag_next_);
  db_->WriteBoundaryConditions(bc_markers(), bc_values());

  // push Dirichlet data into predictor
  applyDirichletBCs(*bc_, *u->Data());

  bool changed(false);
  if (modify_predictor_bc_flux_ ||
      (modify_predictor_first_bc_flux_ &&
       ((S_->Get<int>("cycle", Tags::DEFAULT) == 0) || (S_->Get<int>("cycle", Tags::DEFAULT) == 1)))) {
    changed |= ModifyPredictorFluxBCs_(h,u);
  }

  if (modify_predictor_wc_) {
    changed |= ModifyPredictorWC_(h,u);
  }

  if (modify_predictor_with_consistent_faces_) {
    changed |= ModifyPredictorConsistentFaces_(h,u);
  }
  return changed;
}

bool Richards::ModifyPredictorFluxBCs_(double h, Teuchos::RCP<TreeVector> u)
{
  if (!u->Data()->HasComponent("face")) return false;

  auto& markers = bc_markers();
  auto& values = bc_values();

  Teuchos::OSTab tab = vo_->getOSTab();
  if (vo_->os_OK(Teuchos::VERB_EXTREME))
    *vo_->os() << "  modifications to deal with nonlinearity at flux BCs" << std::endl;

  if (flux_predictor_ == Teuchos::null) {
    flux_predictor_ = Teuchos::rcp(new PredictorDelegateBCFlux(S_, mesh_, matrix_diff_,
            wrms_, &markers, &values));
  }

  UpdatePermeabilityData_(tag_next_);
  Teuchos::RCP<const CompositeVector> rel_perm =
    S_->GetPtr<CompositeVector>(uw_coef_key_, tag_next_);

  matrix_->Init();
  matrix_diff_->SetScalarCoefficient(rel_perm, Teuchos::null);
  Teuchos::RCP<const CompositeVector> rho = S_->GetPtr<CompositeVector>(mass_dens_key_, tag_next_);
  Teuchos::RCP<const CompositeVector> pres = S_->GetPtr<CompositeVector>(key_, tag_next_);
  matrix_diff_->SetDensity(rho);
  matrix_diff_->UpdateMatrices(Teuchos::null, pres.ptr());
  //matrix_diff_->ApplyBCs(true, true, true);

  flux_predictor_->ModifyPredictor(h, u);
  ChangedSolution(); // mark the solution as changed, as modifying with
                      // consistent faces will then get the updated boundary
                      // conditions
  return true;
}

bool Richards::ModifyPredictorConsistentFaces_(double h, Teuchos::RCP<TreeVector> u)
{
  Teuchos::OSTab tab = vo_->getOSTab();
  if (vo_->os_OK(Teuchos::VERB_EXTREME))
    *vo_->os() << "  modifications for consistent face pressures." << std::endl;

  CalculateConsistentFaces(u->Data().ptr());
  return true;
}

bool Richards::ModifyPredictorWC_(double h, Teuchos::RCP<TreeVector> u)
{
  AMANZI_ASSERT(0);
  return false;
}

// void Richards::CalculateConsistentFacesForInfiltration_(
//     const Teuchos::Ptr<CompositeVector>& u) {

//  auto& markers = bc_markers();
//  auto& values = bc_values();

//   if (vo_->os_OK(Teuchos::VERB_EXTREME))
//     *vo_->os() << "  modifications to deal with nonlinearity at flux BCs" << std::endl;

//   if (flux_predictor_ == Teuchos::null) {
//     flux_predictor_ = Teuchos::rcp(new PredictorDelegateBCFlux(S_next_, mesh_, matrix_,
//             wrms_, &markers, &values));
//   }

//   // update boundary conditions
//   bc_pressure_->Compute(S_->get_time(tag_next_));
//   bc_flux_->Compute(S_->get_time(tag_next_));
//   UpdateBoundaryConditions_(S_next_.ptr());

//   bool update = UpdatePermeabilityData_(S_next_.ptr());
//   Teuchos::RCP<const CompositeVector> rel_perm =
//       S_next_->GetPtr<CompositeVector>(uw_coef_key_);
//   matrix_->CreateMFDstiffnessMatrices(rel_perm.ptr());
//   matrix_->CreateMFDrhsVectors();
//   Teuchos::RCP<const CompositeVector> rho = S_next_->GetPtr<CompositeVector>(mass_dens_key_);
//   Teuchos::RCP<const Epetra_Vector> gvec = S_next_->GetConstantVectorData("gravity", Tags::DEFAULT);
//   AddGravityFluxes_(gvec.ptr(), rel_perm.ptr(), rho.ptr(), matrix_.ptr());
//   matrix_->ApplyBoundaryConditions(markers, values);

//   flux_predictor_->ModifyPredictor(u);
// }

void Richards::CalculateConsistentFaces(const Teuchos::Ptr<CompositeVector>& u)
{
  if (!u->HasComponent("face")) return; // not need

  // VerboseObject stuff.
  Teuchos::OSTab tab = vo_->getOSTab();
  if (vo_->os_OK(Teuchos::VERB_EXTREME))
    *vo_->os() << "  Modifying predictor for consistent faces" << std::endl;

  // average cells to faces to give a reasonable initial guess
  u->ScatterMasterToGhosted("cell");
  const Epetra_MultiVector& u_c = *u->ViewComponent("cell",true);
  Epetra_MultiVector& u_f = *u->ViewComponent("face",false);

  int f_owned = u_f.MyLength();
  for (int f=0; f!=f_owned; ++f) {
    AmanziMesh::Entity_ID_List cells;
    mesh_->face_get_cells(f, AmanziMesh::Parallel_type::ALL, &cells);
    int ncells = cells.size();

    double face_value = 0.0;
    for (int n=0; n!=ncells; ++n) {
      face_value += u_c[0][cells[n]];
    }
    u_f[0][f] = face_value / ncells;
  }
  ChangedSolution();

  // Using the old BCs, so should use the old rel perm?
  // update the rel perm according to the scheme of choice
  //  UpdatePermeabilityData_(tag_next_);
  Teuchos::RCP<const CompositeVector> rel_perm =
    S_->GetPtr<CompositeVector>(uw_coef_key_, tag_next_);

  S_->GetEvaluator(mass_dens_key_, tag_next_)
      .Update(*S_, name_);
  Teuchos::RCP<const CompositeVector> rho =
    S_->GetPtr<CompositeVector>(mass_dens_key_, tag_next_);

  // Update the preconditioner with darcy and gravity fluxes
  matrix_->Init();
  matrix_diff_->SetDensity(rho);
  matrix_diff_->SetScalarCoefficient(rel_perm, Teuchos::null);
  matrix_diff_->UpdateMatrices(Teuchos::null, u);
  matrix_diff_->ApplyBCs(true, true, true);

  // derive the consistent faces, involves a solve
  db_->WriteVector(" p_cf guess:", u.ptr(), true);
  matrix_diff_->UpdateConsistentFaces(*u);
  db_->WriteVector(" p_cf soln:", u.ptr(), true);
}

// /* ******************************************************************
// * Clip pressure using pressure threshold.
// ****************************************************************** */
// void Richards::ClipHydrostaticPressure(double pmin, Epetra_MultiVector& p)
// {
//   int ncells_owned = mesh_->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED);
//   for (int c = 0; c < ncells_owned; c++) p[0][c] = std::max(p[0][c], pmin);
// }

// -----------------------------------------------------------------------------
// Check admissibility of the solution guess.
// -----------------------------------------------------------------------------
bool Richards::IsAdmissible(Teuchos::RCP<const TreeVector> up)
{
  Teuchos::OSTab tab = vo_->getOSTab();
  if (vo_->os_OK(Teuchos::VERB_EXTREME))
    *vo_->os() << "  Checking admissibility..." << std::endl;

  // For some reason, wandering PKs break most frequently with an unreasonable
  // pressure.  This simply tries to catch that before it happens.
  Teuchos::RCP<const CompositeVector> pres = up->Data();
  double minT, maxT;

  const Epetra_MultiVector& pres_c = *pres->ViewComponent("cell",false);
  double minT_c(1.e15), maxT_c(-1.e15);
  int min_c(-1), max_c(-1);
  for (int c=0; c!=pres_c.MyLength(); ++c) {
    if (pres_c[0][c] < minT_c) {
      minT_c = pres_c[0][c];
      min_c = c;
    }
    if (pres_c[0][c] > maxT_c) {
      maxT_c = pres_c[0][c];
      max_c = c;
    }
  }

  double minT_f(1.e15), maxT_f(-1.e15);
  int min_f(-1), max_f(-1);
  if (pres->HasComponent("face")) {
    const Epetra_MultiVector& pres_f = *pres->ViewComponent("face",false);
    for (int f=0; f!=pres_f.MyLength(); ++f) {
      if (pres_f[0][f] < minT_f) {
        minT_f = pres_f[0][f];
        min_f = f;
      }
      if (pres_f[0][f] > maxT_f) {
        maxT_f = pres_f[0][f];
        max_f = f;
      }
    }
    minT = std::min(minT_c, minT_f);
    maxT = std::max(maxT_c, maxT_f);

  } else {
    minT = minT_c;
    maxT = maxT_c;
  }

  double minT_l = minT;
  double maxT_l = maxT;
  mesh_->get_comm()->MaxAll(&maxT_l, &maxT, 1);
  mesh_->get_comm()->MinAll(&minT_l, &minT, 1);

  if (vo_->os_OK(Teuchos::VERB_HIGH)) {
    *vo_->os() << "    Admissible p? (min/max): " << minT << ",  " << maxT << std::endl;
  }

  if (minT < -1.e9 || maxT > 1.e8) {
    if (vo_->os_OK(Teuchos::VERB_MEDIUM)) {
      *vo_->os() << " is not admissible, as it is not within bounds of constitutive models:" << std::endl;

      Teuchos::RCP<const Comm_type> comm_p = mesh_->get_comm();
      Teuchos::RCP<const MpiComm_type> mpi_comm_p =
        Teuchos::rcp_dynamic_cast<const MpiComm_type>(comm_p);
      const MPI_Comm& comm = mpi_comm_p->Comm();

      ENorm_t global_minT_c, local_minT_c;
      ENorm_t global_maxT_c, local_maxT_c;

      local_minT_c.value = minT_c;
      local_minT_c.gid = pres_c.Map().GID(min_c);
      local_maxT_c.value = maxT_c;
      local_maxT_c.gid = pres_c.Map().GID(max_c);

      MPI_Allreduce(&local_minT_c, &global_minT_c, 1, MPI_DOUBLE_INT, MPI_MINLOC, comm);
      MPI_Allreduce(&local_maxT_c, &global_maxT_c, 1, MPI_DOUBLE_INT, MPI_MAXLOC, comm);
      *vo_->os() << "   cells (min/max): [" << global_minT_c.gid << "] " << global_minT_c.value
                 << ", [" << global_maxT_c.gid << "] " << global_maxT_c.value << std::endl;

      if (pres->HasComponent("face")) {
        const Epetra_MultiVector& pres_f = *pres->ViewComponent("face",false);
        ENorm_t global_minT_f, local_minT_f;
        ENorm_t global_maxT_f, local_maxT_f;

        local_minT_f.value = minT_f;
        local_minT_f.gid = pres_f.Map().GID(min_f);
        local_maxT_f.value = maxT_f;
        local_maxT_f.gid = pres_f.Map().GID(max_f);

        MPI_Allreduce(&local_minT_f, &global_minT_f, 1, MPI_DOUBLE_INT, MPI_MINLOC, comm);
        MPI_Allreduce(&local_maxT_f, &global_maxT_f, 1, MPI_DOUBLE_INT, MPI_MAXLOC, comm);
        *vo_->os() << "   cells (min/max): [" << global_minT_f.gid << "] " << global_minT_f.value;
        MPI_Allreduce(&local_minT_f, &global_minT_f, 1, MPI_DOUBLE_INT, MPI_MINLOC, comm);
        MPI_Allreduce(&local_maxT_f, &global_maxT_f, 1, MPI_DOUBLE_INT, MPI_MAXLOC, comm);
        *vo_->os() << "   faces (min/max): [" << global_minT_f.gid << "] " << global_minT_f.value
                   << ", [" << global_maxT_f.gid << "] " << global_maxT_f.value << std::endl;
      }
    }
    return false;
  }
  return true;
}

AmanziSolvers::FnBaseDefs::ModifyCorrectionResult
Richards::ModifyCorrection(double h, Teuchos::RCP<const TreeVector> res,
                 Teuchos::RCP<const TreeVector> u,
                 Teuchos::RCP<TreeVector> du)
{
  Teuchos::OSTab tab = vo_->getOSTab();

  // if the primary variable has boundary face, this is for upwinding rel
  // perms and is never actually used.  Make sure it does not go to undefined
  // pressures.
  if (du->Data()->HasComponent("boundary_face"))  {
    du->Data()->ViewComponent("boundary_face")->PutScalar(0.);
  }

  // debugging -- remove me! --etc
  for (CompositeVector::name_iterator comp=du->Data()->begin();
       comp!=du->Data()->end(); ++comp) {
    Epetra_MultiVector& du_c = *du->Data()->ViewComponent(*comp,false);
    double max, l2;
    du_c.NormInf(&max);
    du_c.Norm2(&l2);
    if (vo_->os_OK(Teuchos::VERB_HIGH)) {
      *vo_->os() << "Linf, L2 pressure correction (" << *comp << ") = " << max << ", " << l2 << std::endl;
    }
  }

  // limit by capping corrections when they cross atmospheric pressure
  // (where pressure derivatives are discontinuous)
  int my_limited = 0;
  int n_limited_spurt = 0;
  if (patm_limit_ > 0.) {
    double patm = S_->Get<double>("atmospheric_pressure", Tags::DEFAULT);
    for (CompositeVector::name_iterator comp=du->Data()->begin();
         comp!=du->Data()->end(); ++comp) {
      Epetra_MultiVector& du_c = *du->Data()->ViewComponent(*comp,false);
      const Epetra_MultiVector& u_c = *u->Data()->ViewComponent(*comp,false);

      for (int c=0; c!=du_c.MyLength(); ++c) {
        if ((u_c[0][c] < patm) &&
            (u_c[0][c] - du_c[0][c] > patm + patm_limit_)) {
          du_c[0][c] = u_c[0][c] - (patm + patm_limit_);
          my_limited++;
        } else if ((u_c[0][c] > patm) &&
                   (u_c[0][c] - du_c[0][c] < patm - patm_limit_)) {
          du_c[0][c] = u_c[0][c] - (patm - patm_limit_);
          my_limited++;
        }
      }
    }
    mesh_->get_comm()->MaxAll(&my_limited, &n_limited_spurt, 1);
  }

  if (n_limited_spurt > 0) {
    if (vo_->os_OK(Teuchos::VERB_HIGH)) {
      *vo_->os() << "  limiting the spurt." << std::endl;
    }
  }

  // debugging -- remove me! --etc
  for (CompositeVector::name_iterator comp=du->Data()->begin();
       comp!=du->Data()->end(); ++comp) {
    Epetra_MultiVector& du_c = *du->Data()->ViewComponent(*comp,false);
    double max, l2;
    du_c.NormInf(&max);
    du_c.Norm2(&l2);
    if (vo_->os_OK(Teuchos::VERB_HIGH)) {
      *vo_->os() << "Linf, L2 pressure correction (" << *comp << ") = " << max << ", " << l2 << std::endl;
    }
  }

  // Limit based on a max pressure change
  my_limited = 0;
  int n_limited_change = 0;
  if (p_limit_ >= 0.) {
    for (CompositeVector::name_iterator comp=du->Data()->begin();
         comp!=du->Data()->end(); ++comp) {
      Epetra_MultiVector& du_c = *du->Data()->ViewComponent(*comp,false);

      double max;
      du_c.NormInf(&max);
      if (vo_->os_OK(Teuchos::VERB_HIGH)) {
        *vo_->os() << "Max pressure correction (" << *comp << ") = " << max << std::endl;
      }

      for (int c=0; c!=du_c.MyLength(); ++c) {
        if (std::abs(du_c[0][c]) > p_limit_) {
          du_c[0][c] = ((du_c[0][c] > 0) - (du_c[0][c] < 0)) * p_limit_;
          my_limited++;
        }
      }
    }

    mesh_->get_comm()->MaxAll(&my_limited, &n_limited_change, 1);
  }

  // debugging -- remove me! --etc
  for (CompositeVector::name_iterator comp=du->Data()->begin();
       comp!=du->Data()->end(); ++comp) {
    Epetra_MultiVector& du_c = *du->Data()->ViewComponent(*comp,false);
    double max, l2;
    du_c.NormInf(&max);
    du_c.Norm2(&l2);
    if (vo_->os_OK(Teuchos::VERB_HIGH)) {
      *vo_->os() << "Linf, L2 pressure correction (" << *comp << ") = " << max << ", " << l2 << std::endl;
    }
  }

  if (n_limited_change > 0) {
    if (vo_->os_OK(Teuchos::VERB_HIGH)) {
      *vo_->os() << "  limited by pressure." << std::endl;
    }
  }

  if (n_limited_spurt > 0) {
    return AmanziSolvers::FnBaseDefs::CORRECTION_MODIFIED_LAG_BACKTRACKING;
  } else if (n_limited_change > 0) {
    return AmanziSolvers::FnBaseDefs::CORRECTION_MODIFIED;
  }

  return AmanziSolvers::FnBaseDefs::CORRECTION_NOT_MODIFIED;
}

} // namespace
} // namespace


text:
// Copyright 2018 The Chromium Authors. All rights reserved.
// Use of this source code is governed by a BSD-style license that can be
// found in the LICENSE file.

package org.chromium.chrome.browser.feed;

/**
 * Aggregation point for application lifecycle events that the Feed cares about.
 */
public class FeedAppLifecycle {
    /**
     * Create a FeedAppLifecycle instance. In normal use, this should only be called by {@link
     * FeedAppLifecycleFactory}.
     */
    public FeedAppLifecycle() {}

    /**
     * This is called when a tabbed activity is launched.
     */
    public void onActivityLaunched() {}
}

redpajama_set_name: RedPajamaGithub


text:
var mongoose = require('mongoose'),
  RepositorySchema = require('./repository.schema');

var UserSchema = mongoose.Schema({
  id: String,
  userName: String,
  userAvatar: String,
  repositories: [{ type: mongoose.Schema.Types.ObjectId, ref: 'Repository' }]
});

module.exports = UserSchema;


text:
import { NgModule, CUSTOM_ELEMENTS_SCHEMA }      from '@angular/core';
import { SharedModule }   from '../shared/shared.module';
import { CommonModule } from '@angular/common';

import { HeaderComponent }  from './header.component';
import { FooterComponent }  from './footer.component';
import { HomeComponent } from './home.component';

import { Routes, RouterModule } from '@angular/router';
import { HomeRoutes } from './home.routes';

@NgModule({
  imports:      [
    CommonModule,
    HomeRoutes,
    SharedModule
  ],
  declarations: [ HeaderComponent, HomeComponent, FooterComponent],
  schemas: [CUSTOM_ELEMENTS_SCHEMA]
})

export default class HomeModule { }


text:
// <copyright file="Matcher.cs" company="Hottinger Baldwin Messtechnik GmbH">
//
// SharpJet, a library to communicate with Jet IPC.
//
// The MIT License (MIT)
//
// Copyright (C) Hottinger Baldwin Messtechnik GmbH
//
// Permission is hereby granted, free of charge, to any person obtaining
// a copy of this software and associated documentation files (the
// "Software"), to deal in the Software without restriction, including
// without limitation the rights to use, copy, modify, merge, publish,
// distribute, sublicense, and/or sell copies of the Software, and to
// permit persons to whom the Software is furnished to do so, subject to
// the following conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
// MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
// BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
// ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
// CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
//
// </copyright>

namespace Hbm.Devices.Jet
{
    /// <summary>
    /// This class holds the properties required to fill a path matcher.
    /// </summary>
    public class Matcher
    {
        public string Contains { get; set; }

        public string[] ContainsAllOf { get; set; }

        public string StartsWith { get; set; }

        public string EndsWith { get; set; }

        public string EqualsTo { get; set; }

        public string EqualsNotTo { get; set; }

        public bool CaseInsensitive { get; set; }
    }
}


text:
set(CMAKE_RELATIVE_PATH_TOP_SOURCE "/home/leonidlezner/coffeeio")
set(CMAKE_RELATIVE_PATH_TOP_BINARY "/home/leonidlezner/coffeeio/build")

# Force unix paths in dependencies.
set(CMAKE_FORCE_UNIX_PATHS 1)

# The C and CXX include file regular expressions for this directory.
set(CMAKE_C_INCLUDE_REGEX_SCAN "^.*$")
set(CMAKE_C_INCLUDE_REGEX_COMPLAIN "^$")
set(CMAKE_CXX_INCLUDE_REGEX_SCAN ${CMAKE_C_INCLUDE_REGEX_SCAN})
set(CMAKE_CXX_INCLUDE_REGEX_COMPLAIN ${CMAKE_C_INCLUDE_REGEX_COMPLAIN})


text:
<?php

namespace VIPSoft\CodeCoverageCommon\Driver;

use PHP_CodeCoverage_Driver as DriverInterface;

/**
 * Stub driver
 *
 * @author Anthon Pang <apang@softwaredevelopment.ca>
 */
class Stub implements DriverInterface
{
    private $driver;

    /**
     * Register driver
     *
     * @param DriverInterface $driver
     */
    public function setDriver(DriverInterface $driver)
    {
        $this->driver = $driver;
    }

    /**
     * Get driver
     *
     * @return DriverInterface $driver
     */
    public function getDriver()
    {
        return $this->driver;
    }

    /**
     * {@inheritdoc}
     */
    public function start()
    {
        if ($this->driver) {
            $this->driver->start();
        }
    }

    /**
     * {@inheritdoc}
     */
    public function stop()
    {
        return $this->driver ? $this->driver->stop() : false;
    }
}


text:
package com.adarsh.KeyValueStore.Storage;

public class StorageException extends Exception{
}


text:
<?php

namespace IXF\Util;

class Set
{
    private $_elts;

    public function __construct($members=array())
    {
        $this->_elts = array();
        foreach ($members as $item)
        $this->_elts[$item] = true;
    }

    public function includes($elt)
    {
        return isset($this->_elts[$elt]);
    }

    public function add($elt)
    {
        $this->_elts[$elt] = true;
    }

    public function discard($elt)
    {
        unset($this->_elts[$elt]);
    }

    // TODO: make Set support foreach
    public function toArray()
    {
        return array_keys($this->_elts);
    }
}


text:
export LC_ALL=C.UTF-8

export HOST=s390x-linux-gnu
# The host arch is unknown, so we run the tests through qemu.
# If the host is s390x and wants to run the tests natively, it can set QEMU_USER_CMD to the empty string.
if [ -z ${QEMU_USER_CMD+x} ]; then export QEMU_USER_CMD="${QEMU_USER_CMD:-"qemu-s390x"}"; fi
export PACKAGES="python3-zmq"
if [ -n "$QEMU_USER_CMD" ]; then
  # Likely cross-compiling, so install the needed gcc and qemu-user
  export DPKG_ADD_ARCH="s390x"
  export PACKAGES="$PACKAGES g++-s390x-linux-gnu qemu-user libc6:s390x libstdc++6:s390x libfontconfig1:s390x libxcb1:s390x"
fi
# Use debian to avoid 404 apt errors
export CONTAINER_NAME=ci_s390x
export DOCKER_NAME_TAG="debian:buster"
export RUN_UNIT_TESTS=true
export RUN_FUNCTIONAL_TESTS=true
export GOAL="install"
export BITCOIN_CONFIG="--enable-reduce-exports --with-incompatible-bdb --with-boost-process"


text:
require "vizier/subject.rb"

describe Vizier::Subject do
  before do
    @subject = Vizier::Subject.new
  end

  it "should allow a field to be required" do
    proc do
      @subject.required_fields([:a_field])
    end.should_not raise_error
  end

  it "should allow multiple fields to be required" do
    proc do
      @subject.required_fields([:field1, :field2, :field3])
    end.should_not raise_error
  end
end

describe Vizier::Subject, "when merging with other subjects" do
  before do
    @subject = Vizier::Subject.new
    @subject.required_fields([:one, :two, :three])
    @subject.one = nil
    @subject.two = nil
    @sub_one = Vizier::Subject.new
    @sub_one.required_fields([:two])
  end

  it "should not allow a merge with colliding fields" do
    @sub_one.two = nil
    proc do
      @subject.merge(nil, @sub_one)
    end.should raise_error(Vizier::CommandError)
  end

  it "should allow a merge with colliding unfulfilled requirements" do
    proc do
      @subject.merge(nil, @sub_one)
    end.should_not raise_error
  end

  it "should allow a colliding merge into a context" do
    @sub_one.two = nil
    proc do
      @subject.merge(:sub, @sub_one)
    end.should_not raise_error
  end

  it "should not allow contexts to collide with fields" do
    proc do
      @subject.merge(:three, @sub_one)
    end.should raise_error(Vizier::CommandError)
  end

  it "should not allow fields to collide with contexts" do
    @subject.merge(:sub, @sub_one)

    proc do
      @subject.required_fields(["sub"])
    end.should raise_error(Vizier::CommandError)
  end
end

describe Vizier::Subject, "with contextual sub-subjects" do
  before do
    @subject = Vizier::Subject.new
    @subject.required_fields([:vegetable, :nut])

    @sub_one = Vizier::Subject.new
    @sub_one.required_fields([:fruit, :nut])
    @sub_one.fruit = 2
    @sub_one.nut = "cashew"

    @subject.merge(:one, @sub_one)

    @sub_two = Vizier::Subject.new
    @sub_two.required_fields([:fruit, :nut])
    @sub_two.fruit = 3
    @sub_two.nut = "chestnut"

    @subject.merge(:two, @sub_one)

    @subject.vegetable = 1
    @subject.protect(:one, :nut)
    @subject.nut = "almond"
  end

  it "should project an image without context" do
    image = @subject.get_image([:vegetable, :nut])
    image.vegetable.should eql(1)
    image.nut.should eql("almond")
  end

  it "should not absorb fields from contexts" do
    proc do
      @subject.get_image([:fruit])
    end.should raise_error(Vizier::CommandError)
  end

  it "should project an image within a context" do
    image = @subject.get_image([:fruit, :vegetable, :nut], [:one])
    image.vegetable.should eql(1)
    image.fruit.should eql(2)
    image.nut.should eql("cashew")
  end

  it "should allow assignment of fields within contexts" do
    @subject.one.nut = "pistacio"
    @subject.get_image([:nut], [:one]).nut.should eql("pistacio")
    @subject.get_image([:nut], [:two]).nut.should eql("almond")
    @subject.get_image([:nut]).nut.should eql("almond")
  end

  it "should allow image assignment back to contexts" do
    image = @subject.get_image([:fruit, :vegetable, :nut], [:one])
    image.nut = "pistacio"
    image.vegetable = 17

    @subject.vegetable.should == 17
    @sub_one.nut.should == "pistacio"
    @subject.nut.should == "almond"
    @sub_two.nut.should == "chestnut"
  end
end

describe Vizier::Subject, "with required fields" do
  before do
    @subject = Vizier::Subject.new
    @subject.required_fields([:a_field])
    @subject.required_fields([:another_field, :yaf])
  end

  it "should not verify if fields are not set" do
    proc do
      @subject.verify
    end.should raise_error(RuntimeError)
  end

  it "should not verify if some fields are not set" do
    @subject.a_field = nil
    @subject.another_field = nil
    proc do
      @subject.verify
    end.should raise_error(RuntimeError)
  end

  it "should verify if all fields are set" do
    @subject.a_field = 1
    @subject.another_field = :something
    @subject.yaf = [:some, "more", 'stuff']
    proc do
      @subject.verify
    end.should_not raise_error
  end

  it "should spawn SubjectImage objects that do respond to readers" do
    @subject.a_field = 1
    @subject.another_field = :something
    @subject.yaf = [:some, "more", 'stuff']

    image = @subject.get_image([:a_field, :another_field, :yaf])

    image.a_field.should eql(1)
    image.another_field.should eql(:something)
    image.yaf.should eql([:some, "more", 'stuff'])
  end
end


text:
layout: post
title: "App Cache Tools and Auto-Generators"
description: ""
article:
  written_on: 2011-10-28
  updated_on: 2011-10-28
authors:
  - sethladd
tags:
  - offline
  - appcache
permalink: /2011/10/App-Cache-Tools-and-Auto-Generators
---
As you know, the Application Cache helps you run your web app offline, significantly speeds up startup time, and reduces bandwidth bills.  All good things to the HTML5 developer!

However, for large apps and sites, it can be cumbersome to add each and every file required for the App Cache.  Luckily, there are now tools that can help you auto-generate your App Cache manifest.

Some options include:

"Manifested":http://dregsoft.com/blog/?p=79
"ManifestR":http://westciv.com/tools/manifestR/
"H5BP build script":https://github.com/jamesgpearce/confess

Give them a shot, let us know in the comments how they work for you.


text:
package org.languagetool.rules.ga;

import org.languagetool.AnalyzedSentence;
import org.languagetool.AnalyzedToken;
import org.languagetool.AnalyzedTokenReadings;
import org.languagetool.rules.*;

import java.util.ArrayList;
import java.util.List;
import java.util.ResourceBundle;

import static org.languagetool.rules.ga.DhaNoBeirtData.getNumberReplacements;
import static org.languagetool.rules.ga.DhaNoBeirtData.getDaoine;

public class DhaNoBeirtRule extends Rule {
  public DhaNoBeirtRule(ResourceBundle messages) {
    super.setCategory(Categories.MISC.getCategory(messages));
    setLocQualityIssueType(ITSIssueType.Misspelling);
    addExamplePair(Example.wrong("Tá <marker>dhá</marker> dheartháireacha agam."),
                   Example.fixed("Tá <marker>beirt</marker> dheartháireacha agam."));
  }

  @Override
  public String getId() {
    return "GA_DHA_NO_BEIRT";
  }

  @Override
  public String getDescription() {
    return "'dhá' nó 'beirt'";
  }

  @Override
  public RuleMatch[] match(AnalyzedSentence sentence) {
    List<RuleMatch> ruleMatches = new ArrayList<>();
    AnalyzedTokenReadings[] tokens = sentence.getTokensWithoutWhitespace();
    int markDeag = 1;
    int prevTokenIndex = 0;
    String replacement = null;
    String msg = null;
    for (int i = 1; i < tokens.length; i++) {  // ignoring token 0, i.e., SENT_START
      if (isNumber(tokens[i]) && (i < tokens.length - 2 && isPerson(tokens[i + 1]))) {
        if ("dhá".equalsIgnoreCase(tokens[i].getToken()) && (i < tokens.length - 2)) {
          for (int j = i + 2; j < tokens.length; j++) {
            if ("déag".equalsIgnoreCase(tokens[j].getToken())) {
              markDeag = j;
              replacement = "dháréag";
              msg = "Ba chóir duit <suggestion>" + replacement + "</suggestion> a scríobh";
              RuleMatch match = new RuleMatch(
                this, sentence, tokens[prevTokenIndex+1].getStartPos(), tokens[prevTokenIndex+1].getEndPos(), msg, "Uimhir phearsanta");
              ruleMatches.add(match);
              msg = "Ba chóir duit \"déag\" a scriosadh.";
              RuleMatch match2 = new RuleMatch(
                this, sentence, tokens[markDeag].getStartPos(), tokens[markDeag].getEndPos(), msg, "Uimhir phearsanta");
              ruleMatches.add(match2);
              msg = null;
            }
          }
        }
        if (replacement == null) {
          replacement = getNumberReplacements().get(tokens[i].getToken());
          if (msg == null) {
            msg = "Ba chóir duit <suggestion>" + replacement + " " + tokens[i + 1].getToken() + "</suggestion> a scríobh";
          }
        }
      }
      if (msg != null) {
        RuleMatch match = new RuleMatch(
          this, sentence, tokens[prevTokenIndex+1].getStartPos(), tokens[prevTokenIndex+1 ].getEndPos(), msg, "Uimhir phearsanta");
        ruleMatches.add(match);
        msg = null;
      }
      prevTokenIndex = i;
    }
    return toRuleMatchArray(ruleMatches);
  }

  private boolean isNumber(AnalyzedTokenReadings tok) {
    for (String num : getNumberReplacements().keySet()) {
      if (num.equalsIgnoreCase(tok.getToken())) {
        return true;
      }
    }
    return false;
  }
  private boolean isPerson(AnalyzedTokenReadings tok) {
    if (getDaoine().contains(tok.getToken().toLowerCase())) {
      return true;
    }
    for (AnalyzedToken reading : tok.getReadings()) {
      if (getDaoine().contains(reading.getLemma())) {
        return true;
      }
    }
    return false;
  }
}


text:
package com.datatorrent.stram;

import com.datatorrent.api.Component;
import com.datatorrent.api.Component.ComponentComplementPair;
import com.datatorrent.api.Context;

/**
 *
ComponentContextPair class.

 *
 * @since 0.3.2
 */
public class ComponentContextPair<COMPONENT extends Component<?>, CONTEXT extends Context> extends ComponentComplementPair<COMPONENT, CONTEXT>
{
  public final CONTEXT context;

  public ComponentContextPair(COMPONENT component, CONTEXT context)
  {
    super(component);
    this.context = context;
  }

  @Override
  public CONTEXT getComplement()
  {
    return context;
  }

  @Override
  public String toString()
  {
    return "ComponentContextPair{component=" + component + ", context=" + context + '}';
  }

}


text:
package com.s4game.server.stage.model.core.element.impl;

import com.s4game.server.stage.model.core.stage.IStage;
import com.s4game.server.stage.model.core.stage.IStageElement;
import com.s4game.server.stage.model.core.stage.Point;

public abstract class AbsElement implements IStageElement {

    private IStage stage;

    private Point position;
    private Point bornPosition;
    private String id;
    private String name;
    private String teamId;
    private String camp;

    public AbsElement(String id, String name, String teamId, String camp) {
        this.id = id;
        this.name = name;
        this.teamId = teamId;
        this.camp = camp;
    }

    public String getId() {
        return this.id;
    }

    public IStage getStage() {
        return this.stage;
    }

    public void setStage(IStage stage) {
        this.stage = stage;
    }

    public void setPosition(int x, int y) {
        this.position = new Point(x, y);
    }

    public Point getPosition() {
        return this.position;
    }

    public Point getBornPosition() {
        return this.bornPosition;
    }

    public void setBornPosition(int x, int y) {
        this.bornPosition = new Point(x, y);
    }

    public String getName() {
        return this.name;
    }

    public String getCamp() {
        return this.camp;
    }

    public String getTeamId() {
        return this.teamId;
    }
}


text:
package it.av.youeat.service;

import it.av.youeat.ocm.model.ActivityEaterRelation;
import it.av.youeat.ocm.model.Eater;

import java.util.List;

/**
 * Operation to manage activities on user relations
 *
 * @author Alessandro Vincelli
 *
 */
public interface ActivityRelationService {

    /**
     * Find relation activities for the given user
     *
     * @param eater
     * @return activities for the given user
     */
    List<ActivityEaterRelation> findByEater(Eater eater);

    /**
     * Find relation activities for the given user
     *
     * @param eater
     * @param firstResult first result
     * @param maxResults number of results
     * @return activities for the given user
     */
    List<ActivityEaterRelation> findByEater(Eater eater, int firstResult, int maxResults);

    /**
     * Find relation activities for the friends given user
     *
     * @param ofUser
     * @return activities for the friends of the given user
     */
    List<ActivityEaterRelation> findByEaterFriend(Eater ofUser);

    /**
     * Find relation activities for the friends given user
     *
     * @param ofUser
     * @param firstResult first result
     * @param maxResults number of results
     * @return activities for the friends of the given user
     */
    List<ActivityEaterRelation> findByEaterFriend(Eater eater, int firstResult, int maxResults);

    /**
     * Find relation activities for the friends given user and for the user himself
     *
     * @param ofUser
     * @param firstResult first result
     * @param maxResults number of results
     * @return activities for the friends of the given user
     */
    List<ActivityEaterRelation> findByEaterFriendAndEater(Eater eater, int firstResult, int maxResults);

    /**
     * Find relation activities for the given users
     *
     * @param eaters
     * @param firstResult first result
     * @param maxResults number of results
     * @return activities for the given user
     */
    List<ActivityEaterRelation> findByEaters(List<Eater> eaters, int firstResult, int maxResults);

    /**
     * store the activity in the database
     *
     * @param activityRistorante the activity to save
     * @return the just saved activity
     */
    ActivityEaterRelation save(ActivityEaterRelation activityRistorante);

    /**
     * remove all the activities related to the given users
     *
     * @param eater
     */
    void removeByEater(Eater eater);
}


text:
package com.intellij.openapi.fileEditor.impl;

import com.intellij.openapi.extensions.ExtensionPointName;
import com.intellij.openapi.vfs.VirtualFile;
import org.jetbrains.annotations.NotNull;

public interface NonProjectFileWritingAccessExtension {
  ExtensionPointName<NonProjectFileWritingAccessExtension> EP_NAME =
    ExtensionPointName.create("com.intellij.nonProjectFileWritingAccessExtension");

  /**
   * @return true if the file should not be protected from accidental writing. false to use default logic.
   */
  default boolean isWritable(@NotNull VirtualFile file) {
    return false;
  }

  /**
   * @return true if the file should be protected from accidental writing. false to use default logic.
   */
  default boolean isNotWritable(@NotNull VirtualFile file) {
    return false;
  }
}


text:
(function($) {
  $.fn.ytchromeless = function(options){

    //Initial configuration
    var config = {
      videoWidth  : '640',
      videoHeight : '360',
      videoIdBase : 'ytplayer',
      params : {
        allowScriptAccess: 'always',
        wmode: 'transparent'
      }
    };

    return this.each(function(i) {

      // initial var setup

        var o    = $.extend(config, options),

            // set jQuery objects
            $link      = $(this),

            // set variables
            url        = $link.attr('href'),
            videoId    = $link.attr('id') || o.videoIdBase + i,
            ytVideoId  = url.substr(31),

            // new DOM elements
            $video     = $link.wrap( '

' ).parent(),
            $controls  = $('

' ).appendTo( $video ),
            $toReplace = $('

').prependTo( $video ).attr('id', videoId),
            $bar,
            $indicator,
            $loaded,
            $mute,
            $play,
            $seek,

            // set up the special player object
            player;

        // bind public methods upfront
        $video.bind({

          // playing, pausing, muting,
          'togglePlay' : function(){ $video.togglePlay(); },
          'play'       : function(){ $video.play(); },
          'pause'      : function(){ $video.pause(); },
          'toggleMute' : function(){ $video.toggleMute(); },
          'mute'       : function(){ $video.mute(); },
          'unMute'     : function(){ $video.unMute(); },
          'seek'       : function(){ $video.seek(); },

          // initializing and revising the player
          'update'     : function(){ $video.update(); },
          'cue'        : function(){ player.cueVideoById( ytVideoId ); }

        });

      // control methods

        // function fired when the play/pause button is hit
        $video.togglePlay = function() {
          if( $play.hasClass('playing') ) {
            $video.trigger('pause');
          } else {
            $video.trigger('play');
          }
          return false;
        };

        // play the video
        $video.play = function() {
          player.playVideo();
          $play.removeClass('paused').addClass('playing').attr('title','Pause');
        };

        // pause
        $video.pause = function() {
          player.pauseVideo();
          $play.removeClass('playing').addClass('paused').attr('title','Play');
        };

        // function fired when the mute/unmute button is hit
        $video.toggleMute = function() {
          if( $mute.hasClass('muted') ) {
            $video.trigger('unMute');
          } else {
            $video.trigger('mute');
          }
          return false;
        };

        // mute the video
        $video.mute = function() {
          player.mute();
          $mute.addClass('muted').attr('title','Un-Mute');
        };

        // unmute
        $video.unMute = function() {
          player.unMute();
          $mute.removeClass('muted').attr('title','Mute');
        };

        //Seek to a position in the video
        $video.seek = function(seekPosition) {
          var seekToPosition = Math.round(player.getDuration() * seekPosition);
          player.seekTo(seekToPosition, false);
        };

      // player init and update methods

        //Update the video status
        $video.update = function() {

          if( player && player.getDuration ) {

            if( player.getPlayerState() === 1 ) {
              $video.play();
            } else if ( player.getPlayerState() === 0 ) {
              $video.pause();
            }

            if( player.getVideoLoadedFraction() > 0) {
              var loadedAmount = Math.round(Math.min(player.getVideoLoadedFraction(), 1) * 100);
              $loaded.css( 'width', loadedAmount + '%' );
            }

            if( player.getCurrentTime() > 0 ) {
              var videoPosition = ( player.getCurrentTime() / player.getDuration() ) * 100;
              $indicator.css( 'left', videoPosition + '%' );
            }

          }

        };

        // the youtube movie calls this method when it loads
        // DO NOT CHANGE THIS METHOD'S NAME
        onYouTubePlayerReady = function( videoId ) {

          var $videoRef = $( document.getElementById( videoId ) ).parent();

          setInterval(function(){
            $videoRef.trigger('update');
          }, 250);

          $videoRef.trigger('cue');

        };

      // init methods

        // the embed!
        $video.init = function() {

          swfobject.embedSWF(
            'http://www.youtube.com/apiplayer?version=3&enablejsapi=1&playerapiid=' + videoId,
            videoId,
            o.videoWidth,
            o.videoHeight,
            '8',
            null,
            null,
            o.params,
            { id: videoId },
            function(){
              player = document.getElementById( videoId );
            }
          );

          $video.addControls();

        };

        // add controls
        $video.addControls = function() {

          //Play and pause button
          $play = $('', {
                    href: '#',
                    'class': 'play-pause',
                    text: 'Play/Pause',
                    title: 'Play',
                    click: function() {
                      $video.trigger('togglePlay');
                      return false;
                    }
                  }).appendTo( $controls );

          //Play and pause button
          $mute = $('', {
                    href: '#',
                    'class': 'volume',
                    text: 'Volume',
                    title: 'Mute',
                    click: function() {
                      $video.trigger('toggleMute');
                      return false;
                    }
                  }).appendTo( $controls );

          //View on YouTube
          $link
            .addClass('view-youtube')
            .attr('title', 'View on YouTube')
            .html('Play/View on YouTube')
            .appendTo( $controls );

          //Play and pause button
          $seek = $('
', {
                    'class': 'status',
                    click: function(e) {
                      var skipTo      = e.pageX - $seek.offset().left,
                          statusWidth = $seek.width();
                      $video.seek( skipTo / statusWidth );
                    }
                  }).appendTo( $controls );

          $bar       = $('

').appendTo($seek);
          $loaded    = $('

').appendTo($bar);
          $indicator = $('').appendTo($bar);

        };

        $video.init();

    });

  };

})(jQuery);


text:
<?php

include("class_scraper.php");

if(!isset($_GET['q'])){
  exit(0);
}

session_start();

$school_road = new scraper($_GET['q']);
$school_road->fetch_auctions();
$school_road->scrape_all_auctions();
$school_road->save_session();
$school_road->dump_auctions();

session_write_close();

 ?>


text:
package io.pivotal.fe.demos.citiesui.model;

public class FormInput {
	private String name;
	private Long size = 20L;

	public FormInput() {

	}
	public FormInput(String name, Long size) {
		this.name = name;
		this.size = size;
	}
	public String getName() {
		return name;
	}
	public void setName(String name) {
		this.name = name;
	}
	public Long getSize() {
		return size;
	}
	public void setSize(Long size) {
		this.size = size;
	}

}


text:
package pl.finder.elmer.serialization;

import com.google.gson.Gson;
import com.google.gson.GsonBuilder;

import lombok.AccessLevel;
import lombok.NoArgsConstructor;

@NoArgsConstructor(access = AccessLevel.PRIVATE)
final class DefaultGson {

    static Gson create() {
        final Gson gson = new GsonBuilder()
                .setDateFormat("yyyy-MM-dd'T'HH:mm:ss.SSSXXX")
                .create();
        return gson;
    }
}


text:
Ander Hidalgo - Jon Ander Vicente

Para crear el Modelo MER relacional hemos utilizado el Oracle Datamodeler Versión 4.1.1.888.

Para crear el Modelo Lógico hemos utilizado el Oracle SQLDeveloper Versión 4.1.1.19.

Para crear el Diagrama de clases hemos utilizado el Visual Paradigm 13.0. Para poder además visualizar el diagrama en cualquier ordernador lo hemos exportado en formato jpeg.

Para crear la aplicación de java hemos utilizado el IDE Netbeans 8.1.

La Base de Datos utilizada es Oracle Database Enterprise Edition Release 12.1.0.2.0.

//
//
// Despliegue de la aplicación
//
//

En SQLDeveloper abriremos y ejecutaremos el modelo lógico para crear primero la base de datos de oracle sobre la que trabajaremos.

Para hacer funcionar la aplicación importaremos el proyecto en Netbeans, al importar nos encontrara varios problemas ya que no encuentra las librerias utilizadas.

Debemos de acordarnos también de modificar el codigo del GenericoBD en el modelo segun nuestras necesidades, ya que en este caso tenemos la conexion con nuestra base de datos que es daw10.

Los campos que se modifican en el codigo son estos dos:
  String user = "daw10";
  String pass = "daw10";

Las librerias las adjuntamos en el repositorio y son Datechooser.jar (para el datachoosercombo de fecha) y ojdbc6.jar (para la conexión de java con la base de datos oracle).

Una vez resuelto los problemas con las librerías ejecutaremos el programa en Netbeans. Esto nos llevara a la primera ventana que es Login. Tendremos que introducir un usuario y contraseña metido en la base de datos. Estos usuarios seran de administración o logística.

Para hacer una prueba hay creado un usuario admin con contraseña admin. Este usuario de prueba es de tipo administración, asi que es capaz de acceder a todo.

También se puede entrar por ejemplo con el usuario aritzGM con contraseña 1A que sería de administración o con amaiaGP con contraseña 2B que es de logística.

Para cualquier otra duda el enlace a nuestro blog semanal esta en el repositorio.


text:
package org.apache.felix.ipojo.runtime.core.test.components.temporal;

import org.apache.felix.ipojo.annotations.Component;
import org.apache.felix.ipojo.handler.temporal.Requires;
import org.apache.felix.ipojo.runtime.core.test.services.FooService;

@Component
public class TemporalWithDI {

    @Requires(onTimeout="org.apache.felix.ipojo.runtime.core.test.components.ProvidesSimple")
    private FooService fs;

}


text:
require 'rails_helper'

RSpec.describe "Materials", type: :request do
  describe "GET /materials" do
    it "works! (now write some real specs)" do
      get materials_path
      expect(response).to have_http_status(200)
    end
  end
end