updates and bug fixes -- richards now working for steadystate and tra…

…nsient fv

src/pks/flow/constitutive_relations/wrm/relative_permeability_model.hh

@@ -42,30 +42,43 @@ class RelativePermeabilityModel {
     my_key_ = Keys::cleanPListName(plist);
     auto domain = Keys::getDomain(my_key_);
     s_key_ = Keys::readKey(plist, domain, "saturation", "saturation_liquid");
+    dens_key_ = Keys::readKey(plist, domain, "density", "molar_density_liquid");
+    visc_key_ = Keys::readKey(plist, domain, "viscosity", "viscosity");
+
+    rescaling_ = 1.0 / plist.get<double>("permeability rescaling", 1.0);
   }
 
   void setViews(const std::vector<cView_type>& deps,
                 const std::vector<View_type>& res,
                 const State& s) {
     res_ = res[0];
     s_ = deps[0];
+    dens_ = deps[1];
+    visc_ = deps[2];
   }
 
   KeyVector getMyKeys() const { return { my_key_ }; }
-  KeyVector getDependencies() const { return { s_key_ }; }
+  KeyVector getDependencies() const { return { s_key_, dens_key_, visc_key_ }; }
 
   KOKKOS_INLINE_FUNCTION void operator()(const int i) const {
-    res_(i,0) = model_.k_relative(s_(i,0));
+    res_(i,0) = rescaling_ * model_.k_relative(s_(i,0)) * dens_(i,0) / visc_(i,0);
   }
 
   KOKKOS_INLINE_FUNCTION void operator()(Deriv<0>, const int i) const {
-    res_(i,0) = model_.d_k_relative(s_(i,0));
+    res_(i,0) = rescaling_ * model_.d_k_relative(s_(i,0)) * dens_(i,0) / visc_(i,0);
+  }
+  KOKKOS_INLINE_FUNCTION void operator()(Deriv<1>, const int i) const {
+    res_(i,0) = rescaling_ * model_.k_relative(s_(i,0)) / visc_(i,0);
+  }
+  KOKKOS_INLINE_FUNCTION void operator()(Deriv<2>, const int i) const {
+    res_(i,0) = - rescaling_ * dens_(i,0) * model_.k_relative(s_(i,0)) / pow(visc_(i,0), 2);
   }
 
  private:
   View_type res_;
-  cView_type s_;
-  Key my_key_, s_key_;
+  cView_type s_, dens_, visc_;
+  Key my_key_, s_key_, dens_key_, visc_key_;
+  double rescaling_;
 
   WRM_type model_;
 };

src/pks/flow/richards.hh

