Skip to content
Snippets Groups Projects
Commit 9fffe1dd authored by HBShaoUFZ's avatar HBShaoUFZ Committed by Dmitri Naumov
Browse files

[PL] Cleanup BHE process implementation. 

Avoid duplicates.
Remove unused code parts.
Update source code docu.
Use STL algorithms.
Stronger assertions and improve error messages.
etc.
parent d2bccf9f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 516 additions and 659 deletions
Applications/ApplicationsLib/ProjectData.cpp
+ 15
− 8
View file @ 9fffe1dd
......@@ -377,6 +377,13 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
#ifdef OGS_BUILD_PROCESS_HEATTRANSPORTBHE
if (type == "HEAT_TRANSPORT_BHE")
{
if (_mesh_vec[0]->getDimension() != 3)
{
OGS_FATAL(
"HEAT_TRANSPORT_BHE can only work with a 3-dimentional "
"mesh! ");
}
process =
ProcessLib::HeatTransportBHE::createHeatTransportBHEProcess(
*_mesh_vec[0], std::move(jacobian_assembler),
......@@ -469,17 +476,17 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
{
case 2:
process =
ProcessLib::PhaseField::createPhaseFieldProcess<
2>(*_mesh_vec[0], std::move(jacobian_assembler),
_process_variables, _parameters,
integration_order, process_config);
ProcessLib::PhaseField::createPhaseFieldProcess<2>(
*_mesh_vec[0], std::move(jacobian_assembler),
_process_variables, _parameters, integration_order,
process_config);
break;
case 3:
process =
ProcessLib::PhaseField::createPhaseFieldProcess<
3>(*_mesh_vec[0], std::move(jacobian_assembler),
_process_variables, _parameters,
integration_order, process_config);
ProcessLib::PhaseField::createPhaseFieldProcess<3>(
*_mesh_vec[0], std::move(jacobian_assembler),
_process_variables, _parameters, integration_order,
process_config);
break;
}
}
......
ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.cpp
+ 7
− 11
View file @ 9fffe1dd
......@@ -17,8 +17,7 @@
namespace ProcessLib
{
BHEInflowDirichletBoundaryCondition::BHEInflowDirichletBoundaryCondition(
GlobalIndexType global_idx_T_in_top,
GlobalIndexType global_idx_T_out_top,
std::pair<GlobalIndexType, GlobalIndexType>&& in_out_global_indices,
MeshLib::Mesh const& bc_mesh,
std::vector<MeshLib::Node*> const& vec_inflow_bc_nodes,
int const variable_id,
......@@ -28,8 +27,7 @@ BHEInflowDirichletBoundaryCondition::BHEInflowDirichletBoundaryCondition(
: _bc_mesh(bc_mesh), _pt_bhe(pt_bhe)
{
DBUG(
"Found %d nodes for BHE Inflow Dirichlet BCs for the variable %d "
"and "
"Found %d nodes for BHE Inflow Dirichlet BCs for the variable %d and "
"component %d",
vec_inflow_bc_nodes.size(), variable_id, component_id);
......@@ -48,8 +46,8 @@ BHEInflowDirichletBoundaryCondition::BHEInflowDirichletBoundaryCondition(
_bc_values.values.reserve(bc_mesh_subset.getNumberOfNodes());
// that might be slow, but only done once
const auto g_idx_T_in = global_idx_T_in_top;
const auto g_idx_T_out = global_idx_T_out_top;
const auto g_idx_T_in = in_out_global_indices.first;
const auto g_idx_T_out = in_out_global_indices.second;
if (g_idx_T_in >= 0 && g_idx_T_out >= 0)
{
......@@ -99,7 +97,7 @@ void BHEInflowDirichletBoundaryCondition::preTimestep(const double /*t*/,
std::unique_ptr<BHEInflowDirichletBoundaryCondition>
createBHEInflowDirichletBoundaryCondition(
GlobalIndexType global_idx_T_in_top, GlobalIndexType global_idx_T_out_top,
std::pair<GlobalIndexType, GlobalIndexType>&& in_out_global_indices,
MeshLib::Mesh const& bc_mesh,
std::vector<MeshLib::Node*> const& vec_inflow_bc_nodes,
int const variable_id, int const component_id,
......@@ -108,10 +106,8 @@ createBHEInflowDirichletBoundaryCondition(
{
DBUG("Constructing BHEInflowDirichletBoundaryCondition.");
//! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__Dirichlet__parameter}
return std::make_unique<BHEInflowDirichletBoundaryCondition>(
global_idx_T_in_top, global_idx_T_out_top, bc_mesh, vec_inflow_bc_nodes,
std::move(in_out_global_indices), bc_mesh, vec_inflow_bc_nodes,
variable_id, component_id, pt_bhe);
}
} // namespace ProcessLib
} // namespace ProcessLib
\ No newline at end of file
ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.h
+ 3
− 3
View file @ 9fffe1dd
......@@ -21,8 +21,7 @@ class BHEInflowDirichletBoundaryCondition final : public BoundaryCondition
{
public:
BHEInflowDirichletBoundaryCondition(
GlobalIndexType global_idx_T_in_top,
GlobalIndexType global_idx_T_out_top,
std::pair<GlobalIndexType, GlobalIndexType>&& in_out_global_indices,
MeshLib::Mesh const& bc_mesh,
std::vector<MeshLib::Node*> const& vec_inflow_bc_nodes,
int const variable_id,
......@@ -37,6 +36,7 @@ public:
void preTimestep(const double t, const GlobalVector& x) override;
private:
// TODO (haibing) re-organize as the bottom BC data structure
MeshLib::Mesh const& _bc_mesh;
/// Stores the results of the outflow temperatures per boundary node.
......@@ -51,7 +51,7 @@ private:
std::unique_ptr<BHEInflowDirichletBoundaryCondition>
createBHEInflowDirichletBoundaryCondition(
GlobalIndexType global_idx_T_in_top, GlobalIndexType global_idx_T_out_top,
std::pair<GlobalIndexType, GlobalIndexType>&& in_out_global_indices,
MeshLib::Mesh const& bc_mesh,
std::vector<MeshLib::Node*> const& vec_outflow_bc_nodes,
int const variable_id, int const component_id,
......
ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h
+ 77
− 155
View file @ 9fffe1dd
......@@ -35,6 +35,11 @@
#include "ProcessLib/Utils/ProcessUtils.h"
#include "boost/math/constants/constants.hpp"
#include "BoreholeGeometry.h"
#include "GroutParameters.h"
#include "PipeParameters.h"
#include "RefrigerantParameters.h"
namespace ProcessLib
{
namespace HeatTransportBHE
......@@ -76,127 +81,6 @@ enum class BHE_DISCHARGE_TYPE
class BHEAbstract
{
public:
struct BoreholeGeometry
{
/**
* length/depth of the BHE
* unit is m
*/
double length;
/**
* diameter of the BHE
* unit is m
*/
double diameter;
};
struct PipeParameters
{
/**
* radius of the pipline inner side
* unit is m
*/
double r_inner;
/**
* radius of the pipline outer side
* unit is m
*/
double r_outer;
/**
* pipe-in wall thickness
* unit is m
*/
double b_in;
/**
* pipe-out wall thickness
* unit is m
*/
double b_out;
/**
* thermal conductivity of the pipe wall
* unit is kg m sec^-3 K^-1
*/
double lambda_p;
/**
* thermal conductivity of the inner pipe wall
* unit is kg m sec^-3 K^-1
*/
double lambda_p_i;
/**
* thermal conductivity of the outer pipe wall
* unit is kg m sec^-3 K^-1
*/
double lambda_p_o;
};
struct RefrigerantParameters
{
/**
* dynamics viscosity of the refrigerant
* unit is kg m-1 sec-1
*/
double mu_r;
/**
* density of the refrigerant
* unit is kg m-3
*/
double rho_r;
/**
* thermal conductivity of the refrigerant
* unit is kg m sec^-3 K^-1
*/
double lambda_r;
/**
* specific heat capacity of the refrigerant
* unit is m^2 sec^-2 K^-1
*/
double heat_cap_r;
/**
* longitudinal dispersivity of the
* referigerant flow in the pipeline
*/
double alpha_L;
};
struct GroutParameters
{
/**
* density of the grout
* unit is kg m-3
*/
double rho_g;
/**
* porosity of the grout
* unit is [-]
*/
double porosity_g;
/**
* specific heat capacity of the grout
* unit is m^2 sec^-2 K^-1
*/
double heat_cap_g;
/**
* thermal conductivity of the grout
* unit is kg m sec^-3 K^-1
*/
double lambda_g;
};
struct ExternallyDefinedRaRb
{
/**
......@@ -300,15 +184,7 @@ public:
* initialization calcultion,
* need to be overwritten.
*/
virtual void initialize()
{
calcPipeFlowVelocity();
calcRenoldsNum();
calcPrandtlNum();
calcNusseltNum();
calcThermalResistances();
calcHeatTransferCoefficients();
};
virtual void initialize() = 0;
/**
* update all parameters based on the new flow rate
......@@ -317,12 +193,7 @@ public:
virtual void updateFlowRate(double new_flow_rate)
{
Q_r = new_flow_rate;
calcPipeFlowVelocity();
calcRenoldsNum();
calcPrandtlNum();
calcNusseltNum();
calcThermalResistances();
calcHeatTransferCoefficients();
initialize();
};
virtual void updateFlowRateFromCurve(double current_time) = 0;
......@@ -334,28 +205,61 @@ public:
virtual void calcThermalResistances() = 0;
/**
* Nusselt number calculation,
* need to be overwritten.
*/
virtual void calcNusseltNum() = 0;
/**
* Renolds number calculation,
* need to be overwritten.
* flow velocity inside the pipeline
*/
virtual void calcRenoldsNum() = 0;
double calcPipeFlowVelocity(double const& flow_rate,
double const& pipe_diameter)
{
return 4.0 * flow_rate / (PI * pipe_diameter * pipe_diameter);
}
/**
* Prandtl number calculation,
* need to be overwritten.
*/
virtual void calcPrandtlNum() = 0;
double calcPrandtlNumber(double const& viscosity,
double const& heat_capacity,
double const& heat_conductivity)
{
return viscosity * heat_capacity / heat_conductivity;
}
/**
* flow velocity inside the pipeline
* need to be overwritten.
*/
virtual void calcPipeFlowVelocity() = 0;
double calcRenoldsNumber(double const velocity_norm,
double const pipe_diameter,
double const viscosity,
double const density)
{
return velocity_norm * pipe_diameter / (viscosity / density);
};
double calcNusseltNumber(double const pipe_diameter,
double const pipe_length)
{
double tmp_Nu(0.0);
double gamma(0.0), xi(0.0);
if (Re < 2300.0)
{
tmp_Nu = 4.364;
}
else if (Re >= 2300.0 && Re < 10000.0)
{
gamma = (Re - 2300) / (10000 - 2300);
tmp_Nu = (1.0 - gamma) * 4.364;
tmp_Nu +=
gamma *
((0.0308 / 8.0 * 1.0e4 * Pr) /
(1.0 + 12.7 * std::sqrt(0.0308 / 8.0) *
(std::pow(Pr, 2.0 / 3.0) - 1.0)) *
(1.0 + std::pow(pipe_diameter / pipe_length, 2.0 / 3.0)));
}
else if (Re > 10000.0)
{
xi = pow(1.8 * std::log10(Re) - 1.5, -2.0);
tmp_Nu = (xi / 8.0 * Re * Re) /
(1.0 + 12.7 * std::sqrt(xi / 8.0) *
(std::pow(Pr, 2.0 / 3.0) - 1.0)) *
(1.0 + std::pow(pipe_diameter / pipe_length, 2.0 / 3.0));
}
return tmp_Nu;
};
/**
* heat transfer coefficient,
......@@ -409,6 +313,9 @@ public:
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
R_s_matrix) const = 0;
virtual std::vector<std::pair<int, int>> const&
inflowOutflowBcComponentIds() const = 0;
public:
/**
* name of the borehole heat exchanger
......@@ -530,6 +437,21 @@ protected:
double Q_r;
const double PI;
/**
* Reynolds number
*/
double Re;
/**
* Prandtl number
*/
double Pr;
/**
* pipe distance
*/
double omega;
};
} // end of namespace BHE
} // end of namespace HeatTransportBHE
......
ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
+ 27
− 81
View file @ 9fffe1dd
......@@ -11,6 +11,33 @@
using namespace ProcessLib::HeatTransportBHE::BHE;
void BHE_1U::initialize()
{
double tmp_u = calcPipeFlowVelocity(Q_r, 2.0 * pipe_param.r_inner);
_u(0) = tmp_u;
_u(1) = tmp_u;
// calculate Renolds number
Re = calcRenoldsNumber(std::abs(_u(0)),
2.0 * pipe_param.r_inner,
refrigerant_param.mu_r,
refrigerant_param.rho_r);
// calculate Prandtl number
Pr = calcPrandtlNumber(refrigerant_param.mu_r,
refrigerant_param.heat_cap_r,
refrigerant_param.lambda_r);
// calculate Nusselt number
double tmp_Nu =
calcNusseltNumber(2.0 * pipe_param.r_inner, borehole_geometry.length);
_Nu(0) = tmp_Nu;
_Nu(1) = tmp_Nu;
calcThermalResistances();
calcHeatTransferCoefficients();
}
/**
* calculate thermal resistance
*/
......@@ -153,74 +180,6 @@ void BHE_1U::calcThermalResistances()
// -------------------------------------------------------------------------
}
/**
* Nusselt number calculation
*/
void BHE_1U::calcNusseltNum()
{
// see Eq. 32 in Diersch_2011_CG
double tmp_Nu = 0.0;
double gamma, xi;
double d;
double const& L = borehole_geometry.length;
d = 2.0 * pipe_param.r_inner;
if (_Re < 2300.0)
{
tmp_Nu = 4.364;
}
else if (_Re >= 2300.0 && _Re < 10000.0)
{
gamma = (_Re - 2300) / (10000 - 2300);
tmp_Nu = (1.0 - gamma) * 4.364;
tmp_Nu += gamma * ((0.0308 / 8.0 * 1.0e4 * _Pr) /
(1.0 + 12.7 * std::sqrt(0.0308 / 8.0) *
(std::pow(_Pr, 2.0 / 3.0) - 1.0)) *
(1.0 + std::pow(d / L, 2.0 / 3.0)));
}
else if (_Re > 10000.0)
{
xi = pow(1.8 * std::log10(_Re) - 1.5, -2.0);
tmp_Nu = (xi / 8.0 * _Re * _Pr) /
(1.0 + 12.7 * std::sqrt(xi / 8.0) *
(std::pow(_Pr, 2.0 / 3.0) - 1.0)) *
(1.0 + std::pow(d / L, 2.0 / 3.0));
}
_Nu(0) = tmp_Nu;
_Nu(1) = tmp_Nu;
}
/**
* Renolds number calculation
*/
void BHE_1U::calcRenoldsNum()
{
double u_norm, d;
double const& mu_r = refrigerant_param.mu_r;
double const& rho_r = refrigerant_param.rho_r;
u_norm = std::abs(_u(0));
d = 2.0 * pipe_param.r_inner; // inner diameter of the pipeline
_Re = u_norm * d / (mu_r / rho_r);
}
/**
* Prandtl number calculation
*/
void BHE_1U::calcPrandtlNum()
{
double const& mu_r = refrigerant_param.mu_r;
double const& heat_cap_r = refrigerant_param.heat_cap_r;
double const& lambda_r = refrigerant_param.lambda_r;
_Pr = mu_r * heat_cap_r / lambda_r;
}
/**
* calculate heat transfer coefficient
*/
......@@ -232,19 +191,6 @@ void BHE_1U::calcHeatTransferCoefficients()
_PHI_gs = 1.0 / _R_gs;
}
/**
* flow velocity inside the pipeline
*/
void BHE_1U::calcPipeFlowVelocity()
{
double tmp_u;
tmp_u = Q_r / (PI * pipe_param.r_inner * pipe_param.r_inner);
_u(0) = tmp_u;
_u(1) = tmp_u;
}
double BHE_1U::getMassCoeff(std::size_t idx_unknown) const
{
double const& rho_r = refrigerant_param.rho_r;
......
ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
+ 16
− 40
View file @ 9fffe1dd
......@@ -76,6 +76,7 @@ public:
{
_u = Eigen::Vector2d::Zero();
_Nu = Eigen::Vector2d::Zero();
// 1U type of BHE has 2 pipes
Q_r = my_Qr;
......@@ -162,6 +163,8 @@ public:
*/
std::size_t getNumUnknowns() const { return 4; }
void initialize();
void updateFlowRateFromCurve(double current_time)
{
if (use_flowrate_curve)
......@@ -177,26 +180,6 @@ public:
*/
void calcThermalResistances();
/**
* Nusselt number calculation
*/
void calcNusseltNum();
/**
* Renolds number calculation
*/
void calcRenoldsNum();
/**
* Prandtl number calculation
*/
void calcPrandtlNum();
/**
* flow velocity inside the pipeline
*/
void calcPipeFlowVelocity();
/**
* calculate heat transfer coefficient
*/
......@@ -247,6 +230,12 @@ public:
*/
double getTinByTout(double T_out, double current_time);
std::vector<std::pair<int, int>> const& inflowOutflowBcComponentIds()
const override
{
return _inflow_outflow_bc_component_ids;
}
/**
* required by eigen library,
* to make sure the dynamically allocated class has
......@@ -255,6 +244,9 @@ public:
EIGEN_MAKE_ALIGNED_OPERATOR_NEW
private:
std::vector<std::pair<int, int>> const _inflow_outflow_bc_component_ids = {
{0, 1}};
/**
* thermal resistances
*/
......@@ -296,37 +288,21 @@ private:
double _PHI_fig, _PHI_fog, _PHI_gg, _PHI_gs;
/**
* Reynolds number
* specific exchange surfaces S
*/
double _Re;
double S_i, S_o, S_g1, S_gs;
/**
* Prandtl number
* cross section area
*/
double _Pr;
double CSA_i, CSA_o, CSA_g1, CSA_g2;
/**
* Nusselt number
*/
Eigen::Vector2d _Nu;
/**
* flow velocity inside the pipeline
*/
Eigen::Vector2d _u;
/**
* pipe distance
*/
double omega;
/**
* specific exchange surfaces S
*/
double S_i, S_o, S_g1, S_gs;
/**
* cross section area
*/
double CSA_i, CSA_o, CSA_g1, CSA_g2;
};
} // namespace BHE
} // namespace HeatTransportBHE
......
ProcessLib/HeatTransportBHE/BHE/BoreholeGeometry.h 0 → 100644
+ 36
− 0
View file @ 9fffe1dd
/**
* \file
*
* \copyright
* Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*/
#pragma once
namespace ProcessLib
{
namespace HeatTransportBHE
{
namespace BHE
{
struct BoreholeGeometry
{
/**
* length/depth of the BHE
* unit is m
*/
double const length;
/**
* diameter of the BHE
* unit is m
*/
double const diameter;
};
} // namespace BHE
} // namespace HeatTransportBHE
} // namespace ProcessLib
ProcessLib/HeatTransportBHE/BHE/CreateBHE1U.cpp
+ 1
− 1
View file @ 9fffe1dd
......@@ -246,4 +246,4 @@ BHE::BHE_1U* CreateBHE1U(
}
} // namespace BHE
} // namespace HeatTransportBHE
} // namespace ProcessLib
} // namespace ProcessLib
\ No newline at end of file
ProcessLib/HeatTransportBHE/BHE/GroutParameters.h 0 → 100644
+ 48
− 0
View file @ 9fffe1dd
/**
* \file
*
* \copyright
* Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*/
#pragma once
namespace ProcessLib
{
namespace HeatTransportBHE
{
namespace BHE
{
struct GroutParameters
{
/**
* density of the grout
* unit is kg m-3
*/
double const rho_g;
/**
* porosity of the grout
* unit is [-]
*/
double const porosity_g;
/**
* specific heat capacity of the grout
* unit is m^2 sec^-2 K^-1
*/
double const heat_cap_g;
/**
* thermal conductivity of the grout
* unit is kg m sec^-3 K^-1
*/
double const lambda_g;
};
} // namespace BHE
} // namespace HeatTransportBHE
} // namespace ProcessLib
ProcessLib/HeatTransportBHE/BHE/MeshUtils.cpp
+ 39
− 41
View file @ 9fffe1dd
......@@ -18,42 +18,40 @@ namespace ProcessLib
{
namespace HeatTransportBHE
{
void getBHEDataInMesh(
MeshLib::Mesh const& mesh,
std::vector<MeshLib::Element*>& vec_soil_elements,
std::vector<int>& vec_BHE_mat_IDs,
std::vector<std::vector<MeshLib::Element*>>& vec_BHE_elements,
std::vector<MeshLib::Node*>& vec_soil_nodes,
std::vector<std::vector<MeshLib::Node*>>& vec_BHE_nodes)
BHEMeshData getBHEDataInMesh(MeshLib::Mesh const& mesh)
{
// get vectors of matrix elements and BHE elements
vec_soil_elements.reserve(mesh.getNumberOfElements());
BHEMeshData bheMeshData;
// partition all the mesh elements, and copy them into
// two seperate vectors, one with only matrix elements
// and the other only BHE elements
bheMeshData.soil_elements.reserve(mesh.getNumberOfElements());
std::vector<MeshLib::Element*> all_BHE_elements;
for (MeshLib::Element* e : mesh.getElements())
{
// As for the first step, all elements are counted as a soil element
// first. Those elements connected with a BHE will picked up and
// reorganized into a seperate vector at the end of the function.
if (e->getDimension() == mesh.getDimension())
vec_soil_elements.push_back(e);
else if (e->getDimension() == (unsigned int)1)
all_BHE_elements.push_back(e);
}
auto& all_mesh_elements = mesh.getElements();
std::partition_copy(
std::begin(all_mesh_elements),
std::end(all_mesh_elements),
std::back_inserter(all_BHE_elements),
std::back_inserter(bheMeshData.soil_elements),
[](MeshLib::Element* e) { return e->getDimension() == 1; });
// get BHE material IDs
auto opt_material_ids(
mesh.getProperties().getPropertyVector<int>("MaterialIDs"));
auto opt_material_ids = MeshLib::materialIDs(mesh);
if (opt_material_ids == nullptr)
{
OGS_FATAL("Not able to get material IDs! ");
}
for (MeshLib::Element* e : all_BHE_elements)
vec_BHE_mat_IDs.push_back((*opt_material_ids)[e->getID()]);
BaseLib::makeVectorUnique(vec_BHE_mat_IDs);
DBUG("-> found %d BHE material groups", vec_BHE_mat_IDs.size());
bheMeshData.BHE_mat_IDs.push_back((*opt_material_ids)[e->getID()]);
BaseLib::makeVectorUnique(bheMeshData.BHE_mat_IDs);
DBUG("-> found %d BHE material groups", bheMeshData.BHE_mat_IDs.size());
// create a vector of BHE elements for each group
vec_BHE_elements.resize(vec_BHE_mat_IDs.size());
for (unsigned bhe_id = 0; bhe_id < vec_BHE_mat_IDs.size(); bhe_id++)
bheMeshData.BHE_elements.resize(bheMeshData.BHE_mat_IDs.size());
for (unsigned bhe_id = 0; bhe_id < bheMeshData.BHE_mat_IDs.size(); bhe_id++)
{
const auto bhe_mat_id = vec_BHE_mat_IDs[bhe_id];
std::vector<MeshLib::Element*>& vec_elements = vec_BHE_elements[bhe_id];
const auto bhe_mat_id = bheMeshData.BHE_mat_IDs[bhe_id];
std::vector<MeshLib::Element*>& vec_elements =
bheMeshData.BHE_elements[bhe_id];
std::copy_if(all_BHE_elements.begin(), all_BHE_elements.end(),
std::back_inserter(vec_elements),
[&](MeshLib::Element* e) {
......@@ -63,11 +61,11 @@ void getBHEDataInMesh(
}
// get a vector of BHE nodes
vec_BHE_nodes.resize(vec_BHE_mat_IDs.size());
for (unsigned bhe_id = 0; bhe_id < vec_BHE_mat_IDs.size(); bhe_id++)
bheMeshData.BHE_nodes.resize(bheMeshData.BHE_mat_IDs.size());
for (unsigned bhe_id = 0; bhe_id < bheMeshData.BHE_mat_IDs.size(); bhe_id++)
{
std::vector<MeshLib::Node*>& vec_nodes = vec_BHE_nodes[bhe_id];
for (MeshLib::Element* e : vec_BHE_elements[bhe_id])
std::vector<MeshLib::Node*>& vec_nodes = bheMeshData.BHE_nodes[bhe_id];
for (MeshLib::Element* e : bheMeshData.BHE_elements[bhe_id])
{
for (unsigned i = 0; i < e->getNumberOfNodes(); i++)
{
......@@ -82,20 +80,20 @@ void getBHEDataInMesh(
}
// get all the nodes into the pure soil nodes vector
vec_soil_nodes.reserve(mesh.getNumberOfNodes());
bheMeshData.soil_nodes.reserve(mesh.getNumberOfNodes());
for (MeshLib::Node* n : mesh.getNodes())
{
// All elements are counted as a soil element
vec_soil_nodes.push_back(n);
bheMeshData.soil_nodes.push_back(n);
}
// final count of 3 types of elements
// They are
// (i) soil,
// (ii) BHE
// finalLy counting two types of elements
// They are (i) soil, and (ii) BHE type of elements
DBUG("-> found total %d soil elements and %d BHE elements",
vec_soil_elements.size(),
all_BHE_elements.size());
bheMeshData.soil_elements.size(),
bheMeshData.BHE_elements.size());
return bheMeshData;
}
} // end of namespace HeatTransportBHE
} // namespace ProcessLib
} // namespace ProcessLib
\ No newline at end of file
ProcessLib/HeatTransportBHE/BHE/MeshUtils.h
+ 10
− 10
View file @ 9fffe1dd
......@@ -21,22 +21,22 @@ namespace ProcessLib
{
namespace HeatTransportBHE
{
struct BHEMeshData
{
std::vector<MeshLib::Element*> soil_elements;
std::vector<int> BHE_mat_IDs;
std::vector<std::vector<MeshLib::Element*>> BHE_elements;
std::vector<MeshLib::Node*> soil_nodes;
std::vector<std::vector<MeshLib::Node*>> BHE_nodes;
};
/**
* get data about fracture and matrix elements/nodes from a mesh
*
* @param mesh A mesh which includes BHE elements, i.e. 1-dimensional
* elements. It is assumed that elements forming a BHE have a distinct
* material ID.
* @param vec_soil_elements a vector of soil elements
* @param vec_BHE_elements a vector of BHE elements (grouped by BHE IDs)
* @param vec_BHE_nodes a vector of BHE nodes (grouped by BHE IDs)
*/
void getBHEDataInMesh(
MeshLib::Mesh const& mesh,
std::vector<MeshLib::Element*>& vec_soil_elements,
std::vector<int>& vec_BHE_mat_IDs,
std::vector<std::vector<MeshLib::Element*>>& vec_BHE_elements,
std::vector<MeshLib::Node*>& vec_pure_soil_nodes,
std::vector<std::vector<MeshLib::Node*>>& vec_BHE_nodes);
BHEMeshData getBHEDataInMesh(MeshLib::Mesh const& mesh);
} // end of namespace HeatTransportBHE
} // namespace ProcessLib
ProcessLib/HeatTransportBHE/BHE/PipeParameters.h 0 → 100644
+ 67
− 0
View file @ 9fffe1dd
/**
* \file
*
* \copyright
* Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*/
#pragma once
namespace ProcessLib
{
namespace HeatTransportBHE
{
namespace BHE
{
struct PipeParameters
{
/**
* radius of the pipline inner side
* unit is m
*/
double const r_inner;
/**
* radius of the pipline outer side
* unit is m
*/
double const r_outer;
/**
* pipe-in wall thickness
* unit is m
*/
double const b_in;
/**
* pipe-out wall thickness
* unit is m
*/
double const b_out;
/**
* thermal conductivity of the pipe wall
* unit is kg m sec^-3 K^-1
*/
double const lambda_p;
/**
* thermal conductivity of the inner pipe wall
* unit is kg m sec^-3 K^-1
*/
double const lambda_p_i;
/**
* thermal conductivity of the outer pipe wall
* unit is kg m sec^-3 K^-1
*/
double const lambda_p_o;
};
} // namespace BHE
} // namespace HeatTransportBHE
} // namespace ProcessLib
ProcessLib/HeatTransportBHE/BHE/RefrigerantParameters.h 0 → 100644
+ 54
− 0
View file @ 9fffe1dd
/**
* \file
*
* \copyright
* Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*/
#pragma once
namespace ProcessLib
{
namespace HeatTransportBHE
{
namespace BHE
{
struct RefrigerantParameters
{
/**
* dynamics viscosity of the refrigerant
* unit is kg m-1 sec-1
*/
double const mu_r;
/**
* density of the refrigerant
* unit is kg m-3
*/
double const rho_r;
/**
* thermal conductivity of the refrigerant
* unit is kg m sec^-3 K^-1
*/
double const lambda_r;
/**
* specific heat capacity of the refrigerant
* unit is m^2 sec^-2 K^-1
*/
double const heat_cap_r;
/**
* longitudinal dispersivity of the
* referigerant flow in the pipeline
*/
double const alpha_L;
};
} // namespace BHE
} // namespace HeatTransportBHE
} // namespace ProcessLib
ProcessLib/HeatTransportBHE/CMakeLists.txt
+ 1
− 1
View file @ 9fffe1dd
......@@ -5,4 +5,4 @@ APPEND_SOURCE_FILES(SOURCES LocalAssemblers)
add_library(HeatTransportBHE ${SOURCES})
target_link_libraries(HeatTransportBHE PUBLIC ProcessLib)
# include(Tests.cmake)
include(Tests.cmake)
ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp
+ 74
− 221
View file @ 9fffe1dd
......@@ -37,39 +37,36 @@ HeatTransportBHEProcess::HeatTransportBHEProcess(
: Process(mesh, std::move(jacobian_assembler), parameters,
integration_order, std::move(process_variables),
std::move(secondary_variables), std::move(named_function_caller)),
_process_data(std::move(process_data))
_process_data(std::move(process_data)),
_bheMeshData(getBHEDataInMesh(mesh))
{
getBHEDataInMesh(mesh,
_vec_pure_soil_elements,
_vec_BHE_mat_IDs,
_vec_BHE_elements,
_vec_pure_soil_nodes,
_vec_BHE_nodes);
if (_vec_BHE_mat_IDs.size() != _process_data._vec_BHE_property.size())
if (_bheMeshData.BHE_mat_IDs.size() !=
_process_data._vec_BHE_property.size())
{
OGS_FATAL(
"The number of the given BHE properties (%d) are not "
"consistent"
" with the number of BHE groups in a mesh (%d).",
_process_data._vec_BHE_property.size(),
_vec_BHE_mat_IDs.size());
_bheMeshData.BHE_mat_IDs.size());
}
auto material_ids = MeshLib::materialIDs(mesh);
if (material_ids == nullptr)
{
OGS_FATAL("Not able to get material IDs! ");
}
// create a map from a material ID to a BHE ID
auto max_BHE_mat_id =
std::max_element(_vec_BHE_mat_IDs.begin(), _vec_BHE_mat_IDs.end());
const std::size_t nBHEmatIDs = *max_BHE_mat_id + 1;
_process_data._map_materialID_to_BHE_ID.resize(nBHEmatIDs);
for (std::size_t i = 0; i < _vec_BHE_mat_IDs.size(); i++)
std::max_element(material_ids->begin(), material_ids->end());
_process_data._map_materialID_to_BHE_ID.resize(*max_BHE_mat_id + 1);
for (std::size_t i = 0; i < _bheMeshData.BHE_mat_IDs.size(); i++)
{
// by default, it is assumed that the soil compartment takes material ID
// 0 and the BHE take the successive material group.
_process_data._map_materialID_to_BHE_ID[_vec_BHE_mat_IDs[i]] = i;
// fill in the map structure
_process_data._map_materialID_to_BHE_ID[_bheMeshData.BHE_mat_IDs[i]] =
i;
}
MeshLib::PropertyVector<int> const* material_ids(
mesh.getProperties().getPropertyVector<int>("MaterialIDs"));
_process_data._mesh_prop_materialIDs = material_ids;
}
......@@ -84,16 +81,16 @@ void HeatTransportBHEProcess::constructDofTable()
//------------------------------------------------------------
// first all the soil nodes
_mesh_subset_pure_soil_nodes =
std::make_unique<MeshLib::MeshSubset>(_mesh, _vec_pure_soil_nodes);
std::make_unique<MeshLib::MeshSubset>(_mesh, _bheMeshData.soil_nodes);
std::vector<int> vec_n_BHE_unknowns;
vec_n_BHE_unknowns.reserve(_vec_BHE_nodes.size());
vec_n_BHE_unknowns.reserve(_bheMeshData.BHE_nodes.size());
// the BHE nodes need to be cherry-picked from the vector
for (unsigned i = 0; i < _vec_BHE_nodes.size(); i++)
for (unsigned i = 0; i < _bheMeshData.BHE_nodes.size(); i++)
{
_mesh_subset_BHE_nodes.push_back(
std::make_unique<MeshLib::MeshSubset const>(_mesh,
_vec_BHE_nodes[i]));
std::make_unique<MeshLib::MeshSubset const>(
_mesh, _bheMeshData.BHE_nodes[i]));
int n_BHE_unknowns =
_process_data._vec_BHE_property[i]->getNumUnknowns();
vec_n_BHE_unknowns.emplace_back(n_BHE_unknowns);
......@@ -122,7 +119,7 @@ void HeatTransportBHEProcess::constructDofTable()
// temperature.
vec_n_components.push_back(1);
// now the BHE subsets
for (std::size_t i = 0; i < _vec_BHE_mat_IDs.size(); i++)
for (std::size_t i = 0; i < _bheMeshData.BHE_mat_IDs.size(); i++)
{
// Here the number of components equals to
// the number of unknowns on the BHE
......@@ -132,9 +129,9 @@ void HeatTransportBHEProcess::constructDofTable()
std::vector<std::vector<MeshLib::Element*> const*> vec_var_elements;
// vec_var_elements.push_back(&_vec_pure_soil_elements);
vec_var_elements.push_back(&(_mesh.getElements()));
for (std::size_t i = 0; i < _vec_BHE_elements.size(); i++)
for (std::size_t i = 0; i < _bheMeshData.BHE_elements.size(); i++)
{
vec_var_elements.push_back(&_vec_BHE_elements[i]);
vec_var_elements.push_back(&_bheMeshData.BHE_elements[i]);
}
_local_to_global_index_map =
......@@ -153,12 +150,10 @@ void HeatTransportBHEProcess::initializeConcreteProcess(
MeshLib::Mesh const& mesh,
unsigned const integration_order)
{
const int process_id = 0;
// this process can only run with 3-dimensional mesh
assert(mesh.getDimension() == 3);
ProcessLib::HeatTransportBHE::createLocalAssemblers<
3, /*mesh.getDimension(),*/
HeatTransportBHELocalAssemblerSoil, HeatTransportBHELocalAssemblerBHE>(
3 /* mesh dimension */, HeatTransportBHELocalAssemblerSoil,
HeatTransportBHELocalAssemblerBHE>(
mesh.getElements(), dof_table, _local_assemblers,
_process_data._vec_ele_connected_BHE_IDs,
_process_data._vec_BHE_property, mesh.isAxiallySymmetric(),
......@@ -169,10 +164,13 @@ void HeatTransportBHEProcess::initializeConcreteProcess(
// and one BC at the bottom.
std::vector<std::unique_ptr<BoundaryCondition>> bc_collections =
createBHEBoundaryConditionTopBottom();
const int process_id = 0;
auto& current_process_BCs = _boundary_conditions[process_id];
auto const bc_collections_size = bc_collections.size();
for (unsigned i = 0; i < bc_collections_size; i++)
current_process_BCs.addCreatedBC(std::move(bc_collections[i]));
for (auto& bc_collection : bc_collections)
{
current_process_BCs.addCreatedBC(std::move(bc_collection));
}
}
void HeatTransportBHEProcess::assembleConcreteProcess(const double t,
......@@ -192,19 +190,12 @@ void HeatTransportBHEProcess::assembleConcreteProcess(const double t,
}
void HeatTransportBHEProcess::assembleWithJacobianConcreteProcess(
const double t, GlobalVector const& x, GlobalVector const& xdot,
const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
GlobalVector& b, GlobalMatrix& Jac)
const double /*t*/, GlobalVector const& /*x*/, GlobalVector const& /*xdot*/,
const double /*dxdot_dx*/, const double /*dx_dx*/, GlobalMatrix& /*M*/,
GlobalMatrix& /*K*/, GlobalVector& /*b*/, GlobalMatrix& /*Jac*/)
{
DBUG("AssembleWithJacobian HeatTransportBHE process.");
std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
dof_table = {std::ref(*_local_to_global_index_map)};
// Call global assembler for each local assembly item.
GlobalExecutor::executeMemberDereferenced(
_global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
_local_assemblers, dof_table, t, x, xdot, dxdot_dx, dx_dx, M, K, b, Jac,
_coupled_solutions);
OGS_FATAL(
"HeatTransportBHE: analytical Jacobian assembly is not implemented");
}
void HeatTransportBHEProcess::computeSecondaryVariableConcrete(
......@@ -227,189 +218,51 @@ HeatTransportBHEProcess::createBHEBoundaryConditionTopBottom()
int const n_BHEs = static_cast<int>(_process_data._vec_BHE_property.size());
// bulk dof table
NumLib::LocalToGlobalIndexMap dof_table_bulk = *_local_to_global_index_map;
std::vector<MeshLib::Node*> bc_top_nodes;
std::vector<MeshLib::Node*> bc_bottom_nodes;
// for each BHE
for (int bhe_i = 0; bhe_i < n_BHEs; bhe_i++)
{
// get the BHE type
auto const& bhe_typ = _process_data._vec_BHE_property[bhe_i]->bhe_type;
auto const& bhe_nodes = _bheMeshData.BHE_nodes[bhe_i];
// find the variable ID
// the soil temperature is 0-th variable
// the BHE temperature is therefore bhe_i + 1
const int variable_id = bhe_i + 1;
// find the node in mesh that are at the top
auto const n_bhe_nodes = _vec_BHE_nodes[bhe_i].size();
unsigned int idx_top = 0;
unsigned int idx_bottom = _vec_BHE_nodes[bhe_i].size() - 1;
double top_z_coord = _vec_BHE_nodes[bhe_i][idx_top]->getCoords()[2];
double bottom_z_coord =
_vec_BHE_nodes[bhe_i][idx_bottom]->getCoords()[2];
for (unsigned bhe_node_i = 0; bhe_node_i < n_bhe_nodes; bhe_node_i++)
{
if (_vec_BHE_nodes[bhe_i][bhe_node_i]->getCoords()[2] >=
top_z_coord)
idx_top = bhe_node_i;
if (_vec_BHE_nodes[bhe_i][bhe_node_i]->getCoords()[2] <=
bottom_z_coord)
idx_bottom = bhe_node_i;
}
bc_top_nodes.clear();
bc_top_nodes.emplace_back(_vec_BHE_nodes[bhe_i][idx_top]);
bc_bottom_nodes.clear();
bc_bottom_nodes.emplace_back(_vec_BHE_nodes[bhe_i][idx_bottom]);
MeshLib::MeshSubset bc_mesh_subset_top{_mesh, bc_top_nodes};
MeshLib::Mesh const& bc_mesh_top = bc_mesh_subset_top.getMesh();
MeshLib::MeshSubset bc_mesh_subset_bottom{_mesh, bc_bottom_nodes};
MeshLib::Mesh const& bc_mesh_bottom = bc_mesh_subset_bottom.getMesh();
auto get_global_bhe_bc_index_top = [&](std::size_t const mesh_id,
int const component_id) {
return dof_table_bulk.getGlobalIndex(
{mesh_id, MeshLib::MeshItemType::Node,
bc_top_nodes[0]->getID()},
variable_id, component_id);
};
auto get_global_bhe_bc_index_bottom = [&](std::size_t const mesh_id,
int const component_id) {
return dof_table_bulk.getGlobalIndex(
{mesh_id, MeshLib::MeshItemType::Node,
bc_bottom_nodes[0]->getID()},
variable_id, component_id);
};
// the create_BC function will be repeatedly used
auto create_BC_top_inflow = [&](std::size_t const mesh_id,
int const component_id_in,
int const component_id_out) {
auto const global_index_in =
get_global_bhe_bc_index_top(mesh_id, component_id_in);
auto const global_index_out =
get_global_bhe_bc_index_top(mesh_id, component_id_out);
return ProcessLib::createBHEInflowDirichletBoundaryCondition(
global_index_in, global_index_out, _mesh, bc_top_nodes,
variable_id, component_id_in,
_process_data._vec_BHE_property[bhe_i]);
};
auto create_BC_bottom_outflow = [&](std::size_t const mesh_id,
int const component_id_in,
int const component_id_out) {
auto const global_index_in =
get_global_bhe_bc_index_bottom(mesh_id, component_id_in);
auto const global_index_out =
get_global_bhe_bc_index_bottom(mesh_id, component_id_out);
return ProcessLib::createBHEBottomDirichletBoundaryCondition(
global_index_in, global_index_out, _mesh, bc_top_nodes,
variable_id, component_id_in);
};
// depending on the BHE type
switch (bhe_typ)
// Bottom and top nodes w.r.t. the z coordinate.
auto const bottom_top_nodes = std::minmax_element(
begin(bhe_nodes), end(bhe_nodes),
[&](auto const& a, auto const& b) {
return a->getCoords()[2] < b->getCoords()[2];
});
auto const bc_bottom_node_id = (*bottom_top_nodes.first)->getID();
MeshLib::Node* const bc_top_node = *bottom_top_nodes.second;
auto get_global_bhe_bc_indices =
[&](std::size_t const node_id,
std::pair<int, int> const& in_out_component_id) {
return std::make_pair(
_local_to_global_index_map->getGlobalIndex(
{_mesh.getID(), MeshLib::MeshItemType::Node, node_id},
variable_id, in_out_component_id.first),
_local_to_global_index_map->getGlobalIndex(
{_mesh.getID(), MeshLib::MeshItemType::Node, node_id},
variable_id, in_out_component_id.second));
};
for (auto const& in_out_component_id :
_process_data._vec_BHE_property[bhe_i]
->inflowOutflowBcComponentIds())
{
case ProcessLib::HeatTransportBHE::BHE::BHE_TYPE::TYPE_2U:
{
unsigned const component_id_T_in_1 = 0;
unsigned const component_id_T_out_1 = 2;
unsigned const component_id_T_out_2 = 3;
unsigned const component_id_T_in_2 = 1;
// there are 2 BCs on the top node
std::unique_ptr<BHEInflowDirichletBoundaryCondition> bc_top_1 =
create_BC_top_inflow(bc_mesh_top.getID(),
component_id_T_in_1,
component_id_T_out_1);
std::unique_ptr<BHEInflowDirichletBoundaryCondition> bc_top_2 =
create_BC_top_inflow(bc_mesh_top.getID(),
component_id_T_in_2,
component_id_T_out_2);
// there are also 2 BCs on the bottom node
std::unique_ptr<BHEBottomDirichletBoundaryCondition>
bc_bottom_1 = create_BC_bottom_outflow(
bc_mesh_bottom.getID(), component_id_T_in_1,
component_id_T_out_1);
std::unique_ptr<BHEBottomDirichletBoundaryCondition>
bc_bottom_2 = create_BC_bottom_outflow(
bc_mesh_bottom.getID(), component_id_T_in_2,
component_id_T_out_2);
// add bc_top and bc_bottom to the vector
bcs.push_back(std::move(bc_top_1));
bcs.push_back(std::move(bc_top_2));
bcs.push_back(std::move(bc_bottom_1));
bcs.push_back(std::move(bc_bottom_2));
}
break;
case ProcessLib::HeatTransportBHE::BHE::BHE_TYPE::TYPE_1U:
{
unsigned const component_id_T_in = 0;
unsigned const component_id_T_out = 1;
// there is one BC on the top node
std::unique_ptr<BHEInflowDirichletBoundaryCondition> bc_top =
create_BC_top_inflow(bc_mesh_top.getID(), component_id_T_in,
component_id_T_out);
// there is also 1 BC on the bottom node
std::unique_ptr<BHEBottomDirichletBoundaryCondition> bc_bottom =
create_BC_bottom_outflow(bc_mesh_bottom.getID(),
component_id_T_in,
component_id_T_out);
// add bc_top and bc_bottom to the vector
bcs.push_back(std::move(bc_top));
bcs.push_back(std::move(bc_bottom));
}
break;
case ProcessLib::HeatTransportBHE::BHE::BHE_TYPE::TYPE_CXC:
{
unsigned const component_id_T_in = 1;
unsigned const component_id_T_out = 0;
// there is one BC on the top node
std::unique_ptr<BHEInflowDirichletBoundaryCondition> bc_top =
create_BC_top_inflow(bc_mesh_top.getID(), component_id_T_in,
component_id_T_out);
// there is also 1 BC on the bottom node
std::unique_ptr<BHEBottomDirichletBoundaryCondition> bc_bottom =
create_BC_bottom_outflow(bc_mesh_bottom.getID(),
component_id_T_in,
component_id_T_out);
// add bc_top and bc_bottom to the vector
bcs.push_back(std::move(bc_top));
bcs.push_back(std::move(bc_bottom));
}
break;
case ProcessLib::HeatTransportBHE::BHE::BHE_TYPE::TYPE_CXA:
{
unsigned const component_id_T_in = 0;
unsigned const component_id_T_out = 1;
// there is one BC on the top node
std::unique_ptr<BHEInflowDirichletBoundaryCondition> bc_top =
create_BC_top_inflow(bc_mesh_top.getID(), component_id_T_in,
component_id_T_out);
// there is also 1 BC on the bottom node
std::unique_ptr<BHEBottomDirichletBoundaryCondition> bc_bottom =
create_BC_bottom_outflow(bc_mesh_bottom.getID(),
component_id_T_in,
component_id_T_out);
// add bc_top and bc_bottom to the vector
bcs.push_back(std::move(bc_top));
bcs.push_back(std::move(bc_bottom));
}
break;
default:
OGS_FATAL("WRONG BHE TYPE");
// Top, inflow.
bcs.push_back(ProcessLib::createBHEInflowDirichletBoundaryCondition(
get_global_bhe_bc_indices(bc_top_node->getID(),
in_out_component_id),
_mesh, {bc_top_node}, variable_id, in_out_component_id.first,
_process_data._vec_BHE_property[bhe_i]));
// Bottom, outflow.
bcs.push_back(ProcessLib::createBHEBottomDirichletBoundaryCondition(
get_global_bhe_bc_indices(bc_bottom_node_id,
in_out_component_id)));
}
}
......
ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.h
+ 8
− 22
View file @ 9fffe1dd
......@@ -10,6 +10,7 @@
#pragma once
#include "HeatTransportBHEProcessData.h"
#include "ProcessLib/HeatTransportBHE/BHE/MeshUtils.h"
#include "ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHEProcessAssemblerInterface.h"
#include "ProcessLib/Process.h"
......@@ -17,6 +18,8 @@ namespace ProcessLib
{
namespace HeatTransportBHE
{
struct BHEMeshData;
class HeatTransportBHEProcess final : public Process
{
public:
......@@ -64,35 +67,18 @@ private:
std::vector<std::unique_ptr<HeatTransportBHELocalAssemblerInterface>>
_local_assemblers;
/**
* These are the elements that are representing BHEs
*/
std::vector<std::vector<MeshLib::Element*>> _vec_BHE_elements;
/**
* These are the elements that are not connected with any BHE
*/
std::vector<MeshLib::Element*> _vec_pure_soil_elements;
/**
* These are the soil nodes that are not connected with a BHE
*/
std::vector<MeshLib::Node*> _vec_pure_soil_nodes;
/**
* Mesh nodes that are located on any BHE
* ordered according to each BHE
*/
std::vector<std::vector<MeshLib::Node*>> _vec_BHE_nodes;
std::vector<std::unique_ptr<MeshLib::MeshSubset const>>
_mesh_subset_BHE_nodes;
std::vector<std::unique_ptr<MeshLib::MeshSubset const>>
_mesh_subset_BHE_soil_nodes;
std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_pure_soil_nodes;
std::unique_ptr<MeshLib::MeshSubset const>
_mesh_subset_soil_nodes_connected_with_BHE;
std::vector<int> _vec_BHE_mat_IDs;
const BHEMeshData _bheMeshData;
};
} // namespace HeatTransportBHE
} // namespace ProcessLib
ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h
+ 0
− 11
View file @ 9fffe1dd
......@@ -42,7 +42,6 @@ public:
HeatTransportBHELocalAssemblerBHE(
MeshLib::Element const& e,
std::size_t const local_matrix_size,
std::vector<unsigned> const& dofIndex_to_localIndex,
bool const is_axially_symmetric,
unsigned const integration_order,
......@@ -53,16 +52,6 @@ public:
std::vector<double>& /*local_K_data*/,
std::vector<double>& /*local_b_data*/) override;
void assembleWithJacobian(double const /*t*/,
Eigen::VectorXd const& /*local_u*/,
Eigen::VectorXd& /*local_b*/,
Eigen::MatrixXd& /*local_J*/) override
{
OGS_FATAL(
"HeatTransportBHELocalAssembler: assembly with jacobian is not "
"implemented.");
}
void preTimestepConcrete(std::vector<double> const& /*local_x*/,
double const /*t*/,
double const /*delta_t*/) override
......
ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE_impl.h
+ 2
− 5
View file @ 9fffe1dd
......@@ -25,7 +25,6 @@ template <typename ShapeFunction, typename IntegrationMethod, int BHE_Dim>
HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod, BHE_Dim>::
HeatTransportBHELocalAssemblerBHE(
MeshLib::Element const& e,
std::size_t const /*local_matrix_size*/,
std::vector<unsigned> const& dofIndex_to_localIndex,
bool const is_axially_symmetric,
unsigned const integration_order,
......@@ -62,9 +61,8 @@ HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod, BHE_Dim>::
{
x_position.setIntegrationPoint(ip);
IntegrationPointDataBHE<ShapeMatricesType> int_Point_Data_BHE(
*(_process_data._vec_BHE_property[BHE_id]));
_ip_data.emplace_back(int_Point_Data_BHE);
// create the class IntegrationPointDataBHE in place
_ip_data.emplace_back(*(_process_data._vec_BHE_property[BHE_id]));
auto const& sm = shape_matrices[ip];
auto& ip_data = _ip_data[ip];
ip_data.integration_weight =
......@@ -109,7 +107,6 @@ HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod, BHE_Dim>::
_process_data._vec_BHE_property[BHE_id]->setRMatrices(
idx_bhe_unknowns, ShapeFunction::NPOINTS, matBHE_loc_R, _R_matrix,
_R_pi_s_matrix, _R_s_matrix);
} // end of loop over BHE unknowns
// debugging
......
ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil.h
+ 2
− 28
View file @ 9fffe1dd
......@@ -26,8 +26,6 @@ namespace ProcessLib
{
namespace HeatTransportBHE
{
const unsigned NUM_NODAL_DOF_SOIL = 1;
template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
class HeatTransportBHELocalAssemblerSoil
: public HeatTransportBHELocalAssemblerInterface
......@@ -45,7 +43,6 @@ public:
HeatTransportBHELocalAssemblerSoil(
MeshLib::Element const& e,
std::size_t const local_matrix_size,
std::vector<unsigned> const& dofIndex_to_localIndex,
bool is_axially_symmetric,
unsigned const integration_order,
......@@ -56,29 +53,6 @@ public:
std::vector<double>& /*local_K_data*/,
std::vector<double>& /*local_b_data*/) override;
void assembleWithJacobian(double const /*t*/,
std::vector<double> const& /*local_x*/,
std::vector<double> const& /*local_xdot*/,
const double /*dxdot_dx*/, const double /*dx_dx*/,
std::vector<double>& /*local_M_data*/,
std::vector<double>& /*local_K_data*/,
std::vector<double>& /*local_b_data*/,
std::vector<double>& /*local_Jac_data*/) override
{
OGS_FATAL(
"HeatTransportBHELocalAssemblerMatrix: assembly with jacobian is "
"not "
"implemented.");
}
void preTimestepConcrete(std::vector<double> const& /*local_x*/,
double const /*t*/,
double const /*delta_t*/) override
{
}
void postTimestepConcrete(std::vector<double> const& /*local_x*/) override;
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
const unsigned integration_point) const override
{
......@@ -92,8 +66,8 @@ private:
HeatTransportBHEProcessData& _process_data;
std::vector<
IntegrationPointDataSoil<ShapeMatricesType>,
Eigen::aligned_allocator<IntegrationPointDataSoil<ShapeMatricesType>>>
IntegrationPointDataSoil<NodalMatrixType>,
Eigen::aligned_allocator<IntegrationPointDataSoil<NodalMatrixType>>>
_ip_data;
IntegrationMethod const _integration_method;
......
ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil_impl.h
+ 29
− 21
View file @ 9fffe1dd
......@@ -37,7 +37,6 @@ HeatTransportBHELocalAssemblerSoil<ShapeFunction,
GlobalDim>::
HeatTransportBHELocalAssemblerSoil(
MeshLib::Element const& e,
std::size_t const /*local_matrix_size*/,
std::vector<unsigned> const& dofIndex_to_localIndex,
bool const is_axially_symmetric,
unsigned const integration_order,
......@@ -60,6 +59,26 @@ HeatTransportBHELocalAssemblerSoil<ShapeFunction,
IntegrationMethod,
GlobalDim>(
e, is_axially_symmetric, _integration_method);
SpatialPosition x_position;
x_position.setElementID(element_id);
// ip data initialization
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
x_position.setIntegrationPoint(ip);
// create the class IntegrationPointDataBHE in place
_ip_data.emplace_back();
auto const& sm = _shape_matrices[ip];
auto& ip_data = _ip_data[ip];
double const w = _integration_method.getWeightedPoint(ip).getWeight() *
sm.integralMeasure * sm.detJ;
ip_data.NTN_product_times_w = sm.N.transpose() * sm.N * w;
ip_data.dNdxTdNdx_product_times_w = sm.dNdx.transpose() * sm.dNdx * w;
_secondary_data.N[ip] = sm.N;
}
}
template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
......@@ -72,14 +91,13 @@ void HeatTransportBHELocalAssemblerSoil<
std::vector<double>& local_K_data,
std::vector<double>& /*local_b_data*/)
{
auto const local_matrix_size = local_x.size();
assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF_SOIL);
assert(local_x.size() == ShapeFunction::NPOINTS);
(void)local_x; // Avoid unused arg warning.
auto local_M = MathLib::createZeroedMatrix<NodalMatrixType>(
local_M_data, local_matrix_size, local_matrix_size);
local_M_data, ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
auto local_K = MathLib::createZeroedMatrix<NodalMatrixType>(
local_K_data, local_matrix_size, local_matrix_size);
local_K_data, ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
......@@ -90,8 +108,9 @@ void HeatTransportBHELocalAssemblerSoil<
for (unsigned ip = 0; ip < n_integration_points; ip++)
{
pos.setIntegrationPoint(ip);
auto const& sm = _shape_matrices[ip];
auto const& wp = _integration_method.getWeightedPoint(ip);
auto& ip_data = _ip_data[ip];
auto const& M_ip = ip_data.NTN_product_times_w;
auto const& K_ip = ip_data.dNdxTdNdx_product_times_w;
// auto const k_f = _process_data.thermal_conductivity_fluid(t, pos)[0];
// auto const k_g = _process_data.thermal_conductivity_gas(t, pos)[0];
......@@ -110,13 +129,10 @@ void HeatTransportBHELocalAssemblerSoil<
// for now only using the solid phase parameters
// assemble Conductance matrix
local_K.noalias() += sm.dNdx.transpose() * k_s * sm.dNdx * sm.detJ *
wp.getWeight() * sm.integralMeasure;
local_K.noalias() += K_ip * k_s;
// assemble Mass matrix
local_M.noalias() += sm.N.transpose() * density_s * heat_capacity_s *
sm.N * sm.detJ * wp.getWeight() *
sm.integralMeasure;
local_M.noalias() += M_ip * density_s * heat_capacity_s;
}
// debugging
......@@ -125,13 +141,5 @@ void HeatTransportBHELocalAssemblerSoil<
// std::cout << local_K.format(CleanFmt) << sep;
// std::cout << local_M.format(CleanFmt) << sep;
}
template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
void HeatTransportBHELocalAssemblerSoil<
ShapeFunction,
IntegrationMethod,
GlobalDim>::postTimestepConcrete(std::vector<double> const& /*local_x*/)
{
}
} // namespace HeatTransportBHE
} // namespace ProcessLib
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment