From d1f2c586e8129f0dcc6693b8ff599d73247d2c95 Mon Sep 17 00:00:00 2001
From: Julie PATELA <julie.patela.ocre@cea.fr>
Date: Mon, 3 Aug 2020 17:05:48 +0200
Subject: [PATCH] WIP debugging diamond scheme for heat equation [ci skip]
---
.../algorithms/HeatDiamondAlgorithm.cpp | 137 ++++++++++++++----
1 file changed, 111 insertions(+), 26 deletions(-)
diff --git a/src/language/algorithms/HeatDiamondAlgorithm.cpp b/src/language/algorithms/HeatDiamondAlgorithm.cpp
index 3a2424705..09c6a1328 100644
--- a/src/language/algorithms/HeatDiamondAlgorithm.cpp
+++ b/src/language/algorithms/HeatDiamondAlgorithm.cpp
@@ -34,7 +34,38 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
std::shared_ptr m_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
if constexpr (Dimension > 1) {
+ // if constexpr (Dimension == 2) {
+ // NodeValue<TinyVector<Dimension>> new_coords{m_mesh->connectivity()};
+
+ // new_coords[NodeId{0}] = TinyVector<2>{-1, -1};
+ // new_coords[NodeId{1}] = TinyVector<2>{-1, 0};
+ // new_coords[NodeId{2}] = TinyVector<2>{-1, 1};
+ // new_coords[NodeId{3}] = TinyVector<2>{0, -1};
+ // new_coords[NodeId{4}] = TinyVector<2>{0.2, 0.2};
+ // new_coords[NodeId{5}] = TinyVector<2>{0, 1};
+ // new_coords[NodeId{6}] = TinyVector<2>{1, -1};
+ // new_coords[NodeId{7}] = TinyVector<2>{1, 0};
+ // new_coords[NodeId{8}] = TinyVector<2>{1, 1};
+
+ // m_mesh = std::make_shared<MeshType>(m_mesh->shared_connectivity(), new_coords);
+ // }
+
+ // for (NodeId node_id = 0; node_id < m_mesh->numberOfNodes(); ++node_id) {
+ // std::cout << "x(" << node_id << ") = " << m_mesh->xr()[node_id] << '\n';
+ // }
+
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+ // {
+ // const auto& cell_to_face_matrix = m_mesh->connectivity().cellToFaceMatrix();
+ // CellId cell0 = CellId{0};
+
+ // const auto& ll = mesh_data.ll();
+
+ // for (size_t i_face = 0; i_face < cell_to_face_matrix[cell0].size(); ++i_face) {
+ // FaceId face_id = cell_to_face_matrix[cell0][i_face];
+ // std::cout << "ll[" << face_id << "]=" << ll[face_id] << '\n';
+ // }
+ // }
CellValue<double> Tj =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(T_id, mesh_data.xj());
@@ -188,9 +219,9 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
for (size_t j_cell = 0; j_cell < primal_face_to_cell.size(); ++j_cell) {
const CellId cell_id2 = primal_face_to_cell[j_cell];
if (i_cell == j_cell) {
- S(cell_id1, cell_id2) -= beta_l;
- } else {
S(cell_id1, cell_id2) += beta_l;
+ } else {
+ S(cell_id1, cell_id2) -= beta_l;
}
}
}
@@ -209,19 +240,35 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
NodeValue<const NodeId> dual_node_primal_node_id = [=]() {
CellValue<NodeId> cell_ignored_id{m_mesh->connectivity()};
- cell_ignored_id.fill(NodeId{std::numeric_limits<NodeId>::max()});
+ cell_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
+
NodeValue<NodeId> node_primal_id{m_mesh->connectivity()};
parallel_for(
- m_mesh->numberOfCells(), PUGS_LAMBDA(NodeId node_id) { node_primal_id[node_id] = node_id; });
+ m_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { node_primal_id[node_id] = node_id; });
NodeValue<NodeId> computed_dual_node_primal_node_id{diamond_mesh->connectivity()};
mapper->toDualNode(node_primal_id, cell_ignored_id, computed_dual_node_primal_node_id);
+ for (NodeId node_id = 0; node_id < diamond_mesh->numberOfNodes(); ++node_id) {
+ std::cout << node_id << " -> " << computed_dual_node_primal_node_id[node_id] << '\n';
+ }
+
return computed_dual_node_primal_node_id;
}();
+ // for (CellId cell_id = 0; cell_id < m_mesh->numberOfCells(); ++cell_id) {
+ // std::cout << "xj[" << cell_id << "] = " << xj[cell_id] << "\n";
+ // }
+ // for (CellId dual_cell_id = 0; dual_cell_id < diamond_mesh->numberOfCells(); ++dual_cell_id) {
+ // std::cout << "Vl[" << dual_cell_id << "] = " << dual_Vj[dual_cell_id] << "\n";
+ // }
+ // for (FaceId face_id = 0; face_id < m_mesh->numberOfFaces(); ++face_id) {
+ // std::cout << "alpha_l[" << face_id << "] = " << alpha_l[face_id] << "\n";
+ // }
+ // std::cout << "S(0,0) = " << S(0, 0) << '\n';
+
const auto& dual_Cjr = diamond_mesh_data.Cjr();
const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
@@ -231,19 +278,30 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
const auto& face_local_numbers_in_their_cells = m_mesh->connectivity().faceLocalNumbersInTheirCells();
const auto& primal_node_to_cell_matrix = m_mesh->connectivity().nodeToCellMatrix();
+ const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
+
for (FaceId face_id = 0; face_id < m_mesh->numberOfFaces(); ++face_id) {
+ std::cout << "*** face_id -> " << face_id << "\n";
+ // std::cout << " x_l -> " << xl[face_id] << "\n";
if (face_to_cell_matrix[face_id].size() > 1) {
const double alpha_face_id = alpha_l[face_id];
+ // std::cout << " alpha_face_id -> " << alpha_face_id << "\n";
for (size_t i_face_cell = 0; i_face_cell < face_to_cell_matrix[face_id].size(); ++i_face_cell) {
CellId i_id = face_to_cell_matrix[face_id][i_face_cell];
-
+ // std::cout << " cell_id_i_ -> " << i_id << "\n";
+ // std::cout << " x_i -> " << xj[i_id] << "\n";
const bool is_face_reversed_for_cell_i =
primal_face_cell_is_reversed(i_id, face_local_numbers_in_their_cells(face_id, i_face_cell));
+ std::cout << " i_is_reversed -> " << is_face_reversed_for_cell_i << "\n";
+
+ // for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
for (size_t j_face_cell = 0; j_face_cell < face_to_cell_matrix[face_id].size(); ++j_face_cell) {
CellId j_id = face_to_cell_matrix[face_id][j_face_cell];
CellId dual_cell_id = face_dual_cell_id[face_id];
+ // std::cout << " cell_id_j -> " << j_id << "\n";
+ // std::cout << " x_j -> " << xj[j_id] << "\n";
for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
const TinyVector<Dimension> nil = [&] {
@@ -253,9 +311,12 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
return primal_nlr(face_id, i_node);
}
}();
+ std::cout << " nil -> " << nil << "[" << i_id << face_id << "]"
+ << "\n";
NodeId face_node_id = primal_face_to_node_matrix[face_id][i_node];
const auto& primal_node_to_cell = primal_node_to_cell_matrix[face_node_id];
+ std::cout << " x_r -> " << xr[face_node_id] << "\n";
if (not primal_node_is_on_boundary[face_node_id]) {
const double weight_rj = [&] {
for (size_t i_cell = 0; i_cell < primal_node_to_cell.size(); ++i_cell) {
@@ -266,11 +327,14 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
Assert(false, "could not determine the weight");
return std::numeric_limits<double>::signaling_NaN();
}();
- for (size_t i_dual_node = 0; i_dual_node < dual_Cjr.numberOfSubValues(dual_cell_id); ++i_dual_node) {
+ for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size();
+ ++i_dual_node) {
const NodeId dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
- if (dual_node_id == face_node_id) {
+ if (dual_node_primal_node_id[dual_node_id] == face_node_id) {
const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
- S(i_id, j_id) += weight_rj * alpha_face_id * (nil, Clr);
+ std::cout << " Clr -> " << Clr << "\n";
+ std::cout << " w_rj -> " << weight_rj << "\n";
+ S(i_id, j_id) -= weight_rj * alpha_face_id * (nil, Clr);
}
}
} // else {
@@ -281,7 +345,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
}
}
}
-
+ // std::exit(0);
CellValue<double> fj =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id, mesh_data.xj());
@@ -299,36 +363,45 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
b[cell_id] = fj[cell_id] * primal_Vj[cell_id];
const auto& primal_cell_to_face = primal_cell_to_face_matrix[cell_id];
- for (size_t face_id = 0; face_id < primal_cell_to_face.size(); ++face_id) { // Num local
- FaceId i_cell_face = primal_cell_to_face_matrix[cell_id][face_id]; // Num global
+ for (size_t i_cell_face = 0; i_cell_face < primal_cell_to_face.size(); ++i_cell_face) {
+ FaceId face_id = primal_cell_to_face_matrix[cell_id][i_cell_face];
const bool is_face_reversed_for_cell_i = [=] {
- for (size_t i_face = 0; i_face < primal_cell_to_face.size(); ++i_face) { // Num local
- FaceId primal_face_id = primal_cell_to_face[i_face]; // Num global
- if (primal_face_id == i_cell_face) {
+ for (size_t i_face = 0; i_face < primal_cell_to_face.size(); ++i_face) {
+ FaceId primal_face_id = primal_cell_to_face[i_face];
+ if (primal_face_id == face_id) {
return primal_face_cell_is_reversed(cell_id, i_face);
}
}
+ Assert(false, "cannot find cell's face");
+ return false;
}();
- const auto& primal_face_to_node = primal_face_to_node_matrix[i_cell_face];
- for (size_t i_node = 0; i_node < primal_face_to_node.size(); ++i_node) { // Num local
+ const auto& primal_face_to_node = primal_face_to_node_matrix[face_id];
+ for (size_t i_node = 0; i_node < primal_face_to_node.size(); ++i_node) {
const TinyVector<Dimension> nil = [=] {
if (is_face_reversed_for_cell_i) {
- return -primal_nlr(i_cell_face, i_node);
+ return -primal_nlr(face_id, i_node);
} else {
- return primal_nlr(i_cell_face, i_node);
+ return primal_nlr(face_id, i_node);
}
}();
- CellId dual_cell_id = face_dual_cell_id[i_cell_face]; // Num global dual de la face
+ CellId dual_cell_id = face_dual_cell_id[face_id];
- NodeId face_node_id = primal_face_to_node[i_node]; // Num global
+ NodeId face_node_id = primal_face_to_node[i_node];
if (primal_node_is_on_boundary[face_node_id]) {
- for (size_t i_dual_node = 0; i_dual_node < dual_Cjr.numberOfSubValues(dual_cell_id); ++i_dual_node) {
+ for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size(); ++i_dual_node) {
const NodeId dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
- if (dual_node_id == face_node_id) {
+ if (dual_node_primal_node_id[dual_node_id] == face_node_id) {
const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
- b[cell_id] -= alpha_l[i_cell_face] * gr[face_node_id] * (nil, Clr);
+ b[cell_id] += alpha_l[face_id] * gr[face_node_id] * (nil, Clr);
+ if (cell_id == 0) {
+ std::cout << "***node_id" << face_node_id << "\n";
+ std::cout << " *** alpha_l " << alpha_l[face_id] << "\n";
+ std::cout << " *** Clr " << Clr << "\n";
+ std::cout << " *** gr " << gr[face_node_id] << "\n";
+ std::cout << " *** nil " << nil << "\n";
+ }
}
}
}
@@ -337,8 +410,20 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
}
Vector<double> T{m_mesh->numberOfCells()};
- T = 0;
- BiCGStab{b, A, T, 1000, 1e-12};
+ T = 1;
+
+ for (size_t i = 0; i < b.size(); ++i) {
+ std::cout << "b[" << i << "]=" << b[i] << '\n';
+ }
+
+ Vector AT = A * T;
+
+ for (size_t i = 0; i < b.size(); ++i) {
+ std::cout << "AT[" << i << "]=" << AT[i] << '\n';
+ }
+
+ // PCG{b, A, A, T, 1000, 1e-12};
+ BiCGStab{b, A, T, 1000, 1e-15};
CellValue<double> Temperature{m_mesh->connectivity()};
@@ -355,7 +440,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
{
VTKWriter vtk_writer("Temperature_" + std::to_string(Dimension), 0.01);
- vtk_writer.write(m_mesh, {NamedItemValue{"T", Temperature}}, 0,
+ vtk_writer.write(m_mesh, {NamedItemValue{"T", Temperature}, NamedItemValue{"Texact", Tj}}, 0,
true); // forces last output
}
}
--
GitLab