Skip to content
Snippets Groups Projects
Commit d2972e3d authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Add a bunch of small optimizations

parent 2814bf23
No related branches found
No related tags found
No related merge requests found
......@@ -7,6 +7,8 @@
#include <utils/PugsUtils.hpp>
#include <utils/RandomEngine.hpp>
#include <utils/Timer.hpp>
int
main(int argc, char* argv[])
{
......@@ -19,8 +21,12 @@ main(int argc, char* argv[])
DualConnectivityManager::create();
DualMeshManager::create();
Timer elapsed_time;
parser(filename);
elapsed_time.pause();
DualMeshManager::destroy();
DualConnectivityManager::destroy();
MeshDataManager::destroy();
......@@ -30,5 +36,7 @@ main(int argc, char* argv[])
finalize();
std::cout << "-- Elapsed time: " << elapsed_time << '\n';
return 0;
}
......@@ -23,6 +23,7 @@ Timer solver_t;
Timer get_A_t;
Timer getF_t;
Timer getGradF_t;
Timer HJ_A_t;
static int count_getA = 0;
static int count_getF = 0;
......@@ -108,6 +109,10 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
CellValue<const double> m_inv_gamma;
CellValue<const double> m_g_1_exp_S_Cv_inv_g;
CellValue<const double> m_tau_iter;
CellValue<const double> m_tau;
CellValue<const double> m_Mj;
NodeValuePerCell<const Rdxd> m_Ajr;
NodeValue<const Rdxd> m_inv_Ar;
......@@ -527,14 +532,10 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
Vector<double>
_getGradJ(const Vector<double>& Un,
const Vector<double>& Uk,
const DiscreteScalarFunction& rho,
const DiscreteScalarFunction& p,
const DiscreteScalarFunction& pi,
const DiscreteVectorFunction& u,
const DiscreteScalarFunction& gamma,
const DiscreteScalarFunction& Cv,
const DiscreteScalarFunction& entropy,
const double dt) const
const double dt)
{
count_getF++;
static bool getF_is_started = false;
......@@ -544,7 +545,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
}
getF_t.start();
CellValue<const double> tau_iter = [&]() {
m_tau_iter = [&]() {
CellValue<double> computed_tau_iter(m_mesh.connectivity());
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
......@@ -558,16 +559,6 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
getF_t.pause();
const CellValue<const double> Mj = [&]() {
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(m_mesh);
const CellValue<const double>& Vj = mesh_data.Vj();
CellValue<double> computed_Mj(m_mesh.connectivity());
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { computed_Mj[j] = rho[j] * Vj[j]; });
return computed_Mj;
}();
Vector<double> computed_gradJ{(Dimension + 1) * m_number_of_implicit_cells};
computed_gradJ = zero;
const auto& node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
......@@ -620,13 +611,14 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
}
}
const double inv_dt = 1 / dt;
for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
if (m_is_implicit_cell[cell_id]) {
size_t j_index_p = mapP(cell_id);
computed_gradJ[j_index_p] += (Mj[cell_id] / dt) * (tau_iter[cell_id] - (1. / rho[cell_id]));
computed_gradJ[j_index_p] += (inv_dt * m_Mj[cell_id]) * (m_tau_iter[cell_id] - m_tau[cell_id]);
for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
size_t j_index_u = mapU(i_dimension, cell_id);
computed_gradJ[j_index_u] += (Mj[cell_id] / dt) * (Uk[j_index_u] - Un[j_index_u]);
computed_gradJ[j_index_u] += (inv_dt * m_Mj[cell_id]) * (Uk[j_index_u] - Un[j_index_u]);
}
}
}
......@@ -815,16 +807,12 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
_getF(const CRSMatrix<double, int>& A,
const Vector<double>& Un,
const Vector<double>& Uk,
const DiscreteScalarFunction& rho,
const DiscreteScalarFunction& p,
const DiscreteScalarFunction& pi,
const DiscreteVectorFunction& u,
const DiscreteScalarFunction& gamma,
const DiscreteScalarFunction& Cv,
const DiscreteScalarFunction& entropy,
const double dt) const
const double dt)
{
Vector<double> gradJ = this->_getGradJ(Un, Uk, rho, p, pi, u, gamma, Cv, entropy, dt);
Vector<double> gradJ = this->_getGradJ(Un, Uk, p, pi, u, dt);
// std::cout << "gradJ" << gradJ << '\n';
// CRSMatrix A_crs{A.getCRSMatrix()};
......@@ -839,25 +827,8 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
}
CRSMatrixDescriptor<double>
_getHessianJ(const DiscreteScalarFunction& rho,
const double dt,
const DiscreteScalarFunction& gamma,
const DiscreteScalarFunction& entropy,
const DiscreteScalarFunction& pi,
const DiscreteScalarFunction& Cv,
const Vector<double>& Uk)
_getHessianJ(const double dt, const DiscreteScalarFunction& pi, const Vector<double>& Uk)
{
const CellValue<const double> Mj = [&]() {
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(m_mesh);
const CellValue<const double>& Vj = mesh_data.Vj();
CellValue<double> computed_Mj(m_mesh.connectivity());
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { computed_Mj[j] = rho[j] * Vj[j]; });
return computed_Mj;
}();
auto is_boundary_node = m_mesh.connectivity().isBoundaryNode();
Array<int> non_zeros{(Dimension + 1) * m_number_of_implicit_cells};
......@@ -915,22 +886,25 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
getGradF_is_started = true;
getGradF_t.stop();
}
getGradF_t.start();
const double inv_dt = 1 / dt;
for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
if (m_is_implicit_cell[cell_id]) {
const int i_index_p = mapP(cell_id);
Hess_J(i_index_p, i_index_p) = (Mj[cell_id] / dt) * m_inv_gamma[cell_id] * m_g_1_exp_S_Cv_inv_g[cell_id] *
std::pow(-Uk[i_index_p] + pi[cell_id], -1 - m_inv_gamma[cell_id]);
Hess_J(i_index_p, i_index_p) =
(inv_dt * m_Mj[cell_id]) * m_inv_gamma[cell_id] * m_tau_iter[cell_id] / (-Uk[i_index_p] + pi[cell_id]);
for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
const int i_index_u = mapU(i_dimension, cell_id);
Hess_J(i_index_u, i_index_u) = (Mj[cell_id] / dt);
Hess_J(i_index_u, i_index_u) = inv_dt * m_Mj[cell_id];
}
}
}
getGradF_t.start();
for (NodeId node_id = 0; node_id < m_mesh.numberOfNodes(); ++node_id) {
if (not is_boundary_node[node_id]) {
const auto& inv_Ar = m_inv_Ar[node_id];
const auto& node_cell = node_to_cell_matrix[node_id];
const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
......@@ -939,13 +913,15 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
if (m_is_implicit_cell[id_cell_i]) {
const size_t node_nb_in_i = node_local_number_in_its_cells[cell_i];
const int i_index_p = mapP(id_cell_i);
const auto invArCir = inv_Ar * Cjr(id_cell_i, node_nb_in_i);
for (size_t cell_j = 0; cell_j < node_cell.size(); ++cell_j) {
const CellId id_cell_j = node_cell[cell_j];
if (m_is_implicit_cell[id_cell_j]) {
const size_t node_nb_in_j = node_local_number_in_its_cells[cell_j];
const int j_index_p = mapP(id_cell_j);
Hess_J(i_index_p, j_index_p) +=
dot(Cjr(id_cell_i, node_nb_in_i), m_inv_Ar[node_id] * Cjr(id_cell_j, node_nb_in_j));
Hess_J(i_index_p, j_index_p) += dot(invArCir, Cjr(id_cell_j, node_nb_in_j));
}
}
}
......@@ -955,9 +931,11 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
if (m_is_implicit_cell[id_cell_i]) {
const size_t node_nb_in_i = node_local_number_in_its_cells[cell_i];
const auto Ajr = m_Ajr(id_cell_i, node_nb_in_i);
for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
const int i_index_u = mapU(i_dimension, id_cell_i);
Hess_J(i_index_u, i_index_u) += m_Ajr(id_cell_i, node_nb_in_i)(i_dimension, i_dimension);
Hess_J(i_index_u, i_index_u) += Ajr(i_dimension, i_dimension);
}
}
}
......@@ -967,18 +945,21 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
if (m_is_implicit_cell[id_cell_i]) {
const size_t node_nb_in_i = node_local_number_in_its_cells[cell_i];
for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
const int i_index_u = mapU(i_dimension, id_cell_i);
const auto Air_invAr = m_Ajr(id_cell_i, node_nb_in_i) * inv_Ar;
for (size_t cell_j = 0; cell_j < node_cell.size(); ++cell_j) {
const CellId id_cell_j = node_cell[cell_j];
if (m_is_implicit_cell[id_cell_j]) {
const size_t node_nb_in_j = node_local_number_in_its_cells[cell_j];
const auto Air_invAr_Ajr = Air_invAr * m_Ajr(id_cell_j, node_nb_in_j);
for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
const int i_index_u = mapU(i_dimension, id_cell_i);
for (size_t j_dimension = 0; j_dimension < Dimension; ++j_dimension) {
const int j_index_u = mapU(j_dimension, id_cell_j);
Hess_J(i_index_u, j_index_u) += -(m_Ajr(id_cell_i, node_nb_in_i) * m_inv_Ar[node_id] *
m_Ajr(id_cell_j, node_nb_in_j))(i_dimension, j_dimension);
Hess_J(i_index_u, j_index_u) -= Air_invAr_Ajr(i_dimension, j_dimension);
}
}
}
......@@ -987,6 +968,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
}
}
}
getGradF_t.pause();
for (const auto& boundary_condition : m_boundary_condition_list) {
std::visit(
......@@ -1159,28 +1141,33 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
},
boundary_condition);
}
getGradF_t.pause();
return Hess_J;
}
CRSMatrix<double>
_getGradientF(const CRSMatrix<double, int>& A,
const Vector<double>& Uk,
const DiscreteScalarFunction& rho,
const DiscreteScalarFunction& gamma,
const DiscreteScalarFunction& pi,
const DiscreteScalarFunction& Cv,
const DiscreteScalarFunction& entropy,
const double dt)
{
CRSMatrixDescriptor<double> Hess_J = this->_getHessianJ(rho, dt, gamma, entropy, pi, Cv, Uk);
static bool HJ_A_is_started = false;
if (not HJ_A_is_started) {
HJ_A_is_started = true;
HJ_A_t.stop();
}
CRSMatrix Hess_J_crs = Hess_J.getCRSMatrix();
HJ_A_t.start();
CRSMatrixDescriptor<double> Hess_J = this->_getHessianJ(dt, pi, Uk);
// std::cout << "hessJ" << Hess_J_crs << '\n';
// std::exit(0);
// CRSMatrix A_crs{A.getCRSMatrix()};
CRSMatrix Hess_J_crs = Hess_J.getCRSMatrix();
CRSMatrix gradient_f = Hess_J_crs - A;
HJ_A_t.pause();
return gradient_f;
}
......@@ -1264,14 +1251,9 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
}
void
_computeGradJU_AU(const DiscreteScalarFunction& rho,
const DiscreteScalarFunction& c,
const DiscreteVectorFunction& u,
_computeGradJU_AU(const DiscreteVectorFunction& u,
const DiscreteScalarFunction& p,
const DiscreteScalarFunction& pi,
const DiscreteScalarFunction& gamma,
const DiscreteScalarFunction& Cv,
const DiscreteScalarFunction& entropy,
const double& dt,
const NodeValuePerCell<const Rd>& Cjr_half)
{
......@@ -1369,9 +1351,9 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
// std::cout << "iteration=" << nb_iter << '\n';
nb_iter++;
Vector<double> f = this->_getF(A, Un, Uk, rho, p, pi, u, gamma, Cv, entropy, dt);
Vector<double> f = this->_getF(A, Un, Uk, p, pi, u, dt);
// std::exit(0);
CRSMatrix<double> gradient_f = this->_getGradientF(A, Uk, rho, gamma, pi, Cv, entropy, dt);
CRSMatrix<double> gradient_f = this->_getGradientF(A, Uk, pi, dt);
Vector<double> sol{Un.size()};
static bool solver_is_init = false;
......@@ -1510,6 +1492,23 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
CellValue<double> new_E = copy(E.cellValues());
std::shared_ptr<const MeshType> new_mesh;
m_Mj = [&]() {
const CellValue<const double>& Vj = mesh_data.Vj();
CellValue<double> computed_Mj(m_mesh.connectivity());
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { computed_Mj[j] = rho[j] * Vj[j]; });
return computed_Mj;
}();
m_tau = [&]() {
CellValue<double> computed_tau(m_mesh.connectivity());
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { computed_tau[j] = 1 / rho[j]; });
return computed_tau;
}();
m_Ajr = this->_computeAjr(rho, c);
// for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
......@@ -1552,7 +1551,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
new_u = copy(u.cellValues());
new_E = copy(E.cellValues());
_computeGradJU_AU(rho, c, u, p, pi, gamma, Cv, entropy, dt, Cjr_half);
_computeGradJU_AU(u, p, pi, dt, Cjr_half);
// u_j+1/2
NodeValue<const Rd> ur = [&]() {
......@@ -1718,7 +1717,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
momentum_fluxes += Fjr(j, R);
energy_fluxes += dot(Fjr(j, R), ur[r]);
}
const double dt_over_Mj = dt / (rho[j] * Vj[j]);
const double dt_over_Mj = dt / m_Mj[j];
new_u[j] -= dt_over_Mj * momentum_fluxes;
new_E[j] -= dt_over_Mj * energy_fluxes;
});
......@@ -1750,17 +1749,17 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
const NodeId r = cell_nodes[R];
new_tau_fluxes += dot(Cjr_half(j, R), ur[r]);
}
new_tau[j] = (1. / rho[j]) + (dt / (rho[j] * Vj[j])) * new_tau_fluxes;
new_tau[j] = m_tau[j] + (dt / m_Mj[j]) * new_tau_fluxes;
// std::cout << "new_tau(" << j << ")=" << new_tau[j] << '\n';
}
double sum_tau_rho = 0;
// double sum_tau_rho = 0;
max_tau_error = 0;
for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
// const auto& cell_nodes = cell_to_node_matrix[j];
// V_j^n+1/tau_j^n+1 - M_j
sum_tau_rho += std::abs(new_tau[j] - 1. / new_rho[j]);
// sum_tau_rho += std::abs(new_tau[j] - 1. / new_rho[j]);
max_tau_error = std::max(max_tau_error, std::abs(new_tau[j] - 1. / new_rho[j]));
}
// std::cout << "sum_tau_rho =" << sum_tau_rho << '\n';
......@@ -1783,7 +1782,8 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
implicit_t.pause();
std::cout << "getA=" << get_A_t << "(" << count_getA << ")"
<< " getF=" << getF_t << "(" << count_getF << ")"
<< " getGradF=" << getGradF_t << "(" << count_getGradF << ")" << '\n';
<< " getGradF=" << getGradF_t << "(" << count_getGradF << ")"
<< " HJ_A=" << HJ_A_t << '\n';
std::cout << "solver=" << solver_t << " total= " << implicit_t << '\n';
return {new_mesh, std::make_shared<DiscreteScalarFunction>(new_mesh, new_rho),
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment