diff --git a/dart/api/mphys.py b/dart/api/mphys.py
index 9a121c55c3bb356b2daaa132036250e5f970d1cf..d23c81e436631b6d30747762b29c13479a8742a6 100644
--- a/dart/api/mphys.py
+++ b/dart/api/mphys.py
@@ -218,7 +218,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')
+        self.options.declare('rerr', desc='raise AnalysisError when solver not fully converged')
 
     def setup(self):
         self.sol = self.options['sol']
@@ -239,7 +239,7 @@ class DartSolver(om.ImplicitComponent):
         # Compute
         self.sol.reset() # resets ICs (slows down convergence, but increases robustness for aero-structural)
         status = self.sol.run()
-        if (status == 1 and self.options['raiseError']) or status == 2: # max number of iterations exceeded or NaN in solution
+        if (status == 1 and self.options['rerr']) 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
@@ -389,36 +389,43 @@ class DartWriter(om.ExplicitComponent):
 
     Attributes
     ----------
-    iscn : int
-        ID of scenario (default: 0)
     sol : dart.Newton object
         Direct Newton sovler
-    adj : dart.Adjoint object
-        Adjoint solver
+    wrt : dart.MshExport object
+        Mesh and solution data writer
+    sname : string
+        Name of scenario
+    sall : bool
+        Whether to save results at each iteration to different files
+    it : int
+        Optimization iteration count
     """
     def initialize(self):
         self.options.declare('sol', desc='direct solver', recordable=False)
-        self.options.declare('wrtr', desc='data writer', recordable=False)
-        self.options.declare('iscn', default=0, desc='ID of the scenario')
-        self.options.declare('morph', default=False, desc='whether the mesh can deform or not')
+        self.options.declare('wrt', desc='data writer', recordable=False)
+        self.options.declare('snam', desc='name of the scenario')
+        self.options.declare('sall', desc='whether to save results at each iteration to different files')
 
     def setup(self):
-        self.iscn = self.options['iscn']
-        self.morph = self.options['morph']
         self.sol = self.options['sol']
-        self.wrtr = self.options['wrtr']
+        self.wrt = self.options['wrt']
+        self.snam = self.options['snam'] if self.options['snam'] is not None else 'default'
+        self.sall = self.options['sall']
+        self.it = 0 # optimization iteration count
 
     def compute(self, inputs, outputs):
         """Write to disk
         """
+        # Create suffix
+        suffix = '_{:s}'.format(self.snam)
+        if self.sall:
+            suffix += '_{:04d}'.format(self.it)
+        self.it += 1
         # Write mesh to disk (.msh only)
-        if self.morph and self.wrtr.type() == 1:
-            if self.iscn == 0:
-                self.wrtr.save(self.sol.pbl.msh.name)
-            else:
-                self.wrtr.save(self.sol.pbl.msh.name + '_{0:d}'.format(self.iscn))
+        if self.wrt.type() == 1:
+            self.wrt.save(self.sol.pbl.msh.name + suffix)
         # Write solution to disk
-        self.sol.save(self.wrtr, self.iscn)
+        self.sol.save(self.wrt, suffix)
 
 # Aerodynamic coupling group
 class DartGroup(om.Group):
@@ -437,17 +444,17 @@ 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')
+        self.options.declare('rerr', desc='raise AnalysisError when solver not fully converged')
 
     def setup(self):
         # Components
         if self.options['mrf'] is None:
-            self.add_subsystem('morpher', DartDummyMorpher(dim = self.options['sol'].pbl.nDim, msh = self.options['sol'].pbl.msh), promotes_inputs=['x_aero'], promotes_outputs=['xv'])
+            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'], raiseError = self.options['raiseError']), promotes_inputs=['aoa', 'xv'], promotes_outputs=['phi'])
+            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'], rerr=self.options['rerr']), promotes_inputs=['aoa', 'xv'], promotes_outputs=['phi'])
         if self.options['qinf'] is not None:
-            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'])
+            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'])
 
 # Aerodynamic post-coupling group
 class DartPostGroup(om.Group):
@@ -462,14 +469,14 @@ class DartPostGroup(om.Group):
     def initialize(self):
         self.options.declare('sol', desc='direct solver', recordable=False)
         self.options.declare('adj', desc='adjoint solver', recordable=False)
-        self.options.declare('wrtr', desc='data writer', recordable=False)
-        self.options.declare('iscn', default=0, desc='ID of the scenario')
-        self.options.declare('morph', default=False, desc='whether the mesh can deform or not')
+        self.options.declare('wrt', desc='data writer', recordable=False)
+        self.options.declare('snam', desc='name of the scenario')
+        self.options.declare('sall', desc='whether to save results at each iteration to different files')
 
     def setup(self):
         # Components
-        self.add_subsystem('coeff', DartCoefficients(sol = self.options['sol'], adj = self.options['adj']), promotes_inputs=['aoa', 'xv', 'phi'], promotes_outputs=['cl', 'cd'])
-        self.add_subsystem('writer', DartWriter(sol = self.options['sol'], wrtr = self.options['wrtr'], iscn = self.options['iscn'], morph = self.options['morph']))
+        self.add_subsystem('coeff', DartCoefficients(sol=self.options['sol'], adj=self.options['adj']), promotes_inputs=['aoa', 'xv', 'phi'], promotes_outputs=['cl', 'cd'])
+        self.add_subsystem('writer', DartWriter(sol=self.options['sol'], wrt=self.options['wrt'], snam=self.options['snam'], sall=self.options['sall']))
 
 # Builder
 class DartBuilder(Builder):
@@ -481,8 +488,8 @@ class DartBuilder(Builder):
         Dimension of the problem
     __qinf : float
         Freestream dynamic pressure
-    __wrtr : tbox.MshExport object
-        Mesh writer
+    __wrt : tbox.MshExport object
+        Mesh and solution data writer
     __mrf : tbox.MshDeform object
         Mesh morpher
     __sol : dart.Newton object
@@ -490,7 +497,7 @@ class DartBuilder(Builder):
     __adj : dart.Adjoint object
         Adjoint solver
     """
