diff --git a/asQ/paradiag.py b/asQ/paradiag.py
index 055e67ea..00092503 100644
--- a/asQ/paradiag.py
+++ b/asQ/paradiag.py
@@ -126,6 +126,9 @@ def reset_diagnostics(self):
         self.total_timesteps = 0
         self.total_windows = 0
         self.block_iterations.data()[:] = 0
+        jacobian = self.solver.jacobian
+        if hasattr(jacobian, "pc") and hasattr(jacobian.pc, "block_iterations"):
+            jacobian.pc.block_iterations.data(deepcopy=False)[:] = 0
 
     @profiler()
     def _record_diagnostics(self):
diff --git a/asQ/post.py b/asQ/post.py
index 453e883c..514bada5 100644
--- a/asQ/post.py
+++ b/asQ/post.py
@@ -185,6 +185,7 @@ def write_solver_parameters(sparams, directory=""):
 
     with open(file_name, "w") as f:
         f.write(dumps(sparams, indent=4))
+        f.write("\n")
 
     return
 
@@ -216,7 +217,7 @@ def write_paradiag_setup(pdg, directory=""):
                 f"global_comm.size = {pdg.ensemble.global_comm.size}\n"
             jacobian = pdg.solver.jacobian
             if hasattr(jacobian, "pc") and hasattr(jacobian.pc, "alpha"):
-                info += f"alpha = {pdg.solver.jacobian.pc.alpha}"
+                info += f"alpha = {pdg.solver.jacobian.pc.alpha}\n"
             f.write(info)
     return
 
@@ -254,7 +255,7 @@ def write_aaos_solve_metrics(pdg, directory=""):
                 f"nonlinear iterations per window = {nlits/nw}\n" + \
                 f"linear iterations per window = {lits/nw}\n" + \
                 f"linear iterations per nonlinear iteration = {lits/nlits}\n" + \
-                f"max iterations per block solve = {blits/lits}"
+                f"max iterations per block solve = {blits/lits}\n"
             f.write(info)
     return
 
diff --git a/utils/compressible_flow/compressible_flow.py b/utils/compressible_flow/compressible_flow.py
index aff10d25..59a85708 100644
--- a/utils/compressible_flow/compressible_flow.py
+++ b/utils/compressible_flow/compressible_flow.py
@@ -135,7 +135,7 @@ def hydrostatic_rho(Vv, V2, mesh, thetan, rhon, pi_boundary,
         Pif = pi_formula(rho, thetan, gas)
 
         rhoeqn = gas.cp*(
-            (fd.inner(v, dv) - fd.div(dv*thetan)*Pif)*fd.dx
+            (fd.inner(v, dv) - fd.div(dv*thetan)*Pif)*fd.dx(degree=6)
             + drho*fd.div(thetan*v)*fd.dx
         )
 
diff --git a/utils/shallow_water/miniapp.py b/utils/shallow_water/miniapp.py
index 7be019fe..6b259e2b 100644
--- a/utils/shallow_water/miniapp.py
+++ b/utils/shallow_water/miniapp.py
@@ -104,8 +104,8 @@ def form_mass(u, h, v, q):
         u0 = w0.subfunctions[self.velocity_index]
         h0 = w0.subfunctions[self.depth_index]
 
-        u0.project(velocity_expression(*x))
-        h0.project(depth_expression(*x))
+        u0.project(velocity_expression(*x), quadrature_degree=6)
+        h0.interpolate(depth_expression(*x))
 
         if reference_state:
             self.reference_state = fd.Function(self.function_space())
@@ -235,7 +235,8 @@ def sync_diagnostics(self):
 
         Until this method is called, diagnostic information is not guaranteed to be valid.
         """
-        self.cfl_series.synchronise()
+        if self.record_diagnostics['cfl']:
+            self.cfl_series.synchronise()
 
     def solve(self,
               nwindows=1,
@@ -261,7 +262,8 @@ def postprocess(pdg, w, rhs=None):
             postproc(self, pdg, w)
 
         # extend cfl array
-        self.cfl_series.resize(self.paradiag.total_windows + nwindows)
+        if self.record_diagnostics['cfl']:
+            self.cfl_series.resize(self.paradiag.total_windows + nwindows)
 
         self.paradiag.solve(nwindows,
                             preproc=preprocess,
diff --git a/utils/timing.py b/utils/timing.py
index 97f8bcc5..0843a29c 100644
--- a/utils/timing.py
+++ b/utils/timing.py
@@ -47,13 +47,16 @@ class SolverTimer(Timer):
     '''
     Time multiple solves and print out total/average etc times.
     '''
-    def string(self, timesteps_per_solve=1, ndigits=None):
+    def string(self, timesteps_per_solve=1,
+               total_iterations=1, ndigits=None):
         rnd = lambda x: x if ndigits is None else round(x, ndigits)
         total_time = self.total_time()
         average_time = self.average_time()
         timestep_time = average_time/timesteps_per_solve
+        iteration_time = total_time/total_iterations
         string = ''\
             + f'Total solution time: {rnd(total_time)}\n' \
             + f'Average solve solution time: {rnd(average_time)}\n' \
-            + f'Average timestep solution time: {rnd(timestep_time)}'
+            + f'Average timestep solution time: {rnd(timestep_time)}\n' \
+            + f'Average iteration solution time: {rnd(iteration_time)}'
         return string