Option to time-evolve the electron shape function

Uses an implicit timestep so that the electron shape function can be
advanced with the ion timestep. Does not test the accuracy of the
timestep for the electron shape function, so that the error control on
the timestepper does not try to resolve fast electron dynamics.

Can be used to find steady state, because then dg_e/dt=0 in the steady
state anyway, so there is no error. When the time-dependence is needed,
dg_e/dt should be a mass-ratio-small term, so maybe is OK anyway?

moment_kinetics/src/electron_fluid_equations.jl

@@ -127,6 +127,11 @@ end
 
 function calculate_electron_moments!(scratch, pdf, moments, composition, collisions, r, z,
                                      vpa)
+    if length(scratch.pdf_electron) > 0
+        pdf_electron = scratch.pdf_electron
+    else
+        pdf_electron = pdf.electron.norm
+    end
     calculate_electron_density!(scratch.electron_density, moments.electron.dens_updated,
                                 scratch.density)
     calculate_electron_upar_from_charge_conservation!(
@@ -144,7 +149,7 @@ function calculate_electron_moments!(scratch, pdf, moments, composition, collisi
     end
     update_electron_vth_temperature!(moments, scratch.electron_ppar,
                                      scratch.electron_density, composition)
-    calculate_electron_qpar!(moments.electron, pdf.electron, scratch.electron_ppar,
+    calculate_electron_qpar!(moments.electron, pdf_electron, scratch.electron_ppar,
                              scratch.electron_upar, scratch.upar,
                              collisions.electron_fluid.nu_ei, composition.me_over_mi,
                              composition.electron_physics, vpa)

moment_kinetics/src/electron_kinetic_equation.jl

@@ -859,6 +859,10 @@ function electron_backward_euler_pseudotimestepping!(scratch, pdf, moments, phi,
                 # to arrays, because f_electron_old and electron_ppar_old are captured by
                 # residual_func!() above, so any change in the things they refer to will
                 # cause type instability in residual_func!().
+                f_electron_new = @view new_scratch.pdf_electron[:,:,:,ir]
+                f_electron_old = @view old_scratch.pdf_electron[:,:,:,ir]
+                electron_ppar_new = @view new_scratch.electron_ppar[:,ir]
+                electron_ppar_old = @view old_scratch.electron_ppar[:,ir]
                 begin_z_vperp_vpa_region()
                 @loop_z_vperp_vpa iz ivperp ivpa begin
                     f_electron_new[ivpa,ivperp,iz] = f_electron_old[ivpa,ivperp,iz]
@@ -901,7 +905,7 @@ function electron_backward_euler_pseudotimestepping!(scratch, pdf, moments, phi,
                             upper_vcut_changed[] = 0
                         else
                             upper_vcut_changed[] = 1
-                            precon_upperz_vcut_ind[ir] = new_upperz_vcut_ind[]
+                            precon_upperz_vcut_inds[ir] = new_upperz_vcut_ind[]
                         end
                     end
                     MPI.Bcast!(upper_vcut_changed, comm_inter_block[]; root=n_blocks[]-1)
@@ -1003,8 +1007,6 @@ function electron_backward_euler_pseudotimestepping!(scratch, pdf, moments, phi,
                 end
             end
 
-            reset_nonlinear_per_stage_counters!(nl_solver_params)
-
             t_params.step_counter[] += 1
             if electron_pdf_converged[]
                 break
@@ -1927,8 +1929,11 @@ global_rank[] == 0 && println("recalculating precon")
                                  moments.electron.qpar[:,ir], buffer_1, buffer_2,
                                  buffer_3, buffer_4, z_spectral, z)
         end
+
     end
 
+    reset_nonlinear_per_stage_counters!(nl_solver_params)
+
     return newton_success
 end
 

moment_kinetics/src/initial_conditions.jl

@@ -814,6 +814,14 @@ function initialize_electron_pdf!(scratch, scratch_electron, pdf, moments, field
             pdf.electron.pdf_before_ion_timestep[ivpa,ivperp,iz,ir] =
                 pdf.electron.norm[ivpa,ivperp,iz,ir]
         end
+        if length(scratch[1].pdf_electron) > 0
+            @loop_r_z_vperp_vpa ir iz ivperp ivpa begin
+                for i ∈ 1:length(scratch)
+                    scratch[i].pdf_electron[ivpa,ivperp,iz,ir] =
+                        pdf.electron.norm[ivpa,ivperp,iz,ir]
+                end
+            end
+        end
 
         # No need to do electron I/O (apart from possibly debug I/O) any more, so if
         # adaptive timestep is used, it does not need to adjust to output times.

moment_kinetics/src/input_structs.jl

@@ -79,6 +79,7 @@ struct time_info{Terrorsum <: Real, T_debug_output, T_electron, Trkimp, Timpzero
     maximum_dt::mk_float
     implicit_braginskii_conduction::Bool
     implicit_electron_advance::Bool
+    implicit_electron_time_evolving::Bool
     implicit_ion_advance::Bool
     implicit_vpa_advection::Bool
     implicit_electron_ppar::Bool
@@ -130,6 +131,7 @@ struct advance_info
     continuity::Bool
     force_balance::Bool
     energy::Bool
+    electron_pdf::Bool
     electron_energy::Bool
     electron_conduction::Bool
     neutral_external_source::Bool

moment_kinetics/src/moment_kinetics_input.jl

@@ -162,6 +162,7 @@ function mk_input(input_dict=OptionsDict(); save_inputs_to_txt=false, ignore_MPI
         maximum_dt=Inf,
         implicit_braginskii_conduction=true,
         implicit_electron_advance=false,
+        implicit_electron_time_evolving=false,
         implicit_ion_advance=false,
         implicit_vpa_advection=false,
         implicit_electron_ppar=true,

moment_kinetics/src/moment_kinetics_structs.jl

@@ -9,8 +9,8 @@ using ..type_definitions: mk_float
 
 """
 """
-struct scratch_pdf{n_distribution_ion, n_moment, n_moment_electron,
-                   n_distribution_neutral, n_moment_neutral}
+struct scratch_pdf{n_distribution_ion, n_moment, n_distribution_electron,
+                   n_moment_electron, n_distribution_neutral, n_moment_neutral}
     # ions
     pdf::MPISharedArray{mk_float, n_distribution_ion}
     density::MPISharedArray{mk_float, n_moment}
@@ -19,6 +19,7 @@ struct scratch_pdf{n_distribution_ion, n_moment, n_moment_electron,
     pperp::MPISharedArray{mk_float, n_moment}
     temp_z_s::MPISharedArray{mk_float, n_moment}
     # electrons
+    pdf_electron::MPISharedArray{mk_float, n_distribution_electron}
     electron_density::MPISharedArray{mk_float, n_moment_electron}
     electron_upar::MPISharedArray{mk_float, n_moment_electron}
     electron_ppar::MPISharedArray{mk_float, n_moment_electron}