From 72bd6496dc75cac72d931490f9e2c68d0d7dc325 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Mon, 8 Jul 2024 16:45:07 +0100 Subject: [PATCH] Add option to simulate a 0D2V plasma with collisions and advection terms 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 | 9 ++++ moment_kinetics/src/geo.jl | 31 +++++++++++++- moment_kinetics/src/input_structs.jl | 6 ++- moment_kinetics/src/time_advance.jl | 11 ++--- moment_kinetics/src/vperp_advection.jl | 42 ++++++++++++++++++- 5 files changed, 90 insertions(+), 9 deletions(-) diff --git a/moment_kinetics/ext/manufactured_solns_ext.jl b/moment_kinetics/ext/manufactured_solns_ext.jl index 314853725..e65c964db 100644 --- a/moment_kinetics/ext/manufactured_solns_ext.jl +++ b/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 diff --git a/moment_kinetics/src/geo.jl b/moment_kinetics/src/geo.jl index 75637425a..3635230dc 100644 --- a/moment_kinetics/src/geo.jl +++ b/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 diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl index 634b2fadd..cddb11226 100644 --- a/moment_kinetics/src/input_structs.jl +++ b/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 diff --git a/moment_kinetics/src/time_advance.jl b/moment_kinetics/src/time_advance.jl index 0fbc1d1ac..f9df6e910 100644 --- a/moment_kinetics/src/time_advance.jl +++ b/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,7 +2085,7 @@ 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 @@ -2093,7 +2093,8 @@ function euler_time_advance!(fvec_out, fvec_in, pdf, fields, moments, # 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 diff --git a/moment_kinetics/src/vperp_advection.jl b/moment_kinetics/src/vperp_advection.jl index 16c2697df..2dc6bd5c8 100644 --- a/moment_kinetics/src/vperp_advection.jl +++ b/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