I have issues with nonstationarytransportproblem and nltransienttransportproblem, and I created replacement, just called "TransientTransport".
It was easier than trying to decipher what the other solvers were up to, so that the results could be compared.
Recently, I have fixed up various bugs, trying out the solver on every test.

First, there are tests that use both of them, tmquad.in and tmquadnl.in
These are completely linear models, they use the same time stepping, same initial conditions.. but still yield different results. The nonlinear solver also converges in 1 iteration, but still doesn't reduce to the linear response.
The root cause of this is the inconsistency of boundary conditions, initial conditions time stepping, and time discretization.
This caused me quite a lot of headache, and is the primary reason i opted on created a third solver seperately instead.

My main gripe with the transport solvers: awkward nonlinear setup

When setting up the FE-problem for transport problems, we should always have the interna forces computed as

R(u) = integral B^T * flux(grad(u), u) dV

(i call the unknown u here, i.e. temperature or concentration)
where the flux is defined as a, possible nonlinear, function in terms of the gradient and the field itself.
This is just like we do for the stress.

I see that the nonstationarytransportproblem and nltransienttransportproblem both handle this in a completely different way, (I've seen similar setups as these in transport problems before, so I suspect it's there for historical reasons).
It's just backwards for no tangible benefit whatsoever.
They start with linear problem, and tries to back on nonlinearity, which I find very backwards.
It forces one to add a "bc" RHS vector.
It can't handle models that are nonlinear in the gradient.
It requires around twice as much code compared to just the general nonlinear approach.

There are no benefits gained from this approach.

Inconsistent time discretization

The time discretization scheme should always decide VM_Total and VM_Velocity. This is the only way we have consistency of turning off/on b.c.s etc.
So, b.c.s and i.c.s should apply on a specific target time. If the VM_Velocity is requested, it should not go to the b.c and ask for the analytical dF/dt at the intrinsic time. That is inconsistent, and leads to obscure analysis bugs.
Take for example the heaviside function, which is frequently used. VM_Velocity is always zero, but the value changes.

To give a clear example, in pseudocode

Setup:
Initial condition:  u(t=0) = u0
Boundary condition: u(t>0) = u1     (using heaviside func.)
Using generalized midpoint, alpha, to define intrinsic time.

Currently, we get this:
u = u1
du/dt = 0

We should do this:
u = u0*(1-alpha) + uBC*alpha 
du/dt = (uBC-u0)/deltaT

I.e. we should do the same thing for all dofs, unknown or prescribed, initial condition or not.

Exactly how BC and IC should be applied is up for debate. Maybe the BC should override the IC, but either way, this is how it should work.

TransientTransport
This is the replacement i made for these. The code is much simpler, it's solves all the problems, however, it won't give the same results as the other solvers, because I believe those to be inconsistent and, well, wrong. So I choose not to do the same.
With initial time set in the old solvers, and with properly defined IC and BC, they all produce the same results.

The discrepencies come from
1. Slightly different time stepping setup: Starting from 0->deltaT, instead of going -deltaT->o for the first step. I think this is more sensible.

2. The other solvers are inconsistent in time discretization: the time vector set at the "state" vector of the models  (when calling updateInternalState)  are set to the values at the target time not the intrinsic time, i.e. it doesn't respect "alpha".
I consider this a bug, but it's one that can't be fixed easily. Which I why i made a new version.

3. Again, more inconsistencies with time discretization; The flux reported in the output from one time step is incorrect in nonstationary transport. It also doesn't respect the "alpha" parameter.
This is due to nonstationarytransport trying to handle alpha purely on the global FE-level, which can only ever be done if you assume linear material models. That won't work correctly for anything nonlinear, and gives misleading statevector and output.

4. It's slightly easier to use. Doesn't need "changingProblemSize" parameter. Uses the standard NRSolver input parameters. We could do all the mstep-steps as well.