Commit 1a26db70 authored by Shuang Chen's avatar Shuang Chen
Browse files

remove soil temperature part in the variable x vector during Schur Complement procedure

parent ab4162fd
......@@ -125,12 +125,17 @@ void HeatTransportBHEProcess::constructDofTable()
// With adopting the Schur complement procedure.
if (_process_data.use_schur_complement_method)
{
// Dof of the Soil part.
// Dof of the Soil and BHE part.
std::vector<MeshLib::MeshSubset> soil_mesh_subsets{
*_mesh_subset_soil_nodes};
std::vector<std::vector<MeshLib::Element*> const*>
vec_var_elements_soil;
std::vector<MeshLib::Element*> vec_elements_soil;
std::vector<MeshLib::MeshSubset> BHE_mesh_subsets{};
std::vector<std::vector<MeshLib::Element*> const*> vec_var_elements_BHE;
std::vector<MeshLib::Element*> vec_elements_BHE;
for (auto e : _mesh.getElements())
{
if (e->getDimension() == 3)
......@@ -141,7 +146,34 @@ void HeatTransportBHEProcess::constructDofTable()
std::vector<int> vec_n_components_soil{
1}; // one component for the soil temperature variable.
std::vector<int> vec_n_components_BHE{}; // component for the BHE
// temperature variable.
// the BHE nodes need to be cherry-picked from the vector
for (int i = 0; i < n_BHEs; i++)
{
auto const number_of_unknowns =
visit([](auto const& bhe) { return bhe.number_of_unknowns; },
_process_data._vec_BHE_property[i]);
auto const& bhe_nodes = _bheMeshData.BHE_nodes[i];
auto const& bhe_elements = _bheMeshData.BHE_elements[i];
// All the BHE nodes have additional variables.
_mesh_subset_BHE_nodes.push_back(
std::make_unique<MeshLib::MeshSubset const>(_mesh, bhe_nodes));
std::generate_n(std::back_inserter(BHE_mesh_subsets),
// Here the number of components equals to the
// number of unknowns on the BHE
number_of_unknowns,
[&ms = _mesh_subset_BHE_nodes.back()]()
{ return *ms; });
vec_n_components_BHE.push_back(number_of_unknowns);
vec_var_elements_BHE.push_back(&bhe_elements);
}
// soil part
_local_to_global_index_map =
std::make_unique<NumLib::LocalToGlobalIndexMap>(
std::move(soil_mesh_subsets),
......@@ -152,6 +184,17 @@ void HeatTransportBHEProcess::constructDofTable()
// in case of debugging the dof table, activate the following line
// std::cout << *_local_to_global_index_map << "\n";
// BHE part
_local_to_global_index_map_BHE =
std::make_unique<NumLib::LocalToGlobalIndexMap>(
std::move(BHE_mesh_subsets),
vec_n_components_BHE,
vec_var_elements_BHE,
NumLib::ComponentOrder::BY_COMPONENT);
// in case of debugging the dof table, activate the following line
// std::cout << *_local_to_global_index_map_BHE << "\n";
// Global part, which contains soil part and BHE part.
_local_to_global_index_map_global =
std::make_unique<NumLib::LocalToGlobalIndexMap>(
......@@ -348,7 +391,7 @@ void HeatTransportBHEProcess::SchurComplementConcreteProcess(
size_t size_BHE_var = BHE_var->size();
// TODO: This is potentially problematic.
// The soil temperature from the same
// location is needed.
// location is needed.
double T_ini = x[process_id]->get(0);
for (int i = 0; i < size_BHE_var; i++)
{
......@@ -373,21 +416,8 @@ void HeatTransportBHEProcess::SchurComplementConcreteProcess(
std::size_t n_SOIL_nodes = x[process_id]->size();
// TODO: Why is x_soil_BHE_var needed?
// update the global entire soil and BHE variables vector.
GlobalVector x_soil_BHE_var(total_BHE_SOIL_unknowns_num);
x_soil_BHE_var.setZero();
for (auto i = 0; i < total_BHE_SOIL_unknowns_num; i++)
{
if (i < n_SOIL_nodes)
// TODO: Is this correct?
x_soil_BHE_var[i] = x[process_id]->get(0);
else
{
int pos = i - n_SOIL_nodes;
x_soil_BHE_var[i] = cur_BHE_var[pos];
}
}
auto size_soil_DOFs = n_SOIL_nodes;
auto size_BHE_DOFs = total_BHE_SOIL_unknowns_num - n_SOIL_nodes;
// compute A and Rhs
// compute A
......@@ -399,18 +429,24 @@ void HeatTransportBHEProcess::SchurComplementConcreteProcess(
LinAlg::copy(*M_BHE, *A_entire);
LinAlg::aypx(*A_entire, dxdot_dx, *K_BHE);
// compute Rhs
auto Rhs_entire =
std::make_unique<GlobalVector>(total_BHE_SOIL_unknowns_num);
// compute Rhs for BHE part only
auto Rhs_entire = std::make_unique<GlobalVector>(size_BHE_DOFs);
// only use the BHE part in M_BHE and b_BHE.
GlobalMatrix M_BHE_(size_BHE_DOFs);
M_BHE_.getRawMatrix() = M_BHE->getRawMatrix().block(
n_SOIL_nodes, n_SOIL_nodes, size_BHE_DOFs, size_BHE_DOFs);
GlobalVector b_BHE_(size_BHE_DOFs);
b_BHE_.getRawVector() = b_BHE->getRawVector().tail(size_BHE_DOFs);
std::size_t _tmp_id = 0u;
auto& tmp = NumLib::GlobalVectorProvider::provider.getVector(_tmp_id);
const GlobalVector x_prev = x_soil_BHE_var;
const GlobalVector x_prev = cur_BHE_var;
// getWeightedOldX: y = x_old / delta_t
LinAlg::copy(x_prev, tmp);
LinAlg::scale(tmp, 1.0 / dt);
// rhs = M * weighted_old_x + b
LinAlg::matMultAdd(*M_BHE, tmp, *b_BHE, *Rhs_entire);
LinAlg::matMultAdd(M_BHE_, tmp, b_BHE_, *Rhs_entire);
NumLib::GlobalVectorProvider::provider.releaseVector(tmp);
// output check
......@@ -418,9 +454,6 @@ void HeatTransportBHEProcess::SchurComplementConcreteProcess(
// Rhs_entire.get()->write(".\\output\\Rhs.txt");
// divide A into A_pi, R_pi_s, R_pi and stored.
auto size_soil_DOFs = n_SOIL_nodes;
auto size_BHE_DOFs = total_BHE_SOIL_unknowns_num - n_SOIL_nodes;
// R_pi matrix, to be noticed here the -R_pi is obtained.
GlobalMatrix minus_R_pi(n_SOIL_nodes);
minus_R_pi.getRawMatrix() =
......@@ -445,17 +478,16 @@ void HeatTransportBHEProcess::SchurComplementConcreteProcess(
// BHE boundary condition procedure
auto& bc = _boundary_conditions[process_id];
auto vec_new_bc =
bc.getKnownSolutions_SchurComplement(t, x_soil_BHE_var);
auto vec_new_bc = bc.getKnownSolutions_SchurComplement(t, cur_BHE_var);
int vec_new_bc_inflow_bottomNode_ids = vec_new_bc[0][1].ids[0];
// HS: TODO! from double to int? What's going on here?
// HS: TODO! from double to int? What's going on here?
int vec_new_bc_outflow_bottomNode_ids = vec_new_bc[0][1].values[0];
// Adding additional coorelation for solving the under-determined linear
// equation system.
int id_bottomNode_inflow_pipe_bhe_primary_var_x_pos =
vec_new_bc_inflow_bottomNode_ids - size_soil_DOFs;
vec_new_bc_inflow_bottomNode_ids;
int id_bottomNode_outflow_pipe_bhe_primary_var_x_pos =
vec_new_bc_outflow_bottomNode_ids - size_soil_DOFs;
vec_new_bc_outflow_bottomNode_ids;
// Adding an additional row
GlobalVector A_pi_plus(size_BHE_DOFs);
......@@ -482,7 +514,7 @@ void HeatTransportBHEProcess::SchurComplementConcreteProcess(
B_pi_new.setZero();
for (std::size_t i = 0; i < size_BHE_DOFs; i++)
{
B_pi_new[i] = Rhs_entire->get(size_soil_DOFs + i);
B_pi_new[i] = Rhs_entire->get(i);
}
// output check
// B_pi_new.write(".\\output\\B_pi_new.txt");
......@@ -493,7 +525,7 @@ void HeatTransportBHEProcess::SchurComplementConcreteProcess(
// The imposed inflow BC at the top of the BHE
double bc_T = vec_new_bc_inflow_topNode_val;
int id_topNode_inflow_pipe_bhe_primary_var_x_pos =
vec_new_bc_inflow_topNode_ids - size_soil_DOFs;
vec_new_bc_inflow_topNode_ids;
int bc_id = id_topNode_inflow_pipe_bhe_primary_var_x_pos;
// rhs - col_x
GlobalVector rhs_minus(size_BHE_DOFs + 1);
......@@ -838,9 +870,11 @@ void HeatTransportBHEProcess::createBHEBoundaryConditionTopBottom(
{
if (_process_data.use_schur_complement_method)
{
return _local_to_global_index_map_global->getGlobalIndex(
// no soil temperature, the BHE temperature is therefore bhe_i
auto variable_id_BHE = variable_id - 1;
return _local_to_global_index_map_BHE->getGlobalIndex(
{_mesh.getID(), MeshLib::MeshItemType::Node, node_id},
variable_id, component);
variable_id_BHE, component);
}
else
{
......
......@@ -103,6 +103,9 @@ private:
const BHEMeshData _bheMeshData;
std::unique_ptr<NumLib::LocalToGlobalIndexMap>
_local_to_global_index_map_BHE;
std::unique_ptr<NumLib::LocalToGlobalIndexMap>
_local_to_global_index_map_global;
......
Markdown is supported
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