modified carbon decomp coeff by multiplying cell thickness

amanzi · Jan 8, 2025 · 7daf3d5 · 7daf3d5
1 parent e5becc0
commit 7daf3d5
Show file tree

Hide file tree

Showing 2 changed files with 20 additions and 11 deletions.
diff --git a/src/constitutive_relations/eos/carbon_decomposition_rate_evaluator.cc b/src/constitutive_relations/eos/carbon_decomposition_rate_evaluator.cc
@@ -23,23 +23,21 @@ CarbonDecomposeRateEvaluator::CarbonDecomposeRateEvaluator(Teuchos::ParameterLis
   : EvaluatorSecondaryMonotypeCV(plist)
 {
   Tag tag = my_keys_.front().second;
-  Key domain = Keys::getDomain(my_keys_.front().first); //column domain
+  domain_ = Keys::getDomain(my_keys_.front().first); // column, domain
+  domain_surf_ = Keys::readDomainHint(plist, domain_, "subsurface", "surface"); 
 
-  temp_key_ = Keys::readKey(plist, domain, "temperature", "temperature");
+  temp_key_ = Keys::readKey(plist, domain_, "temperature", "temperature");
   dependencies_.insert(KeyTag{ temp_key_, tag });
 
-  pres_key_ = Keys::readKey(plist, domain, "pressure", "pressure");
+  pres_key_ = Keys::readKey(plist, domain_, "pressure", "pressure");
   dependencies_.insert(KeyTag{ pres_key_, tag });
 
-  por_key_ = Keys::readKey(plist, domain, "porosity", "porosity");
+  por_key_ = Keys::readKey(plist, domain_, "porosity", "porosity");
   dependencies_.insert(KeyTag{ por_key_, tag });
 
-  depth_key_ = Keys::readKey(plist, domain, "depth", "depth");
+  depth_key_ = Keys::readKey(plist, domain_, "depth", "depth");
   dependencies_.insert(KeyTag{ depth_key_, tag });
 
-  cv_key_ =  Keys::readKey(plist, domain, "cell volume","cell_volume");
-  dependencies_.insert(KeyTag{ cv_key_, tag });
-
   q10_ = plist_.get<double>("Q10 [-]", 2.0);
 }
 
@@ -56,19 +54,29 @@ CarbonDecomposeRateEvaluator::Evaluate_(const State& S, const std::vector<Compos
 {
   Tag tag = my_keys_.front().second;
   Epetra_MultiVector& res_c = *result[0]->ViewComponent("cell", false);
+  Epetra_MultiVector& dz_c = *result[0]->ViewComponent("cell", false);
 
   const auto& temp_c = *S.Get<CompositeVector>(temp_key_, tag).ViewComponent("cell", false);
   const auto& pres_c = *S.Get<CompositeVector>(pres_key_, tag).ViewComponent("cell", false);
   const auto& por_c = *S.Get<CompositeVector>(por_key_, tag).ViewComponent("cell", false);
   const auto& depth_c = *S.Get<CompositeVector>(depth_key_, tag).ViewComponent("cell", false);
-  const auto& cv_c = *S.Get<CompositeVector>(cv_key_, tag).ViewComponent("cell", false);
+
+  const auto& mesh = *S.GetMesh(domain_);
+  const auto& mesh_surf = *S.GetMesh(domain_surf_);
+
+  for (int col = 0; col != mesh.columns.num_columns_owned; ++col) {
+    const auto& col_cells = mesh.columns.getCells(col);
+    for (int i = 0; i != col_cells.size(); ++i) {
+      dz_c[0][col_cells[i]] = mesh.getCellVolume(col_cells[i]) / mesh_surf.getCellVolume(col);
+    }
+  }
 
   for (int c = 0; c != res_c.MyLength(); ++c) {
     if (temp_c[0][c] >= 273.15) {
       double f_temp = Func_Temp(temp_c[0][c], q10_);
       double f_depth = Func_Depth(depth_c[0][c]);
       double f_pres_temp = Func_TempPres(temp_c[0][c], pres_c[0][c]);
-      res_c[0][c] = (f_temp * f_depth * f_pres_temp) * (1 - por_c[0][c]) * cv_c[0][c];
+      res_c[0][c] = (f_temp * f_depth * f_pres_temp * dz_c[0][c]) * (1 - por_c[0][c]);
     } else {
       res_c[0][c] = 0.;
     }

diff --git a/src/constitutive_relations/eos/carbon_decomposition_rate_evaluator.hh b/src/constitutive_relations/eos/carbon_decomposition_rate_evaluator.hh
@@ -53,7 +53,8 @@ class CarbonDecomposeRateEvaluator : public EvaluatorSecondaryMonotypeCV {
   Key sat_key_;
   Key por_key_;
   Key depth_key_;
-  Key cv_key_;
+  Key domain_;
+  Key domain_surf_;
 
   double q10_;