From 1f167137d44ed84af686403cff4b83d86a745183 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Sat, 7 Dec 2024 18:02:55 +0000 Subject: [PATCH] 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? --- .../src/electron_fluid_equations.jl | 7 +- .../src/electron_kinetic_equation.jl | 11 +- moment_kinetics/src/initial_conditions.jl | 8 + moment_kinetics/src/input_structs.jl | 2 + moment_kinetics/src/moment_kinetics_input.jl | 1 + .../src/moment_kinetics_structs.jl | 5 +- moment_kinetics/src/time_advance.jl | 150 +++++++++++++++--- 7 files changed, 152 insertions(+), 32 deletions(-) diff --git a/moment_kinetics/src/electron_fluid_equations.jl b/moment_kinetics/src/electron_fluid_equations.jl index bd2f1cdbb..406c778e2 100644 --- a/moment_kinetics/src/electron_fluid_equations.jl +++ b/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) diff --git a/moment_kinetics/src/electron_kinetic_equation.jl b/moment_kinetics/src/electron_kinetic_equation.jl index 6a9ce4379..11888c357 100644 --- a/moment_kinetics/src/electron_kinetic_equation.jl +++ b/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 diff --git a/moment_kinetics/src/initial_conditions.jl b/moment_kinetics/src/initial_conditions.jl index df191e544..1b852d787 100644 --- a/moment_kinetics/src/initial_conditions.jl +++ b/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. diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl index 241dd6958..3b3134db0 100644 --- a/moment_kinetics/src/input_structs.jl +++ b/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 diff --git a/moment_kinetics/src/moment_kinetics_input.jl b/moment_kinetics/src/moment_kinetics_input.jl index aaf124309..25b4c8e36 100644 --- a/moment_kinetics/src/moment_kinetics_input.jl +++ b/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, diff --git a/moment_kinetics/src/moment_kinetics_structs.jl b/moment_kinetics/src/moment_kinetics_structs.jl index 18b5b0189..a20a9ac60 100644 --- a/moment_kinetics/src/moment_kinetics_structs.jl +++ b/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} diff --git a/moment_kinetics/src/time_advance.jl b/moment_kinetics/src/time_advance.jl index 27e0e29ac..efacb7c6e 100644 --- a/moment_kinetics/src/time_advance.jl +++ b/moment_kinetics/src/time_advance.jl @@ -46,7 +46,9 @@ using ..charge_exchange: ion_charge_exchange_collisions_1V!, neutral_charge_exchange_collisions_1V!, ion_charge_exchange_collisions_3V!, neutral_charge_exchange_collisions_3V! -using ..electron_kinetic_equation: update_electron_pdf!, implicit_electron_advance! +using ..electron_kinetic_equation: update_electron_pdf!, implicit_electron_advance!, + electron_backward_euler!, + electron_kinetic_equation_euler_update! using ..ionization: ion_ionization_collisions_1V!, neutral_ionization_collisions_1V!, ion_ionization_collisions_3V!, neutral_ionization_collisions_3V! using ..krook_collisions: krook_collisions! @@ -384,6 +386,7 @@ function setup_time_info(t_input, n_variables, code_time, dt_reload, # Not an IMEX scheme, so cannot have any implicit terms t_input["implicit_braginskii_conduction"] = false t_input["implicit_electron_advance"] = false + t_input["implicit_electron_time_evolving"] = false t_input["implicit_ion_advance"] = false t_input["implicit_vpa_advection"] = false t_input["implicit_electron_ppar"] = false @@ -394,6 +397,7 @@ function setup_time_info(t_input, n_variables, code_time, dt_reload, if composition.electron_physics ∉ (kinetic_electrons, kinetic_electrons_with_temperature_equation) t_input["implicit_electron_advance"] = false + t_input["implicit_electron_time_evolving"] = false t_input["implicit_electron_ppar"] = false end end @@ -436,6 +440,7 @@ function setup_time_info(t_input, n_variables, code_time, dt_reload, elseif electron === false debug_io = nothing implicit_electron_ppar = false + implicit_electron_time_evolving = false electron_preconditioner_type = nothing decrease_dt_iteration_threshold = -1 increase_dt_iteration_threshold = typemax(mk_int) @@ -449,6 +454,8 @@ function setup_time_info(t_input, n_variables, code_time, dt_reload, implicit_electron_ppar = (t_input["implicit_electron_ppar"] !== false) implicit_electron_advance = (t_input["implicit_electron_advance"] !== false) + implicit_electron_time_evolving = (t_input["implicit_electron_time_evolving"] !== false) + electron_precon_types = Dict("lu" => :electron_lu, "adi" => :electron_adi) if implicit_electron_ppar if t_input["implicit_electron_ppar"] === true if block_size[] == 1 @@ -459,7 +466,6 @@ function setup_time_info(t_input, n_variables, code_time, dt_reload, electron_preconditioner_type = Val(:electron_adi) end else - electron_precon_types = Dict("lu" => :electron_lu, "adi" => :electron_adi) if t_input["implicit_electron_ppar"] ∈ keys(electron_precon_types) electron_preconditioner_type = Val(electron_precon_types[t_input["implicit_electron_ppar"]]) else @@ -476,6 +482,26 @@ function setup_time_info(t_input, n_variables, code_time, dt_reload, * "implicit_electron_advance.") end electron_preconditioner_type = Val(:electron_lu) + elseif implicit_electron_time_evolving + if t_input["implicit_electron_time_evolving"] === true + if block_size[] == 1 + # No need to parallelise, so un-split LU solver should be most efficient. + electron_preconditioner_type = Val(:electron_lu) + else + # Want to parallelise preconditioner, so use ADI method. + electron_preconditioner_type = Val(:electron_adi) + end + else + if t_input["implicit_electron_time_evolving"] ∈ keys(electron_precon_types) + electron_preconditioner_type = Val(electron_precon_types[t_input["implicit_electron_time_evolving"]]) + else + precon_keys = collect(keys(electron_precon_types)) + error("Unrecognised option implicit_electron_time_evolving=" + * "\"$(t_input["implicit_electron_time_evolving"])\" which " + * "should be either false/true or a string giving the type of " + * "preconditioner to use - one of $precon_keys.") + end + end else electron_preconditioner_type = Val(:none) end @@ -506,6 +532,7 @@ function setup_time_info(t_input, n_variables, code_time, dt_reload, mk_float(t_input["minimum_dt"]), mk_float(t_input["maximum_dt"]), electron !== nothing && t_input["implicit_braginskii_conduction"], electron !== nothing && implicit_electron_advance, + electron !== nothing && implicit_electron_time_evolving, electron !== nothing && t_input["implicit_ion_advance"], electron !== nothing && t_input["implicit_vpa_advection"], electron !== nothing && implicit_electron_ppar, @@ -720,7 +747,7 @@ function setup_time_advance!(pdf, fields, vz, vr, vzeta, vpa, vperp, z, r, gyrop default_rtol=t_params.rtol / 10.0, default_atol=t_params.atol / 10.0) nl_solver_electron_advance_params = - setup_nonlinear_solve(t_params.implicit_electron_advance || composition.electron_physics ∈ (kinetic_electrons, kinetic_electrons_with_temperature_equation), + setup_nonlinear_solve(t_params.implicit_electron_advance || t_params.implicit_electron_time_evolving || t_params.implicit_electron_ppar, input_dict, (z=z, vperp=vperp, vpa=vpa), (r,); @@ -749,9 +776,10 @@ function setup_time_advance!(pdf, fields, vz, vr, vzeta, vpa, vperp, z, r, gyrop error("Cannot use implicit_ion_advance and implicit_vpa_advection at the same " * "time") end - if t_params.implicit_electron_advance && t_params.implicit_electron_ppar - error("Cannot use implicit_electron_advance and implicit_electron_ppar at the " - * "same time.") + if (t_params.implicit_electron_advance + t_params.implicit_electron_time_evolving + + t_params.implicit_electron_ppar) > 1 + error("Can only use one of implicit_electron_advance, " + * "implicit_electron_time_evolving, and implicit_electron_ppar.") end nl_solver_params = (electron_conduction=electron_conduction_nl_solve_parameters, electron_advance=nl_solver_electron_advance_params, @@ -770,9 +798,11 @@ function setup_time_advance!(pdf, fields, vz, vr, vzeta, vpa, vperp, z, r, gyrop # create an array of structs containing scratch arrays for the pdf and low-order moments # that may be evolved separately via fluid equations - scratch = setup_scratch_arrays(moments, pdf, t_params.n_rk_stages + 1) + scratch = setup_scratch_arrays(moments, pdf, t_params.n_rk_stages + 1, + t_params.implicit_electron_time_evolving) if t_params.rk_coefs_implicit !== nothing - scratch_implicit = setup_scratch_arrays(moments, pdf, t_params.n_rk_stages) + scratch_implicit = setup_scratch_arrays(moments, pdf, t_params.n_rk_stages, + t_params.implicit_electron_time_evolving) else scratch_implicit = nothing end @@ -1169,6 +1199,7 @@ function setup_advance_flags(moments, composition, t_params, collisions, advance_continuity = false advance_force_balance = false advance_energy = false + advance_electron_pdf = false advance_electron_energy = false advance_electron_conduction = false advance_neutral_z_advection = false @@ -1292,7 +1323,9 @@ function setup_advance_flags(moments, composition, t_params, collisions, # moment-kinetically then advance the electron energy equation if composition.electron_physics ∈ (kinetic_electrons, kinetic_electrons_with_temperature_equation) - if !(t_params.implicit_electron_advance || t_params.implicit_electron_ppar) + if !(t_params.implicit_electron_advance || + t_params.implicit_electron_time_evolving || + t_params.implicit_electron_ppar) advance_electron_energy = true advance_electron_conduction = true end @@ -1337,8 +1370,8 @@ function setup_advance_flags(moments, composition, t_params, collisions, advance_external_source, advance_ion_numerical_dissipation, advance_neutral_numerical_dissipation, advance_sources, advance_continuity, advance_force_balance, advance_energy, - advance_electron_energy, advance_electron_conduction, - advance_neutral_external_source, + advance_electron_pdf, advance_electron_energy, + advance_electron_conduction, advance_neutral_external_source, advance_neutral_sources, advance_neutral_continuity, advance_neutral_force_balance, advance_neutral_energy, manufactured_solns_test, r_diffusion, vpa_diffusion, @@ -1377,6 +1410,7 @@ function setup_implicit_advance_flags(moments, composition, t_params, collisions advance_continuity = false advance_force_balance = false advance_energy = false + advance_electron_pdf = false advance_electron_energy = false advance_electron_conduction = false advance_neutral_z_advection = false @@ -1451,10 +1485,14 @@ function setup_implicit_advance_flags(moments, composition, t_params, collisions advance_electron_conduction = true end - if (t_params.implicit_electron_advance || t_params.implicit_electron_ppar) + if (t_params.implicit_electron_advance || t_params.implicit_electron_time_evolving || + t_params.implicit_electron_ppar) advance_electron_energy = true advance_electron_conduction = true end + if t_params.implicit_electron_time_evolving + advance_electron_pdf = true + end manufactured_solns_test = manufactured_solns_input.use_for_advance @@ -1469,11 +1507,12 @@ function setup_implicit_advance_flags(moments, composition, t_params, collisions advance_external_source, advance_ion_numerical_dissipation, advance_neutral_numerical_dissipation, advance_sources, advance_continuity, advance_force_balance, advance_energy, - advance_electron_energy, advance_electron_conduction, - advance_neutral_external_source, advance_neutral_sources, - advance_neutral_continuity, advance_neutral_force_balance, - advance_neutral_energy, manufactured_solns_test, r_diffusion, - vpa_diffusion, vperp_diffusion, vz_diffusion) + advance_electron_pdf, advance_electron_energy, + advance_electron_conduction, advance_neutral_external_source, + advance_neutral_sources, advance_neutral_continuity, + advance_neutral_force_balance, advance_neutral_energy, + manufactured_solns_test, r_diffusion, vpa_diffusion, + vperp_diffusion, vz_diffusion) end function setup_dummy_and_buffer_arrays(nr, nz, nvpa, nvperp, nvz, nvr, nvzeta, @@ -1674,15 +1713,21 @@ end create an array of structs containing scratch arrays for the normalised pdf and low-order moments that may be evolved separately via fluid equations """ -function setup_scratch_arrays(moments, pdf, n) +function setup_scratch_arrays(moments, pdf, n, time_evolve_electrons) # will create n structs, each of which will contain one pdf, # density, parallel flow, parallel pressure, and perpendicular pressure array for ions # (possibly) the same for electrons, and the same for neutrals. The actual array will # be created at the end of the first step of the loop below, once we have a # `scratch_pdf` object of the correct type. - scratch = Vector{scratch_pdf{5,3,2,6,3}}(undef, n) + scratch = Vector{scratch_pdf{5,3,4,2,6,3}}(undef, n) pdf_dims = size(pdf.ion.norm) moment_dims = size(moments.ion.dens) + + if time_evolve_electrons + pdf_electron_dims = size(pdf.electron.norm) + else + pdf_electron_dims = (0,0,0,0) + end moment_electron_dims = size(moments.electron.dens) pdf_neutral_dims = size(pdf.neutral.norm) @@ -1698,6 +1743,7 @@ function setup_scratch_arrays(moments, pdf, n) pperp_array = allocate_shared_float(moment_dims...) temp_array = allocate_shared_float(moment_dims...) + pdf_electron_array = allocate_shared_float(pdf_electron_dims...) density_electron_array = allocate_shared_float(moment_electron_dims...) upar_electron_array = allocate_shared_float(moment_electron_dims...) ppar_electron_array = allocate_shared_float(moment_electron_dims...) @@ -1712,11 +1758,11 @@ function setup_scratch_arrays(moments, pdf, n) scratch[istage] = scratch_pdf(pdf_array, density_array, upar_array, ppar_array, pperp_array, temp_array, - density_electron_array, upar_electron_array, - ppar_electron_array, pperp_electron_array, - temp_electron_array, pdf_neutral_array, - density_neutral_array, uz_neutral_array, - pz_neutral_array) + pdf_electron_array, density_electron_array, + upar_electron_array, ppar_electron_array, + pperp_electron_array, temp_electron_array, + pdf_neutral_array, density_neutral_array, + uz_neutral_array, pz_neutral_array) @serial_region begin scratch[istage].pdf .= pdf.ion.norm scratch[istage].density .= moments.ion.dens @@ -1724,6 +1770,9 @@ function setup_scratch_arrays(moments, pdf, n) scratch[istage].ppar .= moments.ion.ppar scratch[istage].pperp .= moments.ion.pperp + if time_evolve_electrons + scratch[istage].pdf_electron .= pdf.electron.norm + end scratch[istage].electron_density .= moments.electron.dens scratch[istage].electron_upar .= moments.electron.upar scratch[istage].electron_ppar .= moments.electron.ppar @@ -2678,6 +2727,7 @@ appropriate. scratch[t_params.n_rk_stages+1].ppar, scratch[t_params.n_rk_stages+1].pperp, scratch[t_params.n_rk_stages+1].temp_z_s, + scratch[t_params.n_rk_stages+1].pdf_electron, scratch[t_params.n_rk_stages+1].electron_density, scratch[t_params.n_rk_stages+1].electron_upar, scratch[t_params.n_rk_stages+1].electron_ppar, @@ -3011,6 +3061,12 @@ end first_scratch.pperp[iz,ir,is] = moments.ion.pperp[iz,ir,is] end + if length(first_scratch.pdf_electron) > 0 + begin_r_z_vperp_vpa_region() + @loop_r_z_vperp_vpa ir iz ivperp ivpa begin + first_scratch.pdf_electron[ivpa,ivperp,iz,ir] = pdf.electron.norm[ivpa,ivperp,iz,ir] + end + end begin_r_z_region() @loop_r_z ir iz begin first_scratch.electron_density[iz,ir] = moments.electron.dens[iz,ir] @@ -3093,7 +3149,9 @@ end # The result of the implicit solve gives the state vector at 'istage' # which is used as input to the explicit part of the IMEX time step. old_scratch = scratch_implicit[istage] - update_electrons = !(t_params.implicit_electron_advance || t_params.implicit_electron_ppar) + update_electrons = !(t_params.implicit_electron_advance || + t_params.implicit_electron_time_evolving || + t_params.implicit_electron_ppar) bcs_constraints_success = apply_all_bcs_constraints_update_moments!( scratch_implicit[istage], pdf, moments, fields, boundary_distributions, scratch_electron, vz, vr, vzeta, vpa, vperp, @@ -3142,7 +3200,9 @@ end || (istage == n_rk_stages && t_params.implicit_coefficient_is_zero[1]) || t_params.implicit_coefficient_is_zero[istage+1]) update_electrons = (t_params.rk_coefs_implicit === nothing - || !(t_params.implicit_electron_advance || t_params.implicit_electron_ppar)) + || !(t_params.implicit_electron_advance || + t_params.implicit_electron_time_evolving || + t_params.implicit_electron_ppar)) diagnostic_moments = diagnostic_checks && istage == n_rk_stages bcs_constraints_success = apply_all_bcs_constraints_update_moments!( scratch[istage+1], pdf, moments, fields, boundary_distributions, @@ -3233,6 +3293,12 @@ end end # No need to synchronize here as we only change electron quantities and previous # region only changed ion quantities. + if length(final_scratch.pdf_electron) > 0 + begin_r_z_vperp_vpa_region() + @loop_r_z_vperp_vpa ir iz ivperp ivpa begin + pdf.electron.norm[ivpa,ivperp,iz,ir] = final_scratch.pdf_electron[ivpa,ivperp,iz,ir] + end + end begin_r_z_region(no_synchronize=true) @loop_r_z ir iz begin moments.electron.dens[iz,ir] = final_scratch.electron_density[iz,ir] @@ -3545,6 +3611,18 @@ with fvec_in an input and fvec_out the output z_spectral, composition, external_source_settings.neutral, num_diss_params) end + + if advance.electron_pdf + electron_t_params = (dt=Ref(dt),) + for ir ∈ 1:r.n + @views electron_kinetic_equation_euler_update!( + fvec_out.pdf_electron[:,:,:,ir], fvec_out.electron_ppar[:,ir], + fvec_in.pdf_electron[:,:,:,ir], fvec_in.electron_ppar[:,ir], + moments, z, vperp, vpa, z_spectral, vpa_spectral, z_advect, + vpa_advect, scratch_dummy, collisions, composition, + external_source_settings, num_diss_params, electron_t_params, ir) + end + end if advance.electron_energy electron_energy_equation!(fvec_out.electron_ppar, fvec_in.electron_ppar, fvec_in.density, fvec_in.electron_upar, fvec_in.density, @@ -3600,6 +3678,20 @@ end nl_solver_params.electron_advance) success = success && (electron_success == "") + elseif t_params.implicit_electron_time_evolving + begin_serial_region() + @serial_region begin + t_params.electron.dt[] = dt + end + for ir ∈ 1:r.n + electron_success = electron_backward_euler!( + fvec_in, fvec_out, moments, fields.phi, collisions, composition, r, + z, vperp, vpa, z_spectral, vperp_spectral, vpa_spectral, z_advect, + vpa_advect, scratch_dummy, t_params.electron, + external_source_settings, num_diss_params, + nl_solver_params.electron_advance, ir; evolve_ppar=true) + success = success && electron_success + end elseif t_params.implicit_electron_ppar max_electron_pdf_iterations = t_params.electron.max_pseudotimesteps max_electron_sim_time = t_params.electron.max_pseudotime @@ -3931,6 +4023,12 @@ function update_solution_vector!(new_evolved, old_evolved, moments, composition, new_evolved.upar[iz,ir,is] = old_evolved.upar[iz,ir,is] new_evolved.ppar[iz,ir,is] = old_evolved.ppar[iz,ir,is] end + if length(new_evolved.pdf_electron) > 0 + begin_r_z_vperp_vpa_region() + @loop_r_z_vperp_vpa ir iz ivperp ivpa begin + new_evolved.pdf_electron[ivpa,ivperp,iz,ir] = old_evolved.pdf_electron[ivpa,ivperp,iz,ir] + end + end begin_r_z_region() @loop_r_z ir iz begin new_evolved.electron_density[iz,ir] = old_evolved.electron_density[iz,ir]