diff --git a/.github/workflows/format-and-lint.yaml b/.github/workflows/format-and-lint.yaml
new file mode 100644
index 000000000..5607a79c4
--- /dev/null
+++ b/.github/workflows/format-and-lint.yaml
@@ -0,0 +1,12 @@
+name: Format and Lint
+
+on:
+  push:
+    branches: [main]
+    tags:
+      - v*.*.*
+  pull_request:
+
+jobs:
+  pre-commit:
+    uses: mdolab/.github/.github/workflows/format-and-lint.yaml@ruffConfig
diff --git a/adflow/mphys/mphys_adflow.py b/adflow/mphys/mphys_adflow.py
index df61c2141..6c41bdb49 100644
--- a/adflow/mphys/mphys_adflow.py
+++ b/adflow/mphys/mphys_adflow.py
@@ -13,7 +13,7 @@
 
 # Set this to true to print out the name of the function being called and the class it's being called from along with
 # printing messages when node coordinates and states are updated from OpenMDAO inputs and outputs.
-DEBUG_LOGGING = False
+DEBUG_LOGGING = True
 
 
 X_AERO0 = MPhysVariables.Aerodynamics.Surface.COORDINATES_INITIAL
@@ -23,7 +23,6 @@
 F_AERO = MPhysVariables.Aerodynamics.Surface.LOADS
 Q_AERO = MPhysVariables.Aerodynamics.Surface.HEAT_FLOW
 
-
 def print_func_call(component):
     """Prints the name of the class and function being. Useful for debugging when you want to see what order OpenMDAO
     is calling things in.
@@ -499,6 +498,8 @@ def _setup_vectors(self, root_vectors):
         self.set_val("adflow_states", self.solver.getStates())
 
     def apply_nonlinear(self, inputs, outputs, residuals):
+        if DEBUG_LOGGING:
+            print_func_call(self)
         solver = self.solver
         ap = self.ap
         setAeroProblem(solver, ap, self.ap_vars, inputs=inputs, outputs=outputs, print_dict=False)
@@ -506,6 +507,15 @@ def apply_nonlinear(self, inputs, outputs, residuals):
         # flow residuals
         residuals["adflow_states"] = solver.getResidual(ap)
 
+    def printStatesRange(self, states):
+        stateNames = ["density", "velocityX", "velocityY", "velocityZ", "Energy", "Turb"]
+        for ii, stateName in enumerate(stateNames):
+            singleState = states[ii::6]
+            minState = self.comm.allreduce(np.min(singleState), op=MPI.MIN)
+            maxState = self.comm.allreduce(np.max(singleState), op=MPI.MAX)
+            if self.comm.rank == 0:
+                print(f"{stateName}: min = {minState:.11e}, max = {maxState:.11e}")
+
     def solve_nonlinear(self, inputs, outputs):
         if DEBUG_LOGGING:
             print_func_call(self)
@@ -513,10 +523,42 @@ def solve_nonlinear(self, inputs, outputs):
         ap = self.ap
 
         if self._do_solve:
+
             setAeroProblem(solver, ap, self.ap_vars, inputs=inputs, outputs=outputs, print_dict=False)
             ap.solveFailed = False  # might need to clear this out?
             ap.fatalFail = False
 
+            # For some reason, the ADflow solver will randomly NaN in the first iteration of a solve when starting from a state and mesh that it has already converged successfully. I have tried many things to debug this but nothing has worked:
+            # - Explicitly setting state before solve
+            # - compiling with `-O1 -xCORE-AVX2`
+            # - Compiling with `-fp-model=precise`
+            # - Compiling with only `-O2`
+            # - Disabling blockettes
+            # - running with a different number of procs
+            # - Running on `cas_ait` instead of `sky_ele` nodes
+            # - As above but setting slightly perturbed states
+            # - Running on a different mesh
+            # - Running on stampede
+
+            # The thing I have found to work is:
+            # - Set the OpenMDAO inputs and outputs
+            # - Evaluate the residual
+            # - If there is a NaN in the residual, call `resetFlow` and set the same OpenMDAO inputs and outputs again
+            # - Run the solver
+            res = solver.getResidual(ap)
+
+            nanInResidual = self.comm.allreduce(any(np.isnan(res)), op=MPI.LOR)
+            if self.comm.rank == 0:
+                if nanInResidual:
+                    print("NaN present in residual")
+                else:
+                    print("No NaNs in residual")
+
+            if nanInResidual:
+                # Reset the solution then set the inputs and states again, somehow this fixes the issue
+                solver.resetFlow(ap)
+                setAeroProblem(solver, ap, self.ap_vars, inputs=inputs, outputs=outputs, print_dict=False)
+
             # do not write solution files inside the solver loop
             solver(ap, writeSolution=False)
 
@@ -529,6 +571,10 @@ def solve_nonlinear(self, inputs, outputs):
                     print("# Solve Fatal Fail. Analysis Error")
                     print("###############################################################")
 
+                # write the solution so that we can diagnose
+                solver.writeSolution(baseName=fail_name, number=self.solution_counter)
+                self.solution_counter += 1
+
                 raise AnalysisError("ADFLOW Solver Fatal Fail")
 
             if ap.solveFailed:
diff --git a/ruff.toml b/ruff.toml
new file mode 100644
index 000000000..c20358d82
--- /dev/null
+++ b/ruff.toml
@@ -0,0 +1,8 @@
+extend-exclude =[
+    "src/build/fort_depend.py",
+    "src/adjoint/autoEdit/autoEdit*.py",
+    "src_cs/build/complexify.py",
+    "adflow/pyWeightAndBalance.py",
+    "adflow/pyWingCG.py",
+    "adflow/om_*.py",
+]