From 6b9fdeba5d89863e05097d21da8a5a00366b676f Mon Sep 17 00:00:00 2001 From: gcapodag Date: Wed, 3 Apr 2024 12:24:50 -0600 Subject: [PATCH] Add Pipe Drain Evaluator --- src/executables/CMakeLists.txt | 2 + src/executables/main.cc | 1 + .../constitutive_relations/CMakeLists.txt | 207 ++++++++-- .../sources/pipe_drain_evaluator.cc | 360 ++++++++++++++++++ .../sources/pipe_drain_evaluator.hh | 68 ++++ .../sources/pipe_drain_evaluator_reg.hh | 6 + 6 files changed, 606 insertions(+), 38 deletions(-) create mode 100644 src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc create mode 100644 src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh create mode 100644 src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator_reg.hh diff --git a/src/executables/CMakeLists.txt b/src/executables/CMakeLists.txt index 7aea75219..8d2fbafdb 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) @@ -97,6 +98,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 9b4f802a0..e90670617 100644 --- a/src/executables/main.cc +++ b/src/executables/main.cc @@ -34,6 +34,7 @@ // registration files #include "ats_registration_files.hh" +#include "pks_shallow_water_registration.hh" int main(int argc, char* argv[]) diff --git a/src/pks/flow/constitutive_relations/CMakeLists.txt b/src/pks/flow/constitutive_relations/CMakeLists.txt index 0ee2a1947..8d2fbafdb 100644 --- a/src/pks/flow/constitutive_relations/CMakeLists.txt +++ b/src/pks/flow/constitutive_relations/CMakeLists.txt @@ -1,54 +1,185 @@ # -*- mode: cmake -*- + +# # Need this define. Errors from MSTK include files +# # about MPI_COMM_WORLD. --lpritch +# add_definitions("-DMSTK_HAVE_MPI") +# add_definitions("-DDISABLE_PHYSICS") + +# if (WITH_MSTK_2_20rc1_OR_NEWER) +# add_definitions("-DMSTK_2_20rc1_OR_NEWER") +# endif () + +# if (WITH_MSTK_2_21rc1_OR_NEWER) +# add_definitions("-DMSTK_2_21rc1_OR_NEWER") +# endif () + +project(EXECUTABLE) + # # ATS -# Constitutive relations for flow +# Executable # +include_directories(${GEOCHEM_SOURCE_DIR}) +include_directories(${MESH_FACTORY_SOURCE_DIR}) +include_directories(${MESH_LOGICAL_SOURCE_DIR}) +include_directories(${MESH_EXTRACTED_SOURCE_DIR}) +include_directories(${MESH_MSTK_SOURCE_DIR}) +include_directories(${CHEMPK_SOURCE_DIR}) +include_directories(${MPC_TREE_SOURCE_DIR}) +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) +include_directories(${ATS_SOURCE_DIR}/src/constitutive_relations/surface_subsurface_fluxes) +include_directories(${ATS_SOURCE_DIR}/src/constitutive_relations/generic_evaluators) +include_directories(${ATS_SOURCE_DIR}/src/constitutive_relations/column_integrators) +include_directories(${ATS_SOURCE_DIR}/src/pks) +include_directories(${ATS_SOURCE_DIR}/src/pks/mpc) +include_directories(${ATS_SOURCE_DIR}/src/pks/energy) +include_directories(${ATS_SOURCE_DIR}/src/pks/flow) +include_directories(${ATS_SOURCE_DIR}/src/pks/deformation) +include_directories(${ATS_SOURCE_DIR}/src/pks/transport) +include_directories(${ATS_SOURCE_DIR}/src/operators/upwinding) +include_directories(${ATS_SOURCE_DIR}/src/operators/advection) +include_directories(${ATS_SOURCE_DIR}/src/operators/deformation) -# collect all sources - -list(APPEND subdirs elevation overland_conductivity porosity sources water_content wrm) +include_directories(${AMANZI_BINARY_DIR}) # required to pick up amanzi_version.hh +include_directories(${ATS_BINARY_DIR}) -set(ats_flow_relations_src_files "") -set(ats_flow_relations_inc_files "") +include_evaluators_directories(LISTNAME REGISTER_AMANZI_STATE_EVALUATORS_INCLUDES) +include_evaluators_directories(LISTNAME ATS_RELATIONS_REG_INCLUDES) +include_evaluators_directories(LISTNAME ATS_TRANSPORT_REG_INCLUDES) +include_evaluators_directories(LISTNAME ATS_ENERGY_PKS_REG_INCLUDES) +include_evaluators_directories(LISTNAME ATS_ENERGY_RELATIONS_REG_INCLUDES) +include_evaluators_directories(LISTNAME ATS_FLOW_PKS_REG_INCLUDES) +include_evaluators_directories(LISTNAME ATS_FLOW_RELATIONS_REG_INCLUDES) +include_evaluators_directories(LISTNAME ATS_DEFORMATION_REG_INCLUDES) +include_evaluators_directories(LISTNAME ATS_SURFACE_BALANCE_REG_INCLUDES) +include_evaluators_directories(LISTNAME ATS_BGC_REG_INCLUDES) +include_evaluators_directories(LISTNAME ATS_MPC_REG_INCLUDES) +include_evaluators_directories(LISTNAME SED_TRANSPORT_REG_INCLUDES) -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_src_files + ats_mesh_factory.cc + coordinator.cc + ats_driver.cc + ) - file(GLOB subdir_incs "./${lcv}/*.hh") - set(ats_flow_relations_inc_files ${ats_flow_relations_inc_files} ${subdir_incs}) +set(ats_inc_files + ats_mesh_factory.hh + coordinator.hh + ats_driver.hh + ) - file(GLOB registrations "./${lcv}/*_reg.hh" ) - foreach(reg_lcv IN LISTS registrations) - register_abs_evaluator_with_factory(HEADERFILE ${reg_lcv} LISTNAME ATS_FLOW_RELATIONS_REG) - endforeach(reg_lcv) +set(amanzi_link_libs + operators + pks + geochemrxns + geochembase + geochemutil + geochemsolvers + state + whetstone + time_integration + solvers + dbg + data_structures + mesh + mesh_functions + functions + geometry + mesh_factory + output + mesh_mstk + mesh_logical + chemistry_pk + #mpc_tree + shallow_water + transport + ) + +set(ats_link_libs + ats_operators + ats_generic_evals + ats_column_integrator + ats_surf_subsurf + ats_eos + ats_pks + ats_transport + #ats_sed_transport + ats_energy + ats_energy_relations + ats_flow + ats_flow_relations + ats_deform + ats_bgc + ats_surface_balance + ats_mpc + ats_mpc_relations + ) -endforeach(lcv) -set(ats_flow_relations_link_libs +# note, we can be inclusive here, because if they aren't enabled, +# these won't be defined and will result in empty strings. +set(tpl_link_libs + ${ALQUIMIA_LIBRARIES} + ${PFLOTRAN_LIBRARIES} ${Teuchos_LIBRARIES} ${Epetra_LIBRARIES} - error_handling - atk - mesh - data_structures - whetstone - solvers - state + ${Boost_LIBRARIES} + ${PETSC_LIBRARIES} + ${MSTK_LIBRARIES} + ${SILO_LIBRARIES} + ${HYPRE_LIBRARIES} + ${HDF5_LIBRARIES} + ${CLM_LIBRARIES} ) -# make the library -add_amanzi_library(ats_flow_relations - SOURCE ${ats_flow_relations_src_files} - HEADERS ${ats_flow_relations_inc_files} - LINK_LIBS ${ats_flow_relations_link_libs}) - -generate_evaluators_registration_header( - HEADERFILE ats_flow_relations_registration.hh - LISTNAME ATS_FLOW_RELATIONS_REG - INSTALL True - ) - +add_amanzi_library(ats_executable + SOURCE ${ats_src_files} + HEADERS ${ats_inc_files} + LINK_LIBS ${amanzi_link_libs} ${tpl_link_libs}) +if (APPLE AND BUILD_SHARED_LIBS) + set_target_properties(ats_executable PROPERTIES LINK_FLAGS "-Wl,-undefined,dynamic_lookup") +endif() + +if (BUILD_TESTS) + + # Add UnitTest includes + include_directories(${UnitTest_INCLUDE_DIRS}) + include_directories(${EXECUTABLE_SOURCE_DIR}) + + # Copy test subdirectory for out of source builds + if (NOT ("${EXECUTABLE_SOURCE_DIR}" STREQUAL "${EXECUTABLE_BINARY_DIR}")) + file(GLOB DataFiles "${EXECUTABLE_SOURCE_DIR}/test/*.xml" + "${EXECUTABLE_SOURCE_DIR}/test/*.exo" + "${EXECUTABLE_SOURCE_DIR}/test/*.h5" + "${EXECUTABLE_SOURCE_DIR}/test/*.bin") + file(COPY ${DataFiles} DESTINATION ${EXECUTABLE_BINARY_DIR}/test/) + endif() + + # test for ats mesh factory + add_amanzi_test(executable_mesh_factory executable_mesh_factory + KIND int + SOURCE test/Main.cc test/executable_mesh_factory.cc + LINK_LIBS ats_executable ${ats_link_libs} ${UnitTest_LIBRARIES} ${NOX_LIBRARIES} ${HDF5_LIBRARIES}) + add_amanzi_test(executable_mesh_factory_np2 executable_mesh_factory NPROCS 2 KIND uint) + add_amanzi_test(executable_mesh_factory_np4 executable_mesh_factory NPROCS 4 KIND uint) + + # test for coupled water preconditioners + add_amanzi_test(executable_coupled_water executable_coupled_water + KIND int + SOURCE test/Main.cc test/executable_coupled_water.cc + LINK_LIBS ats_executable ${ats_link_libs} ${UnitTest_LIBRARIES} ${NOX_LIBRARIES} ${HDF5_LIBRARIES}) + +endif() + +add_amanzi_executable(ats + SOURCE main.cc + LINK_LIBS ats_executable ${fates_link_libs} ${tpl_link_libs} ${ats_link_libs} ${amanzi_link_libs} + OUTPUT_NAME ats + OUTPUT_DIRECTORY ${ATS_BINARY_DIR}) 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..1c15944ed --- /dev/null +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.cc @@ -0,0 +1,360 @@ +/* -*- mode: c++; indent-tabs-mode: nil -*- */ + +/* + 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 (pressure_head_key_) + 3) pipe drain length that is determined by the bathymetry difference + + Authors: Giacomo Capodaglio (gcapodaglio@lanl.gov) + Naren Vohra (vohra@lanl.gov) +*/ + +#include "pipe_drain_evaluator.hh" +#include "Geometry.hh" + +namespace Amanzi { +namespace Flow { + + +PipeDrainEvaluator::PipeDrainEvaluator(Teuchos::ParameterList& plist) + : EvaluatorSecondaryMonotypeCV(plist) +{ + manhole_radius_ = plist_.get("manhole radius", 0.24); + 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("surface domain name", "surface"); + pipe_domain_name_ = plist.get("pipe domain name", "pipe"); + + Tag tag = my_keys_.front().second; + + // my dependencies + surface_depth_key_ = Keys::readKey(plist_, sw_domain_name_, "ponded depth", "ponded_depth"); + dependencies_.insert(KeyTag{ surface_depth_key_, tag }); + + // 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"); + dependencies_.insert(KeyTag{ pressure_head_key_, tag }); + } + + // figure out if SW or pipe is calling the evaluator + auto domain = Keys::getDomain(my_keys_.front().first); + + if (domain == pipe_domain_name_) { + pipe_flag_ = true; + sw_flag_ = false; + + mask_key_ = Keys::readKey(plist_, pipe_domain_name_, "manhole locations", "manhole_locations"); + dependencies_.insert(KeyTag{ mask_key_, tag }); + + sink_source_coeff_ = -1.0; + } else if (domain == sw_domain_name_) { + pipe_flag_ = false; + sw_flag_ = true; + + mask_key_ = Keys::readKey(plist_, sw_domain_name_, "manhole locations", "manhole_locations"); + dependencies_.insert(KeyTag{ mask_key_, tag }); + + sink_source_coeff_ = 1.0; + } else { + std::cout << "Unknown domain in pipe drain evaluator" << std::endl; + AMANZI_ASSERT(0); + } +} + + +Teuchos::RCP +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; + + 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); + + // 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->getCellCentroid(c_pipe); + for (int c_sw = 0; c_sw < ncells_sw; ++c_sw) { + std::vector coords; + auto cnodes = surface_mesh->getCellNodes(c_sw); + for (int node_sw = 0; node_sw < cnodes.size(); node_sw++) { + coords.push_back(surface_mesh->getNodeCoordinate(node_sw)); + } + + if (AmanziGeometry::point_in_polygon(xc_pipe, coords) == true) { + pipe_map_[c_pipe] = c_sw; + break; + } + } + } + } + pipe_map_created_ = true; + } + + // Pipe map from SW -> pipe domain + if (sw_flag_ == true) { + 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->getCellCentroid(c_sw); + for (int c_pipe = 0; c_pipe < ncells_pipe; ++c_pipe) { + std::vector coords; + auto cnodes = pipe_mesh->getCellNodes(c_pipe); + for (int node_pipe = 0; node_pipe < cnodes.size(); node_pipe++) { + coords.push_back(pipe_mesh->getNodeCoordinate(node_pipe)); + } + + if (AmanziGeometry::point_in_polygon(xc_sw, coords) == true) { + sw_map_[c_sw] = c_pipe; + break; + } + } + } + } + sw_map_created_ = true; + } +} + +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_) { + 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 { + 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, 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); + + // 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); + + 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_; + double sqrtTwoG = sqrt(2.0 * g); + + int ncells = res.MyLength(); + + // generate cell maps + if (pipe_map_created_ == false) { + CreateCellMap(S); + pipe_map_created_ = true; + } + if (sw_map_created_ == false) { + CreateCellMap(S); + sw_map_created_ = true; + } + + 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; + 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]; + + // throw error if bathymetrys are not physical + if (drain_length < 0.0) { + std::cout << "Pipe drain length negative; bathymetrys not physical" << std::endl; + AMANZI_ASSERT(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_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][c_sw] + drain_length - pressHead[0][c_pipe]); + } 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_; + } + } 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_; + } + } +} + +void +PipeDrainEvaluator::EvaluatePartialDerivative_(const State& S, + const Key& wrt_key, + const Tag& wrt_tag, + 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); + + 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_; + double sqrtTwoG = sqrt(2.0 * g); + + 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); + + 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) { + // 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_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) { + // 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 { + 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]); + } + } else { + AMANZI_ASSERT(0); + } + } +} + + +} // namespace Flow +} // namespace Amanzi diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh new file mode 100644 index 000000000..e9141af2d --- /dev/null +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator.hh @@ -0,0 +1,68 @@ +/* -*- 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 "EvaluatorSecondary.hh" +#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; + 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_; + Key sw_domain_name_, pipe_domain_name_; + Key manhole_map_key_; + Key surface_bathymetry_key_, pipe_bathymetry_key_; + + double manhole_radius_; + // 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 + 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_ = false, sw_map_created_ = false; + + std::vector pipe_map_, sw_map_; + + private: + static Utils::RegisteredFactory reg_; +}; + +} // namespace Flow +} // namespace Amanzi + +#endif diff --git a/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator_reg.hh b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator_reg.hh new file mode 100644 index 000000000..a4182e842 --- /dev/null +++ b/src/pks/flow/constitutive_relations/sources/pipe_drain_evaluator_reg.hh @@ -0,0 +1,6 @@ +#include "pipe_drain_evaluator.hh" +namespace Amanzi { +namespace Flow { +Utils::RegisteredFactory PipeDrainEvaluator::reg_("pipe drain"); +} +} // namespace Amanzi