From 5ba134b14938d6a8f4559936cdc55d45877d8112 Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Fri, 22 Nov 2024 13:54:00 +0100
Subject: [PATCH] Add tbb loop in Adjoint. Add doc. Refactoring.
---
dart/api/core.py | 6 +++---
dart/src/wAdjoint.cpp | 27 +++++++++++++--------------
dart/src/wAdjoint.h | 3 ++-
dart/src/wBlowing.cpp | 19 -------------------
dart/src/wBlowing.h | 4 ++--
dart/src/wBlowingResidual.cpp | 21 +++++++++++++++------
dart/src/wSolver.h | 2 +-
7 files changed, 36 insertions(+), 46 deletions(-)
diff --git a/dart/api/core.py b/dart/api/core.py
index 6476b4a..7cb61f7 100644
--- a/dart/api/core.py
+++ b/dart/api/core.py
@@ -219,12 +219,12 @@ def init_dart(cfg, scenario='aerodynamic', task='analysis', viscous=False):
linsol = tbox.Pardiso()
elif cfg['LSolver'] == 'MUMPS':
linsol = tbox.Mumps()
+ elif cfg['LSolver'] == 'SparseLU':
+ linsol = tbox.SparseLu()
elif cfg['LSolver'] == 'GMRES':
linsol = tbox.Gmres(cfg.get('G_fill', 2), 1e-6, cfg.get('G_restart', 50), cfg.get('G_tol', 1e-5), 200)
- elif cfg['LSolver'] == 'SparseLu':
- linsol = tbox.SparseLu()
else:
- raise RuntimeError('Available linear solvers: PARDISO, MUMPS or GMRES, but ' + cfg['LSolver'] + ' was given!\n')
+ raise RuntimeError('Available linear solvers: PARDISO, MUMPS, SparseLU or GMRES, but ' + cfg['LSolver'] + ' was given!\n')
# initialize the nonlinear (outer) solver
_sol = dart.Newton(_pbl, linsol, rTol=cfg['Rel_tol'], aTol=cfg['Abs_tol'], mIt=cfg['Max_it'], nthrds=nthrd, vrb=verb)
diff --git a/dart/src/wAdjoint.cpp b/dart/src/wAdjoint.cpp
index 367921e..030bbbb 100644
--- a/dart/src/wAdjoint.cpp
+++ b/dart/src/wAdjoint.cpp
@@ -564,25 +564,24 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
// Farfield B.C. are not taken into account because they have no influence on the solution
// Blowing boundary
- if (sol->pbl->bBCs.size() > 0)
+ for (auto blow : sol->pbl->bBCs)
{
- for (auto b : sol->pbl->bBCs)
- {
- for (auto p : b->uE)
+ tbb::parallel_for_each(blow->uE.begin(), blow->uE.end(), [&](std::pair<Element *, F0ElC *> p) {
+ // Build elementary matrices
+ Eigen::MatrixXd A = BlowingResidual::buildGradientMesh(*p.first, sol->phi, *p.second, sol->pbl->nDim);
+ // Assembly
+ tbb::spin_mutex::scoped_lock lock(mutex);
+ for (size_t i = 0; i < p.first->nodes.size(); ++i)
{
- Eigen::MatrixXd A = BlowingResidual::buildGradientMesh(*p.first, sol->phi, *p.second, sol->pbl->nDim);
- for (size_t i = 0; i < p.first->nodes.size(); ++i)
+ size_t rowi = p.first->nodes[i]->row;
+ for (size_t j = 0; j < p.first->nodes.size(); ++j)
{
- size_t rowi = p.first->nodes[i]->row;
- for (size_t j = 0; j < p.first->nodes.size(); ++j)
- {
- int rowj = p.first->nodes[j]->row;
- for (int n = 0; n < sol->pbl->nDim; ++n)
- T.push_back(Eigen::Triplet<double>(sol->rows[rowi], sol->pbl->nDim * rowj + n, A(i, sol->pbl->nDim * j + n)));
- }
+ int rowj = p.first->nodes[j]->row;
+ for (int n = 0; n < sol->pbl->nDim; ++n)
+ T.push_back(Eigen::Triplet<double>(sol->rows[rowi], sol->pbl->nDim * rowj + n, A(i, sol->pbl->nDim * j + n)));
}
}
- }
+ });
}
// Wake
diff --git a/dart/src/wAdjoint.h b/dart/src/wAdjoint.h
index 79f4eee..5c1e1cd 100644
--- a/dart/src/wAdjoint.h
+++ b/dart/src/wAdjoint.h
@@ -87,12 +87,13 @@ public:
std::vector<double> computeGradientCoefficientsAoa();
std::vector<std::vector<double>> computeGradientCoefficientsMesh();
+#ifndef SWIG
+ // Getters for partials (required by VII-BLASTER)
Eigen::SparseMatrix<double, Eigen::RowMajor> getRu_U() const { return dRu_U; }
Eigen::VectorXd getRu_A() const { return dRu_A; }
Eigen::SparseMatrix<double, Eigen::RowMajor> getRu_X() const { return dRu_X; }
Eigen::SparseMatrix<double, Eigen::RowMajor> getRx_X() const { return dRx_X; }
-#ifndef SWIG
virtual void write(std::ostream &out) const override;
#endif
diff --git a/dart/src/wBlowing.cpp b/dart/src/wBlowing.cpp
index ef2c86a..c4163f3 100644
--- a/dart/src/wBlowing.cpp
+++ b/dart/src/wBlowing.cpp
@@ -59,25 +59,6 @@ void Blowing::create()
}
}
-/**
- * @brief Set blowing velocity
- */
-void Blowing::setU(size_t i, double ue)
-{
- uE[i].second->update(ue);
-}
-
-/**
- * @brief Get blowing velocity
- */
-double Blowing::getU(size_t i) const
-{
- // Return a constant value
- tbox::Element *e = uE[i].first;
- std::vector<double> phi(0, 0.);
- return uE[i].second->eval(*e, phi, 0);
-}
-
void Blowing::write(std::ostream &out) const
{
out << " Blowing boundary condition on " << *tag << std::endl;
diff --git a/dart/src/wBlowing.h b/dart/src/wBlowing.h
index 2a0cd7d..93b27c7 100644
--- a/dart/src/wBlowing.h
+++ b/dart/src/wBlowing.h
@@ -41,8 +41,8 @@ public:
Blowing(std::shared_ptr<tbox::MshData> _msh, std::string const &names);
virtual ~Blowing();
- void setU(size_t i, double ue);
- double getU(size_t i) const;
+ void setU(size_t i, double ue) { uE[i].second->update(ue); }
+ double getU(size_t i) const { return uE[i].second->eval(*uE[i].first, std::vector<double>(), 0); }
#ifndef SWIG
virtual void write(std::ostream &out) const override;
diff --git a/dart/src/wBlowingResidual.cpp b/dart/src/wBlowingResidual.cpp
index b8b4ea3..0c07b67 100644
--- a/dart/src/wBlowingResidual.cpp
+++ b/dart/src/wBlowingResidual.cpp
@@ -27,7 +27,7 @@ using namespace dart;
/**
* @brief Build the residual vector, on one boundary element
*
- * b = \int psi * vn dV
+ * b = \int psi * vn dS
*/
Eigen::VectorXd BlowingResidual::build(Element const &e, std::vector<double> const &phi, F0El const &f)
{
@@ -47,11 +47,16 @@ Eigen::VectorXd BlowingResidual::build(Element const &e, std::vector<double> con
/**
* @brief Build the gradient of the residual vector with respect to the blowing velocity, on one blowing boundary element
+ *
+ * b = \int psi dS
*/
Eigen::VectorXd BlowingResidual::buildGradientBlowing(tbox::Element const &e)
{
+ // Get pre-computed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
+
+ // Build velocity gradient
Eigen::VectorXd b = Eigen::VectorXd::Zero(e.nodes.size());
for (size_t k = 0; k < gauss.getN(); ++k)
b += cache.getSf(k) * gauss.getW(k) * e.getDetJ(k);
@@ -60,19 +65,23 @@ Eigen::VectorXd BlowingResidual::buildGradientBlowing(tbox::Element const &e)
/**
* @brief Build the gradient of the residual vector with respect to the mesh, on one blowing boundary element
+ *
+ * A = \int psi * vn ddS
*/
Eigen::MatrixXd BlowingResidual::buildGradientMesh(tbox::Element const &e, std::vector<double> const &phi, F0El const &f, int const nDim)
{
- Eigen::MatrixXd A = Eigen::MatrixXd::Zero(e.nodes.size(), nDim * e.nodes.size());
+ // Get pre-computed values
Cache &cache = e.getVCache();
Gauss &gauss = cache.getVGauss();
- // Blowing velocity
+
+ // Build velocity gradient
double vn = f.eval(e, phi, 0);
+ Eigen::MatrixXd A = Eigen::MatrixXd::Zero(e.nodes.size(), nDim * e.nodes.size());
for (size_t k = 0; k < gauss.getN(); ++k)
{
- std::vector<double> gradJ = e.getGradDetJ(k);
- for (size_t i = 0; i < gradJ.size(); i++)
- A.col(i) += cache.getSf(k) * vn * gauss.getW(k) * gradJ[i];
+ std::vector<double> dDetJ = e.getGradDetJ(k);
+ for (size_t i = 0; i < dDetJ.size(); i++)
+ A.col(i) += cache.getSf(k) * vn * gauss.getW(k) * dDetJ[i];
}
return A;
}
diff --git a/dart/src/wSolver.h b/dart/src/wSolver.h
index 62a2968..a988bca 100644
--- a/dart/src/wSolver.h
+++ b/dart/src/wSolver.h
@@ -80,7 +80,7 @@ public:
virtual Status run();
void save(tbox::MshExport *mshWriter, std::string const &suffix = "");
- int getRow(int i) const { return rows[i]; };
+ int getRow(size_t i) const { return rows[i]; };
protected:
void computeFlow();
--
GitLab