Add PI_controller for temperature, update PI_controller for density, and add coll_krook heat flux closure

examples/gk-ions/2D-periodic-gk.toml

@@ -45,7 +45,7 @@ discretization = "chebyshev_pseudospectral"
 n_ion_species = 1
 n_neutral_species = 0
 electron_physics = "boltzmann_electron_response"
-gyrokinetic_ions = true
+ion_physics = "gyrokinetic_ions"
 T_e = 1.0
 T_wall = 1.0
 

moment_kinetics/debug_test/gyroaverage_inputs.jl

@@ -6,7 +6,7 @@ test_input = OptionsDict(
     "output" => OptionsDict("run_name" => "gyroaverage"),
     "composition" => OptionsDict("n_ion_species" => 1,
                           "n_neutral_species" => 0,
-                          "gyrokinetic_ions" => true,
+                          "ion_physics" => "gyrokinetic_ions",
                           "T_e" => 1.0,
                           "T_wall" => 1.0),
     "evolve_moments" => OptionsDict("density" => false,

moment_kinetics/src/collision_frequencies.jl

@@ -0,0 +1,118 @@
+module collision_frequencies
+
+export get_collision_frequency_ii, get_collision_frequency_ee, get_collision_frequency_ei
+
+using ..reference_parameters: get_reference_collision_frequency_ii,
+                              get_reference_collision_frequency_ee,
+                              get_reference_collision_frequency_ei
+using ..reference_parameters: setup_reference_parameters
+
+"""
+    get_collision_frequency_ii(collisions, n, vth)
+
+Calculate the ion-ion collision frequency, depending on the settings/parameters in
+`collisions`, for the given density `n` and thermal speed `vth`.
+
+`n` and `vth` may be scalars or arrays, but should have shapes that can be broadcasted
+together.
+"""
+function get_collision_frequency_ii(collisions, n, vth)
+    # extract krook options from collisions struct
+    colk = collisions.krook
+    nuii0 = colk.nuii0
+    frequency_option = colk.frequency_option
+    if frequency_option ∈ ("reference_parameters", "collisionality_scan")
+        return @. nuii0 * n * vth^(-3)
+    elseif frequency_option == "manual"
+        # Include 0.0*n so that the result gets promoted to an array if n is an array,
+        # which hopefully means this function will have a fixed return type given the
+        # types of the arguments (we don't want to be 'type unstable' for array inputs by
+        # returning a scalar from this branch but an array from the "reference_parameters"
+        # branch).
+        return @. nuii0 + 0.0 * n
+    elseif frequency_option == "none"
+        # Include 0.0*n so that the result gets promoted to an array if n is an array,
+        # which hopefully means this function will have a fixed return type given the
+        # types of the arguments (we don't want to be 'type unstable' for array inputs by
+        # returning a scalar from this branch but an array from the "reference_parameters"
+        # branch).
+        return @. 0.0 * n
+    else
+        error("Unrecognised option [krook_collisions] "
+              * "frequency_option=$(frequency_option)")
+    end
+end
+
+"""
+    get_collision_frequency_ee(collisions, n, vthe)
+
+Calculate the electron-electron collision frequency, depending on the settings/parameters
+in `collisions`, for the given density `n` and electron thermal speed `vthe`.
+
+`n` and `vthe` may be scalars or arrays, but should have shapes that can be broadcasted
+together.
+"""
+function get_collision_frequency_ee(collisions, n, vthe)
+    # extract krook options from collisions struct
+    colk = collisions.krook
+    nuee0 = colk.nuee0
+    frequency_option = colk.frequency_option
+    if frequency_option == "reference_parameters"
+        return @. nuee0 * n * vthe^(-3)
+    elseif frequency_option == "manual"
+        # Include 0.0*n so that the result gets promoted to an array if n is an array,
+        # which hopefully means this function will have a fixed return type given the
+        # types of the arguments (we don't want to be 'type unstable' for array inputs by
+        # returning a scalar from this branch but an array from the "reference_parameters"
+        # branch).
+        return @. nuee0 + 0.0 * n
+    elseif frequency_option == "none"
+        # Include 0.0*n so that the result gets promoted to an array if n is an array,
+        # which hopefully means this function will have a fixed return type given the
+        # types of the arguments (we don't want to be 'type unstable' for array inputs by
+        # returning a scalar from this branch but an array from the "reference_parameters"
+        # branch).
+        return @. 0.0 * n
+    else
+        error("Unrecognised option [krook_collisions] "
+              * "frequency_option=$(frequency_option)")
+    end
+end
+
+"""
+    get_collision_frequency_ei(collisions, n, vthe)
+
+Calculate the electron-electron collision frequency, depending on the settings/parameters
+in `collisions`, for the given density `n` and electron thermal speed `vthe`.
+
+`n` and `vthe` may be scalars or arrays, but should have shapes that can be broadcasted
+together.
+"""
+function get_collision_frequency_ei(collisions, n, vthe)
+    # extract krook options from collisions struct
+    colk = collisions.krook
+    nuei0 = colk.nuei0
+    frequency_option = colk.frequency_option
+    if frequency_option == "reference_parameters"
+        return @. nuei0 * n * vthe^(-3)
+    elseif frequency_option == "manual"
+        # Include 0.0*n so that the result gets promoted to an array if n is an array,
+        # which hopefully means this function will have a fixed return type given the
+        # types of the arguments (we don't want to be 'type unstable' for array inputs by
+        # returning a scalar from this branch but an array from the "reference_parameters"
+        # branch).
+        return @. nuei0 + 0.0 * n
+    elseif frequency_option == "none"
+        # Include 0.0*n so that the result gets promoted to an array if n is an array,
+        # which hopefully means this function will have a fixed return type given the
+        # types of the arguments (we don't want to be 'type unstable' for array inputs by
+        # returning a scalar from this branch but an array from the "reference_parameters"
+        # branch).
+        return @. 0.0 * n
+    else
+        error("Unrecognised option [krook_collisions] "
+              * "frequency_option=$(frequency_option)")
+    end
+end
+
+end

moment_kinetics/src/electron_kinetic_equation.jl

@@ -38,8 +38,9 @@ using ..em_fields: update_phi!
 using ..external_sources: total_external_electron_sources!,
                           add_total_external_electron_source_to_Jacobian!
 using ..file_io: get_electron_io_info, write_electron_state, finish_electron_io
-using ..krook_collisions: electron_krook_collisions!, get_collision_frequency_ee,
-                          get_collision_frequency_ei,
+using ..collision_frequencies: get_collision_frequency_ee,
+                                 get_collision_frequency_ei
+using ..krook_collisions: electron_krook_collisions!,
                           add_electron_krook_collisions_to_Jacobian!
 using ..timer_utils
 using ..moment_constraints: hard_force_moment_constraints!,

moment_kinetics/src/em_fields.jl

@@ -158,7 +158,7 @@ function update_phi!(fields, fvec, vperp, z, r, composition, collisions, moments
     end
 
     # get gyroaveraged field arrays for distribution function advance
-    gkions = composition.gyrokinetic_ions
+    gkions = composition.ion_physics == gyrokinetic_ions
     if gkions
         gyroaverage_field!(fields.gphi,fields.phi,gyroavs,vperp,z,r,composition)
         gyroaverage_field!(fields.gEz,fields.Ez,gyroavs,vperp,z,r,composition)

moment_kinetics/src/external_sources.jl

@@ -74,6 +74,9 @@ function setup_external_sources!(input_dict, r, z, electron_physics)
                      PI_density_target_z_profile="constant",
                      PI_density_target_z_width=1.0,
                      PI_density_target_z_relative_minimum=0.0,
+                     PI_temperature_controller_P=0.0,
+                     PI_temperature_controller_I=0.0,
+                     PI_temperature_target_amplitude=1.0,
                      recycling_controller_fraction=0.0,
                     )
 
@@ -104,6 +107,10 @@ function setup_external_sources!(input_dict, r, z, electron_physics)
             PI_density_target_ir = nothing
             PI_density_target_iz = nothing
             PI_density_target_rank = nothing
+            PI_temperature_target = nothing
+            PI_temperature_target_ir = nothing
+            PI_temperature_target_iz = nothing
+            PI_temperature_target_rank = nothing
         elseif input["source_type"] == "density_midpoint_control"
             PI_density_target = input["PI_density_target_amplitude"]
 
@@ -143,23 +150,76 @@ function setup_external_sources!(input_dict, r, z, electron_physics)
             else
                 PI_density_target_rank = nothing
             end
+            PI_temperature_target = nothing
+            PI_temperature_target_ir = nothing
+            PI_temperature_target_iz = nothing
+            PI_temperature_target_rank = nothing
+        elseif input["source_type"] == "temperature_midpoint_control"
+            PI_temperature_target = input["PI_temperature_target_amplitude"]
+            PI_density_target = nothing
+            PI_density_target_ir = nothing
+            PI_density_target_iz = nothing
+            PI_density_target_rank = nothing
+
+            if comm_block[] != MPI.COMM_NULL
+                PI_controller_amplitude = allocate_shared_float(1)
+                controller_source_profile = allocate_shared_float(z.n, r.n)
+            else
+                PI_controller_amplitude = allocate_float(1)
+                controller_source_profile = allocate_float(z.n, r.n)
+            end
+            for ir ∈ 1:r.n, iz ∈ 1:z.n
+                controller_source_profile[iz,ir] = r_amplitude[ir] * z_amplitude[iz]
+            end
+
+            # Find the indices, and process rank of the point at r=0, z=0.
+            # The result of findfirst() will be `nothing` if the point was not found.
+            PI_temperature_target_ir = findfirst(x->abs(x)<1.e-14, r.grid)
+            PI_temperature_target_iz = findfirst(x->abs(x)<1.e-14, z.grid)
+            if block_rank[] == 0
+                # Only need to do communications from the root process of each
+                # shared-memory block
+                if PI_temperature_target_ir !== nothing && PI_temperature_target_iz !== nothing
+                    PI_temperature_target_rank = iblock_index[]
+                else
+                    PI_temperature_target_rank = 0
+                end
+                if comm_inter_block[] != MPI.COMM_NULL
+                    PI_temperature_target_rank = MPI.Allreduce(PI_temperature_target_rank, +,
+                                                           comm_inter_block[])
+                end
+                if PI_temperature_target_rank == 0 && iblock_index[] == 0 &&
+                        (PI_temperature_target_ir === nothing ||
+                         PI_temperature_target_iz === nothing)
+                    error("No grid point with r=0 and z=0 was found for the "
+                          * "'temperature_midpoint' controller.")
+                end
+            else
+                PI_temperature_target_rank = nothing
+            end
         elseif input["source_type"] ∈ ("Maxwellian", "energy", "alphas", "alphas-with-losses", "beam", "beam-with-losses")
             PI_density_target = nothing
             PI_controller_amplitude = nothing
             controller_source_profile = nothing
             PI_density_target_ir = nothing
             PI_density_target_iz = nothing
             PI_density_target_rank = nothing
+            PI_temperature_target = nothing
+            PI_temperature_target_ir = nothing
+            PI_temperature_target_iz = nothing
+            PI_temperature_target_rank = nothing
         else
             error("Unrecognised ion source_type=$(input["source_type"])."
                   * "Possible values are: Maxwellian, density_profile_control, "
-                  * "density_midpoint_control, energy, alphas, alphas-with-losses, "
-                  * "beam, beam-with-losses")
+                  * "density_midpoint_control, temperature_midpoint_control, energy, "
+                  * "alphas, alphas-with-losses, beam, beam-with-losses")
         end
         return ion_source_data(; Dict(Symbol(k)=>v for (k,v) ∈ input)..., r_amplitude,
                 z_amplitude=z_amplitude, PI_density_target=PI_density_target,
                 PI_controller_amplitude, controller_source_profile,
-                PI_density_target_ir, PI_density_target_iz, PI_density_target_rank)
+                PI_density_target_ir, PI_density_target_iz, PI_density_target_rank,
+                PI_temperature_target, PI_temperature_target_ir, PI_temperature_target_iz,
+                PI_temperature_target_rank)
     end
 
     function get_settings_neutrals(source_index, active_flag)
