Add option to simulate a 0D2V plasma with collisions and advection te…

…rms in vpa and vperp, for testing of the compatibility of the pseudospectral and the finite-element methods for advection and collisons, respectively. The option is chosen by selecting the [geometry] option = "0D-Spitzer-test", and picking fixed constant values of Ez, dBdr, Er, and dBdr. The vperp_advection!() routine must update the z_advect and r_advect structs if nr == 1 and nz == 1 and the "0D-Spitzer-test" option is being used, leading to some duplication of code. This could be avoided if the speed updates happened outside the function where the advection terms are computed.

moment_kinetics/ext/manufactured_solns_ext.jl

@@ -94,6 +94,15 @@ using IfElse
             end
             dBdz = expand_derivatives(Dz(Bmag))
             jacobian = 1.0
+        elseif option == "0D-Spitzer-test"
+            Bmag = 1.0
+            dBdz = geometry_input_data.dBdz_constant
+            dBdr = geometry_input_data.dBdr_constant
+            bzed = pitch
+            bzeta = sqrt(1 - bzed^2)
+            Bzed = Bmag*bzed
+            Bzeta = Bmag*bzeta
+            jacobian = 1.0
         else
             input_option_error("$option", option)
         end

moment_kinetics/src/geo.jl

@@ -87,7 +87,11 @@ function setup_geometry_input(toml_input::Dict, reference_rhostar)
         # constant for testing nonzero Er when nr = 1
         Er_constant = 0.0,
         # constant for testing nonzero Ez when nz = 1
-        Ez_constant = 0.0)
+        Ez_constant = 0.0,
+        # constant for testing nonzero dBdz when nz = 1
+        dBdz_constant = 0.0,
+        # constant for testing nonzero dBdr when nr = 1
+        dBdr_constant = 0.0)
 
     input = Dict(Symbol(k)=>v for (k,v) in input_section)
     #println(input)
@@ -198,6 +202,31 @@ function init_magnetic_geometry(geometry_input_data::geometry_input,z,r)
                 gbdriftz[iz,ir] = cvdriftz[iz,ir]
             end
         end
+    elseif option == "0D-Spitzer-test"
+     # a 0D configuration with certain geometrical factors
+     # set to be constants to enable testing of velocity
+     # space operators such as mirror or vperp advection terms
+     pitch = geometry_input_data.pitch
+     dBdz_constant = geometry_input_data.dBdz_constant
+     dBdr_constant = geometry_input_data.dBdr_constant
+     B0 = 1.0 # chose reference field strength to be Bzeta at r = 1
+     for ir in 1:nr
+            for iz in 1:nz
+                Bmag[iz,ir] = B0
+                bzed[iz,ir] = pitch
+                bzeta[iz,ir] = sqrt(1 - pitch^2)
+                Bzed[iz,ir] = bzed[iz,ir]*Bmag[iz,ir]
+                Bzeta[iz,ir] = bzeta[iz,ir]*Bmag[iz,ir]
+                dBdz[iz,ir] = dBdz_constant
+                dBdr[iz,ir] = dBdr_constant
+                jacobian[iz,ir] = 1.0
+
+                cvdriftr[iz,ir] = 0.0
+                cvdriftz[iz,ir] = 0.0
+                gbdriftr[iz,ir] = 0.0
+                gbdriftz[iz,ir] = 0.0
+            end
+        end
     else 
         input_option_error("$option", option)
     end

moment_kinetics/src/input_structs.jl

@@ -417,7 +417,11 @@ Base.@kwdef struct geometry_input
     # constant for testing nonzero Er when nr = 1
     Er_constant::mk_float
     # constant for testing nonzero Ez when nz = 1
-    Ez_constant::mk_float  
+    Ez_constant::mk_float
+    # constant for testing nonzero dBdz when nz = 1
+    dBdz_constant::mk_float
+    # constant for testing nonzero dBdr when nr = 1
+    dBdr_constant::mk_float
 end
 
 @enum binary_format_type hdf5 netcdf

moment_kinetics/src/time_advance.jl

