From 473495d72dd17017387ae0de3fe5d7377f45047c Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Mon, 21 Nov 2022 19:57:31 -0700 Subject: [PATCH 01/58] Adding evaluator for pipe drainage network Work in progress --- .../flow/source_terms/pipe_drain_evaluator.cc | 131 ++++++++++++++++++ .../flow/source_terms/pipe_drain_evaluator.hh | 59 ++++++++ .../source_terms/pipe_drain_evaluator_reg.hh | 6 + 3 files changed, 196 insertions(+) create mode 100644 src/pks/flow/source_terms/pipe_drain_evaluator.cc create mode 100644 src/pks/flow/source_terms/pipe_drain_evaluator.hh create mode 100644 src/pks/flow/source_terms/pipe_drain_evaluator_reg.hh diff --git a/src/pks/flow/source_terms/pipe_drain_evaluator.cc b/src/pks/flow/source_terms/pipe_drain_evaluator.cc new file mode 100644 index 000000000..830fd3c15 --- /dev/null +++ b/src/pks/flow/source_terms/pipe_drain_evaluator.cc @@ -0,0 +1,131 @@ +/* -*- mode: c++; indent-tabs-mode: nil -*- */ + +/* + Evaluator for water drainage through a pipe network. + This depends on: + 1) flow depth above surface elevation + 2) hydraulic head in the pipe network + + Authors: Giacomo Capodaglio (gcapodaglio@lanl.gov) +*/ + +#include "pipe_drain_evaluator.hh" + +namespace Amanzi { +namespace Flow { + + +PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : + EvaluatorSecondaryMonotypeCV(plist) +{ + M_ = plist_.get("molar mass", 0.0180153); + bar_ = plist_.get("allow negative water content", false); + rollover_ = plist_.get("water content rollover", 0.); + + manhole_radius_ = plist_.get("manhole radius", 0.24); + energy_losses_coeff_ = plist.get("energy losses coeff", 0.1); + H_max_ = plist.get("pipe max height", 0.1); + gravity_ = plist.get("gravity", 9.81);//TODO this needs to be converted to an array + H_ = plist.get("bed flume height", 0.478); + + Key domain_name = Keys::getDomain(my_keys_.front().first); + Tag tag = my_keys_.front().second; + + // TODO + // my dependencies + //pres_key_ = Keys::readKey(plist_, domain_name, "pressure", "pressure"); + //dependencies_.insert(KeyTag{pres_key_, tag}); + //cv_key_ = Keys::readKey(plist_, domain_name, "cell volume", "cell_volume"); + //dependencies_.insert(KeyTag{cv_key_, tag}); +} + + +Teuchos::RCP +PipeDrainEvaluator::Clone() const { + return Teuchos::rcp(new PipeDrainEvaluator(*this)); +} + + +void PipeDrainEvaluator::Evaluate_(const State& S, + const std::vector& result) +{ + Tag tag = my_keys_.front().second; + Epetra_MultiVector& res = *result[0]->ViewComponent("cell",false); + const Epetra_MultiVector& pres = *S.GetPtr(pres_key_, tag) + ->ViewComponent("cell",false); + + const Epetra_MultiVector& cv = *S.GetPtr(cv_key_, tag) + ->ViewComponent("cell",false); + + const double& p_atm = S.Get("atmospheric_pressure", Tags::DEFAULT); + const auto& gravity = S.Get("gravity", Tags::DEFAULT); + double gz = -gravity[2]; // check this + + int ncells = res.MyLength(); + if (bar_) { + for (int c=0; c!=ncells; ++c) { + res[0][c] = cv[0][c] * (pres[0][c] - p_atm) / (gz * M_); + } + } else if (rollover_ > 0.) { + for (int c=0; c!=ncells; ++c) { + double dp = pres[0][c] - p_atm; + double dp_eff = dp < 0. ? 0. : + dp < rollover_ ? + dp*dp/(2*rollover_) : + dp - rollover_/2.; + res[0][c] = cv[0][c] * dp_eff / (gz * M_); + } + } else { + for (int c=0; c!=ncells; ++c) { + res[0][c] = pres[0][c] < p_atm ? 0. : + cv[0][c] * (pres[0][c] - p_atm) / (gz * M_); + } + } +} + + +void OverlandPressureWaterContentEvaluator::EvaluatePartialDerivative_(const State& S, + const Key& wrt_key, const Tag& wrt_tag, + const std::vector& result) +{ + Tag tag = my_keys_.front().second; + AMANZI_ASSERT(wrt_key == pres_key_); + + Epetra_MultiVector& res = *result[0]->ViewComponent("cell",false); + const Epetra_MultiVector& pres = *S.GetPtr(pres_key_, tag) + ->ViewComponent("cell",false); + + const Epetra_MultiVector& cv = *S.GetPtr(cv_key_, tag) + ->ViewComponent("cell",false); + + const double& p_atm = S.Get("atmospheric_pressure", Tags::DEFAULT); + const auto& gravity = S.Get("gravity", Tags::DEFAULT); + double gz = -gravity[2]; // check this + + if (wrt_key == pres_key_) { + int ncells = res.MyLength(); + if (bar_) { + for (int c=0; c!=ncells; ++c) { + res[0][c] = cv[0][c] / (gz * M_); + } + } else if (rollover_ > 0.) { + for (int c=0; c!=ncells; ++c) { + double dp = pres[0][c] - p_atm; + double ddp_eff = dp < 0. ? 0. : + dp < rollover_ ? dp/rollover_ : 1.; + res[0][c] = cv[0][c] * ddp_eff / (gz * M_); + } + } else { + for (int c=0; c!=ncells; ++c) { + res[0][c] = pres[0][c] < p_atm ? 0. : + cv[0][c] / (gz * M_); + } + } + } else { + res.PutScalar(0.); + } +} + + +} //namespace +} //namespace diff --git a/src/pks/flow/source_terms/pipe_drain_evaluator.hh b/src/pks/flow/source_terms/pipe_drain_evaluator.hh new file mode 100644 index 000000000..ff630360b --- /dev/null +++ b/src/pks/flow/source_terms/pipe_drain_evaluator.hh @@ -0,0 +1,59 @@ +/* -*- mode: c++; indent-tabs-mode: nil -*- */ + +/* + Evaluator for water drainage through a pipe network. + This depends on: + 1) flow depth above surface elevation + 2) hydraulic head in the pipe network + + Authors: Giacomo Capodaglio (gcapodaglio@lanl.gov) +*/ + +#ifndef AMANZI_FLOW_RELATIONS_PIPE_DRAIN_EVALUATOR_ +#define AMANZI_FLOW_RELATIONS_PIPE_DRAIN_EVALUATOR_ + +#include "EvaluatorSecondaryMonotype.hh" +#include "Factory.hh" + +namespace Amanzi { +namespace Flow { + +class PipeDrainEvaluator : public EvaluatorSecondaryMonotypeCV { + + public: + // constructor format for all derived classes + explicit + PipeDrainEvaluator(Teuchos::ParameterList& plist); + PipeDrainEvaluator(const PipeDrainEvaluator& other) = default; + + virtual Teuchos::RCP Clone() const override; + + protected: + void InitializeFromPlist_(); + + // Required methods from EvaluatorSecondaryMonotypeCV + virtual void Evaluate_(const State& S, + const std::vector& result) override; + virtual void EvaluatePartialDerivative_(const State& S, + const Key& wrt_key, const Tag& wrt_tag, + const std::vector& result) override; + + protected: + //What are the names for the flow depth above surface and hydraulic head in pipe network? + //Key pres_key_, cv_key_; + + double manhole_radius_; + double energy_losses_coeff_; //at manhole + double H_max_; //max height of pipe (considering a rectangular cross section) + double gravity_; //how to say this is an array? std vectors? + double H_; //surface elevation relative to the pipe flow hydraulic head + + private: + static Utils::RegisteredFactory reg_; + +}; + +} //namespace +} //namespace + +#endif diff --git a/src/pks/flow/source_terms/pipe_drain_evaluator_reg.hh b/src/pks/flow/source_terms/pipe_drain_evaluator_reg.hh new file mode 100644 index 000000000..91fbb2e7b --- /dev/null +++ b/src/pks/flow/source_terms/pipe_drain_evaluator_reg.hh @@ -0,0 +1,6 @@ +#include "pipe_drain_evaluator.hh" +namespace Amanzi { +namespace Flow { +Utils::RegisteredFactory PipeDrainEvaluator::reg_("pipe drain"); +} +} From 58f7a6a13ba1267ea3ece333c14f68114a03c956 Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Mon, 21 Nov 2022 20:29:07 -0700 Subject: [PATCH 02/58] Working on pipe drain evaluator --- .../flow/source_terms/pipe_drain_evaluator.cc | 42 ++++++++----------- .../flow/source_terms/pipe_drain_evaluator.hh | 2 +- 2 files changed, 19 insertions(+), 25 deletions(-) diff --git a/src/pks/flow/source_terms/pipe_drain_evaluator.cc b/src/pks/flow/source_terms/pipe_drain_evaluator.cc index 830fd3c15..d6337a3cc 100644 --- a/src/pks/flow/source_terms/pipe_drain_evaluator.cc +++ b/src/pks/flow/source_terms/pipe_drain_evaluator.cc @@ -18,15 +18,12 @@ namespace Flow { PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : EvaluatorSecondaryMonotypeCV(plist) { - M_ = plist_.get("molar mass", 0.0180153); - bar_ = plist_.get("allow negative water content", false); - rollover_ = plist_.get("water content rollover", 0.); manhole_radius_ = plist_.get("manhole radius", 0.24); energy_losses_coeff_ = plist.get("energy losses coeff", 0.1); H_max_ = plist.get("pipe max height", 0.1); - gravity_ = plist.get("gravity", 9.81);//TODO this needs to be converted to an array H_ = plist.get("bed flume height", 0.478); +// need to get an input array with info on where the manhole is Key domain_name = Keys::getDomain(my_keys_.front().first); Tag tag = my_keys_.front().second; @@ -51,39 +48,36 @@ void PipeDrainEvaluator::Evaluate_(const State& S, { Tag tag = my_keys_.front().second; Epetra_MultiVector& res = *result[0]->ViewComponent("cell",false); + +//change these two to what we need, let's call then hp and h for now +////// const Epetra_MultiVector& pres = *S.GetPtr(pres_key_, tag) ->ViewComponent("cell",false); const Epetra_MultiVector& cv = *S.GetPtr(cv_key_, tag) ->ViewComponent("cell",false); +///// - const double& p_atm = S.Get("atmospheric_pressure", Tags::DEFAULT); const auto& gravity = S.Get("gravity", Tags::DEFAULT); double gz = -gravity[2]; // check this + //do we also have pi somewhere? let's assume we do + //let's assume res for us is also defined on the cells int ncells = res.MyLength(); - if (bar_) { - for (int c=0; c!=ncells; ++c) { - res[0][c] = cv[0][c] * (pres[0][c] - p_atm) / (gz * M_); - } - } else if (rollover_ > 0.) { - for (int c=0; c!=ncells; ++c) { - double dp = pres[0][c] - p_atm; - double dp_eff = dp < 0. ? 0. : - dp < rollover_ ? - dp*dp/(2*rollover_) : - dp - rollover_/2.; - res[0][c] = cv[0][c] * dp_eff / (gz * M_); - } - } else { - for (int c=0; c!=ncells; ++c) { - res[0][c] = pres[0][c] < p_atm ? 0. : - cv[0][c] * (pres[0][c] - p_atm) / (gz * M_); - } +// we need a characteristic function to apply the source term ONLY where the manhole is + for (int c=0; c!=ncells; ++c) { + if (hp < H_) { + res[0][c] = - 4.0 / 3.0 * energy_losses_coeff_ * pi * manhole_radius * sqrt(2.0 * gz) * pow(h[c],3.0/2.0); + } else if (H_ < hp[c] && hp < (H_ + h[c]) ){ + res[0][c] = - energy_losses_coeff_ * pi * manhole_radius * manhole_radius * sqrt(2.0 * gz) * sqrt(h[c] + H_ - hp[c]); + } else if (hp > (H_ + h[c])) { + res[0][c] = energy_losses_coeff_ * pi * manhole_radius * manhole_radius * sqrt(2.0 * gz) * sqrt(hp[c] - H_ - h[c]); + } } -} +} +// do we need a derivative? void OverlandPressureWaterContentEvaluator::EvaluatePartialDerivative_(const State& S, const Key& wrt_key, const Tag& wrt_tag, const std::vector& result) diff --git a/src/pks/flow/source_terms/pipe_drain_evaluator.hh b/src/pks/flow/source_terms/pipe_drain_evaluator.hh index ff630360b..84852276a 100644 --- a/src/pks/flow/source_terms/pipe_drain_evaluator.hh +++ b/src/pks/flow/source_terms/pipe_drain_evaluator.hh @@ -45,8 +45,8 @@ class PipeDrainEvaluator : public EvaluatorSecondaryMonotypeCV { double manhole_radius_; double energy_losses_coeff_; //at manhole double H_max_; //max height of pipe (considering a rectangular cross section) - double gravity_; //how to say this is an array? std vectors? double H_; //surface elevation relative to the pipe flow hydraulic head +//we also need some array of cells that says in what cell the manhole is private: static Utils::RegisteredFactory reg_; From cdd816731dfc075e87393a4541dabad5da0c1d90 Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Mon, 28 Nov 2022 19:02:48 -0700 Subject: [PATCH 03/58] Worked on pipe drain evaluator Need to add partial derivative evaluator --- .../constitutive_relations/CMakeLists.txt | 2 +- .../sources/pipe_drain_evaluator.cc | 123 +++++++++++++++++ .../sources}/pipe_drain_evaluator.hh | 6 +- .../sources}/pipe_drain_evaluator_reg.hh | 0 .../flow/source_terms/pipe_drain_evaluator.cc | 125 ------------------ 5 files changed, 126 insertions(+), 130 deletions(-) create mode 100644 src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc rename src/pks/flow/{source_terms => constitutive_relations/sources}/pipe_drain_evaluator.hh (88%) rename src/pks/flow/{source_terms => constitutive_relations/sources}/pipe_drain_evaluator_reg.hh (100%) delete mode 100644 src/pks/flow/source_terms/pipe_drain_evaluator.cc diff --git a/src/pks/flow/constitutive_relations/CMakeLists.txt b/src/pks/flow/constitutive_relations/CMakeLists.txt index 37910dd55..e531f7887 100644 --- a/src/pks/flow/constitutive_relations/CMakeLists.txt +++ b/src/pks/flow/constitutive_relations/CMakeLists.txt @@ -6,7 +6,7 @@ # collect all sources -list(APPEND subdirs elevation overland_conductivity porosity water_content wrm) +list(APPEND subdirs elevation overland_conductivity porosity water_content wrm sources) set(ats_flow_relations_src_files "") set(ats_flow_relations_inc_files "") diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc new file mode 100644 index 000000000..b808570ca --- /dev/null +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -0,0 +1,123 @@ +/* -*- mode: c++; indent-tabs-mode: nil -*- */ + +/* + Evaluator for water drainage through a pipe network. + This depends on: + 1) flow depth above surface elevation + 2) hydraulic head in the pipe network + + Authors: Giacomo Capodaglio (gcapodaglio@lanl.gov) +*/ + +#include "pipe_drain_evaluator.hh" + +namespace Amanzi { +namespace Flow { + + +PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : + EvaluatorSecondaryMonotypeCV(plist) +{ + + manhole_radius_ = plist_.get("manhole radius", 0.24); + energy_losses_coeff_ = plist.get("energy losses coeff", 0.1); + H_max_ = plist.get("pipe max height", 0.1); + H_ = plist.get("bed flume height", 0.478); + + Key domain_name = Keys::getDomain(my_keys_.front().first); + Tag tag = my_keys_.front().second; + + // TODO + // my dependencies + head_key_ = Keys::readKey(plist_, domain_name, "ponded depth", "ponded_depth"); + dependencies_.insert(KeyTag{head_key_, tag}); + + mark_key_ = Keys::readKey(plist_, domain_name, "manhole locations", "manhole_locations"); + dependencies_.insert(KeyTag{mark_key_, tag}); + + + +} + + +Teuchos::RCP +PipeDrainEvaluator::Clone() const { + return Teuchos::rcp(new PipeDrainEvaluator(*this)); +} + + +void PipeDrainEvaluator::Evaluate_(const State& S, + const std::vector& result) +{ + Tag tag = my_keys_.front().second; + Epetra_MultiVector& res = *result[0]->ViewComponent("cell",false); + + const Epetra_MultiVector& head = *S.GetPtr(head_key_, tag) + ->ViewComponent("cell",false); + + const Epetra_MultiVector& mark = *S.GetPtr(mark_key_, tag) + ->ViewComponent("cell",false); + + const auto& gravity = S.Get("gravity", Tags::DEFAULT); + double gz = -gravity[2]; // check this + double pi = 3.14159265359; + + int ncells = res.MyLength(); + for (int c=0; c!=ncells; ++c) { + //if (hp < H_) { + res[0][c] = - mark[0][c] * 4.0 / 3.0 * energy_losses_coeff_ * pi * manhole_radius_ * sqrt(2.0 * gz) * pow(head[0][c],3.0/2.0); + // } else if (H_ < hp[c] && hp < (H_ + head[0][c]) ){ + //res[0][c] = - energy_losses_coeff_ * pi * manhole_radius_ * manhole_radius_ * sqrt(2.0 * gz) * sqrt(head[0][c] + H_ - hp[c]); + //} else if (hp > (H_ + head[0][c])) { + //res[0][c] = energy_losses_coeff_ * pi * manhole_radius_ * manhole_radius_ * sqrt(2.0 * gz) * sqrt(hp[c] - H_ - head[0][c]); + //} + } + +} + +// do we need a derivative? YES! +void PipeDrainEvaluator::EvaluatePartialDerivative_(const State& S, + const Key& wrt_key, const Tag& wrt_tag, + const std::vector& result) +{ + // Tag tag = my_keys_.front().second; + // AMANZI_ASSERT(wrt_key == pres_key_); + + // Epetra_MultiVector& res = *result[0]->ViewComponent("cell",false); + // const Epetra_MultiVector& pres = *S.GetPtr(pres_key_, tag) + // ->ViewComponent("cell",false); + + // const Epetra_MultiVector& cv = *S.GetPtr(cv_key_, tag) + // ->ViewComponent("cell",false); + + // const double& p_atm = S.Get("atmospheric_pressure", Tags::DEFAULT); + // const auto& gravity = S.Get("gravity", Tags::DEFAULT); + // double gz = -gravity[2]; // check this + + // if (wrt_key == pres_key_) { + // int ncells = res.MyLength(); + // if (bar_) { + // for (int c=0; c!=ncells; ++c) { + // res[0][c] = cv[0][c] / (gz * M_); + // } + // } else if (rollover_ > 0.) { + // for (int c=0; c!=ncells; ++c) { + // double dp = pres[0][c] - p_atm; + // double ddp_eff = dp < 0. ? 0. : + // dp < rollover_ ? dp/rollover_ : 1.; + // res[0][c] = cv[0][c] * ddp_eff / (gz * M_); + // } + // } else { + // for (int c=0; c!=ncells; ++c) { + // res[0][c] = pres[0][c] < p_atm ? 0. : + // cv[0][c] / (gz * M_); + // } + // } + // } else { + // res.PutScalar(0.); + // } +} + + +} //namespace +} //namespace diff --git a/src/pks/flow/source_terms/pipe_drain_evaluator.hh b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh similarity index 88% rename from src/pks/flow/source_terms/pipe_drain_evaluator.hh rename to src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh index 84852276a..d4721e739 100644 --- a/src/pks/flow/source_terms/pipe_drain_evaluator.hh +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh @@ -39,14 +39,12 @@ class PipeDrainEvaluator : public EvaluatorSecondaryMonotypeCV { const std::vector& result) override; protected: - //What are the names for the flow depth above surface and hydraulic head in pipe network? - //Key pres_key_, cv_key_; - + Key head_key_, mark_key_; + double manhole_radius_; double energy_losses_coeff_; //at manhole double H_max_; //max height of pipe (considering a rectangular cross section) double H_; //surface elevation relative to the pipe flow hydraulic head -//we also need some array of cells that says in what cell the manhole is private: static Utils::RegisteredFactory reg_; diff --git a/src/pks/flow/source_terms/pipe_drain_evaluator_reg.hh b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator_reg.hh similarity index 100% rename from src/pks/flow/source_terms/pipe_drain_evaluator_reg.hh rename to src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator_reg.hh diff --git a/src/pks/flow/source_terms/pipe_drain_evaluator.cc b/src/pks/flow/source_terms/pipe_drain_evaluator.cc deleted file mode 100644 index d6337a3cc..000000000 --- a/src/pks/flow/source_terms/pipe_drain_evaluator.cc +++ /dev/null @@ -1,125 +0,0 @@ -/* -*- mode: c++; indent-tabs-mode: nil -*- */ - -/* - Evaluator for water drainage through a pipe network. - This depends on: - 1) flow depth above surface elevation - 2) hydraulic head in the pipe network - - Authors: Giacomo Capodaglio (gcapodaglio@lanl.gov) -*/ - -#include "pipe_drain_evaluator.hh" - -namespace Amanzi { -namespace Flow { - - -PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : - EvaluatorSecondaryMonotypeCV(plist) -{ - - manhole_radius_ = plist_.get("manhole radius", 0.24); - energy_losses_coeff_ = plist.get("energy losses coeff", 0.1); - H_max_ = plist.get("pipe max height", 0.1); - H_ = plist.get("bed flume height", 0.478); -// need to get an input array with info on where the manhole is - - Key domain_name = Keys::getDomain(my_keys_.front().first); - Tag tag = my_keys_.front().second; - - // TODO - // my dependencies - //pres_key_ = Keys::readKey(plist_, domain_name, "pressure", "pressure"); - //dependencies_.insert(KeyTag{pres_key_, tag}); - //cv_key_ = Keys::readKey(plist_, domain_name, "cell volume", "cell_volume"); - //dependencies_.insert(KeyTag{cv_key_, tag}); -} - - -Teuchos::RCP -PipeDrainEvaluator::Clone() const { - return Teuchos::rcp(new PipeDrainEvaluator(*this)); -} - - -void PipeDrainEvaluator::Evaluate_(const State& S, - const std::vector& result) -{ - Tag tag = my_keys_.front().second; - Epetra_MultiVector& res = *result[0]->ViewComponent("cell",false); - -//change these two to what we need, let's call then hp and h for now -////// - const Epetra_MultiVector& pres = *S.GetPtr(pres_key_, tag) - ->ViewComponent("cell",false); - - const Epetra_MultiVector& cv = *S.GetPtr(cv_key_, tag) - ->ViewComponent("cell",false); -///// - - const auto& gravity = S.Get("gravity", Tags::DEFAULT); - double gz = -gravity[2]; // check this - //do we also have pi somewhere? let's assume we do - - //let's assume res for us is also defined on the cells - int ncells = res.MyLength(); -// we need a characteristic function to apply the source term ONLY where the manhole is - for (int c=0; c!=ncells; ++c) { - if (hp < H_) { - res[0][c] = - 4.0 / 3.0 * energy_losses_coeff_ * pi * manhole_radius * sqrt(2.0 * gz) * pow(h[c],3.0/2.0); - } else if (H_ < hp[c] && hp < (H_ + h[c]) ){ - res[0][c] = - energy_losses_coeff_ * pi * manhole_radius * manhole_radius * sqrt(2.0 * gz) * sqrt(h[c] + H_ - hp[c]); - } else if (hp > (H_ + h[c])) { - res[0][c] = energy_losses_coeff_ * pi * manhole_radius * manhole_radius * sqrt(2.0 * gz) * sqrt(hp[c] - H_ - h[c]); - } - } - -} - -// do we need a derivative? -void OverlandPressureWaterContentEvaluator::EvaluatePartialDerivative_(const State& S, - const Key& wrt_key, const Tag& wrt_tag, - const std::vector& result) -{ - Tag tag = my_keys_.front().second; - AMANZI_ASSERT(wrt_key == pres_key_); - - Epetra_MultiVector& res = *result[0]->ViewComponent("cell",false); - const Epetra_MultiVector& pres = *S.GetPtr(pres_key_, tag) - ->ViewComponent("cell",false); - - const Epetra_MultiVector& cv = *S.GetPtr(cv_key_, tag) - ->ViewComponent("cell",false); - - const double& p_atm = S.Get("atmospheric_pressure", Tags::DEFAULT); - const auto& gravity = S.Get("gravity", Tags::DEFAULT); - double gz = -gravity[2]; // check this - - if (wrt_key == pres_key_) { - int ncells = res.MyLength(); - if (bar_) { - for (int c=0; c!=ncells; ++c) { - res[0][c] = cv[0][c] / (gz * M_); - } - } else if (rollover_ > 0.) { - for (int c=0; c!=ncells; ++c) { - double dp = pres[0][c] - p_atm; - double ddp_eff = dp < 0. ? 0. : - dp < rollover_ ? dp/rollover_ : 1.; - res[0][c] = cv[0][c] * ddp_eff / (gz * M_); - } - } else { - for (int c=0; c!=ncells; ++c) { - res[0][c] = pres[0][c] < p_atm ? 0. : - cv[0][c] / (gz * M_); - } - } - } else { - res.PutScalar(0.); - } -} - - -} //namespace -} //namespace From 6fc5140a900434bc04c71a51ef120012125b266b Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Wed, 30 Nov 2022 10:25:46 -0700 Subject: [PATCH 04/58] Cleaned up code and added partial derivative --- .../sources/pipe_drain_evaluator.cc | 70 +++++++------------ 1 file changed, 27 insertions(+), 43 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index b808570ca..95c45cef4 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -27,7 +27,6 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : Key domain_name = Keys::getDomain(my_keys_.front().first); Tag tag = my_keys_.front().second; - // TODO // my dependencies head_key_ = Keys::readKey(plist_, domain_name, "ponded depth", "ponded_depth"); dependencies_.insert(KeyTag{head_key_, tag}); @@ -35,8 +34,6 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : mark_key_ = Keys::readKey(plist_, domain_name, "manhole locations", "manhole_locations"); dependencies_.insert(KeyTag{mark_key_, tag}); - - } @@ -61,61 +58,48 @@ void PipeDrainEvaluator::Evaluate_(const State& S, const auto& gravity = S.Get("gravity", Tags::DEFAULT); double gz = -gravity[2]; // check this double pi = 3.14159265359; + double manhole_area_ = pi * manhole_radius_ * manhole_radius_; int ncells = res.MyLength(); for (int c=0; c!=ncells; ++c) { //if (hp < H_) { - res[0][c] = - mark[0][c] * 4.0 / 3.0 * energy_losses_coeff_ * pi * manhole_radius_ * sqrt(2.0 * gz) * pow(head[0][c],3.0/2.0); + res[0][c] = - mark[0][c] * 4.0 / 3.0 * energy_losses_coeff_ * pi * manhole_radius_ * sqrt(2.0 * gz) * pow(head[0][c],3.0/2.0) / manhole_area_; // } else if (H_ < hp[c] && hp < (H_ + head[0][c]) ){ - //res[0][c] = - energy_losses_coeff_ * pi * manhole_radius_ * manhole_radius_ * sqrt(2.0 * gz) * sqrt(head[0][c] + H_ - hp[c]); + //res[0][c] = - energy_losses_coeff_ * manhole_area_ * sqrt(2.0 * gz) * sqrt(head[0][c] + H_ - hp[c]); //} else if (hp > (H_ + head[0][c])) { - //res[0][c] = energy_losses_coeff_ * pi * manhole_radius_ * manhole_radius_ * sqrt(2.0 * gz) * sqrt(hp[c] - H_ - head[0][c]); + //res[0][c] = energy_losses_coeff_ * manhole_area_ * sqrt(2.0 * gz) * sqrt(hp[c] - H_ - head[0][c]); //} } } -// do we need a derivative? YES! void PipeDrainEvaluator::EvaluatePartialDerivative_(const State& S, const Key& wrt_key, const Tag& wrt_tag, const std::vector& result) { - // Tag tag = my_keys_.front().second; - // AMANZI_ASSERT(wrt_key == pres_key_); - - // Epetra_MultiVector& res = *result[0]->ViewComponent("cell",false); - // const Epetra_MultiVector& pres = *S.GetPtr(pres_key_, tag) - // ->ViewComponent("cell",false); - - // const Epetra_MultiVector& cv = *S.GetPtr(cv_key_, tag) - // ->ViewComponent("cell",false); - - // const double& p_atm = S.Get("atmospheric_pressure", Tags::DEFAULT); - // const auto& gravity = S.Get("gravity", Tags::DEFAULT); - // double gz = -gravity[2]; // check this - - // if (wrt_key == pres_key_) { - // int ncells = res.MyLength(); - // if (bar_) { - // for (int c=0; c!=ncells; ++c) { - // res[0][c] = cv[0][c] / (gz * M_); - // } - // } else if (rollover_ > 0.) { - // for (int c=0; c!=ncells; ++c) { - // double dp = pres[0][c] - p_atm; - // double ddp_eff = dp < 0. ? 0. : - // dp < rollover_ ? dp/rollover_ : 1.; - // res[0][c] = cv[0][c] * ddp_eff / (gz * M_); - // } - // } else { - // for (int c=0; c!=ncells; ++c) { - // res[0][c] = pres[0][c] < p_atm ? 0. : - // cv[0][c] / (gz * M_); - // } - // } - // } else { - // res.PutScalar(0.); - // } + Tag tag = my_keys_.front().second; + AMANZI_ASSERT(wrt_key == head_key_); + + Epetra_MultiVector& res = *result[0]->ViewComponent("cell",false); + + const Epetra_MultiVector& head = *S.GetPtr(head_key_, tag) + ->ViewComponent("cell",false); + + const Epetra_MultiVector& mark = *S.GetPtr(mark_key_, tag) + ->ViewComponent("cell",false); + + const auto& gravity = S.Get("gravity", Tags::DEFAULT); + double gz = -gravity[2]; // check this + double pi = 3.14159265359; + double manhole_area_ = pi * manhole_radius_ * manhole_radius_; + + if (wrt_key == head_key_) { + int ncells = res.MyLength(); + for (int c=0; c!=ncells; ++c) { + res[0][c] = - mark[0][c] * 2.0 * energy_losses_coeff_ * pi * manhole_radius_ * sqrt(2.0 * gz) * sqrt(head[0][c]) / manhole_area_; + } + + } } From 7d37625aa9b3212bf2e03f3016c3f1e7cceae5c6 Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Wed, 30 Nov 2022 17:29:59 -0700 Subject: [PATCH 05/58] Commented out variables not needed at the moment --- .../constitutive_relations/sources/pipe_drain_evaluator.cc | 4 ++-- .../constitutive_relations/sources/pipe_drain_evaluator.hh | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index 95c45cef4..771049042 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -21,8 +21,8 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : manhole_radius_ = plist_.get("manhole radius", 0.24); energy_losses_coeff_ = plist.get("energy losses coeff", 0.1); - H_max_ = plist.get("pipe max height", 0.1); - H_ = plist.get("bed flume height", 0.478); + //H_max_ = plist.get("pipe max height", 0.1); + //H_ = plist.get("bed flume height", 0.478); Key domain_name = Keys::getDomain(my_keys_.front().first); Tag tag = my_keys_.front().second; diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh index d4721e739..95dd91fff 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh @@ -43,8 +43,8 @@ class PipeDrainEvaluator : public EvaluatorSecondaryMonotypeCV { double manhole_radius_; double energy_losses_coeff_; //at manhole - double H_max_; //max height of pipe (considering a rectangular cross section) - double H_; //surface elevation relative to the pipe flow hydraulic head +// double H_max_; //max height of pipe (considering a rectangular cross section) +// double H_; //surface elevation relative to the pipe flow hydraulic head private: static Utils::RegisteredFactory reg_; From 3d205269fe607e97b06d10dc09b92e999fa4aa78 Mon Sep 17 00:00:00 2001 From: Daniil Svyatsky Date: Wed, 14 Dec 2022 17:22:36 -0700 Subject: [PATCH 06/58] Added support for shallow water pk --- docs/documentation/source/ats_demos | 2 +- src/executables/CMakeLists.txt | 2 ++ src/executables/main.cc | 2 ++ testing/ats-regression-tests | 2 +- 4 files changed, 6 insertions(+), 2 deletions(-) diff --git a/docs/documentation/source/ats_demos b/docs/documentation/source/ats_demos index 21a87402b..50790e8d9 160000 --- a/docs/documentation/source/ats_demos +++ b/docs/documentation/source/ats_demos @@ -1 +1 @@ -Subproject commit 21a87402bc42404d6db8f038bb6e44814c203591 +Subproject commit 50790e8d956fc7cf4d0f63021464470ebb427f7f diff --git a/src/executables/CMakeLists.txt b/src/executables/CMakeLists.txt index adf7e704f..40c25b33d 100644 --- a/src/executables/CMakeLists.txt +++ b/src/executables/CMakeLists.txt @@ -30,6 +30,7 @@ include_directories(${DBG_SOURCE_DIR}) include_directories(${TRANSPORT_SOURCE_DIR}) get_property(CHEM_INCLUDES_DIR GLOBAL PROPERTY CHEM_INCLUDES_DIR) include_directories(${CHEM_INCLUDES_DIR}) +include_directories(${SHALLOW_WATER_SOURCE_DIR}) include_directories(${ATS_SOURCE_DIR}/src/constitutive_relations) include_directories(${ATS_SOURCE_DIR}/src/constitutive_relations/eos) @@ -98,6 +99,7 @@ set(amanzi_link_libs mesh_logical chemistry_pk #mpc_tree + shallow_water transport ) diff --git a/src/executables/main.cc b/src/executables/main.cc index 1027d0eda..aa6fbbb89 100644 --- a/src/executables/main.cc +++ b/src/executables/main.cc @@ -23,6 +23,8 @@ // registration files #include "ats_registration_files.hh" +#include "pks_shallow_water_registration.hh" +#include "numerical_flux_registration.hh" // include fenv if it exists diff --git a/testing/ats-regression-tests b/testing/ats-regression-tests index 6e8f21b0b..5d06f6100 160000 --- a/testing/ats-regression-tests +++ b/testing/ats-regression-tests @@ -1 +1 @@ -Subproject commit 6e8f21b0bf704cc612c7eda0413f32aa4f29f3cd +Subproject commit 5d06f61000f5d60787ad978ad2aef8555da95d9e From e248f68bf4334623d80b0fdc2e3c9be2468682e5 Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Fri, 5 May 2023 15:23:51 -0600 Subject: [PATCH 07/58] Fixed gravity in pipe drain evaluator --- .../constitutive_relations/sources/pipe_drain_evaluator.cc | 5 +++-- testing/ats-regression-tests | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index 771049042..c73a25f0d 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -29,6 +29,7 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : // my dependencies head_key_ = Keys::readKey(plist_, domain_name, "ponded depth", "ponded_depth"); + dependencies_.insert(KeyTag{head_key_, tag}); mark_key_ = Keys::readKey(plist_, domain_name, "manhole locations", "manhole_locations"); @@ -56,14 +57,14 @@ void PipeDrainEvaluator::Evaluate_(const State& S, ->ViewComponent("cell",false); const auto& gravity = S.Get("gravity", Tags::DEFAULT); - double gz = -gravity[2]; // check this + double gz = -gravity[1]; double pi = 3.14159265359; double manhole_area_ = pi * manhole_radius_ * manhole_radius_; int ncells = res.MyLength(); for (int c=0; c!=ncells; ++c) { //if (hp < H_) { - res[0][c] = - mark[0][c] * 4.0 / 3.0 * energy_losses_coeff_ * pi * manhole_radius_ * sqrt(2.0 * gz) * pow(head[0][c],3.0/2.0) / manhole_area_; + res[0][c] = - mark[0][c] * 4.0 / 3.0 * energy_losses_coeff_ * pi * manhole_radius_ * sqrt(2.0 * gz) * pow(head[0][c],3.0/2.0) / manhole_area_; // } else if (H_ < hp[c] && hp < (H_ + head[0][c]) ){ //res[0][c] = - energy_losses_coeff_ * manhole_area_ * sqrt(2.0 * gz) * sqrt(head[0][c] + H_ - hp[c]); //} else if (hp > (H_ + head[0][c])) { diff --git a/testing/ats-regression-tests b/testing/ats-regression-tests index 5d06f6100..42f0b6ca8 160000 --- a/testing/ats-regression-tests +++ b/testing/ats-regression-tests @@ -1 +1 @@ -Subproject commit 5d06f61000f5d60787ad978ad2aef8555da95d9e +Subproject commit 42f0b6ca80049b7129733e6c137d8f227fb8ff97 From 022009b3a1e23cc447cd9d04cd38ccffe153f3b6 Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Tue, 30 May 2023 16:42:36 -0600 Subject: [PATCH 08/58] Modified pipe_drain_evaluator.cc --- .../flow/constitutive_relations/sources/pipe_drain_evaluator.cc | 2 +- .../flow/constitutive_relations/sources/pipe_drain_evaluator.hh | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index c73a25f0d..0cdc41a5a 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -64,7 +64,7 @@ void PipeDrainEvaluator::Evaluate_(const State& S, int ncells = res.MyLength(); for (int c=0; c!=ncells; ++c) { //if (hp < H_) { - res[0][c] = - mark[0][c] * 4.0 / 3.0 * energy_losses_coeff_ * pi * manhole_radius_ * sqrt(2.0 * gz) * pow(head[0][c],3.0/2.0) / manhole_area_; + res[0][c] = - mark[0][c] * 4.0 / 3.0 * energy_losses_coeff_ * pi * manhole_radius_ * sqrt(2.0 * gz) * pow(head[0][c],3.0/2.0); // } else if (H_ < hp[c] && hp < (H_ + head[0][c]) ){ //res[0][c] = - energy_losses_coeff_ * manhole_area_ * sqrt(2.0 * gz) * sqrt(head[0][c] + H_ - hp[c]); //} else if (hp > (H_ + head[0][c])) { diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh index 95dd91fff..bfe5e4a5b 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh @@ -43,6 +43,7 @@ class PipeDrainEvaluator : public EvaluatorSecondaryMonotypeCV { double manhole_radius_; double energy_losses_coeff_; //at manhole + // double H_max_; //max height of pipe (considering a rectangular cross section) // double H_; //surface elevation relative to the pipe flow hydraulic head From cafe2329c09065f0e384a94dc6c131d8a204d851 Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Mon, 12 Jun 2023 10:40:14 -0600 Subject: [PATCH 09/58] Modified pipe drain evaluator To account for two way coupling between pipe and overland flow (need to test the new capability) --- .../sources/pipe_drain_evaluator.cc | 97 +++++++++++++------ .../sources/pipe_drain_evaluator.hh | 8 +- 2 files changed, 69 insertions(+), 36 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index 0cdc41a5a..b50e4efa4 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -3,8 +3,8 @@ /* Evaluator for water drainage through a pipe network. This depends on: - 1) flow depth above surface elevation - 2) hydraulic head in the pipe network + 1) flow depth above surface elevation (surface_depth_key_) + 2) hydraulic head in the pipe network (pipe_depth_key_) Authors: Giacomo Capodaglio (gcapodaglio@lanl.gov) */ @@ -20,20 +20,21 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : { manhole_radius_ = plist_.get("manhole radius", 0.24); - energy_losses_coeff_ = plist.get("energy losses coeff", 0.1); - //H_max_ = plist.get("pipe max height", 0.1); - //H_ = plist.get("bed flume height", 0.478); + energ_loss_coeff_ = plist.get("energy losses coeff", 0.1); + drain_length_ = plist.get("drain length", 0.478); Key domain_name = Keys::getDomain(my_keys_.front().first); Tag tag = my_keys_.front().second; // my dependencies - head_key_ = Keys::readKey(plist_, domain_name, "ponded depth", "ponded_depth"); + surface_depth_key_ = Keys::readKey(plist_, domain_name, "ponded depth", "ponded_depth"); + dependencies_.insert(KeyTag{surface_depth_key_, tag}); - dependencies_.insert(KeyTag{head_key_, tag}); + pipe_depth_key_ = Keys::readKey(plist_, domain_name, "water depth", "water_depth"); + dependencies_.insert(KeyTag{pipe_depth_key_, tag}); - mark_key_ = Keys::readKey(plist_, domain_name, "manhole locations", "manhole_locations"); - dependencies_.insert(KeyTag{mark_key_, tag}); + mask_key_ = Keys::readKey(plist_, domain_name, "manhole locations", "manhole_locations"); + dependencies_.insert(KeyTag{mask_key_, tag}); } @@ -50,26 +51,33 @@ void PipeDrainEvaluator::Evaluate_(const State& S, Tag tag = my_keys_.front().second; Epetra_MultiVector& res = *result[0]->ViewComponent("cell",false); - const Epetra_MultiVector& head = *S.GetPtr(head_key_, tag) + const Epetra_MultiVector& srfcDepth = *S.GetPtr(surface_depth_key_, tag) ->ViewComponent("cell",false); - const Epetra_MultiVector& mark = *S.GetPtr(mark_key_, tag) + const Epetra_MultiVector& pipeDepth = *S.GetPtr(pipe_depth_key_, tag) + ->ViewComponent("cell",false); + + const Epetra_MultiVector& mnhMask = *S.GetPtr(mask_key_, tag) ->ViewComponent("cell",false); const auto& gravity = S.Get("gravity", Tags::DEFAULT); - double gz = -gravity[1]; + double g = -gravity[1]; double pi = 3.14159265359; - double manhole_area_ = pi * manhole_radius_ * manhole_radius_; + double mnhArea = pi * manhole_radius_ * manhole_radius_; + double mnhPerimeter = 2.0 * pi * manhole_radius_; + double sqrtTwoG = sqrt(2.0 * g); int ncells = res.MyLength(); for (int c=0; c!=ncells; ++c) { - //if (hp < H_) { - res[0][c] = - mark[0][c] * 4.0 / 3.0 * energy_losses_coeff_ * pi * manhole_radius_ * sqrt(2.0 * gz) * pow(head[0][c],3.0/2.0); - // } else if (H_ < hp[c] && hp < (H_ + head[0][c]) ){ - //res[0][c] = - energy_losses_coeff_ * manhole_area_ * sqrt(2.0 * gz) * sqrt(head[0][c] + H_ - hp[c]); - //} else if (hp > (H_ + head[0][c])) { - //res[0][c] = energy_losses_coeff_ * manhole_area_ * sqrt(2.0 * gz) * sqrt(hp[c] - H_ - head[0][c]); - //} + if (pipeDepth[0][c] < drain_length_) { + res[0][c] = - mnhMask[0][c] * 2.0 / 3.0 * energ_loss_coeff_ * mnhPerimeter * sqrtTwoG * pow(srfcDepth[0][c],3.0/2.0); + } + else if (drain_length_ < pipeDepth[0][c] && pipeDepth[0][c] < (drain_length_ + srfcDepth[0][c]) ){ + res[0][c] = - mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG * sqrt(srfcDepth[0][c] + drain_length_ - pipeDepth[0][c]); + } + else if (pipeDepth[0][c] > (drain_length_ + srfcDepth[0][c])) { + res[0][c] = mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG * sqrt(pipeDepth[0][c] - drain_length_ - srfcDepth[0][c]); + } } } @@ -79,27 +87,54 @@ void PipeDrainEvaluator::EvaluatePartialDerivative_(const State& S, const std::vector& result) { Tag tag = my_keys_.front().second; - AMANZI_ASSERT(wrt_key == head_key_); Epetra_MultiVector& res = *result[0]->ViewComponent("cell",false); - - const Epetra_MultiVector& head = *S.GetPtr(head_key_, tag) - ->ViewComponent("cell",false); - const Epetra_MultiVector& mark = *S.GetPtr(mark_key_, tag) + const Epetra_MultiVector& srfcDepth = *S.GetPtr(surface_depth_key_, tag) + ->ViewComponent("cell",false); + + const Epetra_MultiVector& pipeDepth = *S.GetPtr(pipe_depth_key_, tag) + ->ViewComponent("cell",false); + + const Epetra_MultiVector& mnhMask = *S.GetPtr(mask_key_, tag) ->ViewComponent("cell",false); const auto& gravity = S.Get("gravity", Tags::DEFAULT); - double gz = -gravity[2]; // check this + double g = -gravity[1]; double pi = 3.14159265359; - double manhole_area_ = pi * manhole_radius_ * manhole_radius_; + double mnhArea = pi * manhole_radius_ * manhole_radius_; + double mnhPerimeter = 2.0 * pi * manhole_radius_; + double sqrtTwoG = sqrt(2.0 * g); - if (wrt_key == head_key_) { - int ncells = res.MyLength(); + int ncells = res.MyLength(); + if (wrt_key == surface_depth_key_) { for (int c=0; c!=ncells; ++c) { - res[0][c] = - mark[0][c] * 2.0 * energy_losses_coeff_ * pi * manhole_radius_ * sqrt(2.0 * gz) * sqrt(head[0][c]) / manhole_area_; + if (pipeDepth[0][c] < drain_length_) { + res[0][c] = - mnhMask[0][c] * energ_loss_coeff_ * mnhPerimeter* sqrtTwoG * sqrt(srfcDepth[0][c]); + } + else if (drain_length_ < pipeDepth[0][c] && pipeDepth[0][c] < (drain_length_ + srfcDepth[0][c]) ){ + res[0][c] = - 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(srfcDepth[0][c] + drain_length_ - pipeDepth[0][c]); + } + else if (pipeDepth[0][c] > (drain_length_ + srfcDepth[0][c])) { + res[0][c] = - 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(pipeDepth[0][c] - drain_length_ - srfcDepth[0][c]); + } + } + } + else if (wrt_key == pipe_depth_key_) { + for (int c=0; c!=ncells; ++c) { + if (pipeDepth[0][c] < drain_length_) { + res[0][c] = 0.0; + } + else if (drain_length_ < pipeDepth[0][c] && pipeDepth[0][c] < (drain_length_ + srfcDepth[0][c]) ){ + res[0][c] = 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(srfcDepth[0][c] + drain_length_ - pipeDepth[0][c]); + } + else if (pipeDepth[0][c] > (drain_length_ + srfcDepth[0][c])) { + res[0][c] = 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(pipeDepth[0][c] - drain_length_ - srfcDepth[0][c]); + } } - + } + else { + AMANZI_ASSERT(0); } } diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh index bfe5e4a5b..6054c5928 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh @@ -39,13 +39,11 @@ class PipeDrainEvaluator : public EvaluatorSecondaryMonotypeCV { const std::vector& result) override; protected: - Key head_key_, mark_key_; + Key surface_depth_key_, pipe_depth_key_, mask_key_; double manhole_radius_; - double energy_losses_coeff_; //at manhole - -// double H_max_; //max height of pipe (considering a rectangular cross section) -// double H_; //surface elevation relative to the pipe flow hydraulic head + double energ_loss_coeff_; // energy losses coefficient at manhole + double drain_length_; // drain conduit length private: static Utils::RegisteredFactory reg_; From 5e6e6e09175da7b6de42e56a1223ebc05e32236b Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Mon, 12 Jun 2023 12:27:38 -0600 Subject: [PATCH 10/58] Had to modify a dependency of the pipe drain eval Since it does not depend on the pipe water depth but on the pipe pressure head --- .../sources/pipe_drain_evaluator.cc | 43 +++++++++---------- .../sources/pipe_drain_evaluator.hh | 2 +- 2 files changed, 22 insertions(+), 23 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index b50e4efa4..736440aa6 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -4,7 +4,7 @@ Evaluator for water drainage through a pipe network. This depends on: 1) flow depth above surface elevation (surface_depth_key_) - 2) hydraulic head in the pipe network (pipe_depth_key_) + 2) hydraulic head in the pipe network (pressure_head_key_) Authors: Giacomo Capodaglio (gcapodaglio@lanl.gov) */ @@ -30,11 +30,10 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : surface_depth_key_ = Keys::readKey(plist_, domain_name, "ponded depth", "ponded_depth"); dependencies_.insert(KeyTag{surface_depth_key_, tag}); - pipe_depth_key_ = Keys::readKey(plist_, domain_name, "water depth", "water_depth"); - dependencies_.insert(KeyTag{pipe_depth_key_, tag}); + pressure_head_key_ = Keys::readKey(plist_, domain_name, " pressure head", "pressure_head"); + dependencies_.insert(KeyTag{pressure_head_key_, tag}); mask_key_ = Keys::readKey(plist_, domain_name, "manhole locations", "manhole_locations"); - dependencies_.insert(KeyTag{mask_key_, tag}); } @@ -54,7 +53,7 @@ void PipeDrainEvaluator::Evaluate_(const State& S, const Epetra_MultiVector& srfcDepth = *S.GetPtr(surface_depth_key_, tag) ->ViewComponent("cell",false); - const Epetra_MultiVector& pipeDepth = *S.GetPtr(pipe_depth_key_, tag) + const Epetra_MultiVector& pressHead = *S.GetPtr(pressure_head_key_, tag) ->ViewComponent("cell",false); const Epetra_MultiVector& mnhMask = *S.GetPtr(mask_key_, tag) @@ -69,14 +68,14 @@ void PipeDrainEvaluator::Evaluate_(const State& S, int ncells = res.MyLength(); for (int c=0; c!=ncells; ++c) { - if (pipeDepth[0][c] < drain_length_) { + if (pressHead[0][c] < drain_length_) { res[0][c] = - mnhMask[0][c] * 2.0 / 3.0 * energ_loss_coeff_ * mnhPerimeter * sqrtTwoG * pow(srfcDepth[0][c],3.0/2.0); } - else if (drain_length_ < pipeDepth[0][c] && pipeDepth[0][c] < (drain_length_ + srfcDepth[0][c]) ){ - res[0][c] = - mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG * sqrt(srfcDepth[0][c] + drain_length_ - pipeDepth[0][c]); + else if (drain_length_ < pressHead[0][c] && pressHead[0][c] < (drain_length_ + srfcDepth[0][c]) ){ + res[0][c] = - mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG * sqrt(srfcDepth[0][c] + drain_length_ - pressHead[0][c]); } - else if (pipeDepth[0][c] > (drain_length_ + srfcDepth[0][c])) { - res[0][c] = mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG * sqrt(pipeDepth[0][c] - drain_length_ - srfcDepth[0][c]); + else if (pressHead[0][c] > (drain_length_ + srfcDepth[0][c])) { + res[0][c] = mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG * sqrt(pressHead[0][c] - drain_length_ - srfcDepth[0][c]); } } @@ -93,7 +92,7 @@ void PipeDrainEvaluator::EvaluatePartialDerivative_(const State& S, const Epetra_MultiVector& srfcDepth = *S.GetPtr(surface_depth_key_, tag) ->ViewComponent("cell",false); - const Epetra_MultiVector& pipeDepth = *S.GetPtr(pipe_depth_key_, tag) + const Epetra_MultiVector& pressHead = *S.GetPtr(pressure_head_key_, tag) ->ViewComponent("cell",false); const Epetra_MultiVector& mnhMask = *S.GetPtr(mask_key_, tag) @@ -109,27 +108,27 @@ void PipeDrainEvaluator::EvaluatePartialDerivative_(const State& S, int ncells = res.MyLength(); if (wrt_key == surface_depth_key_) { for (int c=0; c!=ncells; ++c) { - if (pipeDepth[0][c] < drain_length_) { + if (pressHead[0][c] < drain_length_) { res[0][c] = - mnhMask[0][c] * energ_loss_coeff_ * mnhPerimeter* sqrtTwoG * sqrt(srfcDepth[0][c]); } - else if (drain_length_ < pipeDepth[0][c] && pipeDepth[0][c] < (drain_length_ + srfcDepth[0][c]) ){ - res[0][c] = - 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(srfcDepth[0][c] + drain_length_ - pipeDepth[0][c]); + else if (drain_length_ < pressHead[0][c] && pressHead[0][c] < (drain_length_ + srfcDepth[0][c]) ){ + res[0][c] = - 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(srfcDepth[0][c] + drain_length_ - pressHead[0][c]); } - else if (pipeDepth[0][c] > (drain_length_ + srfcDepth[0][c])) { - res[0][c] = - 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(pipeDepth[0][c] - drain_length_ - srfcDepth[0][c]); + else if (pressHead[0][c] > (drain_length_ + srfcDepth[0][c])) { + res[0][c] = - 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(pressHead[0][c] - drain_length_ - srfcDepth[0][c]); } } } - else if (wrt_key == pipe_depth_key_) { + else if (wrt_key == pressure_head_key_) { for (int c=0; c!=ncells; ++c) { - if (pipeDepth[0][c] < drain_length_) { + if (pressHead[0][c] < drain_length_) { res[0][c] = 0.0; } - else if (drain_length_ < pipeDepth[0][c] && pipeDepth[0][c] < (drain_length_ + srfcDepth[0][c]) ){ - res[0][c] = 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(srfcDepth[0][c] + drain_length_ - pipeDepth[0][c]); + else if (drain_length_ < pressHead[0][c] && pressHead[0][c] < (drain_length_ + srfcDepth[0][c]) ){ + res[0][c] = 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(srfcDepth[0][c] + drain_length_ - pressHead[0][c]); } - else if (pipeDepth[0][c] > (drain_length_ + srfcDepth[0][c])) { - res[0][c] = 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(pipeDepth[0][c] - drain_length_ - srfcDepth[0][c]); + else if (pressHead[0][c] > (drain_length_ + srfcDepth[0][c])) { + res[0][c] = 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(pressHead[0][c] - drain_length_ - srfcDepth[0][c]); } } } diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh index 6054c5928..41d08c600 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh @@ -39,7 +39,7 @@ class PipeDrainEvaluator : public EvaluatorSecondaryMonotypeCV { const std::vector& result) override; protected: - Key surface_depth_key_, pipe_depth_key_, mask_key_; + Key surface_depth_key_, pressure_head_key_, mask_key_; double manhole_radius_; double energ_loss_coeff_; // energy losses coefficient at manhole From f8a26a5a585694364453638f240ea1f5befc07fe Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Mon, 12 Jun 2023 16:03:58 -0600 Subject: [PATCH 11/58] Fixed a silly typo --- .../flow/constitutive_relations/sources/pipe_drain_evaluator.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index 736440aa6..57ddcec2e 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -30,7 +30,7 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : surface_depth_key_ = Keys::readKey(plist_, domain_name, "ponded depth", "ponded_depth"); dependencies_.insert(KeyTag{surface_depth_key_, tag}); - pressure_head_key_ = Keys::readKey(plist_, domain_name, " pressure head", "pressure_head"); + pressure_head_key_ = Keys::readKey(plist_, domain_name, "pressure head", "pressure_head"); dependencies_.insert(KeyTag{pressure_head_key_, tag}); mask_key_ = Keys::readKey(plist_, domain_name, "manhole locations", "manhole_locations"); From 8ba1483c0aca59d05182976bd043e760d27b206c Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Thu, 15 Jun 2023 09:54:09 -0600 Subject: [PATCH 12/58] Modified pipe drain evaluator With new options for two way coupling between pipe and sw --- .../sources/pipe_drain_evaluator.cc | 126 +++++++++++------- .../sources/pipe_drain_evaluator.hh | 7 +- 2 files changed, 80 insertions(+), 53 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index 57ddcec2e..1eb5e4e0f 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -22,18 +22,21 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : manhole_radius_ = plist_.get("manhole radius", 0.24); energ_loss_coeff_ = plist.get("energy losses coeff", 0.1); drain_length_ = plist.get("drain length", 0.478); - - Key domain_name = Keys::getDomain(my_keys_.front().first); + sw_domain_name_ = plist.get("sw domain name", "surface"); + pipe_domain_name_ = plist.get("pipe domain name", ""); Tag tag = my_keys_.front().second; // my dependencies - surface_depth_key_ = Keys::readKey(plist_, domain_name, "ponded depth", "ponded_depth"); + surface_depth_key_ = Keys::readKey(plist_, sw_domain_name_, "ponded depth", "ponded_depth"); dependencies_.insert(KeyTag{surface_depth_key_, tag}); - pressure_head_key_ = Keys::readKey(plist_, domain_name, "pressure head", "pressure_head"); - dependencies_.insert(KeyTag{pressure_head_key_, tag}); + if(!pipe_domain_name_.empty()){ + pressure_head_key_ = Keys::readKey(plist_, pipe_domain_name_, "pressure head", "pressure_head"); + dependencies_.insert(KeyTag{pressure_head_key_, tag}); + } - mask_key_ = Keys::readKey(plist_, domain_name, "manhole locations", "manhole_locations"); + mask_key_ = Keys::readKey(plist_, sw_domain_name_, "manhole locations", "manhole_locations"); + dependencies_.insert(KeyTag{mask_key_, tag}); } @@ -51,33 +54,42 @@ void PipeDrainEvaluator::Evaluate_(const State& S, Epetra_MultiVector& res = *result[0]->ViewComponent("cell",false); const Epetra_MultiVector& srfcDepth = *S.GetPtr(surface_depth_key_, tag) - ->ViewComponent("cell",false); - - const Epetra_MultiVector& pressHead = *S.GetPtr(pressure_head_key_, tag) - ->ViewComponent("cell",false); + ->ViewComponent("cell",false); const Epetra_MultiVector& mnhMask = *S.GetPtr(mask_key_, tag) ->ViewComponent("cell",false); const auto& gravity = S.Get("gravity", Tags::DEFAULT); - double g = -gravity[1]; + double g = -gravity[1]; double pi = 3.14159265359; double mnhArea = pi * manhole_radius_ * manhole_radius_; - double mnhPerimeter = 2.0 * pi * manhole_radius_; - double sqrtTwoG = sqrt(2.0 * g); + double mnhPerimeter = 2.0 * pi * manhole_radius_; + double sqrtTwoG = sqrt(2.0 * g); int ncells = res.MyLength(); - for (int c=0; c!=ncells; ++c) { - if (pressHead[0][c] < drain_length_) { - res[0][c] = - mnhMask[0][c] * 2.0 / 3.0 * energ_loss_coeff_ * mnhPerimeter * sqrtTwoG * pow(srfcDepth[0][c],3.0/2.0); - } - else if (drain_length_ < pressHead[0][c] && pressHead[0][c] < (drain_length_ + srfcDepth[0][c]) ){ - res[0][c] = - mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG * sqrt(srfcDepth[0][c] + drain_length_ - pressHead[0][c]); - } - else if (pressHead[0][c] > (drain_length_ + srfcDepth[0][c])) { - res[0][c] = mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG * sqrt(pressHead[0][c] - drain_length_ - srfcDepth[0][c]); - } + + if(!pipe_domain_name_.empty()){ + + const Epetra_MultiVector& pressHead = *S.GetPtr(pressure_head_key_, tag) + ->ViewComponent("cell",false); + + for (int c=0; c!=ncells; ++c) { + if (pressHead[0][c] < drain_length_) { + res[0][c] = - mnhMask[0][c] * 2.0 / 3.0 * energ_loss_coeff_ * mnhPerimeter * sqrtTwoG * pow(srfcDepth[0][c],3.0/2.0); + } + else if (drain_length_ < pressHead[0][c] && pressHead[0][c] < (drain_length_ + srfcDepth[0][c]) ){ + res[0][c] = - mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG * sqrt(srfcDepth[0][c] + drain_length_ - pressHead[0][c]); + } + else if (pressHead[0][c] > (drain_length_ + srfcDepth[0][c])) { + res[0][c] = mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG * sqrt(pressHead[0][c] - drain_length_ - srfcDepth[0][c]); + } + } } + else { + for (int c=0; c!=ncells; ++c) { + res[0][c] = - mnhMask[0][c] * 2.0 / 3.0 * energ_loss_coeff_ * mnhPerimeter * sqrtTwoG * pow(srfcDepth[0][c],3.0/2.0); + } + } } @@ -92,9 +104,6 @@ void PipeDrainEvaluator::EvaluatePartialDerivative_(const State& S, const Epetra_MultiVector& srfcDepth = *S.GetPtr(surface_depth_key_, tag) ->ViewComponent("cell",false); - const Epetra_MultiVector& pressHead = *S.GetPtr(pressure_head_key_, tag) - ->ViewComponent("cell",false); - const Epetra_MultiVector& mnhMask = *S.GetPtr(mask_key_, tag) ->ViewComponent("cell",false); @@ -106,34 +115,51 @@ void PipeDrainEvaluator::EvaluatePartialDerivative_(const State& S, double sqrtTwoG = sqrt(2.0 * g); int ncells = res.MyLength(); - if (wrt_key == surface_depth_key_) { - for (int c=0; c!=ncells; ++c) { - if (pressHead[0][c] < drain_length_) { - res[0][c] = - mnhMask[0][c] * energ_loss_coeff_ * mnhPerimeter* sqrtTwoG * sqrt(srfcDepth[0][c]); - } - else if (drain_length_ < pressHead[0][c] && pressHead[0][c] < (drain_length_ + srfcDepth[0][c]) ){ - res[0][c] = - 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(srfcDepth[0][c] + drain_length_ - pressHead[0][c]); - } - else if (pressHead[0][c] > (drain_length_ + srfcDepth[0][c])) { - res[0][c] = - 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(pressHead[0][c] - drain_length_ - srfcDepth[0][c]); - } - } - } - else if (wrt_key == pressure_head_key_) { - for (int c=0; c!=ncells; ++c) { - if (pressHead[0][c] < drain_length_) { - res[0][c] = 0.0; - } - else if (drain_length_ < pressHead[0][c] && pressHead[0][c] < (drain_length_ + srfcDepth[0][c]) ){ - res[0][c] = 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(srfcDepth[0][c] + drain_length_ - pressHead[0][c]); - } - else if (pressHead[0][c] > (drain_length_ + srfcDepth[0][c])) { - res[0][c] = 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(pressHead[0][c] - drain_length_ - srfcDepth[0][c]); + + if(!pipe_domain_name_.empty()){ + + const Epetra_MultiVector& pressHead = *S.GetPtr(pressure_head_key_, tag) + ->ViewComponent("cell",false); + + if (wrt_key == surface_depth_key_) { + for (int c=0; c!=ncells; ++c) { + if (pressHead[0][c] < drain_length_) { + res[0][c] = - mnhMask[0][c] * energ_loss_coeff_ * mnhPerimeter* sqrtTwoG * sqrt(srfcDepth[0][c]); + } + else if (drain_length_ < pressHead[0][c] && pressHead[0][c] < (drain_length_ + srfcDepth[0][c]) ){ + res[0][c] = - 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(srfcDepth[0][c] + drain_length_ - pressHead[0][c]); + } + else if (pressHead[0][c] > (drain_length_ + srfcDepth[0][c])) { + res[0][c] = - 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(pressHead[0][c] - drain_length_ - srfcDepth[0][c]); + } + } + } + else if (wrt_key == pressure_head_key_) { + for (int c=0; c!=ncells; ++c) { + if (pressHead[0][c] < drain_length_) { + res[0][c] = 0.0; + } + else if (drain_length_ < pressHead[0][c] && pressHead[0][c] < (drain_length_ + srfcDepth[0][c]) ){ + res[0][c] = 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(srfcDepth[0][c] + drain_length_ - pressHead[0][c]); + } + else if (pressHead[0][c] > (drain_length_ + srfcDepth[0][c])) { + res[0][c] = 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(pressHead[0][c] - drain_length_ - srfcDepth[0][c]); + } } } + else { + AMANZI_ASSERT(0); + } } else { - AMANZI_ASSERT(0); + if (wrt_key == surface_depth_key_) { + for (int c=0; c!=ncells; ++c) { + res[0][c] = - mnhMask[0][c] * energ_loss_coeff_ * mnhPerimeter* sqrtTwoG * sqrt(srfcDepth[0][c]); + } + } + else { + AMANZI_ASSERT(0); + } } } diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh index 41d08c600..3d7af6c58 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh @@ -40,10 +40,11 @@ class PipeDrainEvaluator : public EvaluatorSecondaryMonotypeCV { protected: Key surface_depth_key_, pressure_head_key_, mask_key_; + Key sw_domain_name_, pipe_domain_name_; - double manhole_radius_; - double energ_loss_coeff_; // energy losses coefficient at manhole - double drain_length_; // drain conduit length + double manhole_radius_; + double energ_loss_coeff_; // energy losses coefficient at manhole + double drain_length_; // drain conduit length private: static Utils::RegisteredFactory reg_; From ac856b5791c9d9202d4f13e748201684dc7bee45 Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Thu, 24 Aug 2023 10:30:12 -0600 Subject: [PATCH 13/58] Modified energy loss coefficients for pipe drain According to Table 3 in Rubinato et al. "Experimental calibration and validation of sewer/surface flow exchange equations in steady and unsteady flow conditions" --- .../sources/pipe_drain_evaluator.cc | 29 ++++++++++++------- .../sources/pipe_drain_evaluator.hh | 7 ++++- testing/ats-regression-tests | 2 +- 3 files changed, 25 insertions(+), 13 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index 1eb5e4e0f..17a01c659 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -20,7 +20,9 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : { manhole_radius_ = plist_.get("manhole radius", 0.24); - energ_loss_coeff_ = plist.get("energy losses coeff", 0.1); + energ_loss_coeff_weir_ = plist.get("energy losses coeff weir", 0.54); + energ_loss_coeff_subweir_ = plist.get("energy losses coeff submerged weir", 0.056); + energ_loss_coeff_orifice_ = plist.get("energy losses coeff orifice", 0.167); drain_length_ = plist.get("drain length", 0.478); sw_domain_name_ = plist.get("sw domain name", "surface"); pipe_domain_name_ = plist.get("pipe domain name", ""); @@ -75,19 +77,20 @@ void PipeDrainEvaluator::Evaluate_(const State& S, for (int c=0; c!=ncells; ++c) { if (pressHead[0][c] < drain_length_) { - res[0][c] = - mnhMask[0][c] * 2.0 / 3.0 * energ_loss_coeff_ * mnhPerimeter * sqrtTwoG * pow(srfcDepth[0][c],3.0/2.0); + res[0][c] = - mnhMask[0][c] * 2.0 / 3.0 * energ_loss_coeff_weir_ * mnhPerimeter * sqrtTwoG * pow(srfcDepth[0][c],3.0/2.0); } else if (drain_length_ < pressHead[0][c] && pressHead[0][c] < (drain_length_ + srfcDepth[0][c]) ){ - res[0][c] = - mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG * sqrt(srfcDepth[0][c] + drain_length_ - pressHead[0][c]); + res[0][c] = - mnhMask[0][c] * energ_loss_coeff_subweir_ * mnhArea * sqrtTwoG + * sqrt(srfcDepth[0][c] + drain_length_ - pressHead[0][c]); } else if (pressHead[0][c] > (drain_length_ + srfcDepth[0][c])) { - res[0][c] = mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG * sqrt(pressHead[0][c] - drain_length_ - srfcDepth[0][c]); + res[0][c] = mnhMask[0][c] * energ_loss_coeff_orifice_ * mnhArea * sqrtTwoG * sqrt(pressHead[0][c] - drain_length_ - srfcDepth[0][c]); } } } else { for (int c=0; c!=ncells; ++c) { - res[0][c] = - mnhMask[0][c] * 2.0 / 3.0 * energ_loss_coeff_ * mnhPerimeter * sqrtTwoG * pow(srfcDepth[0][c],3.0/2.0); + res[0][c] = - mnhMask[0][c] * 2.0 / 3.0 * energ_loss_coeff_weir_ * mnhPerimeter * sqrtTwoG * pow(srfcDepth[0][c],3.0/2.0); } } @@ -124,13 +127,15 @@ void PipeDrainEvaluator::EvaluatePartialDerivative_(const State& S, if (wrt_key == surface_depth_key_) { for (int c=0; c!=ncells; ++c) { if (pressHead[0][c] < drain_length_) { - res[0][c] = - mnhMask[0][c] * energ_loss_coeff_ * mnhPerimeter* sqrtTwoG * sqrt(srfcDepth[0][c]); + res[0][c] = - mnhMask[0][c] * energ_loss_coeff_weir_ * mnhPerimeter* sqrtTwoG * sqrt(srfcDepth[0][c]); } else if (drain_length_ < pressHead[0][c] && pressHead[0][c] < (drain_length_ + srfcDepth[0][c]) ){ - res[0][c] = - 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(srfcDepth[0][c] + drain_length_ - pressHead[0][c]); + res[0][c] = - 0.5 * mnhMask[0][c] * energ_loss_coeff_subweir_ * mnhArea * sqrtTwoG + / sqrt(srfcDepth[0][c] + drain_length_ - pressHead[0][c]); } else if (pressHead[0][c] > (drain_length_ + srfcDepth[0][c])) { - res[0][c] = - 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(pressHead[0][c] - drain_length_ - srfcDepth[0][c]); + res[0][c] = - 0.5 * mnhMask[0][c] * energ_loss_coeff_orifice_ * mnhArea * sqrtTwoG + / sqrt(pressHead[0][c] - drain_length_ - srfcDepth[0][c]); } } } @@ -140,10 +145,12 @@ void PipeDrainEvaluator::EvaluatePartialDerivative_(const State& S, res[0][c] = 0.0; } else if (drain_length_ < pressHead[0][c] && pressHead[0][c] < (drain_length_ + srfcDepth[0][c]) ){ - res[0][c] = 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(srfcDepth[0][c] + drain_length_ - pressHead[0][c]); + res[0][c] = 0.5 * mnhMask[0][c] * energ_loss_coeff_subweir_ * mnhArea * sqrtTwoG + / sqrt(srfcDepth[0][c] + drain_length_ - pressHead[0][c]); } else if (pressHead[0][c] > (drain_length_ + srfcDepth[0][c])) { - res[0][c] = 0.5 * mnhMask[0][c] * energ_loss_coeff_ * mnhArea * sqrtTwoG / sqrt(pressHead[0][c] - drain_length_ - srfcDepth[0][c]); + res[0][c] = 0.5 * mnhMask[0][c] * energ_loss_coeff_orifice_ * mnhArea * sqrtTwoG + / sqrt(pressHead[0][c] - drain_length_ - srfcDepth[0][c]); } } } @@ -154,7 +161,7 @@ void PipeDrainEvaluator::EvaluatePartialDerivative_(const State& S, else { if (wrt_key == surface_depth_key_) { for (int c=0; c!=ncells; ++c) { - res[0][c] = - mnhMask[0][c] * energ_loss_coeff_ * mnhPerimeter* sqrtTwoG * sqrt(srfcDepth[0][c]); + res[0][c] = - mnhMask[0][c] * energ_loss_coeff_weir_ * mnhPerimeter* sqrtTwoG * sqrt(srfcDepth[0][c]); } } else { diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh index 3d7af6c58..d8d845088 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh @@ -43,7 +43,12 @@ class PipeDrainEvaluator : public EvaluatorSecondaryMonotypeCV { Key sw_domain_name_, pipe_domain_name_; double manhole_radius_; - double energ_loss_coeff_; // energy losses coefficient at manhole + // energy losses coefficients at manhole + // see Table 3 in "Experimental calibration and validation of sewer/surface + // flow exchange equations in steady and unsteady flow conditions" by Rubinato et al. + double energ_loss_coeff_weir_; // weir + double energ_loss_coeff_subweir_; // submerged weir + double energ_loss_coeff_orifice_; // orifice double drain_length_; // drain conduit length private: diff --git a/testing/ats-regression-tests b/testing/ats-regression-tests index 42f0b6ca8..7742df4ce 160000 --- a/testing/ats-regression-tests +++ b/testing/ats-regression-tests @@ -1 +1 @@ -Subproject commit 42f0b6ca80049b7129733e6c137d8f227fb8ff97 +Subproject commit 7742df4ce98333c4d0e975829f7e073a85ec6409 From 29b976f6f307349b2d1a1188aa3a63c7854f3ab6 Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Thu, 21 Sep 2023 10:36:11 -0600 Subject: [PATCH 14/58] Changed def of g in pipe drain evaluator --- .../constitutive_relations/sources/pipe_drain_evaluator.cc | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index 17a01c659..f73e5817d 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -61,8 +61,7 @@ void PipeDrainEvaluator::Evaluate_(const State& S, const Epetra_MultiVector& mnhMask = *S.GetPtr(mask_key_, tag) ->ViewComponent("cell",false); - const auto& gravity = S.Get("gravity", Tags::DEFAULT); - double g = -gravity[1]; + double g = norm(S.Get("gravity")); double pi = 3.14159265359; double mnhArea = pi * manhole_radius_ * manhole_radius_; double mnhPerimeter = 2.0 * pi * manhole_radius_; @@ -110,8 +109,7 @@ void PipeDrainEvaluator::EvaluatePartialDerivative_(const State& S, const Epetra_MultiVector& mnhMask = *S.GetPtr(mask_key_, tag) ->ViewComponent("cell",false); - const auto& gravity = S.Get("gravity", Tags::DEFAULT); - double g = -gravity[1]; + double g = norm(S.Get("gravity")); double pi = 3.14159265359; double mnhArea = pi * manhole_radius_ * manhole_radius_; double mnhPerimeter = 2.0 * pi * manhole_radius_; From bb1c4e680b454a6eb10cc206118d4b3acc66e104 Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Thu, 21 Sep 2023 10:42:51 -0600 Subject: [PATCH 15/58] Fixed typo --- .../constitutive_relations/sources/pipe_drain_evaluator.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index f73e5817d..c9dc7fe4b 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -61,7 +61,7 @@ void PipeDrainEvaluator::Evaluate_(const State& S, const Epetra_MultiVector& mnhMask = *S.GetPtr(mask_key_, tag) ->ViewComponent("cell",false); - double g = norm(S.Get("gravity")); + double g = norm(S.Get("gravity", Tags::DEFAULT)); double pi = 3.14159265359; double mnhArea = pi * manhole_radius_ * manhole_radius_; double mnhPerimeter = 2.0 * pi * manhole_radius_; @@ -109,7 +109,7 @@ void PipeDrainEvaluator::EvaluatePartialDerivative_(const State& S, const Epetra_MultiVector& mnhMask = *S.GetPtr(mask_key_, tag) ->ViewComponent("cell",false); - double g = norm(S.Get("gravity")); + double g = norm(S.Get("gravity", Tags::DEFAULT)); double pi = 3.14159265359; double mnhArea = pi * manhole_radius_ * manhole_radius_; double mnhPerimeter = 2.0 * pi * manhole_radius_; From 8baa654ebe824f2eee860a24bcfc375195d62732 Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Mon, 20 Nov 2023 17:31:28 -0700 Subject: [PATCH 16/58] Brought ats-regressions-tests repo up to speed --- testing/ats-regression-tests | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testing/ats-regression-tests b/testing/ats-regression-tests index 7742df4ce..b78214037 160000 --- a/testing/ats-regression-tests +++ b/testing/ats-regression-tests @@ -1 +1 @@ -Subproject commit 7742df4ce98333c4d0e975829f7e073a85ec6409 +Subproject commit b78214037f3724200085b6bdcf4bbcdf9d55c3d3 From 00c2a7acc1ef90566308e711816d8a8a910dcde5 Mon Sep 17 00:00:00 2001 From: Naren Vohra Date: Fri, 5 Jan 2024 22:05:29 +0000 Subject: [PATCH 17/58] WIP: SW + pipe coupling pipe drain different domains --- .../sources/pipe_drain_evaluator.cc | 33 ++++++++++++++++--- .../sources/pipe_drain_evaluator.hh | 2 ++ 2 files changed, 30 insertions(+), 5 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index c9dc7fe4b..6388cb62a 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -39,6 +39,10 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : mask_key_ = Keys::readKey(plist_, sw_domain_name_, "manhole locations", "manhole_locations"); dependencies_.insert(KeyTag{mask_key_, tag}); + + + manhole_map_key_ = Keys::readKey(plist_, sw_domain_name_, "manhole map", "manhole_map"); + dependencies_.insert(KeyTag{manhole_map_key_, tag}); } @@ -69,21 +73,35 @@ void PipeDrainEvaluator::Evaluate_(const State& S, int ncells = res.MyLength(); + + // manhole cell map + const Epetra_MultiVector& mhmap_c = *S.GetPtr(manhole_map_key_, tag)->ViewComponent("cell",false); + /* + if(manhole_map_key_.empty()){ + for (int c = 0; c < ncells; c++) { + mhmap_c[0][c] = 1.0; + } + } + */ + if(!pipe_domain_name_.empty()){ const Epetra_MultiVector& pressHead = *S.GetPtr(pressure_head_key_, tag) ->ViewComponent("cell",false); for (int c=0; c!=ncells; ++c) { + + int c1 = mhmap_c[0][c]; + if (pressHead[0][c] < drain_length_) { - res[0][c] = - mnhMask[0][c] * 2.0 / 3.0 * energ_loss_coeff_weir_ * mnhPerimeter * sqrtTwoG * pow(srfcDepth[0][c],3.0/2.0); + res[0][c] = - mnhMask[0][c] * 2.0 / 3.0 * energ_loss_coeff_weir_ * mnhPerimeter * sqrtTwoG * pow(srfcDepth[0][c1],3.0/2.0); } - else if (drain_length_ < pressHead[0][c] && pressHead[0][c] < (drain_length_ + srfcDepth[0][c]) ){ + else if (drain_length_ < pressHead[0][c] && pressHead[0][c] < (drain_length_ + srfcDepth[0][c1]) ){ res[0][c] = - mnhMask[0][c] * energ_loss_coeff_subweir_ * mnhArea * sqrtTwoG - * sqrt(srfcDepth[0][c] + drain_length_ - pressHead[0][c]); + * sqrt(srfcDepth[0][c1] + drain_length_ - pressHead[0][c]); } - else if (pressHead[0][c] > (drain_length_ + srfcDepth[0][c])) { - res[0][c] = mnhMask[0][c] * energ_loss_coeff_orifice_ * mnhArea * sqrtTwoG * sqrt(pressHead[0][c] - drain_length_ - srfcDepth[0][c]); + else if (pressHead[0][c] > (drain_length_ + srfcDepth[0][c1])) { + res[0][c] = mnhMask[0][c] * energ_loss_coeff_orifice_ * mnhArea * sqrtTwoG * sqrt(pressHead[0][c] - drain_length_ - srfcDepth[0][c1]); } } } @@ -117,6 +135,7 @@ void PipeDrainEvaluator::EvaluatePartialDerivative_(const State& S, int ncells = res.MyLength(); + if(!pipe_domain_name_.empty()){ const Epetra_MultiVector& pressHead = *S.GetPtr(pressure_head_key_, tag) @@ -168,6 +187,10 @@ void PipeDrainEvaluator::EvaluatePartialDerivative_(const State& S, } } +void PipeDrainEvaluator::isManhole(AmanziGeometry::Point xc) +{ + double x = xc[0], y = xc[1]; +} } //namespace } //namespace diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh index d8d845088..c35feb58c 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh @@ -37,10 +37,12 @@ class PipeDrainEvaluator : public EvaluatorSecondaryMonotypeCV { virtual void EvaluatePartialDerivative_(const State& S, const Key& wrt_key, const Tag& wrt_tag, const std::vector& result) override; + virtual void isManhole(AmanziGeometry::Point xc); protected: Key surface_depth_key_, pressure_head_key_, mask_key_; Key sw_domain_name_, pipe_domain_name_; + Key manhole_map_key_; double manhole_radius_; // energy losses coefficients at manhole From 4f32fbd977ae908d8b9c7a36632c393f0ac7d686 Mon Sep 17 00:00:00 2001 From: Naren Vohra Date: Wed, 10 Jan 2024 16:26:31 +0000 Subject: [PATCH 18/58] pipe drain evaluator changes from Daniil's branch --- .../sources/pipe_drain_evaluator.cc | 33 ++++++++++++++----- .../sources/pipe_drain_evaluator.hh | 3 +- 2 files changed, 26 insertions(+), 10 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index 6388cb62a..c010ca313 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -37,31 +37,46 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : dependencies_.insert(KeyTag{pressure_head_key_, tag}); } - mask_key_ = Keys::readKey(plist_, sw_domain_name_, "manhole locations", "manhole_locations"); + mask_key_ = Keys::readKey(plist_, pipe_domain_name_, "manhole locations", "manhole_locations"); dependencies_.insert(KeyTag{mask_key_, tag}); - + /* manhole_map_key_ = Keys::readKey(plist_, sw_domain_name_, "manhole map", "manhole_map"); dependencies_.insert(KeyTag{manhole_map_key_, tag}); - + */ + } Teuchos::RCP PipeDrainEvaluator::Clone() const { + return Teuchos::rcp(new PipeDrainEvaluator(*this)); } +void PipeDrainEvaluator::EnsureCompatibility_ToDeps_(State& S, const CompositeVectorSpace& fac) +{ + for (const auto& dep : dependencies_) { + auto domain = Keys::getDomain(dep.first); + if (pipe_domain_name_ == domain) { + auto& dep_fac = S.Require(dep.first, dep.second); + dep_fac.Update(fac); + } + } +} void PipeDrainEvaluator::Evaluate_(const State& S, const std::vector& result) { + Tag tag = my_keys_.front().second; Epetra_MultiVector& res = *result[0]->ViewComponent("cell",false); + const Epetra_MultiVector& srfcDepth = *S.GetPtr(surface_depth_key_, tag) ->ViewComponent("cell",false); + const Epetra_MultiVector& mnhMask = *S.GetPtr(mask_key_, tag) ->ViewComponent("cell",false); @@ -75,14 +90,16 @@ void PipeDrainEvaluator::Evaluate_(const State& S, // manhole cell map - const Epetra_MultiVector& mhmap_c = *S.GetPtr(manhole_map_key_, tag)->ViewComponent("cell",false); /* + const Epetra_MultiVector& mhmap_c = *S.GetPtr(manhole_map_key_, tag)->ViewComponent("cell",false); + if(manhole_map_key_.empty()){ for (int c = 0; c < ncells; c++) { mhmap_c[0][c] = 1.0; } } */ + if(!pipe_domain_name_.empty()){ @@ -91,7 +108,8 @@ void PipeDrainEvaluator::Evaluate_(const State& S, for (int c=0; c!=ncells; ++c) { - int c1 = mhmap_c[0][c]; + //int c1 = mhmap_c[0][c]; + int c1 = c; if (pressHead[0][c] < drain_length_) { res[0][c] = - mnhMask[0][c] * 2.0 / 3.0 * energ_loss_coeff_weir_ * mnhPerimeter * sqrtTwoG * pow(srfcDepth[0][c1],3.0/2.0); @@ -117,6 +135,7 @@ void PipeDrainEvaluator::EvaluatePartialDerivative_(const State& S, const Key& wrt_key, const Tag& wrt_tag, const std::vector& result) { + std::cout<<"Good till here? 5"<ViewComponent("cell",false); @@ -187,10 +206,6 @@ void PipeDrainEvaluator::EvaluatePartialDerivative_(const State& S, } } -void PipeDrainEvaluator::isManhole(AmanziGeometry::Point xc) -{ - double x = xc[0], y = xc[1]; -} } //namespace } //namespace diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh index c35feb58c..bd801e06f 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh @@ -12,6 +12,7 @@ #ifndef AMANZI_FLOW_RELATIONS_PIPE_DRAIN_EVALUATOR_ #define AMANZI_FLOW_RELATIONS_PIPE_DRAIN_EVALUATOR_ +#include "EvaluatorSecondary.hh" #include "EvaluatorSecondaryMonotype.hh" #include "Factory.hh" @@ -37,7 +38,7 @@ class PipeDrainEvaluator : public EvaluatorSecondaryMonotypeCV { virtual void EvaluatePartialDerivative_(const State& S, const Key& wrt_key, const Tag& wrt_tag, const std::vector& result) override; - virtual void isManhole(AmanziGeometry::Point xc); + virtual void EnsureCompatibility_ToDeps_(State& S, const CompositeVectorSpace& fac) override; protected: Key surface_depth_key_, pressure_head_key_, mask_key_; From 3fff78c1ad4a02a143158017adf987a31f864340 Mon Sep 17 00:00:00 2001 From: Naren Vohra Date: Wed, 10 Jan 2024 20:36:29 +0000 Subject: [PATCH 19/58] WIP: first results on one sided coupling from SW->pipe --- .../sources/pipe_drain_evaluator.cc | 17 +++++++++++------ .../sources/pipe_drain_evaluator.hh | 1 + 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index c010ca313..e133dbe04 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -26,6 +26,7 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : drain_length_ = plist.get("drain length", 0.478); sw_domain_name_ = plist.get("sw domain name", "surface"); pipe_domain_name_ = plist.get("pipe domain name", ""); + sink_source_coeff_ = plist.get("sink-source coefficient", 1.0); Tag tag = my_keys_.front().second; // my dependencies @@ -40,10 +41,11 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : mask_key_ = Keys::readKey(plist_, pipe_domain_name_, "manhole locations", "manhole_locations"); dependencies_.insert(KeyTag{mask_key_, tag}); - /* - manhole_map_key_ = Keys::readKey(plist_, sw_domain_name_, "manhole map", "manhole_map"); + + manhole_map_key_ = Keys::readKey(plist_, pipe_domain_name_, "manhole map", "manhole_map"); dependencies_.insert(KeyTag{manhole_map_key_, tag}); - */ + + } @@ -90,9 +92,10 @@ void PipeDrainEvaluator::Evaluate_(const State& S, // manhole cell map - /* + const Epetra_MultiVector& mhmap_c = *S.GetPtr(manhole_map_key_, tag)->ViewComponent("cell",false); + /* if(manhole_map_key_.empty()){ for (int c = 0; c < ncells; c++) { mhmap_c[0][c] = 1.0; @@ -108,8 +111,8 @@ void PipeDrainEvaluator::Evaluate_(const State& S, for (int c=0; c!=ncells; ++c) { - //int c1 = mhmap_c[0][c]; - int c1 = c; + int c1 = mhmap_c[0][c]; + //int c1 = c; if (pressHead[0][c] < drain_length_) { res[0][c] = - mnhMask[0][c] * 2.0 / 3.0 * energ_loss_coeff_weir_ * mnhPerimeter * sqrtTwoG * pow(srfcDepth[0][c1],3.0/2.0); @@ -121,11 +124,13 @@ void PipeDrainEvaluator::Evaluate_(const State& S, else if (pressHead[0][c] > (drain_length_ + srfcDepth[0][c1])) { res[0][c] = mnhMask[0][c] * energ_loss_coeff_orifice_ * mnhArea * sqrtTwoG * sqrt(pressHead[0][c] - drain_length_ - srfcDepth[0][c1]); } + res[0][c] *= sink_source_coeff_; } } else { for (int c=0; c!=ncells; ++c) { res[0][c] = - mnhMask[0][c] * 2.0 / 3.0 * energ_loss_coeff_weir_ * mnhPerimeter * sqrtTwoG * pow(srfcDepth[0][c],3.0/2.0); + res[0][c] *= sink_source_coeff_; } } diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh index bd801e06f..37e8a0f45 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh @@ -53,6 +53,7 @@ class PipeDrainEvaluator : public EvaluatorSecondaryMonotypeCV { double energ_loss_coeff_subweir_; // submerged weir double energ_loss_coeff_orifice_; // orifice double drain_length_; // drain conduit length + double sink_source_coeff_; // coefficient that determines sink or source (when using same evaluator file for pipe or surface flow) private: static Utils::RegisteredFactory reg_; From abf199825e35e82c6f42b435e32f9fa909dc095a Mon Sep 17 00:00:00 2001 From: Naren Vohra Date: Thu, 11 Jan 2024 15:26:58 +0000 Subject: [PATCH 20/58] WIP: slight update; SW evaluator still not working but pipe evaluator is working --- .../sources/pipe_drain_evaluator.cc | 54 +++++++++++++------ 1 file changed, 39 insertions(+), 15 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index e133dbe04..332acdc90 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -28,6 +28,7 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : pipe_domain_name_ = plist.get("pipe domain name", ""); sink_source_coeff_ = plist.get("sink-source coefficient", 1.0); Tag tag = my_keys_.front().second; + //std::cout<<"my_keys_.front().first "<(dep.first, dep.second); - dep_fac.Update(fac); - } - } +{ + + if (my_keys_.front().first == "network-source_drain") { + for (const auto& dep : dependencies_) { + auto domain = Keys::getDomain(dep.first); + if (pipe_domain_name_ == domain) { + auto& dep_fac = S.Require(dep.first, dep.second); + dep_fac.Update(fac); + } + } + } + + else if (my_keys_.front().first == "surface-source_drain") { + for (const auto& dep : dependencies_) { + auto domain = Keys::getDomain(dep.first); + if (sw_domain_name_ == domain) { + auto& dep_fac = S.Require(dep.first, dep.second); + dep_fac.Update(fac); + } + } + } } void PipeDrainEvaluator::Evaluate_(const State& S, @@ -92,7 +117,7 @@ void PipeDrainEvaluator::Evaluate_(const State& S, // manhole cell map - + const Epetra_MultiVector& mhmap_c = *S.GetPtr(manhole_map_key_, tag)->ViewComponent("cell",false); /* @@ -140,7 +165,6 @@ void PipeDrainEvaluator::EvaluatePartialDerivative_(const State& S, const Key& wrt_key, const Tag& wrt_tag, const std::vector& result) { - std::cout<<"Good till here? 5"<ViewComponent("cell",false); From 870cab72dfeba96ed6882cd7f393c36885622703 Mon Sep 17 00:00:00 2001 From: Naren Vohra Date: Thu, 11 Jan 2024 22:06:52 +0000 Subject: [PATCH 21/58] WIP: prototype example of two way SW <-> pipe drain coupling working for one manhole; also see working ex coupled_surface_pipe_network_mild_incline.xml in ats-regression-test --- .../sources/pipe_drain_evaluator.cc | 42 ++++++++++--------- .../sources/pipe_drain_evaluator.hh | 1 + 2 files changed, 24 insertions(+), 19 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index 332acdc90..d87d8d103 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -28,7 +28,6 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : pipe_domain_name_ = plist.get("pipe domain name", ""); sink_source_coeff_ = plist.get("sink-source coefficient", 1.0); Tag tag = my_keys_.front().second; - //std::cout<<"my_keys_.front().first "<(surface_depth_key_, tag) ->ViewComponent("cell",false); - - + + const Epetra_MultiVector& mnhMask = *S.GetPtr(mask_key_, tag) ->ViewComponent("cell",false); + double g = norm(S.Get("gravity", Tags::DEFAULT)); double pi = 3.14159265359; @@ -119,35 +125,33 @@ void PipeDrainEvaluator::Evaluate_(const State& S, // manhole cell map const Epetra_MultiVector& mhmap_c = *S.GetPtr(manhole_map_key_, tag)->ViewComponent("cell",false); - - /* - if(manhole_map_key_.empty()){ - for (int c = 0; c < ncells; c++) { - mhmap_c[0][c] = 1.0; - } - } - */ - if(!pipe_domain_name_.empty()){ const Epetra_MultiVector& pressHead = *S.GetPtr(pressure_head_key_, tag) ->ViewComponent("cell",false); - + int c_pipe, c_sw; for (int c=0; c!=ncells; ++c) { int c1 = mhmap_c[0][c]; //int c1 = c; + if (pipe_flag == true) { + c_pipe = c; + c_sw = mhmap_c[0][c]; + } else if (sw_flag == true) { + c_pipe = mhmap_c[0][c]; + c_sw = c; + } - if (pressHead[0][c] < drain_length_) { - res[0][c] = - mnhMask[0][c] * 2.0 / 3.0 * energ_loss_coeff_weir_ * mnhPerimeter * sqrtTwoG * pow(srfcDepth[0][c1],3.0/2.0); + if (pressHead[0][c_pipe] < drain_length_) { + res[0][c] = - mnhMask[0][c] * 2.0 / 3.0 * energ_loss_coeff_weir_ * mnhPerimeter * sqrtTwoG * pow(srfcDepth[0][c_sw],3.0/2.0); } - else if (drain_length_ < pressHead[0][c] && pressHead[0][c] < (drain_length_ + srfcDepth[0][c1]) ){ + else if (drain_length_ < pressHead[0][c_pipe] && pressHead[0][c_pipe] < (drain_length_ + srfcDepth[0][c_sw]) ){ res[0][c] = - mnhMask[0][c] * energ_loss_coeff_subweir_ * mnhArea * sqrtTwoG - * sqrt(srfcDepth[0][c1] + drain_length_ - pressHead[0][c]); + * sqrt(srfcDepth[0][c_sw] + drain_length_ - pressHead[0][c_pipe]); } - else if (pressHead[0][c] > (drain_length_ + srfcDepth[0][c1])) { - res[0][c] = mnhMask[0][c] * energ_loss_coeff_orifice_ * mnhArea * sqrtTwoG * sqrt(pressHead[0][c] - drain_length_ - srfcDepth[0][c1]); + else if (pressHead[0][c_pipe] > (drain_length_ + srfcDepth[0][c_sw])) { + res[0][c] = mnhMask[0][c] * energ_loss_coeff_orifice_ * mnhArea * sqrtTwoG * sqrt(pressHead[0][c_pipe] - drain_length_ - srfcDepth[0][c_sw]); } res[0][c] *= sink_source_coeff_; } diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh index 37e8a0f45..3a7503d89 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh @@ -54,6 +54,7 @@ class PipeDrainEvaluator : public EvaluatorSecondaryMonotypeCV { double energ_loss_coeff_orifice_; // orifice double drain_length_; // drain conduit length double sink_source_coeff_; // coefficient that determines sink or source (when using same evaluator file for pipe or surface flow) + bool pipe_flag, sw_flag; private: static Utils::RegisteredFactory reg_; From 0c476323565d9e3a0cbfe7c7171fad7d67503ad1 Mon Sep 17 00:00:00 2001 From: Naren Vohra Date: Tue, 16 Jan 2024 23:01:57 +0000 Subject: [PATCH 22/58] WIP: first results two way coupling SW<->pipe with MULTIPLE manholes --- .../sources/pipe_drain_evaluator.cc | 62 +++++++++++++++++-- .../sources/pipe_drain_evaluator.hh | 1 + 2 files changed, 59 insertions(+), 4 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index d87d8d103..40d14b2d1 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -7,6 +7,7 @@ 2) hydraulic head in the pipe network (pressure_head_key_) Authors: Giacomo Capodaglio (gcapodaglio@lanl.gov) + Naren Vohra (vohra@lanl.gov) */ #include "pipe_drain_evaluator.hh" @@ -47,8 +48,8 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : manhole_map_key_ = Keys::readKey(plist_, pipe_domain_name_, "manhole map", "manhole_map"); dependencies_.insert(KeyTag{manhole_map_key_, tag}); - } + } else if (my_keys_.front().first == "surface-source_drain") { pipe_flag = false; @@ -60,9 +61,6 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : manhole_map_key_ = Keys::readKey(plist_, sw_domain_name_, "manhole map", "manhole_map"); dependencies_.insert(KeyTag{manhole_map_key_, tag}); } - - - } @@ -73,6 +71,61 @@ PipeDrainEvaluator::Clone() const { return Teuchos::rcp(new PipeDrainEvaluator(*this)); } + +void PipeDrainEvaluator::CreateCellMap(const State& S) +{ + // grab the meshes + Teuchos::RCP pipe_mesh = S.GetMesh(pipe_domain_name_); + Teuchos::RCP surface_mesh = S.GetMesh(sw_domain_name_); + + + // get the manhole locations/ cell maps + Tag tag = my_keys_.front().second; + + Key pipe_mask_key_ = Keys::readKey(plist_, pipe_domain_name_, "manhole locations", "manhole_locations"); + Key sw_mask_key_ = Keys::readKey(plist_, sw_domain_name_, "manhole locations", "manhole_locations"); + + Key pipe_manhole_map_key_ = Keys::readKey(plist_, pipe_domain_name_, "manhole map", "manhole_map"); + Key sw_manhole_map_key_ = Keys::readKey(plist_, sw_domain_name_, "manhole map", "manhole_map"); + + const Epetra_MultiVector& mnhMask_pipe = *S.GetPtr(pipe_mask_key_, tag)->ViewComponent("cell",false); + const Epetra_MultiVector& mnhMask_sw = *S.GetPtr(sw_mask_key_, tag)->ViewComponent("cell",false); + + const Epetra_MultiVector& mnhmap_pipe = *S.GetPtr(pipe_manhole_map_key_, tag)->ViewComponent("cell",false); + const Epetra_MultiVector& mnhmap_sw = *S.GetPtr(sw_manhole_map_key_, tag)->ViewComponent("cell",false); + + + + // loop over mesh cells and create map + + int ncells_pipe = pipe_mesh->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED); + int ncells_sw = surface_mesh->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED); + + for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { + const Amanzi::AmanziGeometry::Point &xc_sw = surface_mesh->cell_centroid(c_sw); + for (int c_pipe = 0; c_pipe < ncells_pipe; ++c_pipe) { + const Amanzi::AmanziGeometry::Point& xc_pipe = pipe_mesh->cell_centroid(c_pipe); + if (mnhMask_pipe[0][c_pipe] == 1 && (std::abs(xc_sw[0] - xc_pipe[0]) < 1.e-12 ) && (std::abs(xc_sw[1] - xc_pipe[1]) < 1.e-12 ) ) { + mnhMask_sw[0][c_sw] = c_pipe; + } + } + } + + + + for (int c_pipe = 0; c_pipe < ncells_pipe; ++c_pipe) { + const Amanzi::AmanziGeometry::Point &xc_pipe = pipe_mesh->cell_centroid(c_pipe); + for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { + const Amanzi::AmanziGeometry::Point& xc_sw = surface_mesh->cell_centroid(c_sw); + if (mnhMask_sw[0][c_sw] == 1 && (std::abs(xc_sw[0] - xc_pipe[0]) < 1.e-12 ) && (std::abs(xc_sw[1] - xc_pipe[1]) < 1.e-12 ) ) { + mnhMask_pipe[0][c_pipe] = c_sw; + } + } + } + +} + + void PipeDrainEvaluator::EnsureCompatibility_ToDeps_(State& S, const CompositeVectorSpace& fac) { @@ -121,6 +174,7 @@ void PipeDrainEvaluator::Evaluate_(const State& S, int ncells = res.MyLength(); + CreateCellMap(S); // manhole cell map diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh index 3a7503d89..cd00fe421 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh @@ -39,6 +39,7 @@ class PipeDrainEvaluator : public EvaluatorSecondaryMonotypeCV { const Key& wrt_key, const Tag& wrt_tag, const std::vector& result) override; virtual void EnsureCompatibility_ToDeps_(State& S, const CompositeVectorSpace& fac) override; + void CreateCellMap(const State& S); protected: Key surface_depth_key_, pressure_head_key_, mask_key_; From 47d83487f4c58c5814285a324f177ddc4b9aaae7 Mon Sep 17 00:00:00 2001 From: Naren Vohra Date: Wed, 17 Jan 2024 02:52:54 +0000 Subject: [PATCH 23/58] WIP: SW<->pipe for multiple manholes --- .../sources/pipe_drain_evaluator.cc | 13 ++++++++++--- .../sources/pipe_drain_evaluator.hh | 1 + 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index 40d14b2d1..a8dd9f504 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -62,6 +62,8 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : dependencies_.insert(KeyTag{manhole_map_key_, tag}); } + cell_map_flag = 0.0; + } @@ -105,7 +107,7 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) const Amanzi::AmanziGeometry::Point &xc_sw = surface_mesh->cell_centroid(c_sw); for (int c_pipe = 0; c_pipe < ncells_pipe; ++c_pipe) { const Amanzi::AmanziGeometry::Point& xc_pipe = pipe_mesh->cell_centroid(c_pipe); - if (mnhMask_pipe[0][c_pipe] == 1 && (std::abs(xc_sw[0] - xc_pipe[0]) < 1.e-12 ) && (std::abs(xc_sw[1] - xc_pipe[1]) < 1.e-12 ) ) { + if (std::abs(mnhMask_sw[0][c_sw] - 1.0) < 1.e-12 && (std::abs(xc_sw[0] - xc_pipe[0]) < 1.e-12 ) && (std::abs(xc_sw[1] - xc_pipe[1]) < 1.e-12 ) ) { mnhMask_sw[0][c_sw] = c_pipe; } } @@ -117,12 +119,14 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) const Amanzi::AmanziGeometry::Point &xc_pipe = pipe_mesh->cell_centroid(c_pipe); for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { const Amanzi::AmanziGeometry::Point& xc_sw = surface_mesh->cell_centroid(c_sw); - if (mnhMask_sw[0][c_sw] == 1 && (std::abs(xc_sw[0] - xc_pipe[0]) < 1.e-12 ) && (std::abs(xc_sw[1] - xc_pipe[1]) < 1.e-12 ) ) { + if (std::abs(mnhMask_sw[0][c_pipe] - 1.0) < 1.e-12 && (std::abs(xc_sw[0] - xc_pipe[0]) < 1.e-12 ) && (std::abs(xc_sw[1] - xc_pipe[1]) < 1.e-12 ) ) { mnhMask_pipe[0][c_pipe] = c_sw; } } } + cell_map_flag += 1.0; + } @@ -148,6 +152,7 @@ void PipeDrainEvaluator::EnsureCompatibility_ToDeps_(State& S, const CompositeVe } } } + } void PipeDrainEvaluator::Evaluate_(const State& S, @@ -174,7 +179,9 @@ void PipeDrainEvaluator::Evaluate_(const State& S, int ncells = res.MyLength(); - CreateCellMap(S); + if (std::abs(cell_map_flag) < 1.0) { + CreateCellMap(S); + } // manhole cell map diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh index cd00fe421..85d290082 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh @@ -56,6 +56,7 @@ class PipeDrainEvaluator : public EvaluatorSecondaryMonotypeCV { double drain_length_; // drain conduit length double sink_source_coeff_; // coefficient that determines sink or source (when using same evaluator file for pipe or surface flow) bool pipe_flag, sw_flag; + double cell_map_flag; private: static Utils::RegisteredFactory reg_; From f8e652672dd854a7bae40d0e30f76d665f860ab1 Mon Sep 17 00:00:00 2001 From: Naren Vohra Date: Wed, 17 Jan 2024 03:29:00 +0000 Subject: [PATCH 24/58] WIP: bug fix --- .../sources/pipe_drain_evaluator.cc | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index a8dd9f504..99c0c5094 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -108,7 +108,7 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) for (int c_pipe = 0; c_pipe < ncells_pipe; ++c_pipe) { const Amanzi::AmanziGeometry::Point& xc_pipe = pipe_mesh->cell_centroid(c_pipe); if (std::abs(mnhMask_sw[0][c_sw] - 1.0) < 1.e-12 && (std::abs(xc_sw[0] - xc_pipe[0]) < 1.e-12 ) && (std::abs(xc_sw[1] - xc_pipe[1]) < 1.e-12 ) ) { - mnhMask_sw[0][c_sw] = c_pipe; + mnhmap_sw[0][c_sw] = c_pipe; } } } @@ -120,13 +120,21 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { const Amanzi::AmanziGeometry::Point& xc_sw = surface_mesh->cell_centroid(c_sw); if (std::abs(mnhMask_sw[0][c_pipe] - 1.0) < 1.e-12 && (std::abs(xc_sw[0] - xc_pipe[0]) < 1.e-12 ) && (std::abs(xc_sw[1] - xc_pipe[1]) < 1.e-12 ) ) { - mnhMask_pipe[0][c_pipe] = c_sw; + mnhmap_pipe[0][c_pipe] = c_sw; } } } cell_map_flag += 1.0; + // output + std::cout<<"pipe map -------------"< Date: Wed, 17 Jan 2024 17:36:42 +0000 Subject: [PATCH 25/58] WIP: new map implementation for SW<->pipe domain; first physically sound results on multiple manholes --- .../sources/pipe_drain_evaluator.cc | 65 ++++++++++++------- .../sources/pipe_drain_evaluator.hh | 2 + 2 files changed, 42 insertions(+), 25 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index 99c0c5094..ed5281f88 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -90,51 +90,57 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) Key pipe_manhole_map_key_ = Keys::readKey(plist_, pipe_domain_name_, "manhole map", "manhole_map"); Key sw_manhole_map_key_ = Keys::readKey(plist_, sw_domain_name_, "manhole map", "manhole_map"); - const Epetra_MultiVector& mnhMask_pipe = *S.GetPtr(pipe_mask_key_, tag)->ViewComponent("cell",false); - const Epetra_MultiVector& mnhMask_sw = *S.GetPtr(sw_mask_key_, tag)->ViewComponent("cell",false); + const Epetra_MultiVector& mnhMask = *S.GetPtr(mask_key_, tag)->ViewComponent("cell",false); - const Epetra_MultiVector& mnhmap_pipe = *S.GetPtr(pipe_manhole_map_key_, tag)->ViewComponent("cell",false); - const Epetra_MultiVector& mnhmap_sw = *S.GetPtr(sw_manhole_map_key_, tag)->ViewComponent("cell",false); - - + const Epetra_MultiVector& mnhmap = *S.GetPtr(manhole_map_key_, tag)->ViewComponent("cell",false); // loop over mesh cells and create map int ncells_pipe = pipe_mesh->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED); int ncells_sw = surface_mesh->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED); - for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { - const Amanzi::AmanziGeometry::Point &xc_sw = surface_mesh->cell_centroid(c_sw); - for (int c_pipe = 0; c_pipe < ncells_pipe; ++c_pipe) { - const Amanzi::AmanziGeometry::Point& xc_pipe = pipe_mesh->cell_centroid(c_pipe); - if (std::abs(mnhMask_sw[0][c_sw] - 1.0) < 1.e-12 && (std::abs(xc_sw[0] - xc_pipe[0]) < 1.e-12 ) && (std::abs(xc_sw[1] - xc_pipe[1]) < 1.e-12 ) ) { - mnhmap_sw[0][c_sw] = c_pipe; - } - } - } - - + // SW map from pipe -> SW domain + if (pipe_flag == true) { + pipe_map.resize(ncells_pipe); for (int c_pipe = 0; c_pipe < ncells_pipe; ++c_pipe) { const Amanzi::AmanziGeometry::Point &xc_pipe = pipe_mesh->cell_centroid(c_pipe); for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { const Amanzi::AmanziGeometry::Point& xc_sw = surface_mesh->cell_centroid(c_sw); - if (std::abs(mnhMask_sw[0][c_pipe] - 1.0) < 1.e-12 && (std::abs(xc_sw[0] - xc_pipe[0]) < 1.e-12 ) && (std::abs(xc_sw[1] - xc_pipe[1]) < 1.e-12 ) ) { - mnhmap_pipe[0][c_pipe] = c_sw; + if ( (std::abs(mnhMask[0][c_pipe] - 1.0) < 1.e-12 ) && (std::abs(xc_sw[0] - xc_pipe[0]) < 1.e-12 ) && (std::abs(xc_sw[1] - xc_pipe[1]) < 1.e-12 ) ) { + mnhmap[0][c_pipe] = c_sw; + pipe_map[c_pipe] = c_sw; + break; } } } - cell_map_flag += 1.0; - // output std::cout<<"pipe map -------------"< pipe domain + if (pipe_flag == false) { + sw_map.resize(ncells_sw); + for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { + const Amanzi::AmanziGeometry::Point &xc_sw = surface_mesh->cell_centroid(c_sw); + for (int c_pipe = 0; c_pipe < ncells_pipe; ++c_pipe) { + const Amanzi::AmanziGeometry::Point& xc_pipe = pipe_mesh->cell_centroid(c_pipe); + if ( (std::abs(mnhMask[0][c_sw] - 1.0) < 1.e-12 ) && (std::abs(xc_sw[0] - xc_pipe[0]) < 1.e-12 ) && (std::abs(xc_sw[1] - xc_pipe[1]) < 1.e-12 ) ) { + mnhmap[0][c_sw] = c_pipe; + sw_map[c_sw] = c_pipe; + break; + } + } + } std::cout<<"sw map -------------"<("sw domain name", "surface"); pipe_domain_name_ = plist.get("pipe domain name", ""); - sink_source_coeff_ = plist.get("sink-source coefficient", 1.0); Tag tag = my_keys_.front().second; // my dependencies @@ -46,24 +45,28 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : mask_key_ = Keys::readKey(plist_, pipe_domain_name_, "manhole locations", "manhole_locations"); dependencies_.insert(KeyTag{mask_key_, tag}); - manhole_map_key_ = Keys::readKey(plist_, pipe_domain_name_, "manhole map", "manhole_map"); - dependencies_.insert(KeyTag{manhole_map_key_, tag}); - + sink_source_coeff_ = -1.0; + std::cout<<"pipe initialize called"<(mask_key_, tag)->ViewComponent("cell",false); - - const Epetra_MultiVector& mnhmap = *S.GetPtr(manhole_map_key_, tag)->ViewComponent("cell",false); - // loop over mesh cells and create map int ncells_pipe = pipe_mesh->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED); @@ -102,45 +96,35 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) // SW map from pipe -> SW domain if (pipe_flag == true) { - pipe_map.resize(ncells_pipe); - for (int c_pipe = 0; c_pipe < ncells_pipe; ++c_pipe) { - const Amanzi::AmanziGeometry::Point &xc_pipe = pipe_mesh->cell_centroid(c_pipe); - for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { - const Amanzi::AmanziGeometry::Point& xc_sw = surface_mesh->cell_centroid(c_sw); - if ( (std::abs(mnhMask[0][c_pipe] - 1.0) < 1.e-12 ) && (std::abs(xc_sw[0] - xc_pipe[0]) < 1.e-12 ) && (std::abs(xc_sw[1] - xc_pipe[1]) < 1.e-12 ) ) { - mnhmap[0][c_pipe] = c_sw; - pipe_map[c_pipe] = c_sw; - break; + pipe_map.resize(ncells_pipe); + for (int c_pipe = 0; c_pipe < ncells_pipe; ++c_pipe) { + const Amanzi::AmanziGeometry::Point &xc_pipe = pipe_mesh->cell_centroid(c_pipe); + for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { + const Amanzi::AmanziGeometry::Point& xc_sw = surface_mesh->cell_centroid(c_sw); + if ( (std::abs(mnhMask[0][c_pipe] - 1.0) < 1.e-12 ) && (std::abs(xc_sw[0] - xc_pipe[0]) < 1.e-12 ) && (std::abs(xc_sw[1] - xc_pipe[1]) < 1.e-12 ) ) { + pipe_map[c_pipe] = c_sw; + break; + } } } - } - - // output - std::cout<<"pipe map -------------"< pipe domain - if (pipe_flag == false) { - sw_map.resize(ncells_sw); - for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { - const Amanzi::AmanziGeometry::Point &xc_sw = surface_mesh->cell_centroid(c_sw); - for (int c_pipe = 0; c_pipe < ncells_pipe; ++c_pipe) { - const Amanzi::AmanziGeometry::Point& xc_pipe = pipe_mesh->cell_centroid(c_pipe); - if ( (std::abs(mnhMask[0][c_sw] - 1.0) < 1.e-12 ) && (std::abs(xc_sw[0] - xc_pipe[0]) < 1.e-12 ) && (std::abs(xc_sw[1] - xc_pipe[1]) < 1.e-12 ) ) { - mnhmap[0][c_sw] = c_pipe; - sw_map[c_sw] = c_pipe; - break; + if (sw_flag == true) { + sw_map.resize(ncells_sw); + for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { + const Amanzi::AmanziGeometry::Point &xc_sw = surface_mesh->cell_centroid(c_sw); + for (int c_pipe = 0; c_pipe < ncells_pipe; ++c_pipe) { + const Amanzi::AmanziGeometry::Point& xc_pipe = pipe_mesh->cell_centroid(c_pipe); + if ( (std::abs(mnhMask[0][c_sw] - 1.0) < 1.e-12 ) && (std::abs(xc_sw[0] - xc_pipe[0]) < 1.e-12 ) && (std::abs(xc_sw[1] - xc_pipe[1]) < 1.e-12 ) ) { + sw_map[c_sw] = c_pipe; + break; + } } - } + } + sw_map_created = true; } - std::cout<<"sw map -------------"<(manhole_map_key_, tag)->ViewComponent("cell",false); if(!pipe_domain_name_.empty()){ const Epetra_MultiVector& pressHead = *S.GetPtr(pressure_head_key_, tag) ->ViewComponent("cell",false); + int c_pipe, c_sw; + for (int c=0; c!=ncells; ++c) { + // use cell map if (pipe_flag == true) { c_pipe = c; - c_sw = mhmap_c[0][c]; c_sw = pipe_map[c]; } else if (sw_flag == true) { - c_pipe = mhmap_c[0][c]; c_sw = c; c_pipe = sw_map[c]; } - // debugging - /* + // debugging (display cell maps) + /* if (std::abs(mnhMask[0][c] - 1.0) < 1.e-12) { std::cout<<"pipe_flag = "<(mask_key_, tag) ->ViewComponent("cell",false); + // surface and pipe bathymetry + const Epetra_MultiVector& srfcB = *S.GetPtr(surface_bathymetry_key_, tag)->ViewComponent("cell", false); + const Epetra_MultiVector& pipeB = *S.GetPtr(pipe_bathymetry_key_, tag)->ViewComponent("cell", false); double g = norm(S.Get("gravity", Tags::DEFAULT)); double pi = 3.14159265359; @@ -205,7 +216,11 @@ void PipeDrainEvaluator::Evaluate_(const State& S, c_sw = c; c_pipe = sw_map[c]; } - + + // pipe drain length using bathymetry + drain_length_ = srfcB[0][c_sw] - pipeB[0][c_pipe]; + std::cout<<"bathymetry surf = "<("energy losses coeff weir", 0.54); energ_loss_coeff_subweir_ = plist.get("energy losses coeff submerged weir", 0.056); energ_loss_coeff_orifice_ = plist.get("energy losses coeff orifice", 0.167); - drain_length_ = plist.get("drain length", 0.478); sw_domain_name_ = plist.get("sw domain name", "surface"); pipe_domain_name_ = plist.get("pipe domain name", ""); Tag tag = my_keys_.front().second; @@ -34,13 +33,8 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : dependencies_.insert(KeyTag{surface_depth_key_, tag}); // bathymetry - surface_bathymetry_key_ = Keys::readKey(plist_, sw_domain_name_, "bathymetry", "bathymetry"); - pipe_bathymetry_key_ = Keys::readKey(plist_, pipe_domain_name_, "bathymetry", "bathymetry"); - - surface_bathymetry_key_ = "surface-bathymetry"; - pipe_bathymetry_key_ = "network-bathymetry"; - - std::cout<<"surface bathymetry = "< pipe_mesh = S.GetMesh(pipe_domain_name_); Teuchos::RCP surface_mesh = S.GetMesh(sw_domain_name_); - // get the manhole locations/ cell maps Tag tag = my_keys_.front().second; @@ -102,7 +94,6 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) int ncells_pipe = pipe_mesh->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED); int ncells_sw = surface_mesh->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED); - // SW map from pipe -> SW domain if (pipe_flag == true) { pipe_map.resize(ncells_pipe); @@ -138,7 +129,6 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) void PipeDrainEvaluator::EnsureCompatibility_ToDeps_(State& S, const CompositeVectorSpace& fac) { - auto domain1 = Keys::getDomain(my_keys_.front().first); if (domain1 == pipe_domain_name_) { for (const auto& dep : dependencies_) { @@ -165,11 +155,9 @@ void PipeDrainEvaluator::EnsureCompatibility_ToDeps_(State& S, const CompositeVe void PipeDrainEvaluator::Evaluate_(const State& S, const std::vector& result) { - Tag tag = my_keys_.front().second; Epetra_MultiVector& res = *result[0]->ViewComponent("cell",false); - const Epetra_MultiVector& srfcDepth = *S.GetPtr(surface_depth_key_, tag) ->ViewComponent("cell",false); @@ -177,7 +165,7 @@ void PipeDrainEvaluator::Evaluate_(const State& S, const Epetra_MultiVector& mnhMask = *S.GetPtr(mask_key_, tag) ->ViewComponent("cell",false); - // surface and pipe bathymetry + // surface and pipe bathymetry needed for pipe drain length const Epetra_MultiVector& srfcB = *S.GetPtr(surface_bathymetry_key_, tag)->ViewComponent("cell", false); const Epetra_MultiVector& pipeB = *S.GetPtr(pipe_bathymetry_key_, tag)->ViewComponent("cell", false); @@ -199,15 +187,14 @@ void PipeDrainEvaluator::Evaluate_(const State& S, sw_map_created = true; } - if(!pipe_domain_name_.empty()){ + double drain_length; + if(!pipe_domain_name_.empty()){ const Epetra_MultiVector& pressHead = *S.GetPtr(pressure_head_key_, tag) ->ViewComponent("cell",false); - int c_pipe, c_sw; for (int c=0; c!=ncells; ++c) { - // use cell map if (pipe_flag == true) { c_pipe = c; @@ -217,19 +204,24 @@ void PipeDrainEvaluator::Evaluate_(const State& S, c_pipe = sw_map[c]; } - // pipe drain length using bathymetry - drain_length_ = srfcB[0][c_sw] - pipeB[0][c_pipe]; - std::cout<<"bathymetry surf = "< (drain_length_ + srfcDepth[0][c_sw])) { - res[0][c] = mnhMask[0][c] * energ_loss_coeff_orifice_ * mnhArea * sqrtTwoG * sqrt(pressHead[0][c_pipe] - drain_length_ - srfcDepth[0][c_sw]); + else if (pressHead[0][c_pipe] > (drain_length + srfcDepth[0][c_sw])) { + res[0][c] = mnhMask[0][c] * energ_loss_coeff_orifice_ * mnhArea * sqrtTwoG * sqrt(pressHead[0][c_pipe] - drain_length - srfcDepth[0][c_sw]); } res[0][c] *= sink_source_coeff_; } @@ -266,39 +258,43 @@ void PipeDrainEvaluator::EvaluatePartialDerivative_(const State& S, int ncells = res.MyLength(); + // surface and pipe bathymetry needed for pipe drain length + const Epetra_MultiVector& srfcB = *S.GetPtr(surface_bathymetry_key_, tag)->ViewComponent("cell", false); + const Epetra_MultiVector& pipeB = *S.GetPtr(pipe_bathymetry_key_, tag)->ViewComponent("cell", false); - if(!pipe_domain_name_.empty()){ + double drain_length = 0.1; // dummy value; FIX ME + if(!pipe_domain_name_.empty()){ const Epetra_MultiVector& pressHead = *S.GetPtr(pressure_head_key_, tag) ->ViewComponent("cell",false); if (wrt_key == surface_depth_key_) { for (int c=0; c!=ncells; ++c) { - if (pressHead[0][c] < drain_length_) { + if (pressHead[0][c] < drain_length) { res[0][c] = - mnhMask[0][c] * energ_loss_coeff_weir_ * mnhPerimeter* sqrtTwoG * sqrt(srfcDepth[0][c]); } - else if (drain_length_ < pressHead[0][c] && pressHead[0][c] < (drain_length_ + srfcDepth[0][c]) ){ + else if (drain_length < pressHead[0][c] && pressHead[0][c] < (drain_length + srfcDepth[0][c]) ){ res[0][c] = - 0.5 * mnhMask[0][c] * energ_loss_coeff_subweir_ * mnhArea * sqrtTwoG - / sqrt(srfcDepth[0][c] + drain_length_ - pressHead[0][c]); + / sqrt(srfcDepth[0][c] + drain_length - pressHead[0][c]); } - else if (pressHead[0][c] > (drain_length_ + srfcDepth[0][c])) { + else if (pressHead[0][c] > (drain_length + srfcDepth[0][c])) { res[0][c] = - 0.5 * mnhMask[0][c] * energ_loss_coeff_orifice_ * mnhArea * sqrtTwoG - / sqrt(pressHead[0][c] - drain_length_ - srfcDepth[0][c]); + / sqrt(pressHead[0][c] - drain_length - srfcDepth[0][c]); } } } else if (wrt_key == pressure_head_key_) { for (int c=0; c!=ncells; ++c) { - if (pressHead[0][c] < drain_length_) { + if (pressHead[0][c] < drain_length) { res[0][c] = 0.0; } - else if (drain_length_ < pressHead[0][c] && pressHead[0][c] < (drain_length_ + srfcDepth[0][c]) ){ + else if (drain_length < pressHead[0][c] && pressHead[0][c] < (drain_length + srfcDepth[0][c]) ){ res[0][c] = 0.5 * mnhMask[0][c] * energ_loss_coeff_subweir_ * mnhArea * sqrtTwoG - / sqrt(srfcDepth[0][c] + drain_length_ - pressHead[0][c]); + / sqrt(srfcDepth[0][c] + drain_length - pressHead[0][c]); } - else if (pressHead[0][c] > (drain_length_ + srfcDepth[0][c])) { + else if (pressHead[0][c] > (drain_length + srfcDepth[0][c])) { res[0][c] = 0.5 * mnhMask[0][c] * energ_loss_coeff_orifice_ * mnhArea * sqrtTwoG - / sqrt(pressHead[0][c] - drain_length_ - srfcDepth[0][c]); + / sqrt(pressHead[0][c] - drain_length - srfcDepth[0][c]); } } } diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh index 9447be1f0..66804c874 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh @@ -54,7 +54,6 @@ class PipeDrainEvaluator : public EvaluatorSecondaryMonotypeCV { double energ_loss_coeff_weir_; // weir double energ_loss_coeff_subweir_; // submerged weir double energ_loss_coeff_orifice_; // orifice - double drain_length_; // drain conduit length double sink_source_coeff_; // coefficient that determines sink or source (when using same evaluator file for pipe or surface flow) bool pipe_flag, sw_flag, pipe_map_created, sw_map_created; double cell_map_flag; From 53edeb003056d5ee4f959df486c94fe9bdd40b39 Mon Sep 17 00:00:00 2001 From: Naren Vohra Date: Wed, 7 Feb 2024 21:32:21 +0000 Subject: [PATCH 30/58] minor changes to address consistency comments --- .../sources/pipe_drain_evaluator.cc | 49 +++++++++---------- .../sources/pipe_drain_evaluator.hh | 5 +- 2 files changed, 26 insertions(+), 28 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index 61b3538be..aae3a5400 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -45,8 +45,8 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : auto domain = Keys::getDomain(my_keys_.front().first); if (domain == pipe_domain_name_) { - pipe_flag = true; - sw_flag = false; + pipe_flag_ = true; + sw_flag_ = false; mask_key_ = Keys::readKey(plist_, pipe_domain_name_, "manhole locations", "manhole_locations"); dependencies_.insert(KeyTag{mask_key_, tag}); @@ -54,8 +54,8 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : sink_source_coeff_ = -1.0; } else { - pipe_flag = false; - sw_flag = true; + pipe_flag_ = false; + sw_flag_ = true; mask_key_ = Keys::readKey(plist_, sw_domain_name_, "manhole locations", "manhole_locations"); dependencies_.insert(KeyTag{mask_key_, tag}); @@ -63,11 +63,11 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : sink_source_coeff_ = 1.0; } // flags to ensure we create the cell maps only once - if (pipe_map_created != true) { - pipe_map_created = false; + if (pipe_map_created_ != true) { + pipe_map_created_ = false; } - if (sw_map_created != true) { - sw_map_created = false; + if (sw_map_created_ != true) { + sw_map_created_ = false; } } @@ -95,35 +95,35 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) int ncells_sw = surface_mesh->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED); // SW map from pipe -> SW domain - if (pipe_flag == true) { - pipe_map.resize(ncells_pipe); + if (pipe_flag_ == true) { + pipe_map_.resize(ncells_pipe); for (int c_pipe = 0; c_pipe < ncells_pipe; ++c_pipe) { const Amanzi::AmanziGeometry::Point &xc_pipe = pipe_mesh->cell_centroid(c_pipe); for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { const Amanzi::AmanziGeometry::Point& xc_sw = surface_mesh->cell_centroid(c_sw); if ( (std::abs(mnhMask[0][c_pipe] - 1.0) < 1.e-12 ) && (std::abs(xc_sw[0] - xc_pipe[0]) < 1.e-12 ) && (std::abs(xc_sw[1] - xc_pipe[1]) < 1.e-12 ) ) { - pipe_map[c_pipe] = c_sw; + pipe_map_[c_pipe] = c_sw; break; } } } - pipe_map_created = true; + pipe_map_created_ = true; } // Pipe map from SW -> pipe domain - if (sw_flag == true) { - sw_map.resize(ncells_sw); + if (sw_flag_ == true) { + sw_map_.resize(ncells_sw); for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { const Amanzi::AmanziGeometry::Point &xc_sw = surface_mesh->cell_centroid(c_sw); for (int c_pipe = 0; c_pipe < ncells_pipe; ++c_pipe) { const Amanzi::AmanziGeometry::Point& xc_pipe = pipe_mesh->cell_centroid(c_pipe); if ( (std::abs(mnhMask[0][c_sw] - 1.0) < 1.e-12 ) && (std::abs(xc_sw[0] - xc_pipe[0]) < 1.e-12 ) && (std::abs(xc_sw[1] - xc_pipe[1]) < 1.e-12 ) ) { - sw_map[c_sw] = c_pipe; + sw_map_[c_sw] = c_pipe; break; } } } - sw_map_created = true; + sw_map_created_ = true; } } @@ -160,7 +160,6 @@ void PipeDrainEvaluator::Evaluate_(const State& S, const Epetra_MultiVector& srfcDepth = *S.GetPtr(surface_depth_key_, tag) ->ViewComponent("cell",false); - const Epetra_MultiVector& mnhMask = *S.GetPtr(mask_key_, tag) ->ViewComponent("cell",false); @@ -178,13 +177,13 @@ void PipeDrainEvaluator::Evaluate_(const State& S, int ncells = res.MyLength(); // generate cell maps - if (pipe_map_created == false){ + if (pipe_map_created_ == false){ CreateCellMap(S); - pipe_map_created = true; + pipe_map_created_ = true; } - if (sw_map_created == false) { + if (sw_map_created_ == false) { CreateCellMap(S); - sw_map_created = true; + sw_map_created_ = true; } double drain_length; @@ -196,12 +195,12 @@ void PipeDrainEvaluator::Evaluate_(const State& S, for (int c=0; c!=ncells; ++c) { // use cell map - if (pipe_flag == true) { + if (pipe_flag_ == true) { c_pipe = c; - c_sw = pipe_map[c]; - } else if (sw_flag == true) { + c_sw = pipe_map_[c]; + } else if (sw_flag_ == true) { c_sw = c; - c_pipe = sw_map[c]; + c_pipe = sw_map_[c]; } // calculate pipe drain length using bathymetry diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh index 66804c874..da07e6092 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh @@ -55,10 +55,9 @@ class PipeDrainEvaluator : public EvaluatorSecondaryMonotypeCV { double energ_loss_coeff_subweir_; // submerged weir double energ_loss_coeff_orifice_; // orifice double sink_source_coeff_; // coefficient that determines sink or source (when using same evaluator file for pipe or surface flow) - bool pipe_flag, sw_flag, pipe_map_created, sw_map_created; - double cell_map_flag; + bool pipe_flag_, sw_flag_, pipe_map_created_, sw_map_created_; - std::vector pipe_map, sw_map; + std::vector pipe_map_, sw_map_; private: static Utils::RegisteredFactory reg_; From 8ed67b99dfd66409808b62dd671aefbc1a433b1d Mon Sep 17 00:00:00 2001 From: Naren Vohra Date: Wed, 7 Feb 2024 22:03:21 +0000 Subject: [PATCH 31/58] updated EvaluatePartialDerivative_ --- .../sources/pipe_drain_evaluator.cc | 91 ++++++++++++------- 1 file changed, 57 insertions(+), 34 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index aae3a5400..f1e3083d5 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -89,8 +89,8 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) Tag tag = my_keys_.front().second; const Epetra_MultiVector& mnhMask = *S.GetPtr(mask_key_, tag)->ViewComponent("cell",false); - // loop over mesh cells and create map + // loop over mesh cells and create map int ncells_pipe = pipe_mesh->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED); int ncells_sw = surface_mesh->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED); @@ -234,7 +234,6 @@ void PipeDrainEvaluator::Evaluate_(const State& S, } -// Needs to be possibly UPDATED void PipeDrainEvaluator::EvaluatePartialDerivative_(const State& S, const Key& wrt_key, const Tag& wrt_tag, const std::vector& result) @@ -261,56 +260,80 @@ void PipeDrainEvaluator::EvaluatePartialDerivative_(const State& S, const Epetra_MultiVector& srfcB = *S.GetPtr(surface_bathymetry_key_, tag)->ViewComponent("cell", false); const Epetra_MultiVector& pipeB = *S.GetPtr(pipe_bathymetry_key_, tag)->ViewComponent("cell", false); - double drain_length = 0.1; // dummy value; FIX ME + double drain_length; if(!pipe_domain_name_.empty()){ const Epetra_MultiVector& pressHead = *S.GetPtr(pressure_head_key_, tag) ->ViewComponent("cell",false); + int c_pipe, c_sw; + if (wrt_key == surface_depth_key_) { - for (int c=0; c!=ncells; ++c) { - if (pressHead[0][c] < drain_length) { - res[0][c] = - mnhMask[0][c] * energ_loss_coeff_weir_ * mnhPerimeter* sqrtTwoG * sqrt(srfcDepth[0][c]); - } - else if (drain_length < pressHead[0][c] && pressHead[0][c] < (drain_length + srfcDepth[0][c]) ){ - res[0][c] = - 0.5 * mnhMask[0][c] * energ_loss_coeff_subweir_ * mnhArea * sqrtTwoG + for (int c=0; c!=ncells; ++c) { + // use cell map + if (pipe_flag_ == true) { + c_pipe = c; + c_sw = pipe_map_[c]; + } else if (sw_flag_ == true) { + c_sw = c; + c_pipe = sw_map_[c]; + } + // calculate pipe drain length using bathymetry + drain_length = srfcB[0][c_sw] - pipeB[0][c_pipe]; + + if (pressHead[0][c_pipe] < drain_length) { + res[0][c] = - mnhMask[0][c] * energ_loss_coeff_weir_ * mnhPerimeter* sqrtTwoG * sqrt(srfcDepth[0][c_sw]); + } + else if (drain_length < pressHead[0][c_pipe] && pressHead[0][c_pipe] < (drain_length + srfcDepth[0][c_sw]) ){ + res[0][c] = - 0.5 * mnhMask[0][c] * energ_loss_coeff_subweir_ * mnhArea * sqrtTwoG / sqrt(srfcDepth[0][c] + drain_length - pressHead[0][c]); - } - else if (pressHead[0][c] > (drain_length + srfcDepth[0][c])) { - res[0][c] = - 0.5 * mnhMask[0][c] * energ_loss_coeff_orifice_ * mnhArea * sqrtTwoG - / sqrt(pressHead[0][c] - drain_length - srfcDepth[0][c]); - } - } + } + else if (pressHead[0][c_pipe] > (drain_length + srfcDepth[0][c_sw])) { + res[0][c] = - 0.5 * mnhMask[0][c] * energ_loss_coeff_orifice_ * mnhArea * sqrtTwoG + / sqrt(pressHead[0][c_pipe] - drain_length - srfcDepth[0][c_sw]); + } + } } else if (wrt_key == pressure_head_key_) { - for (int c=0; c!=ncells; ++c) { - if (pressHead[0][c] < drain_length) { - res[0][c] = 0.0; - } - else if (drain_length < pressHead[0][c] && pressHead[0][c] < (drain_length + srfcDepth[0][c]) ){ - res[0][c] = 0.5 * mnhMask[0][c] * energ_loss_coeff_subweir_ * mnhArea * sqrtTwoG - / sqrt(srfcDepth[0][c] + drain_length - pressHead[0][c]); - } - else if (pressHead[0][c] > (drain_length + srfcDepth[0][c])) { - res[0][c] = 0.5 * mnhMask[0][c] * energ_loss_coeff_orifice_ * mnhArea * sqrtTwoG - / sqrt(pressHead[0][c] - drain_length - srfcDepth[0][c]); - } - } + for (int c=0; c!=ncells; ++c) { + // use cell map + if (pipe_flag_ == true) { + c_pipe = c; + c_sw = pipe_map_[c]; + } else if (sw_flag_ == true) { + c_sw = c; + c_pipe = sw_map_[c]; + } + // calculate pipe drain length using bathymetry + drain_length = srfcB[0][c_sw] - pipeB[0][c_pipe]; + + if (pressHead[0][c_pipe] < drain_length) { + res[0][c] = 0.0; + } + else if (drain_length < pressHead[0][c_pipe] && pressHead[0][c_pipe] < (drain_length + srfcDepth[0][c_sw]) ){ + res[0][c] = 0.5 * mnhMask[0][c] * energ_loss_coeff_subweir_ * mnhArea * sqrtTwoG + / sqrt(srfcDepth[0][c_sw] + drain_length - pressHead[0][c_pipe]); + } + else if (pressHead[0][c_pipe] > (drain_length + srfcDepth[0][c_sw])) { + res[0][c] = 0.5 * mnhMask[0][c] * energ_loss_coeff_orifice_ * mnhArea * sqrtTwoG + / sqrt(pressHead[0][c_pipe] - drain_length - srfcDepth[0][c_sw]); + } + } } else { AMANZI_ASSERT(0); } - } - else { + } + else { if (wrt_key == surface_depth_key_) { - for (int c=0; c!=ncells; ++c) { - res[0][c] = - mnhMask[0][c] * energ_loss_coeff_weir_ * mnhPerimeter* sqrtTwoG * sqrt(srfcDepth[0][c]); - } + for (int c=0; c!=ncells; ++c) { + res[0][c] = - mnhMask[0][c] * energ_loss_coeff_weir_ * mnhPerimeter* sqrtTwoG * sqrt(srfcDepth[0][c]); + } } else { AMANZI_ASSERT(0); } - } + } } From 4d34cf9ce4257ecd7bca22d203d82d66adaa692c Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Wed, 7 Feb 2024 17:24:20 -0700 Subject: [PATCH 32/58] Uopdated testing/ats-regressions-tests --- testing/ats-regression-tests | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testing/ats-regression-tests b/testing/ats-regression-tests index b78214037..4e3497353 160000 --- a/testing/ats-regression-tests +++ b/testing/ats-regression-tests @@ -1 +1 @@ -Subproject commit b78214037f3724200085b6bdcf4bbcdf9d55c3d3 +Subproject commit 4e3497353d39093a63d121b7283eea22427ef8aa From 68e453d3f78ccad4066e338ba2350ebc80935138 Mon Sep 17 00:00:00 2001 From: Naren Vohra Date: Sat, 24 Feb 2024 07:07:49 +0000 Subject: [PATCH 33/58] First results on using different meshes for manhole coupling WITHOUT having the cell centroids coinciding --- .../sources/pipe_drain_evaluator.cc | 37 +++++++++++++------ 1 file changed, 25 insertions(+), 12 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index f1e3083d5..fb57e5c60 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -12,6 +12,7 @@ */ #include "pipe_drain_evaluator.hh" +#include "Geometry.hh" namespace Amanzi { namespace Flow { @@ -98,12 +99,18 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) if (pipe_flag_ == true) { pipe_map_.resize(ncells_pipe); for (int c_pipe = 0; c_pipe < ncells_pipe; ++c_pipe) { - const Amanzi::AmanziGeometry::Point &xc_pipe = pipe_mesh->cell_centroid(c_pipe); - for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { - const Amanzi::AmanziGeometry::Point& xc_sw = surface_mesh->cell_centroid(c_sw); - if ( (std::abs(mnhMask[0][c_pipe] - 1.0) < 1.e-12 ) && (std::abs(xc_sw[0] - xc_pipe[0]) < 1.e-12 ) && (std::abs(xc_sw[1] - xc_pipe[1]) < 1.e-12 ) ) { - pipe_map_[c_pipe] = c_sw; - break; + if (std::abs(mnhMask[0][c_pipe] - 1.0) < 1.e-12) { + const Amanzi::AmanziGeometry::Point &xc_pipe = pipe_mesh->cell_centroid(c_pipe); + for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { + const Amanzi::AmanziGeometry::Point& xc_sw = surface_mesh->cell_centroid(c_sw); + std::vector coords; + surface_mesh->cell_get_coordinates(c_sw, &coords); + + if (AmanziGeometry::point_in_polygon(xc_pipe, coords) == true) { + std::cout<<"Pipe cell = "<("energy losses coeff orifice", 0.167); sw_domain_name_ = plist.get("sw domain name", "surface"); - pipe_domain_name_ = plist.get("pipe domain name", ""); + pipe_domain_name_ = plist.get("pipe domain name", "pipe"); Tag tag = my_keys_.front().second; // my dependencies @@ -54,7 +54,7 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : sink_source_coeff_ = -1.0; } - else { + else if (domain == sw_domain_name_){ pipe_flag_ = false; sw_flag_ = true; @@ -63,6 +63,10 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : sink_source_coeff_ = 1.0; } + else{ + std::cout<<"Unknown domain in pipe drain evaluator"< Date: Thu, 14 Mar 2024 11:43:37 -0600 Subject: [PATCH 37/58] Fixes to comply with master --- .../sources/pipe_drain_evaluator.cc | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index 55bc47ff9..bd08008fc 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -89,19 +89,18 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) const Epetra_MultiVector& mnhMask = *S.GetPtr(mask_key_, tag)->ViewComponent("cell",false); // loop over mesh cells and create map - int ncells_pipe = pipe_mesh->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED); - int ncells_sw = surface_mesh->num_entities(AmanziMesh::CELL, AmanziMesh::Parallel_type::OWNED); + int ncells_pipe = pipe_mesh_->getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::OWNED); + int ncells_sw = surface_mesh->->getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::OWNED); // SW map from pipe -> SW domain if (pipe_flag_ == true) { pipe_map_.resize(ncells_pipe); for (int c_pipe = 0; c_pipe < ncells_pipe; ++c_pipe) { if (std::abs(mnhMask[0][c_pipe] - 1.0) < 1.e-12) { - const Amanzi::AmanziGeometry::Point &xc_pipe = pipe_mesh->cell_centroid(c_pipe); + const Amanzi::AmanziGeometry::Point &xc_pipe = pipe_mesh->getCellCentroid(c_pipe); for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { - const Amanzi::AmanziGeometry::Point& xc_sw = surface_mesh->cell_centroid(c_sw); - std::vector coords; - surface_mesh->cell_get_coordinates(c_sw, &coords); + const Amanzi::AmanziGeometry::Point& xc_sw = surface_mesh->getCellCentroid(c_sw); + auto coords = PKUtils_EntityCoordinates(*c_sw, AmanziMesh::Entity_kind::CELL, *surface_mesh); if (AmanziGeometry::point_in_polygon(xc_pipe, coords) == true) { pipe_map_[c_pipe] = c_sw; @@ -118,11 +117,10 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) sw_map_.resize(ncells_sw); for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { if (std::abs(mnhMask[0][c_sw] - 1.0) < 1.e-12) { - const Amanzi::AmanziGeometry::Point &xc_sw = surface_mesh->cell_centroid(c_sw); + const Amanzi::AmanziGeometry::Point &xc_sw = surface_mesh->getCellCentroid(c_sw); for (int c_pipe = 0; c_pipe < ncells_pipe; ++c_pipe) { - const Amanzi::AmanziGeometry::Point& xc_pipe = pipe_mesh->cell_centroid(c_pipe); - std::vector coords; - pipe_mesh->cell_get_coordinates(c_pipe, &coords); + const Amanzi::AmanziGeometry::Point& xc_pipe = pipe_mesh->getCellCentroid(c_pipe); + auto coords = PKUtils_EntityCoordinates(*c_pipe, AmanziMesh::Entity_kind::CELL, *pipe_mesh); if (AmanziGeometry::point_in_polygon(xc_sw, coords) == true) { sw_map_[c_sw] = c_pipe; From 630834f948d25c6adb3670eff5153153ecf3823c Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Thu, 14 Mar 2024 11:51:36 -0600 Subject: [PATCH 38/58] Making changes to comply with master --- .../sources/pipe_drain_evaluator.cc | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index bd08008fc..0e85b0c92 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -13,6 +13,7 @@ #include "pipe_drain_evaluator.hh" #include "Geometry.hh" +#include "PK_DomainFunctionVolume.hh" namespace Amanzi { namespace Flow { @@ -89,8 +90,8 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) const Epetra_MultiVector& mnhMask = *S.GetPtr(mask_key_, tag)->ViewComponent("cell",false); // loop over mesh cells and create map - int ncells_pipe = pipe_mesh_->getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::OWNED); - int ncells_sw = surface_mesh->->getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::OWNED); + int ncells_pipe = pipe_mesh->getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::OWNED); + int ncells_sw = surface_mesh->getNumEntities(AmanziMesh::Entity_kind::CELL, AmanziMesh::Parallel_kind::OWNED); // SW map from pipe -> SW domain if (pipe_flag_ == true) { @@ -100,7 +101,7 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) const Amanzi::AmanziGeometry::Point &xc_pipe = pipe_mesh->getCellCentroid(c_pipe); for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { const Amanzi::AmanziGeometry::Point& xc_sw = surface_mesh->getCellCentroid(c_sw); - auto coords = PKUtils_EntityCoordinates(*c_sw, AmanziMesh::Entity_kind::CELL, *surface_mesh); + auto coords = PKUtils_EntityCoordinates(c_sw, AmanziMesh::Entity_kind::CELL, *surface_mesh); if (AmanziGeometry::point_in_polygon(xc_pipe, coords) == true) { pipe_map_[c_pipe] = c_sw; @@ -120,7 +121,7 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) const Amanzi::AmanziGeometry::Point &xc_sw = surface_mesh->getCellCentroid(c_sw); for (int c_pipe = 0; c_pipe < ncells_pipe; ++c_pipe) { const Amanzi::AmanziGeometry::Point& xc_pipe = pipe_mesh->getCellCentroid(c_pipe); - auto coords = PKUtils_EntityCoordinates(*c_pipe, AmanziMesh::Entity_kind::CELL, *pipe_mesh); + auto coords = PKUtils_EntityCoordinates(c_pipe, AmanziMesh::Entity_kind::CELL, *pipe_mesh); if (AmanziGeometry::point_in_polygon(xc_sw, coords) == true) { sw_map_[c_sw] = c_pipe; From 1b3f36e877174276e5dc97366739006f14931c40 Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Thu, 14 Mar 2024 11:56:06 -0600 Subject: [PATCH 39/58] Making changes to comply with master --- .../flow/constitutive_relations/sources/pipe_drain_evaluator.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index 0e85b0c92..4c5cc5d3b 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -13,7 +13,7 @@ #include "pipe_drain_evaluator.hh" #include "Geometry.hh" -#include "PK_DomainFunctionVolume.hh" +#include "PK_Utils.hh" namespace Amanzi { namespace Flow { From fa7e35d709a89ca8501de23ebbffbb616e4c341c Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Thu, 14 Mar 2024 12:10:13 -0600 Subject: [PATCH 40/58] Making changes to comply with master --- .../sources/pipe_drain_evaluator.cc | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index 4c5cc5d3b..e7bfe21fb 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -100,8 +100,11 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) if (std::abs(mnhMask[0][c_pipe] - 1.0) < 1.e-12) { const Amanzi::AmanziGeometry::Point &xc_pipe = pipe_mesh->getCellCentroid(c_pipe); for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { - const Amanzi::AmanziGeometry::Point& xc_sw = surface_mesh->getCellCentroid(c_sw); - auto coords = PKUtils_EntityCoordinates(c_sw, AmanziMesh::Entity_kind::CELL, *surface_mesh); + std::vector coords; + auto cnodes = mesh_->getCellNodes(c_sw); + for (int node_sw=0; node_swgetCellCentroid(c_sw); for (int c_pipe = 0; c_pipe < ncells_pipe; ++c_pipe) { - const Amanzi::AmanziGeometry::Point& xc_pipe = pipe_mesh->getCellCentroid(c_pipe); - auto coords = PKUtils_EntityCoordinates(c_pipe, AmanziMesh::Entity_kind::CELL, *pipe_mesh); + std::vector coords; + auto cnodes = mesh_->getCellNodes(c_pipe); + for (int node_pipe=0; node_pipe Date: Thu, 14 Mar 2024 12:12:21 -0600 Subject: [PATCH 41/58] Making changes to comply with master --- .../constitutive_relations/sources/pipe_drain_evaluator.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index e7bfe21fb..e9c79ff12 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -101,7 +101,7 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) const Amanzi::AmanziGeometry::Point &xc_pipe = pipe_mesh->getCellCentroid(c_pipe); for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { std::vector coords; - auto cnodes = mesh_->getCellNodes(c_sw); + auto cnodes = surface_mesh_->getCellNodes(c_sw); for (int node_sw=0; node_swgetCellCentroid(c_sw); for (int c_pipe = 0; c_pipe < ncells_pipe; ++c_pipe) { std::vector coords; - auto cnodes = mesh_->getCellNodes(c_pipe); + auto cnodes = pipe_mesh_->getCellNodes(c_pipe); for (int node_pipe=0; node_pipe Date: Thu, 14 Mar 2024 12:13:45 -0600 Subject: [PATCH 42/58] Making changes to comply with master --- .../constitutive_relations/sources/pipe_drain_evaluator.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index e9c79ff12..df90fedd1 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -101,7 +101,7 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) const Amanzi::AmanziGeometry::Point &xc_pipe = pipe_mesh->getCellCentroid(c_pipe); for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { std::vector coords; - auto cnodes = surface_mesh_->getCellNodes(c_sw); + auto cnodes = surface_mesh->getCellNodes(c_sw); for (int node_sw=0; node_swgetCellCentroid(c_sw); for (int c_pipe = 0; c_pipe < ncells_pipe; ++c_pipe) { std::vector coords; - auto cnodes = pipe_mesh_->getCellNodes(c_pipe); + auto cnodes = pipe_mesh->getCellNodes(c_pipe); for (int node_pipe=0; node_pipe Date: Thu, 14 Mar 2024 13:14:05 -0600 Subject: [PATCH 43/58] Testing something --- src/pks/flow/constitutive_relations/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pks/flow/constitutive_relations/CMakeLists.txt b/src/pks/flow/constitutive_relations/CMakeLists.txt index e531f7887..4ff2a4dea 100644 --- a/src/pks/flow/constitutive_relations/CMakeLists.txt +++ b/src/pks/flow/constitutive_relations/CMakeLists.txt @@ -8,7 +8,7 @@ # collect all sources list(APPEND subdirs elevation overland_conductivity porosity water_content wrm sources) set(ats_flow_relations_src_files "") -set(ats_flow_relations_inc_files "") +set(ats_flow_relations_inc_files "PK_Utils.hh") foreach(lcv IN LISTS subdirs) include_directories(${ATS_SOURCE_DIR}/src/pks/flow/constitutive_relations/${lcv}) From 304ab1eb1b20d22e00db8f9d22ece0337de0067b Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Thu, 14 Mar 2024 13:15:50 -0600 Subject: [PATCH 44/58] Reverted test --- src/pks/flow/constitutive_relations/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pks/flow/constitutive_relations/CMakeLists.txt b/src/pks/flow/constitutive_relations/CMakeLists.txt index 4ff2a4dea..e531f7887 100644 --- a/src/pks/flow/constitutive_relations/CMakeLists.txt +++ b/src/pks/flow/constitutive_relations/CMakeLists.txt @@ -8,7 +8,7 @@ # collect all sources list(APPEND subdirs elevation overland_conductivity porosity water_content wrm sources) set(ats_flow_relations_src_files "") -set(ats_flow_relations_inc_files "PK_Utils.hh") +set(ats_flow_relations_inc_files "") foreach(lcv IN LISTS subdirs) include_directories(${ATS_SOURCE_DIR}/src/pks/flow/constitutive_relations/${lcv}) From a08004cf62ade8844e7872b1cc5acd087842db98 Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Thu, 14 Mar 2024 18:14:37 -0600 Subject: [PATCH 45/58] Added PKs dir linking in CMakeList --- src/pks/flow/constitutive_relations/CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/pks/flow/constitutive_relations/CMakeLists.txt b/src/pks/flow/constitutive_relations/CMakeLists.txt index e531f7887..c6abd6c8a 100644 --- a/src/pks/flow/constitutive_relations/CMakeLists.txt +++ b/src/pks/flow/constitutive_relations/CMakeLists.txt @@ -10,6 +10,8 @@ list(APPEND subdirs elevation overland_conductivity porosity water_content wrm s set(ats_flow_relations_src_files "") set(ats_flow_relations_inc_files "") +include_directories(${PKS_SOURCE_DIR}) + foreach(lcv IN LISTS subdirs) include_directories(${ATS_SOURCE_DIR}/src/pks/flow/constitutive_relations/${lcv}) From c19fcee1e182cd32772de02400f456cdf6d5d2d1 Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Thu, 14 Mar 2024 18:24:00 -0600 Subject: [PATCH 46/58] Added PKS_SOURCE_DIR to CMakeList --- src/pks/flow/constitutive_relations/CMakeLists.txt | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/pks/flow/constitutive_relations/CMakeLists.txt b/src/pks/flow/constitutive_relations/CMakeLists.txt index c6abd6c8a..3c1deba9f 100644 --- a/src/pks/flow/constitutive_relations/CMakeLists.txt +++ b/src/pks/flow/constitutive_relations/CMakeLists.txt @@ -10,8 +10,6 @@ list(APPEND subdirs elevation overland_conductivity porosity water_content wrm s set(ats_flow_relations_src_files "") set(ats_flow_relations_inc_files "") -include_directories(${PKS_SOURCE_DIR}) - foreach(lcv IN LISTS subdirs) include_directories(${ATS_SOURCE_DIR}/src/pks/flow/constitutive_relations/${lcv}) @@ -19,7 +17,7 @@ foreach(lcv IN LISTS subdirs) set(ats_flow_relations_src_files ${ats_flow_relations_src_files} ${subdir_sources}) file(GLOB subdir_incs "./${lcv}/*.hh") - set(ats_flow_relations_inc_files ${ats_flow_relations_inc_files} ${subdir_incs}) + set(ats_flow_relations_inc_files ${ats_flow_relations_inc_files} ${subdir_incs} ${PKS_SOURCE_DIR}) file(GLOB registrations "./${lcv}/*_reg.hh" ) foreach(reg_lcv IN LISTS registrations) From c1485303f03d4bb3d6d23041a0cd555847c74db8 Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Thu, 14 Mar 2024 18:25:47 -0600 Subject: [PATCH 47/58] Testing --- src/pks/flow/constitutive_relations/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pks/flow/constitutive_relations/CMakeLists.txt b/src/pks/flow/constitutive_relations/CMakeLists.txt index 3c1deba9f..02d9d67a9 100644 --- a/src/pks/flow/constitutive_relations/CMakeLists.txt +++ b/src/pks/flow/constitutive_relations/CMakeLists.txt @@ -14,10 +14,10 @@ foreach(lcv IN LISTS subdirs) include_directories(${ATS_SOURCE_DIR}/src/pks/flow/constitutive_relations/${lcv}) file(GLOB subdir_sources "./${lcv}/*.cc") - set(ats_flow_relations_src_files ${ats_flow_relations_src_files} ${subdir_sources}) + set(ats_flow_relations_src_files ${ats_flow_relations_src_files} ${subdir_sources} ${PKS_SOURCE_DIR}) file(GLOB subdir_incs "./${lcv}/*.hh") - set(ats_flow_relations_inc_files ${ats_flow_relations_inc_files} ${subdir_incs} ${PKS_SOURCE_DIR}) + set(ats_flow_relations_inc_files ${ats_flow_relations_inc_files} ${subdir_incs}) file(GLOB registrations "./${lcv}/*_reg.hh" ) foreach(reg_lcv IN LISTS registrations) From a2c9e475063c3ec1ab18d006dc4c9ee84d06cf85 Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Thu, 14 Mar 2024 18:30:47 -0600 Subject: [PATCH 48/58] testing --- src/pks/flow/constitutive_relations/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pks/flow/constitutive_relations/CMakeLists.txt b/src/pks/flow/constitutive_relations/CMakeLists.txt index 02d9d67a9..d7193529f 100644 --- a/src/pks/flow/constitutive_relations/CMakeLists.txt +++ b/src/pks/flow/constitutive_relations/CMakeLists.txt @@ -14,10 +14,10 @@ foreach(lcv IN LISTS subdirs) include_directories(${ATS_SOURCE_DIR}/src/pks/flow/constitutive_relations/${lcv}) file(GLOB subdir_sources "./${lcv}/*.cc") - set(ats_flow_relations_src_files ${ats_flow_relations_src_files} ${subdir_sources} ${PKS_SOURCE_DIR}) + set(ats_flow_relations_src_files ${ats_flow_relations_src_files} ${subdir_sources}) file(GLOB subdir_incs "./${lcv}/*.hh") - set(ats_flow_relations_inc_files ${ats_flow_relations_inc_files} ${subdir_incs}) + set(ats_flow_relations_inc_files ${ats_flow_relations_inc_files} ${subdir_incs} "~/Amanzi/repos/amanzi_master/src/physics/ats/src/pks") file(GLOB registrations "./${lcv}/*_reg.hh" ) foreach(reg_lcv IN LISTS registrations) From b5d4f2408fbf8c41ec78fdacf690e764ae6d0c1e Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Thu, 14 Mar 2024 18:33:56 -0600 Subject: [PATCH 49/58] testing --- src/pks/flow/constitutive_relations/CMakeLists.txt | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/pks/flow/constitutive_relations/CMakeLists.txt b/src/pks/flow/constitutive_relations/CMakeLists.txt index d7193529f..d6aa7bba1 100644 --- a/src/pks/flow/constitutive_relations/CMakeLists.txt +++ b/src/pks/flow/constitutive_relations/CMakeLists.txt @@ -10,6 +10,8 @@ list(APPEND subdirs elevation overland_conductivity porosity water_content wrm s set(ats_flow_relations_src_files "") set(ats_flow_relations_inc_files "") +include_directories(${AMANZI_SOURCE_DIR}/src/pks) + foreach(lcv IN LISTS subdirs) include_directories(${ATS_SOURCE_DIR}/src/pks/flow/constitutive_relations/${lcv}) @@ -17,7 +19,7 @@ foreach(lcv IN LISTS subdirs) set(ats_flow_relations_src_files ${ats_flow_relations_src_files} ${subdir_sources}) file(GLOB subdir_incs "./${lcv}/*.hh") - set(ats_flow_relations_inc_files ${ats_flow_relations_inc_files} ${subdir_incs} "~/Amanzi/repos/amanzi_master/src/physics/ats/src/pks") + set(ats_flow_relations_inc_files ${ats_flow_relations_inc_files} ${subdir_incs}) file(GLOB registrations "./${lcv}/*_reg.hh" ) foreach(reg_lcv IN LISTS registrations) From 5ee3921cc2307bb3b836dd2df897142f281d3014 Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Thu, 14 Mar 2024 18:37:41 -0600 Subject: [PATCH 50/58] testing --- src/pks/flow/constitutive_relations/CMakeLists.txt | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/pks/flow/constitutive_relations/CMakeLists.txt b/src/pks/flow/constitutive_relations/CMakeLists.txt index d6aa7bba1..066cfcb0d 100644 --- a/src/pks/flow/constitutive_relations/CMakeLists.txt +++ b/src/pks/flow/constitutive_relations/CMakeLists.txt @@ -10,8 +10,6 @@ list(APPEND subdirs elevation overland_conductivity porosity water_content wrm s set(ats_flow_relations_src_files "") set(ats_flow_relations_inc_files "") -include_directories(${AMANZI_SOURCE_DIR}/src/pks) - foreach(lcv IN LISTS subdirs) include_directories(${ATS_SOURCE_DIR}/src/pks/flow/constitutive_relations/${lcv}) @@ -19,7 +17,7 @@ foreach(lcv IN LISTS subdirs) set(ats_flow_relations_src_files ${ats_flow_relations_src_files} ${subdir_sources}) file(GLOB subdir_incs "./${lcv}/*.hh") - set(ats_flow_relations_inc_files ${ats_flow_relations_inc_files} ${subdir_incs}) + set(ats_flow_relations_inc_files ${ats_flow_relations_inc_files} ${subdir_incs} ${AMANZI_SOURCE_DIR}/src/pks) file(GLOB registrations "./${lcv}/*_reg.hh" ) foreach(reg_lcv IN LISTS registrations) From 32d981a0b17a3167c470fcba6d50cb1de2fe5b26 Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Thu, 14 Mar 2024 18:38:38 -0600 Subject: [PATCH 51/58] testing --- src/pks/flow/constitutive_relations/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pks/flow/constitutive_relations/CMakeLists.txt b/src/pks/flow/constitutive_relations/CMakeLists.txt index 066cfcb0d..4aafcf8e3 100644 --- a/src/pks/flow/constitutive_relations/CMakeLists.txt +++ b/src/pks/flow/constitutive_relations/CMakeLists.txt @@ -14,10 +14,10 @@ foreach(lcv IN LISTS subdirs) include_directories(${ATS_SOURCE_DIR}/src/pks/flow/constitutive_relations/${lcv}) file(GLOB subdir_sources "./${lcv}/*.cc") - set(ats_flow_relations_src_files ${ats_flow_relations_src_files} ${subdir_sources}) + set(ats_flow_relations_src_files ${ats_flow_relations_src_files} ${subdir_sources} ${AMANZI_SOURCE_DIR}/src/pks) file(GLOB subdir_incs "./${lcv}/*.hh") - set(ats_flow_relations_inc_files ${ats_flow_relations_inc_files} ${subdir_incs} ${AMANZI_SOURCE_DIR}/src/pks) + set(ats_flow_relations_inc_files ${ats_flow_relations_inc_files} ${subdir_incs}) file(GLOB registrations "./${lcv}/*_reg.hh" ) foreach(reg_lcv IN LISTS registrations) From b93c775cc42b4bdc07f6ad580fd82637cd891016 Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Fri, 15 Mar 2024 16:45:54 -0600 Subject: [PATCH 52/58] Finally the linking is OK --- src/pks/flow/constitutive_relations/CMakeLists.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/pks/flow/constitutive_relations/CMakeLists.txt b/src/pks/flow/constitutive_relations/CMakeLists.txt index 4aafcf8e3..a4005067d 100644 --- a/src/pks/flow/constitutive_relations/CMakeLists.txt +++ b/src/pks/flow/constitutive_relations/CMakeLists.txt @@ -12,9 +12,10 @@ set(ats_flow_relations_inc_files "") foreach(lcv IN LISTS subdirs) include_directories(${ATS_SOURCE_DIR}/src/pks/flow/constitutive_relations/${lcv}) + include_directories(${AMANZI_SRC_DIR}/src/pks) file(GLOB subdir_sources "./${lcv}/*.cc") - set(ats_flow_relations_src_files ${ats_flow_relations_src_files} ${subdir_sources} ${AMANZI_SOURCE_DIR}/src/pks) + set(ats_flow_relations_src_files ${ats_flow_relations_src_files} ${subdir_sources}) file(GLOB subdir_incs "./${lcv}/*.hh") set(ats_flow_relations_inc_files ${ats_flow_relations_inc_files} ${subdir_incs}) From 854c00984ebef81b649b1419f8153c9e452c7991 Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Fri, 15 Mar 2024 17:24:41 -0600 Subject: [PATCH 53/58] Let's see this time --- src/pks/flow/constitutive_relations/CMakeLists.txt | 1 - .../constitutive_relations/sources/pipe_drain_evaluator.cc | 5 ++--- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/pks/flow/constitutive_relations/CMakeLists.txt b/src/pks/flow/constitutive_relations/CMakeLists.txt index a4005067d..baf2451ca 100644 --- a/src/pks/flow/constitutive_relations/CMakeLists.txt +++ b/src/pks/flow/constitutive_relations/CMakeLists.txt @@ -12,7 +12,6 @@ set(ats_flow_relations_inc_files "") foreach(lcv IN LISTS subdirs) include_directories(${ATS_SOURCE_DIR}/src/pks/flow/constitutive_relations/${lcv}) - include_directories(${AMANZI_SRC_DIR}/src/pks) file(GLOB subdir_sources "./${lcv}/*.cc") set(ats_flow_relations_src_files ${ats_flow_relations_src_files} ${subdir_sources}) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index df90fedd1..74471e08f 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -13,7 +13,6 @@ #include "pipe_drain_evaluator.hh" #include "Geometry.hh" -#include "PK_Utils.hh" namespace Amanzi { namespace Flow { @@ -103,7 +102,7 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) std::vector coords; auto cnodes = surface_mesh->getCellNodes(c_sw); for (int node_sw=0; node_sw coords; auto cnodes = pipe_mesh->getCellNodes(c_pipe); for (int node_pipe=0; node_pipe Date: Fri, 15 Mar 2024 18:39:51 -0600 Subject: [PATCH 54/58] Now --- .../constitutive_relations/sources/pipe_drain_evaluator.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index 74471e08f..074b97585 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -102,7 +102,7 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) std::vector coords; auto cnodes = surface_mesh->getCellNodes(c_sw); for (int node_sw=0; node_sw coords; auto cnodes = pipe_mesh->getCellNodes(c_pipe); for (int node_pipe=0; node_pipe Date: Fri, 15 Mar 2024 19:38:59 -0600 Subject: [PATCH 55/58] Now --- .../constitutive_relations/sources/pipe_drain_evaluator.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index 074b97585..4c17a519b 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -102,7 +102,7 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) std::vector coords; auto cnodes = surface_mesh->getCellNodes(c_sw); for (int node_sw=0; node_swgetNodeCoordinate(node_sw); } if (AmanziGeometry::point_in_polygon(xc_pipe, coords) == true) { @@ -125,7 +125,7 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) std::vector coords; auto cnodes = pipe_mesh->getCellNodes(c_pipe); for (int node_pipe=0; node_pipegetNodeCoordinate(node_pipe); } if (AmanziGeometry::point_in_polygon(xc_sw, coords) == true) { From 4ce88e5e5cbf45f1d91a0de720e309cc1fefbc01 Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Fri, 15 Mar 2024 20:18:45 -0600 Subject: [PATCH 56/58] Now --- .../constitutive_relations/sources/pipe_drain_evaluator.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index 4c17a519b..ef8d35ebc 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -102,7 +102,7 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) std::vector coords; auto cnodes = surface_mesh->getCellNodes(c_sw); for (int node_sw=0; node_swgetNodeCoordinate(node_sw); + coords.push_back(surface_mesh->getNodeCoordinate(node_sw)); } if (AmanziGeometry::point_in_polygon(xc_pipe, coords) == true) { @@ -125,7 +125,7 @@ void PipeDrainEvaluator::CreateCellMap(const State& S) std::vector coords; auto cnodes = pipe_mesh->getCellNodes(c_pipe); for (int node_pipe=0; node_pipegetNodeCoordinate(node_pipe); + coords.push_back(pipe_mesh->getNodeCoordinate(node_pipe)); } if (AmanziGeometry::point_in_polygon(xc_sw, coords) == true) { From 8d02cf77b5339c45ffdaad1cfab4419dd9ce012d Mon Sep 17 00:00:00 2001 From: Giacomo Capodaglio Date: Sat, 16 Mar 2024 00:27:22 -0600 Subject: [PATCH 57/58] ATS builds --- src/executables/main.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/src/executables/main.cc b/src/executables/main.cc index 9479d8385..e90670617 100644 --- a/src/executables/main.cc +++ b/src/executables/main.cc @@ -35,7 +35,6 @@ // registration files #include "ats_registration_files.hh" #include "pks_shallow_water_registration.hh" -#include "numerical_flux_registration.hh" int main(int argc, char* argv[]) From 447181e49d741be6b0aa99a7f37a9407de38ee29 Mon Sep 17 00:00:00 2001 From: Daniil Svyatsky Date: Tue, 19 Mar 2024 20:21:32 -0600 Subject: [PATCH 58/58] modified reading of keys in pipe_drain_evaluator.cc --- .../sources/pipe_drain_evaluator.cc | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc index 55bc47ff9..38dc56d64 100644 --- a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -25,8 +25,10 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : energ_loss_coeff_weir_ = plist.get("energy losses coeff weir", 0.54); energ_loss_coeff_subweir_ = plist.get("energy losses coeff submerged weir", 0.056); energ_loss_coeff_orifice_ = plist.get("energy losses coeff orifice", 0.167); - sw_domain_name_ = plist.get("sw domain name", "surface"); + + sw_domain_name_ = plist.get("surface domain name", "surface"); pipe_domain_name_ = plist.get("pipe domain name", "pipe"); + Tag tag = my_keys_.front().second; // my dependencies @@ -34,8 +36,8 @@ PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) : dependencies_.insert(KeyTag{surface_depth_key_, tag}); // bathymetry - surface_bathymetry_key_ = Keys::getKey(sw_domain_name_, "bathymetry"); - pipe_bathymetry_key_ = Keys::getKey(pipe_domain_name_, "bathymetry"); + surface_bathymetry_key_ = Keys::readKey(plist_, sw_domain_name_, "surface bathymetry", "bathymetry"); + pipe_bathymetry_key_ = Keys::readKey(plist_, pipe_domain_name_, "pipe bathymetry", "bathymetry"); if(!pipe_domain_name_.empty()){ pressure_head_key_ = Keys::readKey(plist_, pipe_domain_name_, "pressure head", "pressure_head");