Recent changes for tile implementation (#236)

* Improved model for disctributed tiles (linear + nonlinear terms)
* Defined mass density for distributed tiles

---------

Co-authored-by: saubhagya-gatech <saubhagya.rathore@outlook.com>
Co-authored-by: Ethan Coon <coonet@ornl.gov>

src/pks/flow/constitutive_relations/sources/distributed_tiles_evaluator.cc

@@ -28,8 +28,10 @@ DistributedTilesRateEvaluator::DistributedTilesRateEvaluator(Teuchos::ParameterL
 
   // there is no real good way to guess, so we just force the user to provide
   // if they want to override the defaults.
+
   acc_sources_key_ = Keys::readKey(plist, domain_, "accumulated source", "accumulated_source");
   dist_sources_key_ = Keys::readKey(plist, domain_, "distributed source", "distributed_source");
+
   if ((akey != acc_sources_key_) && (akey != dist_sources_key_)) {
     Errors::Message msg;
     msg << "DistributedTilesRateEvaluator: key requested \"" << akey
@@ -52,13 +54,22 @@ DistributedTilesRateEvaluator::DistributedTilesRateEvaluator(Teuchos::ParameterL
   }
 
   pres_key_ = Keys::readKey(plist, domain_, "pressure", "pressure");
+  mol_dens_key_ = Keys::readKey(plist, domain_, "molar density", "molar_density_liquid");
+  mass_dens_key_ = Keys::readKey(plist, domain_, "mass density", "mass_density_liquid");
+  visc_key_ = Keys::readKey(plist, domain_, "viscosity liquid", "viscosity_liquid");
   dependencies_.insert(KeyTag{ pres_key_, tag });
-  mol_dens_key_ = Keys::readKey(plist, domain_, "molar density liquid", "molar_density_liquid");
   dependencies_.insert(KeyTag{ mol_dens_key_, tag });
+  dependencies_.insert(KeyTag{ visc_key_, tag });
+  dependencies_.insert(KeyTag{ mass_dens_key_, tag });
+
 
   // other parameters
   num_ditches_ = plist.get<int>("number of ditches");
-  k_ = plist.get<double>("tile permeability [m^2]");
+  ka_ = plist.get<double>("permeability above [m^2]");
+  kb_ = plist.get<double>("permeability below [m^2]");
+  d_ = plist.get<double>("equivalent distance to impermeable boundary [m]");
+  L_ = plist.get<double>("drain spacing [m]");
+  th_ = plist.get<double>("drain layer thickness [m]", -1.0);
   num_components_ = plist.get<int>("number of components", 1);
   p_enter_ = plist.get<double>("entering pressure [Pa]", 101325);
 }
@@ -67,11 +78,22 @@ void
 DistributedTilesRateEvaluator::Update_(State& S)
 {
   auto tag = my_keys_.front().second;
+  Key akey = my_keys_.front().first;
+  Key domain_name = Keys::getDomain(akey);
+  auto mesh = S.GetMesh(domain_name);
 
   double dt = S.Get<double>("dt", tag);
 
+
+  int z_index = mesh->getSpaceDimension() - 1;
+  const auto& gravity = S.Get<AmanziGeometry::Point>("gravity", Tags::DEFAULT);
+  double gr = -gravity[z_index];
+
   const auto& pres = *S.Get<CompositeVector>(pres_key_, tag).ViewComponent("cell", false);
   const auto& dens = *S.Get<CompositeVector>(mol_dens_key_, tag).ViewComponent("cell", false);
+  const auto& mass_dens = *S.Get<CompositeVector>(mass_dens_key_, tag).ViewComponent("cell", false);
+  const auto& visc = *S.Get<CompositeVector>(visc_key_, tag).ViewComponent("cell", false);
+
   const auto& cv =
     *S.Get<CompositeVector>(Keys::getKey(domain_, "cell_volume"), tag).ViewComponent("cell", false);
   const auto& sub_marks = *S.Get<CompositeVector>(catch_id_key_, tag).ViewComponent("cell", false);
@@ -102,7 +124,26 @@ DistributedTilesRateEvaluator::Update_(State& S)
 
     for (AmanziMesh::Entity_ID c = 0; c != ncells; ++c) {
       if (sub_marks[0][c] > 0) {
-        double val = std::min(p_enter_ - pres[0][c], 0.0) * dens[0][c] * k_;
+        if (th_ < 0) {
+          const auto& [faces, dirs] = mesh->getCellFacesAndDirections(c);
+          double zmax = -1e+98;
+          double zmin = 1e+98;
+          for (auto f : faces) {
+            auto xf = mesh->getFaceCentroid(f);
+            zmax = std::max(zmax, xf[z_index]);
+            zmin = std::min(zmin, xf[z_index]);
+          }
+          th_ = zmax - zmin;
+          AMANZI_ASSERT(th_ > 1e-12);
+        }
+
+
+        double diff_p = std::min(p_enter_ - pres[0][c], 0.0); //represents h_average difference
+        double val = diff_p * dens[0][c] * (12. * kb_ * d_) /
+                     (L_ * L_ * visc[0][c] * th_); //linear term [mol/(m^3 s)]
+        val -= 9. * diff_p * diff_p * dens[0][c] *
+               (ka_ / (visc[0][c] * mass_dens[0][c] * L_ * L_ * gr *
+                       th_)); // nonlinear term [mol/(m^3 s)]
 
         if (!factor_key_.empty()) {
           for (int i = 0; i < num_components_; ++i) {

src/pks/flow/constitutive_relations/sources/distributed_tiles_evaluator.hh

@@ -71,11 +71,11 @@ class DistributedTilesRateEvaluator : public EvaluatorSecondary {
   // my dependencies
   Key catch_id_key_;
   Key pres_key_;
-  Key mol_dens_key_;
+  Key mol_dens_key_, mass_dens_key_, visc_key_;
   Key factor_key_;
 
   double p_enter_;
-  double k_;
+  double ka_, kb_, d_, L_, th_;
   int num_ditches_;
   int num_components_;
 

src/pks/flow/constitutive_relations/sources/surface_distributed_tiles_evaluator.cc

@@ -46,8 +46,7 @@ SurfDistributedTilesRateEvaluator::SurfDistributedTilesRateEvaluator(Teuchos::Pa
 
   auto domain_subsurf = Keys::readDomainHint(plist, domain_, "surface", "domain");
   acc_sources_key_ =
-    Keys::readKey(plist, domain_subsurf, "accumulated source", "accumulated_source");
-  dependencies_.insert(KeyTag{ acc_sources_key_, tag });
+    Keys::readKey(plist, domain_subsurf, "accumulated source", "accumulated_tile_sources");
 
   cv_key_ = Keys::readKey(plist, domain_, "cell volume", "cell_volume");
   dependencies_.insert(KeyTag{ cv_key_, tag });