Porting of commit cd07a26. Refactoring of the update_derived_moments …

…function to use get_density, get_upar and get_ppar functions. In the split moment operation, the moments passed in to these functions along with the distribution are set to the enforced values of (density,upar) = (1.0,0.0) appropriately. This branch now supports the test_velocity_integrals.jl script.

src/velocity_moments.jl

@@ -22,6 +22,12 @@ export update_neutral_pr!
 export update_neutral_pzeta!
 export update_neutral_qz!
 
+# for testing 
+export get_density
+export get_upar
+export get_ppar
+export get_pperp
+
 using ..type_definitions: mk_float
 using ..array_allocation: allocate_shared_float, allocate_bool, allocate_float
 using ..calculus: integral
@@ -431,13 +437,16 @@ function update_density_species!(dens, ff, vpa, vperp, z, r)
     @loop_r_z ir iz begin
         # When evolve_density = false, the evolved pdf is the 'true' pdf, and the vpa
         # coordinate is (dz/dt) / c_s.
-        # Integrating calculates n_s / N_e = (1/√π)∫d(vpa/c_s) (√π f_s c_s / N_e)
-        dens[iz,ir] = integrate_over_vspace(@view(ff[:,:,iz,ir]), 
-            vpa.grid, 0, vpa.wgts, vperp.grid, 0, vperp.wgts)
+        dens[iz,ir] = get_density(@view(ff[:,:,iz,ir]), vpa, vperp)
     end
     return nothing
 end
 
+function get_density(ff, vpa, vperp)
+    # Integrating calculates n_s / N_e = (0/√π)∫d(vpa/c_s) (√π f_s c_s / N_e)
+    return integrate_over_vspace(@view(ff[:,:]), vpa.grid, 0, vpa.wgts, vperp.grid, 0, vperp.wgts)
+end
+
 """
 NB: if this function is called and if upar_updated is false, then
 the incoming pdf is the un-normalized pdf that satisfies int dv pdf = density
@@ -475,35 +484,42 @@ function update_upar_species!(upar, density, ppar, ff, vpa, vperp, z, r, evolve_
         # Integrating calculates
         # (upar_s / vth_s) = (1/√π)∫d(vpa/vth_s) * (vpa/vth_s) * (√π f_s vth_s / n_s)
         # so convert from upar_s / vth_s to upar_s / c_s
+        # we set the input density to get_upar = 1.0 as the normalised distribution has density of 1.0
         @loop_r_z ir iz begin
             vth = sqrt(2.0*ppar[iz,ir]/density[iz,ir])
-            upar[iz,ir] = integrate_over_vspace(@view(ff[:,:,iz,ir]),
-                              vpa.grid, 1, vpa.wgts, vperp.grid, 0, vperp.wgts) * vth
+            upar[iz,ir] = vth*get_upar(@view(ff[:,:,iz,ir]), vpa, vperp, 1.0)
         end
     elseif evolve_density
         # corresponds to case where only the density is evolved separately from the
         # normalised pdf, given by g_s = (√π f_s c_s / n_s); the vpa coordinate is
         # (dz/dt) / c_s.
         # Integrating calculates
         # (upar_s / c_s) = (1/√π)∫d(vpa/c_s) * (vpa/c_s) * (√π f_s c_s / n_s)
+        # we set the input density to get_upar = 1.0 as the normalised distribution has density of 1.0
         @loop_r_z ir iz begin
-            upar[iz,ir] = integrate_over_vspace(@view(ff[:,:,iz,ir]),
-                              vpa.grid, 1, vpa.wgts, vperp.grid, 0, vperp.wgts)
+            upar[iz,ir] = get_upar(@view(ff[:,:,iz,ir]), vpa, vperp, 1.0)
         end
     else
         # When evolve_density = false, the evolved pdf is the 'true' pdf,
         # and the vpa coordinate is (dz/dt) / c_s.
         # Integrating calculates
         # (n_s / N_e) * (upar_s / c_s) = (1/√π)∫d(vpa/c_s) * (vpa/c_s) * (√π f_s c_s / N_e)
         @loop_r_z ir iz begin
-            upar[iz,ir] = integrate_over_vspace(@view(ff[:,:,iz,ir]),
-                              vpa.grid, 1, vpa.wgts, vperp.grid, 0, vperp.wgts) /
-                          density[iz,ir]
+            upar[iz,ir] = get_upar(@view(ff[:,:,iz,ir]), vpa, vperp, density[iz,ir])
         end
     end
     return nothing
 end
 
+function get_upar(ff, vpa, vperp, density)
+    # Integrating calculates
+    # (n_s / N_e) * (upar_s / c_s) = (1/√π)∫d(vpa/c_s) * (vpa/c_s) * (√π f_s c_s / N_e)
+    # so we divide by the density of f_s
+    upar = integrate_over_vspace(@view(ff[:,:]), vpa.grid, 1, vpa.wgts, vperp.grid, 0, vperp.wgts)
+    upar /= density
+    return upar
+end
+
 """
 NB: if this function is called and if ppar_updated is false, then
 the incoming pdf is the un-normalized pdf that satisfies int dv pdf = density
