Fix electron pressure evolution

.github/workflows/parallel_test.yml

@@ -10,7 +10,7 @@ permissions:
 jobs:
   test-ubuntu:
     runs-on: ubuntu-latest
-    timeout-minutes: 150
+    timeout-minutes: 180
 
     steps:
       - uses: actions/checkout@v4

examples/kinetic-electrons/nonuniform_periodic_split3_kinetic-coarse-tails-compressed-z4-PareschiRusso2222.toml

@@ -0,0 +1,152 @@
+[evolve_moments]
+parallel_pressure = true
+density = true
+moments_conservation = true
+parallel_flow = true
+
+[reactions]
+electron_ionization_frequency = 0.0
+ionization_frequency = 0.5
+charge_exchange_frequency = 0.75
+
+[r]
+ngrid = 1
+nelement = 1
+
+[z]
+ngrid = 5
+discretization = "gausslegendre_pseudospectral"
+nelement = 32
+#nelement_local = 4
+bc = "periodic"
+element_spacing_option = "compressed_4"
+
+[vpa]
+ngrid = 6
+discretization = "gausslegendre_pseudospectral"
+nelement = 31
+L = 30.0
+element_spacing_option = "coarse_tails"
+bc = "zero"
+
+[vz]
+ngrid = 6
+discretization = "gausslegendre_pseudospectral"
+nelement = 31
+L = 36.0
+element_spacing_option = "coarse_tails"
+bc = "zero"
+
+[composition]
+T_e = 0.2
+electron_physics = "kinetic_electrons"
+n_ion_species = 1
+n_neutral_species = 1
+
+[ion_species_1]
+initial_temperature = 0.2
+initial_density = 1.0
+
+[z_IC_ion_species_1]
+initialization_option = "sinusoid"
+density_amplitude = 0.0 #0.2
+temperature_amplitude = 0.3
+density_phase = 0.0
+upar_amplitude = 0.0 #0.1
+temperature_phase = 1.0
+upar_phase = 2.0
+
+[vpa_IC_ion_species_1]
+initialization_option = "gaussian"
+density_amplitude = 1.0
+temperature_amplitude = 0.0
+density_phase = 0.0
+upar_amplitude = 0.0
+temperature_phase = 0.0
+upar_phase = 0.0
+
+[neutral_species_1]
+initial_temperature = 0.2
+initial_density = 1.0
+
+[z_IC_neutral_species_1]
+initialization_option = "sinusoid"
+temperature_amplitude = 0.0
+density_amplitude = 0.0
+density_phase = 0.0
+upar_amplitude = 0.0
+temperature_phase = 0.0
+upar_phase = 0.0
+
+[vz_IC_neutral_species_1]
+initialization_option = "gaussian"
+density_amplitude = 1.0
+temperature_amplitude = 0.0
+density_phase = 0.0
+upar_amplitude = 0.0
+temperature_phase = 0.0
+upar_phase = 0.0
+
+[krook_collisions]
+use_krook = true
+
+[timestepping]
+type = "PareschiRusso2(2,2,2)"
+implicit_electron_advance = false
+implicit_electron_ppar = true
+implicit_ion_advance = false
+implicit_vpa_advection = false
+nstep = 100000
+dt = 1.0e-5
+nwrite = 1000
+nwrite_dfns = 1000
+steady_state_residual = true
+converged_residual_value = 1.0e-3
+
+#write_after_fixed_step_count = true
+#nstep = 1
+#nwrite = 1
+#nwrite_dfns = 1
+
+[electron_timestepping]
+nstep = 5000000
+#nstep = 1
+dt = 2.0e-8
+#maximum_dt = 1.0e-8
+nwrite = 10 #10000
+nwrite_dfns = 10 #100000
+#type = "SSPRK4"
+type = "Fekete4(3)"
+rtol = 1.0e-3
+atol = 1.0e-14
+minimum_dt = 1.0e-9
+decrease_dt_iteration_threshold = 100
+increase_dt_iteration_threshold = 20
+cap_factor_ion_dt = 5.0
+initialization_residual_value = 2.5
+converged_residual_value = 1.0e-2
+
+#debug_io = 1
+
+[nonlinear_solver]
+nonlinear_max_iterations = 100
+rtol = 1.0e-6 #1.0e-8
+atol = 1.0e-14 #1.0e-16
+linear_restart = 5
+preconditioner_update_interval = 100
+
+[ion_numerical_dissipation]
+#vpa_dissipation_coefficient = 1.0e-1
+#vpa_dissipation_coefficient = 1.0e-2
+#vpa_dissipation_coefficient = 1.0e-3
+force_minimum_pdf_value = 0.0
+
+[electron_numerical_dissipation]
+#vpa_dissipation_coefficient = 2.0
+force_minimum_pdf_value = 0.0
+
+[neutral_numerical_dissipation]
+#vz_dissipation_coefficient = 1.0e-1
+#vz_dissipation_coefficient = 1.0e-2
+#vz_dissipation_coefficient = 1.0e-3
+force_minimum_pdf_value = 0.0

