Initial phi field not written out correctly when force_maxwell_reinit = .false.

[Stephen Biggs-Fox]

I am interested in phi_t(t, kx, ky, theta, ri). However, I noticed that phi_t(t = 0, ...) = 0 everywhere in all of my data files, which seems wrong given that my input file says to initialize with noise. After some experimentation and debugging, it appears that the problem occurs only when force_maxwell_reinit = .false. as follows:

• phi is initialized as per the specified ginit_option in the init_g module. This is used to initialize g, which is stored in dist_fn_arrays::g and dist_fn_arrays::gnew. However, the initialized phi is not stored but rather is just a local variable within the relevant ginit procedure so deallocated when the procedure returns.
• When force_maxwell_reinit = .true., the initialized g is later used to back-calculate phi and the correct values are written to NetCDF. But when force_maxwell_reinit = .false., the back-calculation is not done and phi is written out as zeros.

This feels like a bug. Regardless of whether the user wants to use force_maxwell_reinit or not, the initial conditions should always be written out correctly.

As a work around for now, I can re-run initial runs with force_maxwell_reinit = .true. for nstep = 0 just to get the initial phi field (and moments thereof). However, this is pretty annoying as I will have to run loads of extra (albeit cheap) jobs and combine data from multiple output files.

NB: I have tested this on next but it probably affects releases from 8.0.4 onward as that was the first release to include the recent force_maxwell_reinit bug fix - before the fix, gs2 ignored the input file value and just always had force_maxwell_reinit = .true.

P.S. For convenience, I have attached zero_phi_init.in which is a fairly minimal working example to reproduce this issue. I have set nstep = 0 and nwrite = 1 to ensure that we write out the initialized state only (you can try nstep = 1 too to see what happens after one timestep), ntheta = 9 and nperiod = 1 just so that it is a small problem that runs quickly, clean_init = chop_side = .false. just to make the correct output clearer and write_phi_over_time = .true. to get the problematic variable phi_t. Obviously, just switch between force_maxwell_reinit = .true. and force_maxwell_reinit = .false. to see correct / incorrect output in the phi_tvariable respectively.

[David Dickinson]

Yes this seems like a bug. We should always call the routine to calculate the initial fields whatever force_maxwell_reinit says – that flag is only supposed to apply when reloading from a restart file.

[David Dickinson]

So looking at gs2_init I think I can see the issue. The routine set_initial_field_and_dist_fn_values is designed to set the initial fields and distribution function, but it has been written with time step changes in mind (i.e. being able to read values from restart files etc.) as well. There are two if blocks, one for dealing with setting g and the other for setting the fields. The fields block has a path through which is a no-op if force_maxwell_reinit = .false. and current%initval_ov%override .and. current%initval_ov%in_memory is false, such that the fields will be left empty.

I think that rather than doing if (force_maxwell_reinit) then set the fields from g we should do if (force_maxwell_reinit .or. .not. current%initval_ov%override) then set the fields from g (There’s a neater implementation that should do the same thing which I’m testing).

I’m also not convinced that if we’re changing the timestep and using in_memory = .false. this does the correct thing, at least I can’t see where the fields are set in this path.

[David Dickinson]

I need to do some more testing but hopefully branch bugfix/fix_issue_141_force_maxwell_false_leading_to_empty_initial_fields should address this issue.

The fields are restored in gs2_restore_many if changing the timestep and using in_memory = .false..

‌

[Joseph Parker]

Is this fully resolved by PR #402?

[Stephen Biggs-Fox]

Maybe… I’ve done some testing and the results are not what I expected but might still be OK. I am discussing them with David on Friday so I will report back after that…

[David Dickinson]

Should be fixed in release 8.1