CreateHeatTransportBHEProcess.cpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#include "CreateHeatTransportBHEProcess.h"
 
#include <pybind11/pybind11.h>
 
#include <algorithm>
#include <cmath>
#include <map>
#include <ranges>
#include <vector>
 
#include "BHE/BHETypes.h"
#include "BHE/CreateBHE1PType.h"
#include "BHE/CreateBHECoaxial.h"
#include "BHE/CreateBHEUType.h"
#include "HeatTransportBHEProcess.h"
#include "HeatTransportBHEProcessData.h"
#include "MaterialLib/Utils/MediaCreation.h"
#include "ParameterLib/Utils.h"
#include "ProcessLib/HeatTransportBHE/BHE/MeshUtils.h"
#include "ProcessLib/Output/CreateSecondaryVariables.h"
 
namespace ProcessLib
{
namespace HeatTransportBHE
{
using BHECreatorFunc = std::function<BHE::BHETypes(
    BaseLib::ConfigTree const&,
    std::vector<std::unique_ptr<ParameterLib::ParameterBase>>&,
    std::map<std::string,
             std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&,
    std::vector<MeshLib::Node*> const&)>;
 
std::map<std::string_view, BHECreatorFunc> bheCreators = {
    {"1U",
     [](auto& config, auto& parameters, auto& curves, auto& bhe_nodes)
     {
         return BHE::BHE_1U(BHE::createBHEUType<BHE::BHE_1U>(
             config, parameters, curves, bhe_nodes));
     }},
    {"2U",
     [](auto& config, auto& parameters, auto& curves, auto& bhe_nodes)
     {
         return BHE::BHE_2U(BHE::createBHEUType<BHE::BHE_2U>(
             config, parameters, curves, bhe_nodes));
     }},
    {"CXA",
     [](auto& config, auto& parameters, auto& curves, auto& bhe_nodes)
     {
         return BHE::BHE_CXA(BHE::createBHECoaxial<BHE::BHE_CXA>(
             config, parameters, curves, bhe_nodes));
     }},
    {"CXC",
     [](auto& config, auto& parameters, auto& curves, auto& bhe_nodes)
     {
         return BHE::BHE_CXC(BHE::createBHECoaxial<BHE::BHE_CXC>(
             config, parameters, curves, bhe_nodes));
     }},
    {"1P", [](auto& config, auto& parameters, auto& curves, auto& bhe_nodes)
     {
         return BHE::BHE_1P(BHE::createBHE1PType<BHE::BHE_1P>(
             config, parameters, curves, bhe_nodes));
     }}};
 
void createAndInsertBHE(
    const std::string& bhe_type,
    const std::vector<int>& bhe_ids_of_this_bhe,
    const BaseLib::ConfigTree& bhe_config,
    std::vector<std::unique_ptr<ParameterLib::ParameterBase>>& parameters,
    std::map<std::string,
             std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
        curves,
    BHEMeshData const& bhe_mesh_data,
    std::map<int, BHE::BHETypes>& bhes_map)
{
    auto bhe_creator_it = bheCreators.find(bhe_type);
    if (bhe_creator_it == bheCreators.end())
    {
        OGS_FATAL("Unknown BHE type: {:s}", bhe_type);
    }
    for (auto const& id : bhe_ids_of_this_bhe)
    {
        if (id < 0 || id >= static_cast<int>(bhe_mesh_data.BHE_nodes.size()))
        {
            OGS_FATAL(
                "BHE id {:d} is out of range. The mesh contains {:d} "
                "BHE(s) (valid ids: 0 to {:d}).",
                id, bhe_mesh_data.BHE_nodes.size(),
                static_cast<int>(bhe_mesh_data.BHE_nodes.size()) - 1);
        }
        std::pair<std::map<int, BHE::BHETypes>::iterator, bool> result;
        if (id == bhe_ids_of_this_bhe[0])
        {
            auto const& bhe_nodes =
                bhe_mesh_data.BHE_topology_ordered_nodes[id];
            result = bhes_map.try_emplace(
                id, bhe_creator_it->second(bhe_config, parameters, curves,
                                           bhe_nodes));
        }
        else
        {
            // ConfigTree enforces single-read semantics, so we cannot
            // re-parse the same config block.  Grouped BHEs share all
            // configuration except the borehole geometry, which is rebuilt
            // per-ID from each BHE's own node set.
            auto const& first_bhe =
                bhes_map.find(bhe_ids_of_this_bhe[0])->second;
            auto const& bhe_nodes =
                bhe_mesh_data.BHE_topology_ordered_nodes[id];
            result = bhes_map.try_emplace(
                id,
                std::visit(
                    [&](auto const& bhe) -> BHE::BHETypes
                    {
                        return bhe.withGeometry(
                            bhe.borehole_geometry.rebuildForNodes(bhe_nodes));
                    },
                    first_bhe));
        }
        if (!result.second)
        {
            OGS_FATAL(
                "BHE with id '{:d}' is already present in the list! Check for "
                "duplicate definitions of BHE ids.",
                id);
        }
    }
}
 
std::unique_ptr<Process> createHeatTransportBHEProcess(
    std::string const& name,
    MeshLib::Mesh& mesh,
    std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
    std::vector<ProcessVariable> const& variables,
    std::vector<std::unique_ptr<ParameterLib::ParameterBase>>& parameters,
    unsigned const integration_order,
    BaseLib::ConfigTree const& config,
    std::map<std::string,
             std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
        curves,
    std::map<int, std::shared_ptr<MaterialPropertyLib::Medium>> const& media)
{
    config.checkConfigParameter("type", "HEAT_TRANSPORT_BHE");
 
    DBUG("Create HeatTransportBHE Process.");


    auto const pv_config = config.getConfigSubtree("process_variables");
    std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
        process_variables;
 
    // reading primary variables for each
    // BHE----------------------------------------------------------
    auto range =
        pv_config.getConfigParameterList<std::string>("process_variable");
    std::vector<std::reference_wrapper<ProcessVariable>> per_process_variables;
 
    for (std::string const& pv_name : range)
    {
        if (pv_name != "temperature_soil" &&
            pv_name.find("temperature_BHE") == std::string::npos)
        {
            OGS_FATAL(
                "Found a process variable name '{}'. It should be "
                "'temperature_soil' or 'temperature_BHE_X'",
                pv_name);
        }
        auto variable = std::find_if(variables.cbegin(), variables.cend(),
                                     [&pv_name](ProcessVariable const& v)
                                     { return v.getName() == pv_name; });
 
        if (variable == variables.end())
        {
            OGS_FATAL(
                "Could not find process variable '{:s}' in the provided "
                "variables list for config tag <{:s}>.",
                pv_name, "process_variable");
        }
        DBUG("Found process variable '{:s}' for config tag <{:s}>.",
             variable->getName(), "process_variable");
 
        per_process_variables.emplace_back(
            const_cast<ProcessVariable&>(*variable));
    }
    process_variables.push_back(std::move(per_process_variables));
    // end of reading primary variables for each
    // BHE----------------------------------------------------------

    // reading BHE parameters --------------------------------------------------
 
    auto bhe_mesh_data = getBHEDataInMesh(mesh);
 
    auto const& bhe_configs =
        config.getConfigSubtree("borehole_heat_exchangers");
 
    auto const using_server_communication =
        config.getConfigParameter<bool>("use_server_communication", false);
 
    auto const mass_lumping =
        config.getConfigParameter<bool>("mass_lumping", false);
 
    auto const using_algebraic_bc =
        config.getConfigParameter<bool>("use_algebraic_bc", false);
 
    auto const weighting_factor =
        config.getConfigParameter<float>("weighting_factor", 1000.0);
 
    auto const is_linear =
        config.getConfigParameter<bool>("linear", false);
    if (is_linear)
    {
        if (!using_algebraic_bc)
        {
            WARN(
                "You specified that the process simulated by OGS is linear. "
                "With that optimization the process will be assembled only "
                "once and the non-linear solver will only do iterations per "
                "time step to fulfill the BHE boundary conditions. No other "
                "non-linearities will be resolved and OGS will not detect if "
                "there are any non-linearities. It is your responsibility to "
                "ensure that the assembled equation systems are linear, "
                "indeed! There is no safety net!");
        }
        else
        {
            WARN(
                "You specified that the process simulated by OGS is linear. "
                "With that optimization the process will be assembled only "
                "once and the non-linear solver will do only one iteration per "
                "time step. No non-linearities will be resolved and OGS will "
                "not detect if there are any non-linearities. It is your "
                "responsibility to ensure that the assembled equation systems "
                "are linear, indeed! There is no safety net!");
        }
    }
 
    std::map<int, BHE::BHETypes> bhes_map;
 
    int bhe_iterator = 0;
 
    for (
        auto const& bhe_config :
        bhe_configs.getConfigSubtreeList("borehole_heat_exchanger"))
    {
        auto const bhe_id_string =
            bhe_config.getConfigAttribute<std::string>(
                "id", std::to_string(bhe_iterator));
 
        std::vector<int> bhe_ids_of_this_bhe;
 
        if (bhe_id_string == "*")
        {
            int const size = static_cast<int>(bhe_mesh_data.BHE_mat_IDs.size());
            bhe_ids_of_this_bhe.resize(size);
            std::iota(bhe_ids_of_this_bhe.begin(), bhe_ids_of_this_bhe.end(),
                      0);
        }
        else
        {
            bhe_ids_of_this_bhe =
                MaterialLib::splitMaterialIdString(bhe_id_string);
        }
 
        // read in the parameters
        const std::string bhe_type =
            bhe_config.getConfigParameter<std::string>("type");
 
        createAndInsertBHE(bhe_type, bhe_ids_of_this_bhe, bhe_config,
                           parameters, curves, bhe_mesh_data, bhes_map);
        bhe_iterator++;
    }
 
    if (static_cast<int>(bhes_map.size()) - 1 != bhes_map.rbegin()->first)
    {
        OGS_FATAL(
            "The maximum given BHE id '{:d}' did not match the number of given "
            "BHE definitions '{:d}'. The BHE ids needs to be defined starting "
            "from 0, so the maximum BHE id needs to be number of BHE "
            "definitions minus 1. After all definitions there are no gaps "
            "allowed between the given ids.",
            bhes_map.rbegin()->first, bhes_map.size());
    }
 
    std::vector<BHE::BHETypes> bhes;
    bhes.reserve(bhes_map.size());
    std::ranges::copy(bhes_map | std::views::values, std::back_inserter(bhes));
    bhe_mesh_data.updateElementSectionIndices(bhes);
    //  end of reading BHE parameters
    //  -------------------------------------------
 
    auto media_map =
        MaterialPropertyLib::createMaterialSpatialDistributionMap(media, mesh);
 
    // find if bhe uses python boundary condition
    auto const using_tespy =
        visit([](auto const& bhe) { return bhe.use_python_bcs; }, bhes[0]);

    BHEInflowPythonBoundaryConditionPythonSideInterface* py_object = nullptr;
    // create a pythonBoundaryCondition object
    if (using_tespy || using_server_communication)
    {
        // Evaluate Python code in scope of main module
        pybind11::object scope =
            pybind11::module::import("__main__").attr("__dict__");
 
        if (!scope.contains("bc_bhe"))
            OGS_FATAL(
                "Function 'bc_bhe' is not defined in the python script file, "
                "or there was no python script file specified.");
 
        py_object =
            scope["bc_bhe"]
                .cast<BHEInflowPythonBoundaryConditionPythonSideInterface*>();
 
        if (py_object == nullptr)
            OGS_FATAL(
                "Not able to access the correct bc pointer from python script "
                "file specified.");
 
        // create BHE network dataframe from Python
        py_object->dataframe_network = py_object->initializeDataContainer();
        if (!py_object->isOverriddenEssential())
        {
            DBUG(
                "Method `initializeDataContainer' not overridden in Python "
                "script.");
        }
        // clear ogs bc_node_id memory in dataframe
        std::get<3>(py_object->dataframe_network).clear();  // ogs_bc_node_id
 
        // here calls the tespyHydroSolver to get the pipe flow velocity in bhe
        // network
        /* for 2U type the flowrate initialization process below causes conflict
        // replace the value in flow velocity Matrix _u
        auto const tespy_flow_rate = std::get<4>(py_object->dataframe_network);
        const std::size_t n_bhe = tespy_flow_rate.size();
        if (bhes.size() != n_bhe)
            OGS_FATAL(
                "The number of BHEs defined in OGS and TESPy are not the "
                "same!");
 
        for (std::size_t idx_bhe = 0; idx_bhe < n_bhe; idx_bhe++)
        {
            // the flow_rate in OGS should be updated from the flow_rate
            // computed by TESPy.
            auto update_flow_rate = [&](auto& bhe) {
                bhe.updateHeatTransferCoefficients(tespy_flow_rate[idx_bhe]);
            };
            visit(update_flow_rate, bhes[idx_bhe]);
        }
        */
    }
 
    HeatTransportBHEProcessData process_data(
        std::move(media_map), std::move(bhes), py_object, using_tespy,
        using_server_communication, mass_lumping,
        {using_algebraic_bc, weighting_factor}, is_linear);
 
    SecondaryVariableCollection secondary_variables;
 
    ProcessLib::createSecondaryVariables(config, secondary_variables);
 
    return std::make_unique<HeatTransportBHEProcess>(
        std::move(name), mesh, std::move(jacobian_assembler), parameters,
        integration_order, std::move(process_variables),
        std::move(process_data), std::move(secondary_variables),
        std::move(bhe_mesh_data));
}
}  // namespace HeatTransportBHE
}  // namespace ProcessLib