@@ -383,7 +383,7 @@ class Richards : public PK_PhysicalBDF_Default {
   Teuchos::RCP<Operators::PDE_DiffusionWithGravity> preconditioner_diff_;
   Teuchos::RCP<Operators::PDE_DiffusionWithGravity> face_matrix_diff_;
   Teuchos::RCP<Operators::PDE_Accumulation> preconditioner_acc_;
-  double perm_scale_;
+  // double perm_scale_;
 
   // flag to do jacobian and therefore coef derivs
   bool precon_used_;

src/pks/flow/richards_physics.cc

@@ -44,12 +44,12 @@ Richards::ApplyDiffusion_(const Tag& tag, const Teuchos::Ptr<CompositeVector>& g
 
   Teuchos::RCP<const CompositeVector> pres = S_->GetPtrW<CompositeVector>(key_, tag, name_);
   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, name_);
   matrix_diff_->UpdateFlux(pres.ptr(), flux.ptr());
   changedEvaluatorPrimary(flux_key_, tag, *S_);
-  matrix_diff_->ApplyBCs(true, true, true);
 
   // calculate the residual
   matrix_->ComputeNegativeResidual(*pres, *g);

src/pks/flow/richards_pk.cc

@@ -49,7 +49,7 @@ Richards::Richards(const Comm_ptr_type& comm,
     modify_predictor_first_bc_flux_(false),
     upwind_from_prev_flux_(false),
     clobber_boundary_flux_dir_(false),
-    perm_scale_(1.),
+    // perm_scale_(1.),
     jacobian_(false),
     jacobian_lag_(0),
     iter_(0),
@@ -93,8 +93,8 @@ Richards::ParseParameterList_()
     deform_key_ = Keys::readKey(*plist_, domain_, "deformation indicator", "base_porosity");
 
   // scaling for permeability for better "nondimensionalization"
-  perm_scale_ = plist_->get<double>("permeability rescaling", 1.e7);
-  S_->GetEvaluatorList(coef_key_).set<double>("permeability rescaling", perm_scale_);
+  // perm_scale_ = plist_->get<double>("permeability rescaling", 1.e7);
+  // S_->GetEvaluatorList(coef_key_).set<double>("permeability rescaling", perm_scale_);
 }
 
 // -------------------------------------------------------------
@@ -502,6 +502,8 @@ Richards::Initialize()
   const AmanziGeometry::Point& g = S_->Get<AmanziGeometry::Point>("gravity", Tags::DEFAULT);
   matrix_diff_->SetGravity(g);
   matrix_diff_->SetBCs(bc_, bc_);
+
+  S_->GetEvaluator(perm_key_, Tags::DEFAULT).Update(*S_, name_);
   auto K = S_->GetPtr<TensorVector>(perm_key_, Tags::DEFAULT);
   matrix_diff_->SetTensorCoefficient(K);
 
@@ -758,13 +760,14 @@ Richards::UpdatePermeabilityData_(const Tag& tag)
 
       if (!deform_key_.empty() &&
           S_->GetEvaluator(deform_key_, tag_next_).Update(*S_, name_ + " flux dir")) {
+        S_->GetEvaluator(perm_key_, Tags::DEFAULT).Update(*S_, name_ + " flux_dir");
         auto K = S_->GetPtr<TensorVector>(perm_key_, Tags::DEFAULT);
         face_matrix_diff_->SetTensorCoefficient(K);
       }
       face_matrix_diff_->SetDensity(rho);
       face_matrix_diff_->UpdateMatrices(Teuchos::null, pres.ptr());
-      face_matrix_diff_->UpdateFlux(pres.ptr(), flux_dir.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);

src/pks/flow/richards_ti.cc

@@ -79,9 +79,9 @@ Richards::FunctionalResidual(double t_old,
   vnames.emplace_back("poro");
   vecs.emplace_back(
     S_->GetPtr<CompositeVector>(Keys::getKey(domain_, "porosity"), tag_next_).ptr());
-  vnames.emplace_back("perm_K");
-  vecs.emplace_back(
-    S_->GetPtr<CompositeVector>(Keys::getKey(domain_, "permeability"), tag_next_).ptr());
+  // vnames.emplace_back("perm_K");
+  // vecs.emplace_back(
+  //   S_->GetPtr<CompositeVector>(Keys::getKey(domain_, "permeability"), tag_next_).ptr());
   vnames.emplace_back("k_rel");
   vecs.emplace_back(S_->GetPtr<CompositeVector>(coef_key_, tag_next_).ptr());
   vnames.emplace_back("wind");
@@ -171,6 +171,7 @@ Richards::UpdatePreconditioner(double t, Teuchos::RCP<const TreeVector> up, doub
     } else {
       dkrdp = S_->GetDerivativePtr<CompositeVector>(coef_key_, tag_next_, key_, tag_next_);
     }
+    dkrdp->print(std::cout);
   }
 
   // -- primary term
@@ -181,14 +182,14 @@ Richards::UpdatePreconditioner(double t, Teuchos::RCP<const TreeVector> up, doub
   // -- local matries, primary term
   preconditioner_->Zero();
   preconditioner_diff_->UpdateMatrices(Teuchos::null, up->getData().ptr());
+  preconditioner_diff_->ApplyBCs(true, true, true);
 
   // -- local matries, Jacobian term
   if (jacobian_ && iter_ >= jacobian_lag_) {
     Teuchos::RCP<CompositeVector> flux = S_->GetPtrW<CompositeVector>(flux_key_, tag_next_, name_);
     preconditioner_diff_->UpdateFlux(up->getData().ptr(), flux.ptr());
     preconditioner_diff_->UpdateMatricesNewtonCorrection(flux.ptr(), up->getData().ptr());
   }
-  preconditioner_diff_->ApplyBCs(true, true, true);
 
   // Update the preconditioner with accumulation terms.
   // -- update the accumulation derivatives

src/pks/pk_physical_default.cc

@@ -144,11 +144,16 @@ PK_Physical_Default::Initialize()
 
     // -- Calculate the IC.
     Teuchos::ParameterList ic_plist = plist_->sublist("initial condition");
+    ic_plist.setName(key_);
     record.Initialize(ic_plist);
 
+    // debug
+    db_->WriteVector(key_ + " IC", record.GetPtr<CompositeVector>().ptr(), true);
+
     // communicate just to make sure values are initialized for valgrind's sake
     record.Get<CompositeVector>().scatterMasterToGhosted();
     ChangedSolutionPK(tag_next_);
+
   }
 };