Finite difference method for coupled PDEs: optimizing performance (time step, iterations per step)

I'm solving coupled PDEs using finite difference method: Incompressible Navier-Stokes and the divergence-free induction equation (Maxwell's equations) with non-uniform electrical conductivity. The dimensionless equations are:

$1) \partial_t \mathbf{u} + \nabla \bullet (\mathbf{u}\mathbf{u}) = -\nabla p + Re^{-1}\nabla^2 \mathbf{u} + N \mathbf{j} \times \mathbf{B}, \qquad \mathbf{j} = Re_m^{-1} \nabla \times \mathbf{B}$

$2) \partial_t \mathbf{B} + Re_m^{-1} \nabla \times (\bar{\sigma}^{-1} \nabla \times \mathbf{B}) = \nabla \times (\mathbf{u} \times \mathbf{B})$

$3) \nabla \bullet \mathbf{u} = 0, \qquad \nabla \bullet \mathbf{B} = 0$

Where $Re,N,Re_m$ are the Reynolds, Interaction and magnetic Reynolds numbers. Also, $\bar{\sigma}$ is the dimensionless electrical conductivity.

The induction equation becomes stiff / ill-conditioned as $\bar{\sigma}$ approaches zero. In my case, $\bar{\sigma}$ is very small and, as a result, I resort to treating $Re_m^{-1} \nabla \times (\bar{\sigma}^{-1} \nabla \times \mathbf{B})$ implicitly in time using diagonally Preconditioned Conjugate Gradient Method (PCG). I use the limiting stopping criteria among 3 choices:

• $N_{\text{max-iterations}}$
• $\frac{||Ax_k-b||}{||Ax_0-b||} \le \varepsilon_{\text{relative}}$
• $||Ax_k-b|| \le \varepsilon_{\text{absolute}}$

Here, $A=I + \Delta t Re_m^{-1} \nabla \times (\bar{\sigma}^{-1} \nabla \times)$. $b = \Delta t (\mathbf{B}^n + \nabla \times (\mathbf{u}^n \times \mathbf{B}^n))$. $x_k$ is the $k$th estimated solution of $\mathbf{B}^{n+1}$ using PCG. Finally, $N,\varepsilon$ are the maximum allowable PCG iterations and tolerances (relative to initial guess and absolute respectively).

Experimenting with stopping criteria and different time steps ($\Delta t$), it seems that there is a non-trivial trade-off in picking $\Delta t$ and stopping criteria to optimize performance, which I'll define as

$\mathcal P = \frac{\text{convective time per step}}{\text{wall clock time per step}} = \frac{\text{convective time}}{\text{wall clock time}}$.

Ultimately, I'm trying to reach a steady solution, however, the coupled equations are non-linear, stiff and can (depending on $Re,N,Re_m$) have seemingly long time scales. If I choose $\Delta t$ too large, or perform too few iterations of PCG on Eq. 2), the solution diverges.

Questions:

1) Are there any methods to optimize the choice of stopping criteria or $\Delta t$ in order to optimize $\mathcal P$?

2) Are there any other tricks I can use to reach steady state faster?