@@ -540,40 +556,46 @@ function update_ppar_species!(ppar, density, upar, ff, vpa, vperp, z, r, evolve_
     if evolve_upar
         # this is the case where the parallel flow and density are evolved separately
         # from the normalized pdf, g_s = (√π f_s c_s / n_s); the vpa coordinate is
-        # ((dz/dt) - upar_s) / c_s>
-        # Integrating calculates (p_parallel/m_s n_s c_s^2) = (1/√π)∫d((vpa-upar_s)/c_s) (1/2)*((vpa-upar_s)/c_s)^2 * (√π f_s c_s / n_s)
-        # so convert from p_s / m_s n_s c_s^2 to ppar_s = p_s / m_s N_e c_s^2
+        # ((dz/dt) - upar_s) / c_s> and so we set upar = 0 in the call to get_ppar
+        # because the mean flow of the normalised ff is zero
         @loop_r_z ir iz begin
-            ppar[iz,ir] = integrate_over_vspace(@view(ff[:,:,iz,ir]), vpa.grid, 2, vpa.wgts, vperp.grid, 0, vperp.wgts) *
-                          density[iz,ir]
+            ppar[iz,ir] = density[iz,ir]*get_ppar(@view(ff[:,:,iz,ir]), vpa, vperp, 0.0)
         end
     elseif evolve_density
         # corresponds to case where only the density is evolved separately from the
         # normalised pdf, given by g_s = (√π f_s c_s / n_s); the vpa coordinate is
         # (dz/dt) / c_s.
-        # Integrating calculates
-        # (p_parallel/m_s n_s c_s^2) + (upar_s/c_s)^2 = (1/√π)∫d(vpa/c_s) (vpa/c_s)^2 * (√π f_s c_s / n_s)
-        # so subtract off the mean kinetic energy and multiply by density to get the
-        # internal energy density (aka pressure)
         @loop_r_z ir iz begin
-            ppar[iz,ir] = (integrate_over_vspace(@view(ff[:,:,iz,ir]), vpa.grid, 2, vpa.wgts, vperp.grid, 0, vperp.wgts) -
-                           upar[iz,ir]^2) * density[iz,ir]
+            ppar[iz,ir] = density[iz,ir]*get_ppar(@view(ff[:,:,iz,ir]), vpa, vperp, upar[iz,ir])
         end
     else
         # When evolve_density = false, the evolved pdf is the 'true' pdf,
         # and the vpa coordinate is (dz/dt) / c_s.
-        # Integrating calculates
-        # (p_parallel/m_s N_e c_s^2) + (n_s/N_e)*(upar_s/c_s)^2 = (1/√π)∫d(vpa/c_s) (vpa/c_s)^2 * (√π f_s c_s / N_e)
-        # so subtract off the mean kinetic energy density to get the internal energy
-        # density (aka pressure)
         @loop_r_z ir iz begin
-            ppar[iz,ir] = integrate_over_vspace(@view(ff[:,:,iz,ir]), vpa.grid, 2, vpa.wgts, vperp.grid, 0, vperp.wgts) -
-                          density[iz,ir]*upar[iz,ir]^2
+            ppar[iz,ir] = get_ppar(@view(ff[:,:,iz,ir]), vpa, vperp, upar[iz,ir])
         end
     end
     return nothing
 end
 
+function get_ppar(ff, vpa, vperp, upar)
+    # Integrating calculates
+    # (p_parallel/m_s N_e c_s^2) = (1/√π)∫d(vpa/c_s) ((vpa-upar)/c_s)^2 * (√π f_s c_s / N_e)
+    # the internal energy density (aka pressure of f_s)
+
+    # modify input vpa.grid to account for the mean flow
+    @. vpa.scratch = vpa.grid - upar
+    norm_fac = 1.0 # normalise to m_s N_e c_s^2
+    #norm_fac = 2.0 # normalise to 0.5 m_s N_e c_s^2 = N_e T_s
+    return norm_fac*integrate_over_vspace(@view(ff[:,:]), vpa.scratch, 2, vpa.wgts, vperp.grid, 0, vperp.wgts)
+end
+
+function get_pperp(ff, vpa, vperp)
+    norm_fac = 0.5 # normalise to m_s N_e c_s^2
+    #norm_fac = 1.0 # normalise to 0.5 m_s N_e c_s^2 = N_e T_s
+    return norm_fac*integrate_over_vspace(@view(ff[:,:]), vpa.grid, 0, vpa.wgts, vperp.grid, 2, vperp.wgts)
+end
+
 """
 NB: the incoming pdf is the normalized pdf
 """

test_velocity_integrals.jl

@@ -0,0 +1,112 @@
+using Printf
+using Plots
+using LaTeXStrings
+using MPI
+using Measures
+
+if abspath(PROGRAM_FILE) == @__FILE__
+    using Pkg
+    Pkg.activate(".")
+
+    import moment_kinetics
+	using moment_kinetics.input_structs: grid_input, advection_input
+	using moment_kinetics.coordinates: define_coordinate
+    using moment_kinetics.velocity_moments: get_density, get_upar, get_ppar, get_pperp
+    using moment_kinetics.array_allocation: allocate_float
+
+    # define inputs needed for the test
+	ngrid = 17 #number of points per element 
+	nelement_local = 20 # number of elements per rank
+	nelement_global = nelement_local # total number of elements 
+	Lvpa = 12.0 #physical box size in reference units 
+	Lvperp = 6.0 #physical box size in reference units 
+	bc = "" #not required to take a particular value, not used 
+	# fd_option and adv_input not actually used so given values unimportant
+	discretization = "chebyshev_pseudospectral"
+    fd_option = "fourth_order_centered"
+	adv_input = advection_input("default", 1.0, 0.0, 0.0)
+	nrank = 1
+    irank = 0
+    comm = MPI.COMM_NULL
+	# create the 'input' struct containing input info needed to create a
+	# coordinate
+    vpa_input = grid_input("vpa", ngrid, nelement_global, nelement_local, 
+		nrank, irank, Lvpa, discretization, fd_option, bc, adv_input,comm)
+	vperp_input = grid_input("vperp", ngrid, nelement_global, nelement_local, 
+		nrank, irank, Lvperp, discretization, fd_option, bc, adv_input,comm)
+	# create the coordinate struct 'x'
+	println("made inputs")
+	println("vpa: ngrid: ",ngrid," nelement: ",nelement_local)
+	println("vperp: ngrid: ",ngrid," nelement: ",nelement_local)
+	vpa, vpa_spectral = define_coordinate(vpa_input)
+	vperp, vperp_spectral = define_coordinate(vperp_input)
+
+    dfn = allocate_float(vpa.n,vperp.n)
+
+    function pressure(ppar,pperp)
+        pres = (1.0/3.0)*(ppar + 2.0*pperp) 
+        return pres
+    end
+    # assign a known isotropic Maxwellian distribution in normalised units
+    dens = 3.0/4.0
+    upar = 2.0/3.0
+    ppar = 2.0/3.0
+    pperp = 2.0/3.0
+    pres = pressure(ppar,pperp) 
+    mass = 1.0
+    vth = sqrt(2.0*pres/(dens*mass))
+    for ivperp in 1:vperp.n
+        for ivpa in 1:vpa.n
+            vpa_val = vpa.grid[ivpa]
+            vperp_val = vperp.grid[ivperp]
+            dfn[ivpa,ivperp] = (dens/vth^3)*exp( - ((vpa_val-upar)^2 + vperp_val^2)/vth^2 )
+        end
+    end
+
+    # now check that we can extract the correct moments from the distribution
+
+    dens_test = get_density(dfn,vpa,vperp)
+    upar_test = get_upar(dfn,vpa,vperp,dens_test)
+    ppar_test = get_ppar(dfn,vpa,vperp,upar_test)
+    pperp_test = get_pperp(dfn,vpa,vperp)
+    pres_test = pressure(ppar_test,pperp_test)
+    # output test results 
+    println("")
+    println("Isotropic Maxwellian")
+    println("dens_test: ", dens_test, " dens: ", dens, " error: ", abs(dens_test-dens))
+    println("upar_test: ", upar_test, " upar: ", upar, " error: ", abs(upar_test-upar))
+    println("pres_test: ", pres_test, " pres: ", pres, " error: ", abs(pres_test-pres))
+    println("")
+
+    # assign a known biMaxwellian distribution in normalised units
+    dens = 3.0/4.0
+    upar = 2.0/3.0
+    ppar = 4.0/5.0
+    pperp = 1.0/4.0 
+    mass = 1.0
+    vthpar = sqrt(2.0*ppar/(dens*mass))
+    vthperp = sqrt(2.0*pperp/(dens*mass))
+    for ivperp in 1:vperp.n
+        for ivpa in 1:vpa.n
+            vpa_val = vpa.grid[ivpa]
+            vperp_val = vperp.grid[ivperp]
+            dfn[ivpa,ivperp] = (dens/(vthpar*vthperp^2))*exp( - ((vpa_val-upar)^2)/vthpar^2 - (vperp_val^2)/vthperp^2 )
+        end
+    end
+
+    # now check that we can extract the correct moments from the distribution
+
+    dens_test = get_density(dfn,vpa,vperp)
+    upar_test = get_upar(dfn,vpa,vperp,dens_test)
+    ppar_test = get_ppar(dfn,vpa,vperp,upar_test)
+    pperp_test = get_pperp(dfn,vpa,vperp)
+    # output test results 
+
+    println("")
+    println("biMaxwellian")
+    println("dens_test: ", dens_test, " dens: ", dens, " error: ", abs(dens_test-dens))
+    println("upar_test: ", upar_test, " upar: ", upar, " error: ", abs(upar_test-upar))
+    println("ppar_test: ", ppar_test, " ppar: ", ppar, " error: ", abs(ppar_test-ppar))
+    println("pperp_test: ", pperp_test, " pperp: ", pperp, " error: ", abs(pperp_test-pperp))
+    println("")
+end