@@ -649,8 +649,8 @@ function setup_advance_flags(moments, composition, t_params, collisions,
     # otherwise, check to see if the flags need to be set to true
     if !t_params.split_operators
         # default for non-split operators is to include both vpa and z advection together
-        advance_vpa_advection = vpa.n > 1 && z.n > 1
-        advance_vperp_advection = vperp.n > 1 && z.n > 1
+        advance_vpa_advection = vpa.n > 1 #&& z.n > 1
+        advance_vperp_advection = vperp.n > 1 #&& z.n > 1
         advance_z_advection = z.n > 1
         advance_r_advection = r.n > 1
         if collisions.fkpl.nuii > 0.0 && vperp.n > 1 
@@ -2085,15 +2085,16 @@ function euler_time_advance!(fvec_out, fvec_in, pdf, fields, moments,
     end
 
     # r advection relies on derivatives in z to get ExB
-    if advance.r_advection && r.n > 1
+    if advance.r_advection
         r_advection!(fvec_out.pdf, fvec_in, moments, fields, r_advect, r, z, vperp, vpa,
                      dt, r_spectral, composition, geometry, scratch_dummy)
     end
     # vperp_advection requires information about z and r advection
     # so call vperp_advection! only after z and r advection routines
     if advance.vperp_advection
         vperp_advection!(fvec_out.pdf, fvec_in, vperp_advect, r, z, vperp, vpa,
-                      dt, vperp_spectral, composition, z_advect, r_advect, geometry)
+                      dt, vperp_spectral, composition, z_advect, r_advect, geometry,
+                      moments, fields, t)
     end
 
     if advance.source_terms
@@ -2106,7 +2107,7 @@ function euler_time_advance!(fvec_out, fvec_in, pdf, fields, moments,
             r, z, vzeta, vr, vz, dt, t, z_spectral, composition, scratch_dummy)
     end
 
-    if advance.neutral_r_advection && r.n > 1
+    if advance.neutral_r_advection
         neutral_advection_r!(fvec_out.pdf_neutral, fvec_in, neutral_r_advect,
             r, z, vzeta, vr, vz, dt, r_spectral, composition, geometry, scratch_dummy)
     end

moment_kinetics/src/vperp_advection.jl

@@ -6,13 +6,21 @@ export update_speed_vperp!
 using ..advection: advance_f_local!
 using ..chebyshev: chebyshev_info
 using ..looping
+using ..z_advection: update_speed_z!
+using ..r_advection: update_speed_r!
 
 # do a single stage time advance (potentially as part of a multi-stage RK scheme)
 function vperp_advection!(f_out, fvec_in, vperp_advect, r, z, vperp, vpa,
-                      dt, vperp_spectral, composition, z_advect, r_advect, geometry)
+                      dt, vperp_spectral, composition, z_advect, r_advect, geometry,
+                      moments, fields, t)
+
+    # if appropriate, update z and r speeds
+    update_z_r_speeds!(z_advect, r_advect, fvec_in, moments, fields,
+                            geometry, vpa, vperp, z, r, t)
+
     begin_s_r_z_vpa_region()
     @loop_s is begin
-        # get the updated speed along the r direction using the current f
+        # get the updated speed along the vperp direction using the current f
         @views update_speed_vperp!(vperp_advect[is], vpa, vperp, z, r, z_advect[is], r_advect[is], geometry)
         @loop_r_z_vpa ir iz ivpa begin
             @views advance_f_local!(f_out[ivpa,:,iz,ir,is], fvec_in.pdf[ivpa,:,iz,ir,is],
@@ -58,4 +66,34 @@ function update_speed_vperp!(vperp_advect, vpa, vperp, z, r, z_advect, r_advect,
     return nothing
 end
 
+function update_z_r_speeds!(z_advect, r_advect, fvec_in, moments, fields,
+                            geometry, vpa, vperp, z, r, t)
+    update_z_speed = (z.n == 1 && geometry.input.option == "0D-Spitzer-test")
+    if update_z_speed
+    # make sure z speed is physical despite
+    # 0D nature of simulation
+        begin_s_r_vperp_vpa_region()
+        @loop_s is begin
+            # get the updated speed along the z direction using the current f
+            @views update_speed_z!(z_advect[is], fvec_in.upar[:,:,is],
+                                   moments.ion.vth[:,:,is], moments.evolve_upar,
+                                   moments.evolve_ppar, fields, vpa, vperp, z, r, t, geometry, is)
+        end
+    end
+
+    update_r_speed = (r.n == 1 && geometry.input.option == "0D-Spitzer-test")
+    if update_r_speed
+    # make sure r speed is physical despite
+    # 0D nature of simulation
+        begin_s_z_vperp_vpa_region()
+        @loop_s is begin
+            # get the updated speed along the r direction using the current f
+            update_speed_r!(r_advect[is], fvec_in.upar[:,:,is],
+                               moments.ion.vth[:,:,is], fields, moments.evolve_upar,
+                               moments.evolve_ppar, vpa, vperp, z, r, geometry, is)
+        end
+    end
+    return nothing
+end    
+
 end