diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp
index 06fe956c09ea9e03fddb41a9abf6e5c170e5f4e6..9367249f71b20a049156aee7557d5c35cb48a736 100644
--- a/src/scheme/ScalarNodalScheme.cpp
+++ b/src/scheme/ScalarNodalScheme.cpp
@@ -386,8 +386,20 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
const unsigned int i_node = node_local_number_in_its_cell[i_cell];
const CellId cell_id = node_to_cell[i_cell];
const auto& cell_nodes = cell_to_node_matrix[cell_id];
- const NodeId prev_node_id = cell_to_node_matrix[cell_id][(i_node - 1) % cell_nodes.size()];
- const NodeId next_node_id = cell_to_node_matrix[cell_id][(i_node + 1) % cell_nodes.size()];
+ NodeId i_node_prev;
+ NodeId i_node_next;
+ if (i_node == 0) {
+ i_node_prev = cell_nodes.size() - 1;
+ } else {
+ i_node_prev = i_node - 1;
+ }
+ if (i_node == cell_nodes.size() - 1) {
+ i_node_next = 0;
+ } else {
+ i_node_next = i_node + 1;
+ }
+ const NodeId prev_node_id = cell_to_node_matrix[cell_id][i_node_prev];
+ const NodeId next_node_id = cell_to_node_matrix[cell_id][i_node_next];
if (is_boundary_node[prev_node_id] and is_boundary_node[next_node_id]) {
compute_node_is_corner[node_id] = true;
}
@@ -411,6 +423,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
}
const double norm_exterior_normal = l2Norm(compute_exterior_normal[node_id]);
compute_exterior_normal[node_id] *= 1. / norm_exterior_normal;
+ // std::cout << node_id << " " << compute_exterior_normal[node_id] << "\n";
}
}
return compute_exterior_normal;
@@ -451,7 +464,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
}
}
}
-
} else if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
const auto& face_list = bc.faceList();
const auto& value_list = bc.valueList();
@@ -480,6 +492,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
} else if ((not node_is_neumann[node_id]) && (node_is_dirichlet[node_id])) {
compute_node_boundary_values[node_id] /= sum_mes_l[node_id];
}
+ // std::cout << node_id << " " << compute_node_boundary_values[node_id] << "\n";
}
return compute_node_boundary_values;
}();
@@ -709,15 +722,17 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
}();
const Array<int> non_zeros{mesh->numberOfCells()};
- non_zeros.fill(Dimension); // Modif pour optimiser
- CRSMatrixDescriptor<double> S(mesh->numberOfCells(), mesh->numberOfCells(), non_zeros);
+ non_zeros.fill(4); // Modif pour optimiser
+ CRSMatrixDescriptor<double> S1(mesh->numberOfCells(), mesh->numberOfCells(), non_zeros);
+ CRSMatrixDescriptor<double> S2(mesh->numberOfCells(), mesh->numberOfCells(), non_zeros);
for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
const auto& node_to_cell = node_to_cell_matrix[node_id];
for (size_t i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) {
const CellId cell_id_j = node_to_cell[i_cell_j];
for (size_t i_cell_k = 0; i_cell_k < node_to_cell.size(); ++i_cell_k) {
- const CellId cell_id_k = node_to_cell[i_cell_k];
- S(cell_id_j, cell_id_k) = 0;
+ const CellId cell_id_k = node_to_cell[i_cell_k];
+ S1(cell_id_j, cell_id_k) = 0;
+ S2(cell_id_j, cell_id_k) = 0;
}
}
}
@@ -734,8 +749,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
const CellId cell_id_k = node_to_cell[i_cell_k];
const size_t i_node_k = node_local_number_in_its_cells[i_cell_k];
- S(cell_id_j, cell_id_k) +=
- dt / Vj[cell_id_j] * dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
+ S1(cell_id_j, cell_id_k) +=
+ dt * (1. - cn_coeff / 2.) / Vj[cell_id_j] *
+ dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
+ S2(cell_id_j, cell_id_k) -=
+ dt * (cn_coeff / 2.) / Vj[cell_id_j] *
+ dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
}
}
} else if ((node_is_neumann[node_id]) && (not node_is_corner[node_id]) && (not node_is_dirichlet[node_id])) {
@@ -752,8 +771,13 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
const CellId cell_id_k = node_to_cell[i_cell_k];
const size_t i_node_k = node_local_number_in_its_cells[i_cell_k];
- S(cell_id_j, cell_id_k) +=
- dt / Vj[cell_id_j] *
+ S1(cell_id_j, cell_id_k) +=
+ dt * (1. - cn_coeff / 2.) / Vj[cell_id_j] *
+ dot(kappar_invBetar[node_id] *
+ (Cjr(cell_id_k, i_node_k) - theta(cell_id_k, i_node_k) / sum_theta[node_id] * sum_Cjr[node_id]),
+ Q * Cjr(cell_id_j, i_node_j));
+ S2(cell_id_j, cell_id_k) -=
+ dt * (cn_coeff / 2.) / Vj[cell_id_j] *
dot(kappar_invBetar[node_id] *
(Cjr(cell_id_k, i_node_k) - theta(cell_id_k, i_node_k) / sum_theta[node_id] * sum_Cjr[node_id]),
Q * Cjr(cell_id_j, i_node_j));
@@ -768,8 +792,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
const CellId cell_id_k = node_to_cell[i_cell_k];
const size_t i_node_k = node_local_number_in_its_cells[i_cell_k];
- S(cell_id_j, cell_id_k) +=
- dt / Vj[cell_id_j] * dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
+ S1(cell_id_j, cell_id_k) +=
+ dt * (1. - cn_coeff / 2.) / Vj[cell_id_j] *
+ dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
+ S2(cell_id_j, cell_id_k) -=
+ dt * (cn_coeff / 2.) / Vj[cell_id_j] *
+ dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
}
}
} else if (node_is_dirichlet[node_id] && node_is_corner[node_id]) {
@@ -781,8 +809,11 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
const CellId cell_id_k = node_to_cell[i_cell_k];
const size_t i_node_k = node_local_number_in_its_cells[i_cell_k];
- S(cell_id_j, cell_id_k) +=
- dt / Vj[cell_id_j] *
+ S1(cell_id_j, cell_id_k) +=
+ dt * (1. - cn_coeff / 2.) / Vj[cell_id_j] *
+ dot(corner_kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
+ S2(cell_id_j, cell_id_k) -=
+ dt * (cn_coeff / 2.) / Vj[cell_id_j] *
dot(corner_kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j));
}
}
@@ -790,13 +821,22 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
}
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
- S(cell_id, cell_id) += 1;
+ S1(cell_id, cell_id) += 1;
+ S2(cell_id, cell_id) += 1;
}
- CellValue<const double> fj = f->cellValues();
- CellValue<const double> Ph = P->cellValues();
+ CellValue<const double> fj = f->cellValues();
+ CellValue<const double> fj_prev = f_prev->cellValues();
+ CellValue<const double> Pj = P->cellValues();
- CRSMatrix A{S.getCRSMatrix()};
+ Vector<double> Ph{mesh->numberOfCells()};
+ Ph = zero;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ Ph[cell_id] = Pj[cell_id];
+ };
+
+ CRSMatrix A1{S1.getCRSMatrix()};
+ CRSMatrix A2{S2.getCRSMatrix()};
Vector<double> b{mesh->numberOfCells()};
b = zero;
@@ -804,6 +844,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
b[cell_id] = fj[cell_id];
};
+ Vector<double> b_prev{mesh->numberOfCells()};
+ b_prev = zero;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ b_prev[cell_id] = fj_prev[cell_id];
+ };
+
for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
const auto& node_to_cell = node_to_cell_matrix[node_id];
const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
@@ -821,6 +867,10 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
(1. / (sum_theta[node_id] * l2Norm(node_kappar[node_id] * n)) *
dot(kappar_invBetar[node_id] * sum_Cjr[node_id], Q * Cjr(cell_id, i_node)) +
dot(Cjr(cell_id, i_node), n));
+ b_prev[cell_id] += 1. / Vj[cell_id] * node_boundary_prev_values[node_id] *
+ (1. / (sum_theta[node_id] * l2Norm(node_kappar[node_id] * n)) *
+ dot(kappar_invBetar[node_id] * sum_Cjr[node_id], Q * Cjr(cell_id, i_node)) +
+ dot(Cjr(cell_id, i_node), n));
}
} else if (node_is_corner[node_id]) {
const auto& node_to_face = node_to_face_matrix[node_id];
@@ -831,6 +881,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
sum_mes_l += mes_l[face_id];
}
b[cell_id] += 1. / Vj[cell_id] * 0.5 * node_boundary_values[node_id] * sum_mes_l;
+ b_prev[cell_id] += 1. / Vj[cell_id] * 0.5 * node_boundary_prev_values[node_id] * sum_mes_l;
}
} else if (node_is_dirichlet[node_id]) {
if (not node_is_corner[node_id]) {
@@ -845,6 +896,10 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
b[cell_id_j] += 1. / Vj[cell_id_j] *
dot(node_boundary_values[node_id] * kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
Cjr(cell_id_j, i_node_j));
+ b_prev[cell_id_j] +=
+ 1. / Vj[cell_id_j] *
+ dot(node_boundary_prev_values[node_id] * kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
+ Cjr(cell_id_j, i_node_j));
}
}
@@ -861,6 +916,10 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
1. / Vj[cell_id_j] *
dot(node_boundary_values[node_id] * corner_kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
Cjr(cell_id_j, i_node_j));
+ b_prev[cell_id_j] +=
+ 1. / Vj[cell_id_j] *
+ dot(node_boundary_prev_values[node_id] * corner_kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
+ Cjr(cell_id_j, i_node_j));
}
}
}
@@ -870,14 +929,19 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
Vector<double> T{mesh->numberOfCells()};
T = zero;
+ // Array<double> Ph2 = Ph;
+ Vector<double> A2P = A2 * Ph;
+
Vector<double> B{mesh->numberOfCells()};
parallel_for(
- mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { B[cell_id] = dt * b[cell_id] + Ph[cell_id]; });
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ B[cell_id] = dt * ((1. - cn_coeff / 2.) * b[cell_id] + cn_coeff / 2. * b_prev[cell_id]) + A2P[cell_id];
+ });
// std::cout << "g = " << node_boundary_values << "\n";
LinearSolver solver;
- solver.solveLocalSystem(A, T, B);
+ solver.solveLocalSystem(A1, T, B);
m_solution = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
auto& solution = *m_solution;