diff --git a/dart/api/mphys.py b/dart/api/mphys.py
index 2edba1ebe352dbcaa57bf879dec24690cdaec229..58200c0023f226513f7ecca1b575836a11cd5a01 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 e17c9d1f1d9428cf0f62835a11fc19396cb42556..bea283450ec26deed8a6dabad98b08af6feddddc 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 1d801953e86634863fba0c643f8c2c3c633c975f..79bcd1375a1487f7125872cfbfcd3c47e284405e 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 76cac5358ecb724d91b41afe7777c97382d2f832..bec77d2248501a96b9048b5a1829e34ab328c35b 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 ea218f7d4c02d3dd4e14c3816776fef51755279b..5d567514923b835699a796af894a5b5b32ace767 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 42a4580529f1822b6df2cd16a8719cb2c0d7cc97..0a8a27d10a6b539d07d6066027e33b1a203483a6 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 3b11893a558e00df78cb679bf60e2dc48fc7535e..8cf95c7ff5c61641401fa9d00c3c3dd4007126e0 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 d48fc3a7369c6d56293665e730ec0504abd6b4f7..aab40ad17beb65aaee05537c73ab5bad693dc710 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 2d5f07835f07e2607f88aa357f6a434ee4292b7b..07d43117a1d25ebf38f3f0e24cdbb843fc383762 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 a121eb5f4cbce8333bf005295012b3e72ea27dac..cc656c11f7211fad86603d24ea9b253615b57698 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 445e3a8f9c8b5eeb472f452d6560f73a98a5ffc4..4d8e8c291a79afe0b4ffe3b430e401b2ddb6967e 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))