Address issue #111 by incorporating the true geometric factor in the …

…ExB drift terms. Note that this means that ExB drifts will only be present if bzeta > 0, so the pitch input parameter must be set to less that 1. This changes significantly the usage of the code so far, which assumed that bzeta ~= 1, whilst allowing an input parameter to set bzed.
mabarnes · Dec 6, 2023 · 250a324 · 250a324
1 parent d4ecf2f
commit 250a324
Show file tree

Hide file tree

Showing 3 changed files with 15 additions and 4 deletions.
diff --git a/src/geo.jl b/src/geo.jl
@@ -44,7 +44,8 @@ bzeta::Array{mk_float,2}
 dBdz::Array{mk_float,2}
 # d Bmag d r
 dBdr::Array{mk_float,2}
-
+# jacobian =  grad r x grad z . grad zeta
+jacobian::Array{mk_float,2}
 end
 
 """
@@ -63,6 +64,7 @@ function init_magnetic_geometry(geometry_input_data::geometry_input,z,r)
     bzeta = allocate_float(nz,nr)
     dBdr = allocate_float(nz,nr)
     dBdz = allocate_float(nz,nr)
+    jacobian = allocate_float(nz,nr)
 
     option = geometry_input_data.option
     rhostar = geometry_input_data.rhostar
@@ -77,14 +79,15 @@ function init_magnetic_geometry(geometry_input_data::geometry_input,z,r)
                 Bzeta[iz,ir] = Bmag[iz,ir]*bzeta[iz,ir]
                 dBdr[iz,ir] = 0.0
                 dBdz[iz,ir] = 0.0
+                jacobian[iz,ir] = 1.0
             end
         end
     else 
         input_option_error("$option", option)
     end
 
     geometry = geometric_coefficients(geometry_input_data, rhostar,
-               Bzed,Bzeta,Bmag,bzed,bzeta,dBdz,dBdr)
+               Bzed,Bzeta,Bmag,bzed,bzeta,dBdz,dBdr,jacobian)
     return geometry
 end
 

diff --git a/src/r_advection.jl b/src/r_advection.jl
@@ -95,9 +95,13 @@ function update_speed_r!(advect, upar, vth, fields, evolve_upar, evolve_ppar, vp
     if r.advection.option == "default" && r.n > 1
         Bmag = geometry.Bmag
         ExBfac = 0.5*geometry.rhostar
+        bzeta = geometry.bzeta
+        jacobian = geometry.jacobian
+        geofac = r.scratch
         @inbounds begin
             @loop_z_vperp_vpa iz ivperp ivpa begin
-                @views @. advect.speed[:,ivpa,ivperp,iz] = ExBfac*fields.Ez[iz,:]/Bmag[iz,:]
+                @. geofac = bzeta[iz,:]*jacobian[iz,:]/Bmag[iz,:]
+                @views @. advect.speed[:,ivpa,ivperp,iz] = ExBfac*geofac*fields.Ez[iz,:]
             end
         end
     elseif r.advection.option == "default" && r.n == 1

diff --git a/src/z_advection.jl b/src/z_advection.jl
@@ -94,10 +94,14 @@ function update_speed_z!(advect, upar, vth, evolve_upar, evolve_ppar, fields, vp
         # bzed = B_z/B only used for z.advection.option == "default"
         bzed = geometry.bzed
         Bmag = geometry.Bmag
+        bzeta = geometry.bzeta
+        jacobian = geometry.jacobian
         ExBfac = -0.5*geometry.rhostar
+        geofac = z.scratch
         @inbounds begin
             @loop_r_vperp_vpa ir ivperp ivpa begin
-                @. @views advect.speed[:,ivpa,ivperp,ir] = vpa.grid[ivpa]*bzed[:,ir] + ExBfac*fields.Er[:,ir]/Bmag[:,ir]
+                @. geofac = bzeta[:,ir]*jacobian[:,ir]/Bmag[:,ir]
+                @. @views advect.speed[:,ivpa,ivperp,ir] = vpa.grid[ivpa]*bzed[:,ir] + ExBfac*geofac*fields.Er[:,ir]
             end
             if evolve_ppar
                 @loop_r_vperp_vpa ir ivperp ivpa begin