From b3e3bcc7ca3005815a8f0cd39eeaa32feab13dd1 Mon Sep 17 00:00:00 2001
From: Julie PATELA <julie.patela.ocre@cea.fr>
Date: Tue, 4 Aug 2020 14:55:11 +0200
Subject: [PATCH] Fix diamond scheme for Heat equation with Dirichlet
---
.../algorithms/HeatDiamondAlgorithm.cpp | 148 ++++--------------
1 file changed, 32 insertions(+), 116 deletions(-)
diff --git a/src/language/algorithms/HeatDiamondAlgorithm.cpp b/src/language/algorithms/HeatDiamondAlgorithm.cpp
index 09c6a1328..c329feb02 100644
--- a/src/language/algorithms/HeatDiamondAlgorithm.cpp
+++ b/src/language/algorithms/HeatDiamondAlgorithm.cpp
@@ -34,38 +34,7 @@ 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());
@@ -210,9 +179,8 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
SparseMatrixDescriptor S{m_mesh->numberOfCells()};
for (FaceId face_id = 0; face_id < m_mesh->numberOfFaces(); ++face_id) {
- const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
- const double beta_l = 0.5 * alpha_l[face_id] * mes_l[face_id];
- const NodeValuePerCell<const TinyVector<Dimension>>& Cjr = mesh_data.Cjr();
+ const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
+ const double beta_l = 0.5 * alpha_l[face_id] * mes_l[face_id];
for (size_t i_cell = 0; i_cell < primal_face_to_cell.size(); ++i_cell) {
const CellId cell_id1 = primal_face_to_cell[i_cell];
@@ -251,24 +219,9 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
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();
@@ -278,32 +231,19 @@ 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 i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
+ NodeId node_id = primal_face_to_node_matrix[face_id][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) {
+ if (not primal_node_is_on_boundary[node_id]) {
const TinyVector<Dimension> nil = [&] {
if (is_face_reversed_for_cell_i) {
return -primal_nlr(face_id, i_node);
@@ -311,41 +251,29 @@ 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) {
- if (j_id == primal_node_to_cell[i_cell]) {
- return w_rj(face_node_id, i_cell);
- }
- }
- 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_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_primal_node_id[dual_node_id] == face_node_id) {
- const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
- std::cout << " Clr -> " << Clr << "\n";
- std::cout << " w_rj -> " << weight_rj << "\n";
- S(i_id, j_id) -= weight_rj * alpha_face_id * (nil, Clr);
+
+ CellId dual_cell_id = face_dual_cell_id[face_id];
+
+ 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_primal_node_id[dual_node_id] == node_id) {
+ const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
+
+ const double a_ir = alpha_face_id * (nil, Clr);
+
+ for (size_t j_cell = 0; j_cell < primal_node_to_cell_matrix[node_id].size(); ++j_cell) {
+ CellId j_id = primal_node_to_cell_matrix[node_id][j_cell];
+ S(i_id, j_id) -= w_rj(node_id, j_cell) * a_ir;
}
}
- } // else {
- // std::cout << "*** ignore node " << face_node_id << '\n';
- // }
+ }
}
}
}
}
}
- // std::exit(0);
+
CellValue<double> fj =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id, mesh_data.xj());
@@ -353,14 +281,14 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, m_mesh->xr());
const CellValue<const double> primal_Vj = mesh_data.Vj();
+
CRSMatrix A{S};
Vector<double> b{m_mesh->numberOfCells()};
- // parallel_for(
- // m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { b[cell_id] = fj[cell_id] * primal_Vj[cell_id]; });
const auto& primal_cell_to_face_matrix = m_mesh->connectivity().cellToFaceMatrix();
for (CellId cell_id = 0; cell_id < m_mesh->numberOfCells(); ++cell_id) {
- b[cell_id] = fj[cell_id] * primal_Vj[cell_id];
+ 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 i_cell_face = 0; i_cell_face < primal_cell_to_face.size(); ++i_cell_face) {
@@ -376,7 +304,9 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
Assert(false, "cannot find cell's face");
return false;
}();
+
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) {
@@ -389,19 +319,13 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
CellId dual_cell_id = face_dual_cell_id[face_id];
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_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_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[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";
- }
}
}
}
@@ -410,19 +334,8 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
}
Vector<double> T{m_mesh->numberOfCells()};
- T = 1;
+ T = 0;
- 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()};
@@ -440,7 +353,10 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
{
VTKWriter vtk_writer("Temperature_" + std::to_string(Dimension), 0.01);
- vtk_writer.write(m_mesh, {NamedItemValue{"T", Temperature}, NamedItemValue{"Texact", Tj}}, 0,
+ vtk_writer.write(m_mesh,
+ {NamedItemValue{"T", Temperature}, NamedItemValue{"Vj", primal_Vj},
+ NamedItemValue{"Texact", Tj}},
+ 0,
true); // forces last output
}
}
--
GitLab