Small inconsistency in nonlinear time step change

[David Dickinson]

In nonlinear_terms

    if (istep /= istep_last) then
       g3 = g2
       g2 = g1

       call debug_message(verb, 'nonlinear_terms::add_explicit copied old terms')

       ! if running nonlinearly, then compute the nonlinear term at grid points                                                                       
       ! and store it in g1                                                                                                                           
       if (present(nl)) then
          call debug_message(verb, 'nonlinear_terms::add_explicit calling add_nl')
          if(gryfx_zonal%on) then
            call add_nl_gryfx (g1)
          else
            call add_nl (g1, phi, apar, bpar)
          endif

          if(reset) return !Return if resetting                                                                                                       

          ! takes g1 at grid points and returns 2*g1 at cell centers                                                                                  

          call center (g1)

       else
          g1 = 0.
       end if

:
:
istep_last = istep

The first operation in this section is to "shuffle" the explicit source term history along one time step, we then calculate the current term g1 by calling add_nl. In this routine we check the CFL condition, which can mean we may decide we need to reduce the timestep. Once we've dropped the timestep we try to redo this iteration. The next time we visit the section of code above (as a part of the attempt to redo the iteration) we'll end up shuffling the g1, g2, g3 arrays again -- i.e. what started out as g1 will now be in g3, the first calculation of the nonlinear term for this iteration will be in g2 (without center being applied) and another calculation of the nonlinear term for this iteration will be in g1.

There are two possible solutions

1. Move the "shuffling" to after the return, such that we only do it if not changing the time step. This requires allocating a temporary array to allow the shuffling.
2. Move the istep_last = istep statement to the top of this if block. This means we essentially skip the shuffling and recalculation of g1 after changing the timestep. In this case we must move the call to center above the return statement.

Option 2 seems to be the preferred option as it allows us to avoid recalculating the nonlinear term. (Note I think this should be valid as this is an explicit term and hence should only depend on the values of g and the fields at the start of the time step, which shouldn't have changed). In testing both options give slightly different results, suggesting there is some further difference between the two presumably because the g and fields aren't actually the same the second time the add_explicit_terms routine is called. This difference is not explained by force_maxwell_reinit=.true. (which will change the fields slightly when restoring from a restart file).

[Joseph Parker]

I agree that option 2 is better

[David Dickinson]

Difference between two approaches now rectified, bugfix incoming.

[Colin Malcolm Roach]

I assume g1 contain explicit value of (dg/dt)_NL (am I right?).

1. suppose timestep is LC(F)LC N g1 will only be identical on 2nd call (after reset) if g,chi at time we evaluate g1 are identical to values at previous attempt at this timestep. BUT when timestep resets, linear advance is over different dt and this results in DIFFERENT g,chi for explicit evaluation of g1!
   THEREFORE, while Option (2) sounds more efficient, its actually wrong and that OPTION ! is the correct approach.
   i.e. You can't avoid reevaluating the NL term if the fields and g at that time are different.

2. BUT IF timestep is N LC(F)LC then Option 2 would be fine because g,chi and the corresponding evaluation of g1 would always be identical on a reset. In this case pptions 1 and 2 should be identical.

[David Dickinson]

@colinmroach yes g1 is (dg/dt)_Explicit

The timestep is option (2), i.e. the nonlinear term is the first thing calculated. There's some more information in the related PR #108 -- this has resolved the difference between the two approaches to the bug fix (as noted above). In conclusion I think the bug has been fixed now and it leads to slightly different nonlinear behaviour (to be expected) but qualitatively little change.

[David Dickinson]

Fixed by PR #107