From cf404acd2931278556201f2b6bea865e3dac81aa Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Thu, 24 Mar 2022 17:38:42 +0100
Subject: [PATCH] Use linear solvers from direct methods in Adjoint. Minor
changes.
---
dart/_src/dartw.i | 2 +-
dart/api/core.py | 21 ++++++---------------
dart/api/internal/trim.py | 4 ++--
dart/api/mphys.py | 2 +-
dart/benchmark/onera.py | 6 +-----
dart/default.py | 7 +++----
dart/src/wAdjoint.cpp | 32 +++++++++++++++++++++++---------
dart/src/wAdjoint.h | 9 +++++----
dart/src/wSolver.cpp | 4 ++--
ext/amfe | 2 +-
10 files changed, 45 insertions(+), 44 deletions(-)
diff --git a/dart/_src/dartw.i b/dart/_src/dartw.i
index e16e839..6847f36 100644
--- a/dart/_src/dartw.i
+++ b/dart/_src/dartw.i
@@ -70,8 +70,8 @@ threads="1"
%immutable dart::Body::Cm;
%immutable dart::Body::nLoads;
%include "wBody.h"
-%include "wBlowing.h"
%warnfilter(+509);
+%include "wBlowing.h"
%immutable dart::Problem::msh; // read-only variable (no setter)
%immutable dart::Problem::nDim;
diff --git a/dart/api/core.py b/dart/api/core.py
index e729b8d..69d1751 100644
--- a/dart/api/core.py
+++ b/dart/api/core.py
@@ -146,7 +146,7 @@ def initDart(params, scenario='aerodynamic', task='analysis'):
# Mesh morpher creation
if scenario == 'aerostructural' or task == 'optimization':
- _mrf = tbox.MshDeform(_msh, _dim, nthrds=nthrd)
+ _mrf = tbox.MshDeform(_msh, tbox.Gmres(1, 1e-6, 30, 1e-8), _dim, nthrds=nthrd)
_mrf.setField(params['Fluid'])
_mrf.addFixed(params['Farfield'])
if _dim == 2:
@@ -213,10 +213,10 @@ def initDart(params, scenario='aerodynamic', task='analysis'):
elif params['LSolver'] == 'MUMPS':
linsol = tbox.Mumps()
elif params['LSolver'] == 'GMRES':
- linsol = tbox.Gmres()
- linsol.setFillFactor(params['G_fill']) if 'G_fill' in params else linsol.setFillFactor(2)
- linsol.setRestart(params['G_restart']) if 'G_restart' in params else linsol.setRestart(50)
- linsol.setTolerance(params['G_tol']) if 'G_tol' in params else linsol.setTolerance(1e-5)
+ gfil = params['G_fill'] if 'G_fill' in params else 2
+ grst = params['G_restart'] if 'G_restart' in params else 50
+ gtol = params['G_tol'] if 'G_tol' in params else 1e-5
+ linsol = tbox.Gmres(gfil, 1e-6, grst, gtol)
else:
raise RuntimeError('Available linear solvers: PARDISO, MUMPS or GMRES, but ' + params['LSolver'] + ' was given!\n')
# initialize the nonlinear (outer) solver
@@ -224,16 +224,7 @@ def initDart(params, scenario='aerodynamic', task='analysis'):
# Adjoint (reverse) solver creation
if task == 'optimization':
- if params['LSolver'] == 'PARDISO':
- alinsol = tbox.Pardiso()
- elif params['LSolver'] == 'MUMPS':
- alinsol = tbox.Mumps()
- else:
- alinsol = tbox.Gmres()
- alinsol.setFillFactor(params['G_fill']) if 'G_fill' in params else alinsol.setFillFactor(2)
- alinsol.setRestart(params['G_restart']) if 'G_restart' in params else alinsol.setRestart(50)
- alinsol.setTolerance(1e-12)
- _adj = dart.Adjoint(_sol, _mrf, alinsol)
+ _adj = dart.Adjoint(_sol, _mrf)
else:
_adj = None
diff --git a/dart/api/internal/trim.py b/dart/api/internal/trim.py
index be1c330..74c0025 100644
--- a/dart/api/internal/trim.py
+++ b/dart/api/internal/trim.py
@@ -90,8 +90,8 @@ class Trim:
self.pbl.update(self.alpha)
# run the solver
self.tms["solver"].start()
- success = self.sol.run()
- if not success:
+ status = self.sol.run()
+ if status >= 2:
raise RuntimeError('Flow solver diverged!\n')
self.tms["solver"].stop()
# update slope
diff --git a/dart/api/mphys.py b/dart/api/mphys.py
index 58200c0..56fcf9a 100644
--- a/dart/api/mphys.py
+++ b/dart/api/mphys.py
@@ -372,7 +372,7 @@ class DartCoefficients(om.ExplicitComponent):
"""
# NB: inputs already up-to-date because DartSolver has already been run
# Write to disk
- if self.options['morph'] and str(type(self.options['wrtr'])) == '<class \'tboxw.GmshExport\'>':
+ if self.options['morph'] and self.options['wrtr'].type() == 1:
if self.iscn == 0:
self.options['wrtr'].save(self.sol.pbl.msh.name)
else:
diff --git a/dart/benchmark/onera.py b/dart/benchmark/onera.py
index 75209e2..b0f9c21 100644
--- a/dart/benchmark/onera.py
+++ b/dart/benchmark/onera.py
@@ -36,11 +36,7 @@ def newton(pbl):
try:
newton = dart.Newton(pbl, tbox.Pardiso(), nthrds=k, vrb=2)
except:
- gmres = tbox.Gmres()
- gmres.setFillFactor(2)
- gmres.setDropTol(1e-6)
- gmres.setRestart(50)
- gmres.setTolerance(1e-5)
+ gmres = tbox.Gmres(2, 1e-6, 50, 1e-5)
newton = dart.Newton(pbl, gmres, nthrds=k, vrb=2)
return newton
diff --git a/dart/default.py b/dart/default.py
index 93cb86c..caf91dc 100644
--- a/dart/default.py
+++ b/dart/default.py
@@ -150,7 +150,7 @@ def picard(pbl):
problem formulation
"""
args = parseargs()
- solver = dart.Picard(pbl, tbox.Gmres(), nthrds=args.k, vrb=args.verb+1)
+ solver = dart.Picard(pbl, tbox.Gmres(1, 1e-6, 30, 1e-8), nthrds=args.k, vrb=args.verb+1)
return solver
def newton(pbl):
@@ -187,7 +187,7 @@ def morpher(msh, dim, mov, fxd = ['upstream', 'farfield', 'downstream'], fld = '
name of physical tag contaning the symmetry plane (y)
"""
args = parseargs()
- mshDef = tbox.MshDeform(msh, dim, nthrds=args.k)
+ mshDef = tbox.MshDeform(msh, tbox.Gmres(1, 1e-6, 30, 1e-8), dim, nthrds=args.k)
mshDef.setField(fld)
mshDef.addFixed(fxd)
mshDef.addMoving(mov)
@@ -207,8 +207,7 @@ def adjoint(solver, morpher):
solver : tbox.MshDeform object
mesh morpher
"""
- from tbox.solvers import LinearSolver
- adjoint = dart.Adjoint(solver, morpher, LinearSolver().pardiso())
+ adjoint = dart.Adjoint(solver, morpher)
return adjoint
def initViewer(problem):
diff --git a/dart/src/wAdjoint.cpp b/dart/src/wAdjoint.cpp
index 3e58eae..8b22008 100644
--- a/dart/src/wAdjoint.cpp
+++ b/dart/src/wAdjoint.cpp
@@ -32,6 +32,7 @@
#include "wNode.h"
#include "wElement.h"
#include "wLinearSolver.h"
+#include "wGmres.h"
#include "wMshDeform.h"
#include "wResults.h"
#include "wMshExport.h"
@@ -51,12 +52,24 @@ using namespace dart;
/**
* @brief Initialize the adjoint solver
*/
-Adjoint::Adjoint(std::shared_ptr<Newton> _sol, std::shared_ptr<tbox::MshDeform> _morph, std::shared_ptr<tbox::LinearSolver> _linsol) : sol(_sol), morph(_morph), linsol(_linsol)
+Adjoint::Adjoint(std::shared_ptr<Newton> _sol, std::shared_ptr<tbox::MshDeform> _morph) : sol(_sol), morph(_morph)
{
// Default parameters
nthreads = sol->nthreads;
verbose = sol->verbose;
+ // Linear solvers
+ if (sol->linsol->type() == LSOLTYPE::GMRES_ILUT)
+ {
+ // Use the same GMRES, but with a thighter tolerance
+ std::shared_ptr<Gmres> gmres = std::make_shared<Gmres>(*dynamic_cast<Gmres*>(sol->linsol.get()));
+ gmres->setTolerance(1e-8);
+ alinsol = gmres;
+ }
+ else
+ alinsol = sol->linsol;
+ mlinsol = morph->linsol;
+
// Update element gradients
sol->pbl->initGradElems();
@@ -67,6 +80,13 @@ Adjoint::Adjoint(std::shared_ptr<Newton> _sol, std::shared_ptr<tbox::MshDeform>
dCdMsh.resize(sol->pbl->msh->nodes.size(), Eigen::Vector3d::Zero());
dClAoa = 0;
dCdAoa = 0;
+
+ // Display solver parameters
+ std::cout << "--- Adjoint solver ---\n"
+ << "Aero inner solver: " << *alinsol
+ << "Mesh inner solver: " << *mlinsol
+ << "Number of threads: " << nthreads << "\n"
+ << std::endl;
}
/**
@@ -81,12 +101,6 @@ void Adjoint::run()
// Init
tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
- // Display solver parameters
- std::cout << "--- Adjoint solver ---\n"
- << "Inner solver: " << *linsol
- << "Number of threads: " << nthreads << "\n"
- << std::endl;
-
// Compute partial gradients of flow residuals and solve flow adjoint
tms["0-AdjFlo"].start();
this->linearizeFlow(); // dRu/dU, dRu/dA, dRu/dX
@@ -195,7 +209,7 @@ std::vector<double> Adjoint::solveMesh(std::vector<double> const &dX)
Eigen::Map<Eigen::VectorXd> dR_(dR.data(), dR.size());
Eigen::VectorXd dX_ = Eigen::Map<const Eigen::VectorXd>(dX.data(), dX.size());
// Solve
- linsol->compute(dRx_X.transpose(), Eigen::Map<Eigen::VectorXd>(dX_.data(), dX_.size()), dR_);
+ mlinsol->compute(dRx_X.transpose(), Eigen::Map<Eigen::VectorXd>(dX_.data(), dX_.size()), dR_);
return std::vector<double>(dR.data(), dR.data() + dR.size());
}
@@ -243,7 +257,7 @@ std::vector<double> Adjoint::solveFlow(std::vector<double> const &dU)
Eigen::Map<Eigen::VectorXd> dR_(dR.data(), dR.size());
Eigen::VectorXd dU_ = Eigen::Map<const Eigen::VectorXd>(dU.data(), dU.size());
// Solve
- linsol->compute(dRu_U.transpose(), Eigen::Map<Eigen::VectorXd>(dU_.data(), dU_.size()), dR_);
+ alinsol->compute(dRu_U.transpose(), Eigen::Map<Eigen::VectorXd>(dU_.data(), dU_.size()), dR_);
return std::vector<double>(dR.data(), dR.data() + dR.size());
}
diff --git a/dart/src/wAdjoint.h b/dart/src/wAdjoint.h
index 308a438..a33efb1 100644
--- a/dart/src/wAdjoint.h
+++ b/dart/src/wAdjoint.h
@@ -37,9 +37,10 @@ namespace dart
class DART_API Adjoint : public fwk::wSharedObject
{
public:
- std::shared_ptr<Newton> sol; ///< Newton solver
- std::shared_ptr<tbox::MshDeform> morph; ///< mesh morpher
- std::shared_ptr<tbox::LinearSolver> linsol; ///< linear solver
+ std::shared_ptr<Newton> sol; ///< Newton solver
+ std::shared_ptr<tbox::MshDeform> morph; ///< mesh morpher
+ std::shared_ptr<tbox::LinearSolver> alinsol; ///< linear solver (adjoint aero)
+ std::shared_ptr<tbox::LinearSolver> mlinsol; ///< linear solver (adjoint mesh)
int nthreads; ///< number of threads for the assembly
int verbose; ///< display more info
@@ -61,7 +62,7 @@ private:
fwk::Timers tms; ///< internal timers
public:
- Adjoint(std::shared_ptr<Newton> _sol, std::shared_ptr<tbox::MshDeform> _morph, std::shared_ptr<tbox::LinearSolver> _linsol);
+ Adjoint(std::shared_ptr<Newton> _sol, std::shared_ptr<tbox::MshDeform> _morph);
virtual ~Adjoint() { std::cout << "~Adjoint()\n"; }
void run();
diff --git a/dart/src/wSolver.cpp b/dart/src/wSolver.cpp
index 8cf95c7..68d2b61 100644
--- a/dart/src/wSolver.cpp
+++ b/dart/src/wSolver.cpp
@@ -64,9 +64,9 @@ Solver::Solver(std::shared_ptr<Problem> _pbl,
<< "** _ /_/ /_ ___ | __ _/_ / **\n"
<< "** /_____/ /_/ |_/_/ |_| /_/ **\n"
<< "*******************************************************************************\n"
- << "** Hi! My name is DART v1.1.0-21.09 **\n"
+ << "** Hi! My name is DART v1.2.0-22.04 **\n"
<< "** Adrien Crovato **\n"
- << "** ULiege 2018-2021 **\n"
+ << "** ULiege 2018-2022 **\n"
<< "*******************************************************************************\n"
<< std::endl;
diff --git a/ext/amfe b/ext/amfe
index 054ba00..46eb2e5 160000
--- a/ext/amfe
+++ b/ext/amfe
@@ -1 +1 @@
-Subproject commit 054ba00c327e4f2745ff9c559832247f90c3cec9
+Subproject commit 46eb2e5f08e56260e21fe610bf572a708c1b1179
--
GitLab