diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp index 72936fe9222f2a9ac273f2190f37d79cf1513b39..846a4f3fa98d98ac99838007f16203c3d88cda80 100644 --- a/src/scheme/ScalarNodalScheme.cpp +++ b/src/scheme/ScalarNodalScheme.cpp @@ -586,7 +586,7 @@ 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) += + S(cell_id_j, cell_id_k) -= dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j)); } } @@ -604,10 +604,10 @@ 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) += + S(cell_id_j, cell_id_k) -= 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_k, i_node_j)); + Q * Cjr(cell_id_j, i_node_j)); } } } else if ((node_is_dirichlet[node_id]) && (not node_is_corner[node_id])) { @@ -619,7 +619,7 @@ 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) += + S(cell_id_j, cell_id_k) -= dot(kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j)); } } @@ -632,7 +632,7 @@ 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) += + S(cell_id_j, cell_id_k) -= dot(corner_kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j)); } } @@ -641,7 +641,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand for (CellId cell_id_j = 0; cell_id_j < mesh->numberOfCells(); ++cell_id_j) { for (CellId cell_id_k = 0; cell_id_k < mesh->numberOfCells(); ++cell_id_k) { - S(cell_id_j, cell_id_k) *= dt / Vj[cell_id_j]; + S(cell_id_j, cell_id_k) = -dt / Vj[cell_id_j] * S(cell_id_j, cell_id_k); } }; @@ -720,23 +720,58 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand } } } - }; + } + }; + + // TESTS // + NodeValue<TinyVector<Dimension>> gradPr{mesh->connectivity()}; + NodeValue<double> Pr{mesh->connectivity()}; + for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) { + gradPr[node_id] = zero; + 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); + if (not is_boundary_node[node_id]) { + for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) { + const CellId cell_id = node_to_cell[i_cell]; + const size_t i_node = node_local_number_in_its_cells[i_cell]; + gradPr[node_id] += Ph[cell_id] * Cjr(cell_id, i_node); + } + gradPr[node_id] = -inverse(node_betar[node_id]) * gradPr[node_id]; + std::cout << "node_id : " << node_id << "; gradPr = " << gradPr[node_id] << "\n"; + } else if (is_boundary_node[node_id] and not node_is_corner[node_id]) { + Pr[node_id] = 0; + for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) { + const CellId cell_id = node_to_cell[i_cell]; + const size_t i_node = node_local_number_in_its_cells[i_cell]; + Pr[node_id] += theta(cell_id, i_node) * Ph[cell_id]; + } + Pr[node_id] += node_boundary_values[node_id] / l2Norm(node_kappar[node_id] * exterior_normal[node_id]); + Pr[node_id] *= 1. / sum_theta[node_id]; + std::cout << "node_id :" << node_id << "; Pr = " << Pr[node_id] << "\n"; + } + }; + // END TESTS// - Vector<double> T{mesh->numberOfCells()}; - T = zero; + Vector<double> T{mesh->numberOfCells()}; + T = zero; - Vector<double> B{mesh->numberOfCells()}; - parallel_for( - mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { B[cell_id] = dt * b[cell_id] + Ph[cell_id]; }); + Vector<double> B{mesh->numberOfCells()}; + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { B[cell_id] = dt * b[cell_id] + Ph[cell_id]; }); - LinearSolver solver; - solver.solveLocalSystem(A, T, B); + Vector<double> Pvec{mesh->numberOfCells()}; + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { Pvec[cell_id] = Ph[cell_id]; }); + std::cout << "B = " << B << "\n" + << "AP = " << A * Pvec << "\n"; - m_solution = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh); - auto& solution = *m_solution; - parallel_for( - mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { solution[cell_id] = T[cell_id]; }); - } + LinearSolver solver; + solver.solveLocalSystem(A, T, B); + + m_solution = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh); + auto& solution = *m_solution; + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { solution[cell_id] = T[cell_id]; }); } } };