Addition of "beam" and "beam-with-losses" options to external sources…

… to model injection of particles around specific vpa0, vperp0, with a Gaussian shape function. In the "beam-with-losses" option, a sink is added to permit steady state solutions in a 0D2V simulation.

moment_kinetics/src/external_sources.jl

@@ -48,15 +48,17 @@ function setup_external_sources!(input_dict, r, z)
                      source_n=1.0,
                      source_T=neutrals ? get(input_dict, "T_wall", 1.0) : 1.0,
                      source_v0=0.0, # birth speed for "alphas" option
-                     sink_strength=1.0, # strength of sink in "alphas-with-losses" option
-                     sink_vth=0.0, # thermal speed for sink in "alphas-with-losses" option 
+                     source_vpa0=0.0, # birth vpa for "beam" option
+                     source_vperp0=0.0, # birth vperp for "beam" option
+                     sink_strength=1.0, # strength of sink in "alphas-with-losses" & "beam-with-losses" option
+                     sink_vth=0.0, # thermal speed for sink in "alphas-with-losses" & "beam-with-losses" option 
                      r_profile="constant",
                      r_width=1.0,
                      r_relative_minimum=0.0,
                      z_profile="constant",
                      z_width=1.0,
                      z_relative_minimum=0.0,
-                     source_type="Maxwellian", # "energy", "alphas", "alphas-with-losses"
+                     source_type="Maxwellian", # "energy", "alphas", "alphas-with-losses", "beam", "beam-with-losses"
                      PI_density_controller_P=0.0,
                      PI_density_controller_I=0.0,
                      PI_density_target_amplitude=1.0,
@@ -175,7 +177,7 @@ function setup_external_sources!(input_dict, r, z)
             PI_density_target_ir = nothing
             PI_density_target_iz = nothing
             PI_density_target_rank = nothing
-        elseif input["source_type"] ∈ ("Maxwellian", "energy", "alphas", "alphas-with-losses")
+        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
@@ -564,6 +566,67 @@ function external_ion_source!(pdf, fvec, moments, ion_source_settings, vperp, vp
             error("Unsupported combination in source_type=$(source_type) evolve_density=$(moments.evolve_density), "
                   * "evolve_upar=$(moments.evolve_upar), evolve_ppar=$(moments.evolve_ppar)")
         end
+    elseif source_type == "beam" || source_type == "beam-with-losses"
+        begin_s_r_z_region()
+        source_vpa0 = ion_source_settings.source_vpa0
+        source_vperp0 = ion_source_settings.source_vperp0
+        if !(source_vpa0 > 1.0e-8)
+            error("source_vpa0=$source_vpa0 < 1.0e-8")
+        end
+        if !(source_vperp0 > 1.0e-8)
+            error("source_vperp0=$source_vperp0 < 1.0e-8")
+        end
+        dummy_vpavperp = scratch_dummy.dummy_vpavperp
+        if !moments.evolve_ppar && !moments.evolve_upar && !moments.evolve_density
+            @loop_s_r_z is ir iz begin
+                this_prefactor = dt * source_amplitude[iz,ir]
+                # first assign source to local scratch array
+                @loop_vperp_vpa ivperp ivpa begin
+                    vth0  = sqrt(2.0*source_T) # sqrt(2 T / m), m = mref = 1
+                    v2 = ((vperp_grid[ivperp]-source_vperp0)^2 + (vpa_grid[ivpa]-source_vpa0)^2)/(vth0^2)
+                    dummy_vpavperp[ivpa,ivperp] = (1.0/vth0^3)*exp(-v2)
+                end
+                # get the density for normalisation purposes
+                normfac = get_density(dummy_vpavperp, vpa, vperp)
+                # add the source
+                @loop_vperp_vpa ivperp ivpa begin
+                    # Factor of 1/sqrt(π) (for 1V) or 1/π^(3/2) (for 2V/3V) is absorbed by the
+                    # normalisation of F
+                    pdf[ivpa,ivperp,iz,ir,is] +=
+                        this_prefactor * dummy_vpavperp[ivpa,ivperp] / normfac
+                end
+            end
+
+            if source_type == "beam-with-losses"
+                sink_vth = ion_source_settings.sink_vth
+                sink_strength = ion_source_settings.sink_strength
+                if !(sink_vth > 1.0e-8)
+                   error("sink_vth=$sink_vth < 1.0e-8")
+                end
+                # subtract a sink function representing the loss of slow ash particles
+                @loop_s_r_z is ir iz begin
+                    # first assign sink to local scratch array
+                    @loop_vperp_vpa ivperp ivpa begin
+                        v2 = vperp_grid[ivperp]^2 + vpa_grid[ivpa]^2
+                        fac = 1.0/(sink_vth^2)
+                        dummy_vpavperp[ivpa,ivperp] = (1.0/(sink_vth^3))*exp(-fac*v2)
+                    end
+                    # numerical correction to normalisation
+                    normfac = get_density(dummy_vpavperp, vpa, vperp)
+                    # println("sink norm", normfac)
+                    # add the source
+                    @loop_vperp_vpa ivperp ivpa begin
+                        # Factor of 1/sqrt(π) (for 1V) or 1/π^(3/2) (for 2V/3V) is absorbed by the
+                        # normalisation of F
+                        pdf[ivpa,ivperp,iz,ir,is] -=
+                            dt * sink_strength * dummy_vpavperp[ivpa,ivperp] * pdf[ivpa,ivperp,iz,ir,is] / normfac
+                    end
+                end
+            end
+        else
+            error("Unsupported combination in source_type=$(source_type) evolve_density=$(moments.evolve_density), "
+                  * "evolve_upar=$(moments.evolve_upar), evolve_ppar=$(moments.evolve_ppar)")
+        end
     else
         error("Unsupported source_type=$(source_type) ")
     end
@@ -786,6 +849,10 @@ function external_ion_source_controller!(fvec_in, moments, ion_source_settings,
         # do nothing
     elseif ion_source_settings.source_type == "alphas-with-losses"
         # do nothing
+    elseif ion_source_settings.source_type == "beam"
+        # do nothing
+    elseif ion_source_settings.source_type == "beam-with-losses"
+        # do nothing
     else
         error("Unrecognised source_type=$(ion_source_settings.source_type)")
     end