-    def __init__(self, cfg, scenario='aerodynamic', task='analysis', scenarioId=0, raiseError=False):
+    def __init__(self, cfg, scenario='aerodynamic', task='analysis', saveAll=False, raiseError=False):
         """Instantiate the builder and save the parameters and options
 
         Parameters
@@ -501,8 +508,8 @@ class DartBuilder(Builder):
             Type of scenario (available: aerodynamic, aerostructural) (default: aerodynamic)
         task : string, optional
             Type of task (available: analysis, optimization) (default: analysis)
-        scenarioId : int, optional
-            ID of the scenario, will be used in the output filename
+        saveAll : bool, optional
+            Whether to write results to different files at each iteration (default: False, results file will be overwritten)
         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
@@ -510,8 +517,8 @@ class DartBuilder(Builder):
         self.__cfg = cfg
         self.__scn = scenario
         self.__tsk = task
-        self.__sid = scenarioId
-        self.__raiseError = raiseError
+        self.__sall = saveAll
+        self.__rerr = raiseError
 
     def initialize(self, comm):
         """Instantiate and initialize all the solver components
@@ -523,7 +530,7 @@ class DartBuilder(Builder):
         """
         from .core import initDart
         def _init():
-            self.__dim, self.__qinf, _, self.__wrtr, self.__mrf, _, self.__bnd, self.__sol, self.__adj = initDart(self.__cfg, scenario=self.__scn, task=self.__tsk)
+            self.__dim, self.__qinf, _, self.__wrt, self.__mrf, _, self.__bnd, self.__sol, self.__adj = initDart(self.__cfg, scenario=self.__scn, task=self.__tsk)
         # If n MPI processes, init DART on rank=0,...,n-1 sequentially to avoid issues with mesh IO
         # If mpi4py is not available or only 1 MPI process, init DART as usual
         try:
@@ -545,20 +552,20 @@ class DartBuilder(Builder):
         except:
             _init()
 
-    def get_mesh_coordinate_subsystem(self):
+    def get_mesh_coordinate_subsystem(self, scenario_name=None):
         """Return openMDAO component to get the initial surface mesh coordinates
         """
-        return DartMesh(dim = self.__dim, bnd = self.__bnd)
+        return DartMesh(dim=self.__dim, bnd=self.__bnd)
 
-    def get_coupling_group_subsystem(self):
+    def get_coupling_group_subsystem(self, scenario_name=None):
         """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, raiseError = self.__raiseError)
