From 4d61def3f975e63432210f9d9f0b678128bf2d86 Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Tue, 11 Jan 2022 16:24:05 +0100
Subject: [PATCH] Add proper exit code (status) for DART solver and handle in
mphys interface
---
dart/api/mphys.py | 16 ++++++++++++----
dart/src/wF0El.cpp | 4 ++--
dart/src/wNewton.cpp | 8 ++++----
dart/src/wNewton.h | 2 +-
dart/src/wPicard.cpp | 8 ++++----
dart/src/wPicard.h | 2 +-
dart/src/wSolver.cpp | 2 +-
dart/src/wSolver.h | 12 +++++++++++-
dart/tests/cylinder.py | 6 ++----
dart/tests/rae_25.py | 6 +++---
dart/tests/rae_3.py | 6 +++---
11 files changed, 44 insertions(+), 28 deletions(-)
diff --git a/dart/api/mphys.py b/dart/api/mphys.py
index 2edba1e..58200c0 100644
--- a/dart/api/mphys.py
+++ b/dart/api/mphys.py
@@ -215,6 +215,7 @@ class DartSolver(om.ImplicitComponent):
def initialize(self):
self.options.declare('sol', desc='direct solver', recordable=False)
self.options.declare('adj', desc='adjoint solver', recordable=False)
+ self.options.declare('raiseError', default=False, desc='raise AnalysisError when solver not fully converged')
def setup(self):
self.sol = self.options['sol']
@@ -234,7 +235,9 @@ class DartSolver(om.ImplicitComponent):
self.sol.pbl.update(inputs['aoa'][0])
# Compute
self.sol.reset() # resets ICs (slows down convergence, but increases robustness for aero-structural)
- self.sol.run()
+ status = self.sol.run()
+ if (status == 1 and self.options['raiseError']) or status == 2: # max number of iterations exceeded or NaN in solution
+ raise om.AnalysisError('DART solver failed to fully converge!\n')
# Update outputs
outputs['phi'] = self.sol.phi
@@ -412,6 +415,7 @@ class DartGroup(om.Group):
self.options.declare('sol', desc='direct solver', recordable=False)
self.options.declare('mrf', desc='mesh morpher', recordable=False)
self.options.declare('adj', desc='adjoint solver', recordable=False)
+ self.options.declare('raiseError', default=False, desc='raise AnalysisError when solver not fully converged')
def setup(self):
# Components
@@ -419,7 +423,7 @@ class DartGroup(om.Group):
self.add_subsystem('morpher', DartDummyMorpher(dim = self.options['sol'].pbl.nDim, msh = self.options['sol'].pbl.msh), promotes_inputs=['x_aero'], promotes_outputs=['xv'])
else:
self.add_subsystem('morpher', DartMorpher(dim = self.options['sol'].pbl.nDim, bnd = self.options['bnd'], mrf = self.options['mrf'], adj = self.options['adj']), promotes_inputs=['x_aero'], promotes_outputs=['xv'])
- self.add_subsystem('solver', DartSolver(sol = self.options['sol'], adj = self.options['adj']), promotes_inputs=['aoa', 'xv'], promotes_outputs=['phi'])
+ self.add_subsystem('solver', DartSolver(sol = self.options['sol'], adj = self.options['adj'], raiseError = self.options['raiseError']), promotes_inputs=['aoa', 'xv'], promotes_outputs=['phi'])
if self.options['qinf'] != 0:
self.add_subsystem('loads', DartLoads(qinf = self.options['qinf'], bnd = self.options['bnd'], adj = self.options['adj']), promotes_inputs=['xv', 'phi'], promotes_outputs=['f_aero'])
@@ -442,7 +446,7 @@ class DartBuilder(Builder):
__adj : dart.Adjoint object
Adjoint solver
"""
- def __init__(self, params, scenario='aerodynamic', task='analysis'):
+ def __init__(self, params, scenario='aerodynamic', task='analysis', raiseError=False):
"""Instantiate and initialize all the solver components.
Because we do not use MPI, we do not use the initialize method.
@@ -454,10 +458,14 @@ class DartBuilder(Builder):
Type of scenario (available: aerodynamic, aerostructural) (default: aerodynamic)
task : string, optional
Type of task (available: analysis, optimization) (default: analysis)
+ raiseError : bool, optional
+ Whether to raise an openMDAO.AnalysisError when the solver does not fully converge (default: False)
+ As of openMDAO V3.9.2, AnalysisError are only handled by pyOptSparse
"""
print('*'*31 + 'mphys.DartBuilder' + '*'*31 + '\n')
from .core import initDart
self.__dim, self.__qinf, _, self.__wrtr, self.__mrf, _, self.__bnd, self.__sol, self.__adj = initDart(params, scenario=scenario, task=task)
+ self.__raiseError = raiseError
print('*'*31 + 'mphys.DartBuilder' + '*'*31 + '\n')
def get_mesh_coordinate_subsystem(self):
@@ -468,7 +476,7 @@ class DartBuilder(Builder):
def get_coupling_group_subsystem(self):
"""Return openMDAO group containing the morpher, solver and loads
"""
- return DartGroup(qinf = self.__qinf, bnd = self.__bnd, sol = self.__sol, mrf = self.__mrf, adj = self.__adj)
+ return DartGroup(qinf = self.__qinf, bnd = self.__bnd, sol = self.__sol, mrf = self.__mrf, adj = self.__adj, raiseError = self.__raiseError)
def get_post_coupling_subsystem(self):
"""Return openMDAO component that computes the aero coefficients and writes data to disk
diff --git a/dart/src/wF0El.cpp b/dart/src/wF0El.cpp
index e17c9d1..bea2834 100644
--- a/dart/src/wF0El.cpp
+++ b/dart/src/wF0El.cpp
@@ -49,8 +49,8 @@ double F0ElC::evalGrad(Element const &e, std::vector<double> const &phi, size_t
F0ElRhoBase::F0ElRhoBase(double _mInf, double _mC2) : F0El(), gamma(1.4), mInf(_mInf)
{
uC = sqrt(_mC2 * (1 / (mInf * mInf) + (gamma - 1) / 2) / (1 + (gamma - 1) / 2 * _mC2)); // critical velocity
- rC = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - uC * uC); // r(uC)
- beta = ((gamma - 1) * mInf * mInf * uC) / rC; // dr/du(uC) / -r(uC)
+ rC = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - uC * uC); // r(uC)
+ beta = ((gamma - 1) * mInf * mInf * uC) / rC; // dr/du(uC) / -r(uC)
}
/**
* @brief Evaluate the base of the isentropic density
diff --git a/dart/src/wNewton.cpp b/dart/src/wNewton.cpp
index 1d80195..79bcd13 100644
--- a/dart/src/wNewton.cpp
+++ b/dart/src/wNewton.cpp
@@ -81,7 +81,7 @@ Newton::Newton(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver
*
* Solve the steady transonic Full Potential Equation
*/
-bool Newton::run()
+STATUS Newton::run()
{
// Multithread
tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
@@ -242,7 +242,7 @@ bool Newton::run()
std::cout << "Newton solver CPU" << std::endl
<< tms;
std::cout << std::endl;
- return true;
+ return STATUS::CONVERGED;
}
else if (std::isnan(relRes))
{
@@ -251,7 +251,7 @@ bool Newton::run()
std::cout << "Newton solver CPU" << std::endl
<< tms;
std::cout << std::endl;
- return false;
+ return STATUS::FAILED;
}
else
{
@@ -260,7 +260,7 @@ bool Newton::run()
std::cout << "Newton solver CPU" << std::endl
<< tms;
std::cout << std::endl;
- return true;
+ return STATUS::MAXIT;
}
}
diff --git a/dart/src/wNewton.h b/dart/src/wNewton.h
index 76cac53..bec77d2 100644
--- a/dart/src/wNewton.h
+++ b/dart/src/wNewton.h
@@ -49,7 +49,7 @@ private:
public:
Newton(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver> _linsol, double rTol = 1e-6, double aTol = 1e-8, int mIt = 25, int nthrds = 1, int vrb = 1);
- virtual bool run() override;
+ virtual STATUS run() override;
#ifndef SWIG
virtual void write(std::ostream &out) const override;
diff --git a/dart/src/wPicard.cpp b/dart/src/wPicard.cpp
index ea218f7..5d56751 100644
--- a/dart/src/wPicard.cpp
+++ b/dart/src/wPicard.cpp
@@ -75,7 +75,7 @@ Picard::Picard(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver
*
* Solve the steady subsonic Full Potential Equation
*/
-bool Picard::run()
+STATUS Picard::run()
{
// Multithread
tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
@@ -173,7 +173,7 @@ bool Picard::run()
std::cout << "Picard solver CPU" << std::endl
<< tms;
std::cout << std::endl;
- return true;
+ return STATUS::CONVERGED;
}
else if (std::isnan(relRes))
{
@@ -182,7 +182,7 @@ bool Picard::run()
std::cout << "Picard solver CPU" << std::endl
<< tms;
std::cout << std::endl;
- return false;
+ return STATUS::FAILED;
}
else
{
@@ -191,7 +191,7 @@ bool Picard::run()
std::cout << "Picard solver CPU" << std::endl
<< tms;
std::cout << std::endl;
- return true;
+ return STATUS::MAXIT;
}
}
diff --git a/dart/src/wPicard.h b/dart/src/wPicard.h
index 42a4580..0a8a27d 100644
--- a/dart/src/wPicard.h
+++ b/dart/src/wPicard.h
@@ -40,7 +40,7 @@ public:
public:
Picard(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver> _linsol, double rTol = 1e-6, double aTol = 1e-8, int mIt = 25, int nthrds = 1, int vrb = 1);
- virtual bool run() override;
+ virtual STATUS run() override;
#ifndef SWIG
virtual void write(std::ostream &out) const override;
diff --git a/dart/src/wSolver.cpp b/dart/src/wSolver.cpp
index 3b11893..8cf95c7 100644
--- a/dart/src/wSolver.cpp
+++ b/dart/src/wSolver.cpp
@@ -153,7 +153,7 @@ void Solver::reset()
/**
* @brief Dummy run
*/
-bool Solver::run()
+STATUS Solver::run()
{
throw std::runtime_error("dart::Solver::run not implemented!\n");
}
diff --git a/dart/src/wSolver.h b/dart/src/wSolver.h
index d48fc3a..aab40ad 100644
--- a/dart/src/wSolver.h
+++ b/dart/src/wSolver.h
@@ -29,6 +29,16 @@
namespace dart
{
+/**
+ * @brief Exit code of DART solver
+ */
+enum class STATUS
+{
+ CONVERGED = 0, // fully converged to prescribed (relative or absolute) tolerance
+ MAXIT = 1, // not fully converged (max. number of iterations exceeded)
+ FAILED = 2 // NaN in the solution vector
+};
+
/**
* @brief Base class for full potential solvers
* @authors Adrien Crovato
@@ -67,7 +77,7 @@ public:
virtual ~Solver();
void reset();
- virtual bool run();
+ virtual STATUS run();
void save(tbox::MshExport *mshWriter, int n = 0);
protected:
diff --git a/dart/tests/cylinder.py b/dart/tests/cylinder.py
index 2d5f078..07d4311 100644
--- a/dart/tests/cylinder.py
+++ b/dart/tests/cylinder.py
@@ -66,12 +66,12 @@ def main():
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['picard'].start()
solver0 = floD.picard(pbl)
- cnvrgd0 = solver0.run()
+ solver0.run()
solver0.save(gmshWriter)
tms['picard'].stop()
tms['newton'].start()
solver1 = floD.newton(pbl)
- cnvrgd1 = solver1.run()
+ solver1.run()
solver1.save(gmshWriter)
tms['newton'].stop()
@@ -126,8 +126,6 @@ def main():
# check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
- if not cnvrgd0 or not cnvrgd1:
- raise Exception(ccolors.ANSI_RED + 'Flow solver failed to converge!' + ccolors.ANSI_RESET)
tests = CTests()
tests.add(CTest('Picard: iteration count', solver0.nIt, 12, 1, forceabs=True))
tests.add(CTest('Picard: Neumann B.C. mean error', errNeu0, 0., 1e-2))
diff --git a/dart/tests/rae_25.py b/dart/tests/rae_25.py
index a121eb5..cc656c1 100644
--- a/dart/tests/rae_25.py
+++ b/dart/tests/rae_25.py
@@ -68,7 +68,7 @@ def main():
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver'].start()
solver = floD.newton(pbl)
- cnvrgd = solver.run()
+ status = solver.run()
solver.save(gmshWriter)
tms['solver'].stop()
@@ -113,8 +113,8 @@ def main():
# check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
- if (not cnvrgd):
- raise Exception(ccolors.ANSI_RED + 'Flow solver failed to converge!' + ccolors.ANSI_RESET)
+ if status > 0:
+ raise Exception(ccolors.ANSI_RED + 'Solver failed to converge!' + ccolors.ANSI_RESET)
tests = CTests()
tests.add(CTest('iteration count', solver.nIt, 6, 1, forceabs=True))
tests.add(CTest('CL', solver.Cl, 0.60, 5e-2)) # 2D, 0.62
diff --git a/dart/tests/rae_3.py b/dart/tests/rae_3.py
index 445e3a8..4d8e8c2 100644
--- a/dart/tests/rae_3.py
+++ b/dart/tests/rae_3.py
@@ -93,7 +93,7 @@ def main():
# solve problem
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver'].start()
- cnvrgd = solver.run()
+ status = solver.run()
solver.save(gmshWriter)
tms['solver'].stop()
@@ -138,8 +138,8 @@ def main():
# check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
- if (not cnvrgd):
- raise Exception(ccolors.ANSI_RED + 'Flow solver failed to converge!' + ccolors.ANSI_RESET)
+ if status > 0:
+ raise Exception(ccolors.ANSI_RED + 'Solver failed to converge!' + ccolors.ANSI_RESET)
tests = CTests()
tests.add(CTest('iteration count', solver.nIt, 7, 1, forceabs=True))
tests.add(CTest('CL', solver.Cl, 0.21, 5e-2))
--
GitLab