moment_kinetics/src/electron_kinetic_equation.jl

@@ -1029,7 +1029,7 @@ function electron_backward_euler!(scratch, pdf, moments, phi, collisions, compos
                 end
 
                 if nl_solver_params.solves_since_precon_update[] ≥ nl_solver_params.preconditioner_update_interval
-println("recalculating precon")
+global_rank[] == 0 && println("recalculating precon")
                     nl_solver_params.solves_since_precon_update[] = 0
                     nl_solver_params.precon_dt[] = t_params.dt[]
 
@@ -1335,26 +1335,26 @@ println("recalculating precon")
                     if t_params.dt[] < t_params.previous_dt[]
                         # Had to decrease timestep on the first step to get convergence,
                         # so start next ion timestep with the decreased value.
-                        print("decreasing previous_dt due to failures ", t_params.previous_dt[])
+                        global_rank[] == 0 && print("decreasing previous_dt due to failures ", t_params.previous_dt[])
                         t_params.previous_dt[] = t_params.dt[]
-                        println(" -> ", t_params.previous_dt[])
+                        global_rank[] == 0 && println(" -> ", t_params.previous_dt[])
                     #elseif nl_solver_params.max_linear_iterations_this_step[] > max(0.4 * nl_solver_params.nonlinear_max_iterations, 5)
                     elseif nl_solver_params.max_linear_iterations_this_step[] > t_params.decrease_dt_iteration_threshold
                         # Step succeeded, but took a lot of iterations so decrease initial
                         # step size.
-                        print("decreasing previous_dt due to iteration count ", t_params.previous_dt[])
+                        global_rank[] == 0 && print("decreasing previous_dt due to iteration count ", t_params.previous_dt[])
                         t_params.previous_dt[] /= t_params.max_increase_factor
-                        println(" -> ", t_params.previous_dt[])
+                        global_rank[] == 0 && println(" -> ", t_params.previous_dt[])
                     #elseif nl_solver_params.max_linear_iterations_this_step[] < max(0.1 * nl_solver_params.nonlinear_max_iterations, 2)
                     elseif nl_solver_params.max_linear_iterations_this_step[] < t_params.increase_dt_iteration_threshold && (ion_dt === nothing || t_params.previous_dt[] < t_params.cap_factor_ion_dt * ion_dt)
                         # Only took a few iterations, so increase initial step size.
-                        print("increasing previous_dt due to iteration count ", t_params.previous_dt[])
+                        global_rank[] == 0 && print("increasing previous_dt due to iteration count ", t_params.previous_dt[])
                         if ion_dt === nothing
                             t_params.previous_dt[] *= t_params.max_increase_factor
                         else
                             t_params.previous_dt[] = min(t_params.previous_dt[] * t_params.max_increase_factor, t_params.cap_factor_ion_dt * ion_dt)
                         end
-                        println(" -> ", t_params.previous_dt[])
+                        global_rank[] == 0 && println(" -> ", t_params.previous_dt[])
                     end
                 end
 

moment_kinetics/src/time_advance.jl

@@ -3530,7 +3530,17 @@ end
                                                 nl_solver_params.electron_advance,
                                                 max_electron_pdf_iterations,
                                                 max_electron_sim_time; ion_dt=dt)
+
+        # Update `fvec_out.electron_ppar` with the new electron pressure
+        begin_r_z_region()
+        fvec_out_electron_ppar = fvec_out.electron_ppar
+        moments_electron_ppar = moments.electron.ppar
+        @loop_r_z ir iz begin
+            fvec_out_electron_ppar[iz,ir] = moments_electron_ppar[iz,ir]
+        end
+
         success = (electron_success == "")
+
     elseif advance.electron_conduction
         success = implicit_braginskii_conduction!(fvec_out, fvec_in, moments, z, r, dt,
                                                   z_spectral, composition, collisions,