@@ -328,7 +388,7 @@ function setup_external_sources!(input_dict, r, z, electron_physics)
                      source_strength=ion_settings.source_strength,
                      source_T=ion_settings.source_T,
                     )
-        if ion_settings.source_type != "energy"
+        if ion_settings.source_type != "energy" && ion_settings.source_type != "temperature_midpoint_control"
             # Need to keep same amplitude for ions and electrons so there is no charge
             # source.
             if input["source_strength"] != ion_settings.source_strength
@@ -413,6 +473,9 @@ function get_source_profile(profile_type, width, relative_minimum, coord)
         L = coord.L
         return @. (1.0 - relative_minimum) * exp(-(x+0.5*L) / width) + relative_minimum +
                   (1.0 - relative_minimum) * exp(-(0.5*L-x) / width) + relative_minimum
+    elseif profile_type == "super_gaussian_4"
+        x = coord.grid
+        return @. (1.0 - relative_minimum) * exp(-(x / width)^4) + relative_minimum
     else
         error("Unrecognised source profile type '$profile_type'.")
     end
@@ -447,7 +510,7 @@ function initialize_external_source_amplitude!(moments, external_source_settings
 
     for index ∈ eachindex(ion_source_settings)
         if ion_source_settings[index].active
-            if ion_source_settings[index].source_type == "energy"
+            if ion_source_settings[index].source_type ∈ ("energy", "temperature_midpoint_control")
                 @loop_r_z ir iz begin
                     moments.ion.external_source_amplitude[iz,ir,index] =
                         ion_source_settings[index].source_strength *
@@ -471,7 +534,7 @@ function initialize_external_source_amplitude!(moments, external_source_settings
                 if moments.evolve_ppar
                     @loop_r_z ir iz begin
                         moments.ion.external_source_pressure_amplitude[iz,ir,index] =
-                            (0.5 * ion_source.source_T +
+                            (0.5 * ion_source_settings[index].source_T +
                             moments.ion.upar[iz,ir]^2 - moments.ion.ppar[iz,ir]) *
                             ion_source_settings[index].source_strength *
                             ion_source_settings[index].r_amplitude[ir] *
@@ -665,7 +728,7 @@ function initialize_external_source_controller_integral!(
         for index ∈ eachindex(ion_source_settings)
             if ion_source_settings[index].active && 
                ion_source_settings[index].PI_density_controller_I != 0.0 && 
-               ion_source_settings[index].source_type ∈ ("density_profile_control", "density_midpoint_control")
+               ion_source_settings[index].source_type ∈ ("density_profile_control", "density_midpoint_control", "temperature_midpoint_control")
                 moments.ion.external_source_controller_integral[:, :, index] .= 0.0
             end
         end
@@ -721,7 +784,7 @@ Add external source term to the ion kinetic equation.
     end
     vpa_grid = vpa.grid
     vperp_grid = vperp.grid
-    if source_type in ("Maxwellian","energy")
+    if source_type in ("Maxwellian","energy","density_midpoint_control","density_profile_control","temperature_midpoint_control")
         begin_s_r_z_vperp_region()
         if moments.evolve_ppar && moments.evolve_upar && moments.evolve_density
             vth = moments.ion.vth
@@ -785,7 +848,7 @@ Add external source term to the ion kinetic equation.
                   * "evolve_upar=$(moments.evolve_upar), evolve_ppar=$(moments.evolve_ppar)")
         end
 
-        if source_type == "energy"
+        if source_type ∈ ("energy", "temperature_midpoint_control")
             if moments.evolve_density
                 # Take particles out of pdf so source does not change density
                 @loop_s_r_z_vperp_vpa is ir iz ivperp ivpa begin
@@ -1254,7 +1317,7 @@ source amplitude.
                     dt * ion_source_settings.PI_density_controller_I * n_error
 
                 # Only want a source, so never allow amplitude to be negative
-                amplitude = max(
+                amplitude = max(ion_source_settings.source_strength +
                     ion_source_settings.PI_density_controller_P * n_error +
                     ion_moments.external_source_controller_integral[1,1,index],
                     0)
@@ -1286,6 +1349,65 @@ source amplitude.
                     amplitude * ion_source_settings.controller_source_profile[iz,ir]
             end
         end
+    elseif ion_source_settings.source_type == "temperature_midpoint_control"
+        begin_serial_region()
+        ion_moments.temp .= 2 .* ppar ./ density
+        # controller_amplitude error is a shared memory Vector of length 1
+        controller_amplitude = ion_source_settings.PI_controller_amplitude
+        @serial_region begin
+            if ion_source_settings.PI_temperature_target_ir !== nothing &&
+                    ion_source_settings.PI_temperature_target_iz !== nothing
+                # This process has the target point
+
+                T_mid = ion_moments.temp[ion_source_settings.PI_temperature_target_iz,
+                                ion_source_settings.PI_temperature_target_ir, is]
+                T_error = ion_source_settings.PI_temperature_target - T_mid
+
+                ion_moments.external_source_controller_integral[1,1,index] +=
+                    dt * ion_source_settings.PI_temperature_controller_I * T_error
+
+                # Only want a source, so never allow amplitude to be negative
+                amplitude = max(ion_source_settings.source_strength +
+                    ion_source_settings.PI_temperature_controller_P * T_error +
+                    ion_moments.external_source_controller_integral[1,1,index],
+                    0)
+            else
+                amplitude = nothing
+            end
+            controller_amplitude[1] =
+                MPI.Bcast(amplitude, ion_source_settings.PI_temperature_target_rank,
+                          comm_inter_block[])
+        end
+
+        begin_r_z_region()
+
+        amplitude = controller_amplitude[1]
+        @loop_r_z ir iz begin
+            ion_moments.external_source_amplitude[iz,ir,index] =
+                amplitude * ion_source_settings.controller_source_profile[iz,ir]
+        end
+
+        if moments.evolve_upar
+            @loop_r_z ir iz begin
+                ion_moments.external_source_momentum_amplitude[iz,ir,index] =
+                    - density[iz,ir] * upar[iz,ir] * amplitude *
+                    ion_source_settings.controller_source_profile[iz,ir]
+            end
+        end
+        if moments.evolve_ppar
+            @loop_r_z ir iz begin
+                ion_moments.external_source_pressure_amplitude[iz,ir,index] =
+                    ((0.5 * ion_source_settings.source_T + 2 * upar[iz,ir]^2) *
+                    amplitude) * ion_source_settings.controller_source_profile[iz,ir]
+            end
+        end
+        #if moments.evolve_ppar
+        #    @loop_r_z ir iz begin
+        #        ion_moments.external_source_pressure_amplitude[iz,ir,index] =
+        #            (0.5 * ion_source_settings.source_T + upar[iz,ir]^2 - ppar[iz,ir]) *
+        #            amplitude * ion_source_settings.controller_source_profile[iz,ir]
+        #    end
+        #end
     elseif ion_source_settings.source_type == "density_profile_control"
         begin_r_z_region()