+        return DartGroup(qinf=self.__qinf, bnd=self.__bnd, sol=self.__sol, mrf=self.__mrf, adj=self.__adj, rerr=self.__rerr)
 
-    def get_post_coupling_subsystem(self):
+    def get_post_coupling_subsystem(self, scenario_name=None):
         """Return openMDAO group that computes the aero coefficients and writes data to disk
         """
-        return DartPostGroup(sol = self.__sol, adj = self.__adj, wrtr = self.__wrtr, iscn = self.__sid, morph = False if self.__mrf is None else True)
+        return DartPostGroup(sol=self.__sol, adj=self.__adj, wrt=self.__wrt, snam=scenario_name, sall=self.__sall)
 
     def get_number_of_nodes(self):
         """Return the number of surface nodes
diff --git a/dart/src/wAdjoint.cpp b/dart/src/wAdjoint.cpp
index e48ce29dd3b46ea7bb9b2ea9403c5f0ac405a04f..ac6acba9f47ea3e1c5a88ac8fc1607b52480dd94 100644
--- a/dart/src/wAdjoint.cpp
+++ b/dart/src/wAdjoint.cpp
@@ -724,7 +724,7 @@ void Adjoint::buildGradientCoefficientsMesh(Eigen::RowVectorXd &dCl, Eigen::RowV
 /**
  * @brief Write the results
  */
-void Adjoint::save(MshExport *mshWriter, int n)
+void Adjoint::save(MshExport *mshWriter, std::string const &suffix)
 {
     // Write files
     std::cout << "Saving files "
@@ -740,18 +740,9 @@ void Adjoint::save(MshExport *mshWriter, int n)
     results.vectors_at_nodes["gradClMsh"] = &dClMsh;
     results.vectors_at_nodes["gradCdMsh"] = &dCdMsh;
     // save (whole mesh and bodies)
-    if (n > 0)
-    {
-        mshWriter->save(sol->pbl->msh->name + "_adjoint_" + std::to_string(n), results);
-        for (auto bnd : sol->pbl->bodies)
-            bnd->save(bnd->groups[0]->tag->name + "adjoint_" + std::to_string(n), results);
-    }
-    else
-    {
-        mshWriter->save(sol->pbl->msh->name + "_adjoint", results);
-        for (auto bnd : sol->pbl->bodies)
-            bnd->save(bnd->groups[0]->tag->name + "_adjoint", results);
-    }
+    mshWriter->save(sol->pbl->msh->name + "_adjoint" + suffix, results);
+    for (auto bnd : sol->pbl->bodies)
+        bnd->save(bnd->groups[0]->tag->name + "_adjoint" + suffix, results);
 }
 
 void Adjoint::write(std::ostream &out) const
diff --git a/dart/src/wAdjoint.h b/dart/src/wAdjoint.h
index a33efb11100d5775f064e882e91b57d56ea63de0..b2dc8ffd116cf3415e30bde0fef648e4fac4fb54 100644
--- a/dart/src/wAdjoint.h
+++ b/dart/src/wAdjoint.h
@@ -66,7 +66,7 @@ public:
     virtual ~Adjoint() { std::cout << "~Adjoint()\n"; }
 
     void run();
-    void save(tbox::MshExport *mshWriter, int n = 0);
+    void save(tbox::MshExport *mshWriter, std::string const &suffix = "");
 
     // Partial gradients of the mesh residuals
     void linearizeMesh();
diff --git a/dart/src/wSolver.cpp b/dart/src/wSolver.cpp
index e16f839db1e4b9e7eae155be35a6fc7235fcf2d5..5ad807f17353ee1ec2c850453468d441788caece 100644
--- a/dart/src/wSolver.cpp
+++ b/dart/src/wSolver.cpp
@@ -161,7 +161,7 @@ STATUS Solver::run()
 /**
  * @brief Write the results
  */
-void Solver::save(MshExport *mshWriter, int n)
+void Solver::save(MshExport *mshWriter, std::string const &suffix)
 {
     // Write files
     std::cout << "Saving files "
@@ -180,18 +180,9 @@ void Solver::save(MshExport *mshWriter, int n)
     results.scalars_at_nodes["Mach"] = &M;
     results.scalars_at_nodes["Cp"] = &Cp;
     // save (whole mesh and bodies)
-    if (n > 0)
-    {
-        mshWriter->save(pbl->msh->name + "_" + std::to_string(n), results);
-        for (auto bnd : pbl->bodies)
-            bnd->save(bnd->groups[0]->tag->name + "_" + std::to_string(n), results);
-    }
-    else
-    {
-        mshWriter->save(pbl->msh->name, results);
-        for (auto bnd : pbl->bodies)
-            bnd->save(bnd->groups[0]->tag->name, results);
-    }
+    mshWriter->save(pbl->msh->name + suffix, results);
+    for (auto bnd : pbl->bodies)
+        bnd->save(bnd->groups[0]->tag->name + suffix, results);
 }
 
 /**
diff --git a/dart/src/wSolver.h b/dart/src/wSolver.h
index aab40ad17beb65aaee05537c73ab5bad693dc710..fdf41515f3278f2056eb7bf977f21dc0ecb18e10 100644
--- a/dart/src/wSolver.h
+++ b/dart/src/wSolver.h
@@ -78,7 +78,7 @@ public:
 
     void reset();
     virtual STATUS run();
-    void save(tbox::MshExport *mshWriter, int n = 0);
+    void save(tbox::MshExport *mshWriter, std::string const &suffix = "");
 
 protected:
     void computeFlow();
diff --git a/dart/tests/rae_25.py b/dart/tests/rae_25.py
index cc656c11f7211fad86603d24ea9b253615b57698..1677f290debd46cfa69b2ef138211e4d8a52092d 100644
--- a/dart/tests/rae_25.py
+++ b/dart/tests/rae_25.py
@@ -122,11 +122,11 @@ def main():
     tests.add(CTest('CS', solver.Cs, 0.0000, 1e-3, forceabs=True))
     tests.add(CTest('CM', solver.Cm, -0.113, 5e-2)) # 2D -0.117
     tests.add(CTest('dCl/dAoA', adjoint.dClAoa, 12.7, 5e-2)) # 2D 12.7, FD 12.73 (1e-6)
-    tests.add(CTest('dCd/dAoA', adjoint.dCdAoa, 0.059, 5e-3, forceabs=True)) # 2D 0.051, FD 0.0557 (1e-6)
+    tests.add(CTest('dCd/dAoA', adjoint.dCdAoa, 0.055, 1e-2, forceabs=True)) # 2D 0.051, FD 0.0557 (1e-6)
     tests.add(CTest('dCl_dMsh', dClX, 22.0, 5e-2))
-    tests.add(CTest('dCd_dMsh', dCdX, 1.0, 5e-2))
+    tests.add(CTest('dCd_dMsh', dCdX, 1.0, 1e-1))
     tests.add(CTest('dCl/dAoA (msh)', dClAoA, 12.9, 5e-2)) # 2D 12.9
-    tests.add(CTest('dCd/dAoA (msh)', dCdAoA, 0.058, 5e-3, forceabs=True)) # 2D 0.051
+    tests.add(CTest('dCd/dAoA (msh)', dCdAoA, 0.055, 1e-2, forceabs=True)) # 2D 0.051
     tests.run()
 
     # eof