From 83c93accefdee4be76cc33684f865f17e37a6172 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Fri, 28 Jun 2024 11:32:58 +0100 Subject: [PATCH 01/34] Add vperp.bc option "zero-no-regularity", intended for use with the collision operator only (or any diffusion operator that has no flux at vperp = 0 -- so not the numerical vperp diffusion at present -- where zero boundary conditions are imposed at vperp = L but no condition is imposed at vperp = 0. Initial short tests suggests that Maxwellian relaxation problems might run successfully without the regularity conditions -TBC. --- moment_kinetics/src/boundary_conditions.jl | 6 ++++-- moment_kinetics/src/fokker_planck_calculus.jl | 10 +++++++--- moment_kinetics/src/moment_kinetics_input.jl | 4 +++- 3 files changed, 14 insertions(+), 6 deletions(-) diff --git a/moment_kinetics/src/boundary_conditions.jl b/moment_kinetics/src/boundary_conditions.jl index 61dff1a08..7082c1383 100644 --- a/moment_kinetics/src/boundary_conditions.jl +++ b/moment_kinetics/src/boundary_conditions.jl @@ -972,7 +972,7 @@ function enforce_vperp_boundary_condition!(f::AbstractArray{mk_float,5}, bc, vpe end function enforce_vperp_boundary_condition!(f::AbstractArray{mk_float,4}, bc, vperp, vperp_spectral, vperp_advect, diffusion) - if bc == "zero" + if bc == "zero" || bc == "zero-no-regularity" nvperp = vperp.n ngrid = vperp.ngrid # set zero boundary condition @@ -982,7 +982,7 @@ function enforce_vperp_boundary_condition!(f::AbstractArray{mk_float,4}, bc, vpe end end # set regularity condition d F / d vperp = 0 at vperp = 0 - if vperp.discretization == "gausslegendre_pseudospectral" || vperp.discretization == "chebyshev_pseudospectral" + if bc == "zero" && (vperp.discretization == "gausslegendre_pseudospectral" || vperp.discretization == "chebyshev_pseudospectral") D0 = vperp_spectral.radau.D0 buffer = @view vperp.scratch[1:ngrid-1] @loop_r_z_vpa ir iz ivpa begin @@ -992,6 +992,8 @@ function enforce_vperp_boundary_condition!(f::AbstractArray{mk_float,4}, bc, vpe f[ivpa,1,iz,ir] = -sum(buffer)/D0[1] end end + elseif bc == "zero-no-regularity" + # do nothing else println("vperp.bc=\"$bc\" not supported by discretization " * "$(vperp.discretization)") diff --git a/moment_kinetics/src/fokker_planck_calculus.jl b/moment_kinetics/src/fokker_planck_calculus.jl index 300e0ca6e..7838608fe 100644 --- a/moment_kinetics/src/fokker_planck_calculus.jl +++ b/moment_kinetics/src/fokker_planck_calculus.jl @@ -2297,11 +2297,15 @@ function enforce_vpavperp_BCs!(pdf,vpa,vperp,vpa_spectral,vperp_spectral) # set regularity condition d F / d vperp = 0 at vperp = 0 # adjust F(vperp = 0) so that d F / d vperp = 0 at vperp = 0 begin_anyv_vpa_region() - buffer = @view vperp.scratch[1:ngrid_vperp-1] @loop_vpa ivpa begin pdf[ivpa,nvperp] = 0.0 - @views @. buffer = D0[2:ngrid_vperp] * pdf[ivpa,2:ngrid_vperp] - pdf[ivpa,1] = -sum(buffer)/D0[1] + end + if vperp.bc == "zero" + buffer = @view vperp.scratch[1:ngrid_vperp-1] + @loop_vpa ivpa begin + @views @. buffer = D0[2:ngrid_vperp] * pdf[ivpa,2:ngrid_vperp] + pdf[ivpa,1] = -sum(buffer)/D0[1] + end end end diff --git a/moment_kinetics/src/moment_kinetics_input.jl b/moment_kinetics/src/moment_kinetics_input.jl index 56358a27b..8e0772442 100644 --- a/moment_kinetics/src/moment_kinetics_input.jl +++ b/moment_kinetics/src/moment_kinetics_input.jl @@ -1196,6 +1196,8 @@ function check_coordinate_input(coord, coord_name, io) println(io,">$coord_name.bc = 'constant'. enforcing constant incoming BC in $coord_name.") elseif coord.bc == "zero" println(io,">$coord_name.bc = 'zero'. enforcing zero incoming BC in $coord_name. Enforcing zero at both boundaries if diffusion operator is present.") + elseif coord.bc == "zero-no-regularity" + println(io,">$coord_name.bc = 'zero'. enforcing zero incoming BC in $coord_name. Enforcing zero at both boundaries if diffusion operator is present. Do not enforce dF/dcoord = 0 at origin if coord = vperp.") elseif coord.bc == "zero_gradient" println(io,">$coord_name.bc = 'zero_gradient'. enforcing zero gradients at both limits of $coord_name domain.") elseif coord.bc == "both_zero" @@ -1216,7 +1218,7 @@ function check_coordinate_input(coord, coord_name, io) if coord.bc != "zero" && coord.n_global > 1 && global_rank[] == 0 println("WARNING: regularity condition (df/dvperp=0 at vperp=0) not being " - * "imposed. Collisions or vperp-diffusion will be unstable.") + * "imposed. Collisions or vperp-diffusion may be unstable.") end else println(io,">using ", coord.ngrid, " grid points per $coord_name element on ", From 1d38f20ac6e79d3814f0d7a79c565b7655fae0cf Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Fri, 28 Jun 2024 12:54:15 +0100 Subject: [PATCH 02/34] Attempt to add new test for the fokker planck time evolution, currently failing. --- .../fokker_planck_time_evolution_tests.jl | 117 ++++++++++++++---- 1 file changed, 96 insertions(+), 21 deletions(-) diff --git a/moment_kinetics/test/fokker_planck_time_evolution_tests.jl b/moment_kinetics/test/fokker_planck_time_evolution_tests.jl index 08239c0f7..e11563d24 100644 --- a/moment_kinetics/test/fokker_planck_time_evolution_tests.jl +++ b/moment_kinetics/test/fokker_planck_time_evolution_tests.jl @@ -34,7 +34,7 @@ struct expected_data f_ion::Array{mk_float, 3} # vpa, vperp, time end -const expected = +const expected_base = expected_data( [-3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0], [0.155051025721682, 0.644948974278318, 1.000000000000000, 1.500000000000000, 2.000000000000000, 2.500000000000000, 3.000000000000000], @@ -81,34 +81,97 @@ const expected = 0.005877384425022965 0.00466291282312835 0.0028530945920350282 0.0008130106752860425 2.4399432485774198e-5 -9.743480904212476e-5 0.0; 0.00017108987342944414 7.097261590252536e-5 -7.822316408657793e-5 -0.0001534897546375058 -9.087984332447822e-5 -3.3079379573126077e-5 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0]) +const expected_no_regularity = +expected_data( + [-3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0], + [0.155051025721682, 0.644948974278318, 1.000000000000000, 1.500000000000000, 2.000000000000000, 2.500000000000000, 3.000000000000000], + # Expected phi: + [-1.249777922692213, -1.249777922692476], + # Expected n_ion: + [0.286568430139637, 0.286568430139562], + # Expected upar_ion: + [0.0, 0.0], + # Expected ppar_ion: + [0.182900778868619, 0.157103495982094], + # Expected pperp_ion + [0.143950337301367, 0.156848978744563], + # Expected qpar_ion + [0.0, 0.0], + # Expected v_t_ion + [1.046547864034549, 1.046547864034539], + # Expected dSdt + [0.000000000000000, -0.000000000409521], + # Expected f_ion: + [0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 ; + 0.000709315517525 0.000479326544463 0.000267291885601 0.000076580407386 0.000013307679383 0.000001402619088 0.000000000000000 ; + 0.006729798254906 0.004547723633218 0.002535994801781 0.000726574675523 0.000126259746577 0.000013307679383 0.000000000000000 ; + 0.038727315040126 0.026170342585036 0.014593642470202 0.004181148571375 0.000726574675400 0.000076580407373 0.000000000000000 ; + 0.135165490857824 0.091339334982523 0.050934510844468 0.014592981682667 0.002535879973697 0.000267279782810 0.000000000000000 ; + 0.262731408806302 0.177543188036297 0.099005269067555 0.028365484502288 0.004929182099845 0.000519531971048 0.000000000000000 ; + 0.000976245858025 0.000659707199564 0.000367879441171 0.000105399224562 0.000018315638889 0.000001930454136 0.000000000000000 ; + 0.262731408806302 0.177543188036297 0.099005269067555 0.028365484502288 0.004929182099845 0.000519531971048 0.000000000000000 ; + 0.135165490857824 0.091339334982523 0.050934510844468 0.014592981682667 0.002535879973697 0.000267279782810 0.000000000000000 ; + 0.038727315040126 0.026170342585036 0.014593642470202 0.004181148571375 0.000726574675400 0.000076580407373 0.000000000000000 ; + 0.006729798254906 0.004547723633218 0.002535994801781 0.000726574675523 0.000126259746577 0.000013307679383 0.000000000000000 ; + 0.000709315517525 0.000479326544463 0.000267291885601 0.000076580407386 0.000013307679383 0.000001402619088 0.000000000000000 ; + 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 ;;; + 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 ; + 0.000474751734681 0.000190566474175 -0.000068453615894 -0.000220474853312 -0.000166467621731 -0.000070269966866 0.000000000000000 ; + 0.008012848857986 0.005724808854515 0.003340228975343 0.000928890127606 -0.000038607624651 -0.000175709447738 0.000000000000000 ; + 0.034159249923911 0.024646971077049 0.014803095427209 0.004994733815267 0.000871813908241 -0.000279350037847 0.000000000000000 ; + 0.095764438891876 0.069477778980504 0.042139603659514 0.014773029429542 0.003247498909427 -0.000201036843458 0.000000000000000 ; + 0.177293203999877 0.129346984538368 0.079203821037167 0.028272781052345 0.006554590511338 0.000118708200311 0.000000000000000 ; + 0.216600640826244 0.158225686082234 0.097107273632873 0.034812857798807 0.008154333928791 0.000303091030245 0.000000000000000 ; + 0.177293203999877 0.129346984538368 0.079203821037167 0.028272781052345 0.006554590511338 0.000118708200311 0.000000000000000 ; + 0.095764438891876 0.069477778980504 0.042139603659514 0.014773029429542 0.003247498909427 -0.000201036843458 0.000000000000000 ; + 0.034159249923911 0.024646971077049 0.014803095427209 0.004994733815267 0.000871813908241 -0.000279350037847 0.000000000000000 ; + 0.008012848857986 0.005724808854515 0.003340228975343 0.000928890127606 -0.000038607624651 -0.000175709447738 0.000000000000000 ; + 0.000474751734681 0.000190566474175 -0.000068453615894 -0.000220474853312 -0.000166467621731 -0.000070269966866 0.000000000000000 ; + 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000]) ########################################################################################### # to modify the test, with a new expected f, print the new f using the following commands # in an interative Julia REPL. The path is the path to the .dfns file. ########################################################################################## """ +using moment_kinetics.load_data: open_readonly_output_file, load_pdf_data, load_ion_moments_data, load_fields_data fid = open_readonly_output_file(path, "dfns") f_ion_vpavperpzrst = load_pdf_data(fid) f_ion = f_ion_vpavperpzrst[:,:,1,1,1,:] ntind = 2 nvpa = 13 #subject to grid choices nvperp = 7 #subject to grid choices +# pdf +for k in 1:ntind + for i in 1:nvpa-1 + for j in 1:nvperp-1 + @printf("%.15f ", f_ion[i,j,k]) + end + @printf("%.15f ", f_ion[i,nvperp,k]) + print(";\n") + end + for j in 1:nvperp-1 + @printf("%.15f ", f_ion[nvpa,j,k]) + end + @printf("%.15f ", f_ion[nvpa,nvperp,k]) + if k < ntind + print(";;;\n") + end +end +# a moment +n_ion_zrst, upar_ion_zrst, ppar_ion_zrst, pperp_ion_zrst, qpar_ion_zrst, v_t_ion_zrst, dSdt_zrst = load_ion_moments_data(fid,extended_moments=true) for k in 1:ntind - for j in 1:nvperp-1 - for i in 1:nvpa-1 - @printf("%.15f ", f_ion[i,j,k]) - print("; ") - end - @printf("%.15f ", f_ion[nvpa,j,k]) - print(";;\n") - end - for i in 1:nvpa-1 - @printf("%.15f ", f_ion[i,nvperp,k]) - print("; ") - end - @printf("%.15f ", f_ion[nvpa,nvperp,k]) - if k < ntind - print(";;;\n") - end + @printf("%.15f", n_ion_zrst[1,1,1,k]) + if k < ntind + print(", ") + end +end +# a field +phi_zrt, Er_zrt, Ez_zrt = load_fields_data(fid) +for k in 1:ntind + @printf("%.15f", phi_zrt[1,1,k]) + if k < ntind + print(", ") + end end """ # default inputs for tests @@ -175,19 +238,21 @@ test_input_gauss_legendre = Dict("run_name" => "gausslegendre_pseudospectral", Run a sound-wave test for a single set of parameters """ # Note 'name' should not be shared by any two tests in this file -function run_test(test_input, rtol, atol, upar_rtol=nothing; args...) +function run_test(test_input, run_name, vperp_bc, expected, rtol, atol, upar_rtol=nothing; args...) # by passing keyword arguments to run_test, args becomes a Dict which can be used to # update the default inputs # Make a copy to make sure nothing modifies the input Dicts defined in this test # script. test_input = deepcopy(test_input) - + test_input["vperp_bc"] = vperp_bc + if upar_rtol === nothing upar_rtol = rtol end # Convert keyword arguments to a unique name + test_input["run_name"] = run_name name = test_input["run_name"] if length(args) > 0 name = string(name, "_", (string(k, "-", v, "_") for (k, v) in args)...) @@ -207,7 +272,7 @@ function run_test(test_input, rtol, atol, upar_rtol=nothing; args...) #input = test_input input["run_name"] = name - + print(input) # Suppress console output while running quietoutput() do # run simulation @@ -296,6 +361,8 @@ function run_test(test_input, rtol, atol, upar_rtol=nothing; args...) ###################################### @test isapprox(expected.n_ion[tind], n_ion[tind], atol=atol) + println(expected.n_ion[tind]) + println(n_ion[tind]) @test isapprox(expected.upar_ion[tind], upar_ion[tind], atol=atol) @test isapprox(expected.ppar_ion[tind], ppar_ion[tind], atol=atol) @test isapprox(expected.pperp_ion[tind], pperp_ion[tind], atol=atol) @@ -325,7 +392,15 @@ function runtests() # GaussLegendre pseudospectral # Benchmark data is taken from this run (GaussLegendre) @testset "Gauss Legendre base" begin - run_test(test_input_gauss_legendre, 1.e-14, 1.0e-14 ) + run_name = "gausslegendre_pseudospectral" + vperp_bc = "zero" + run_test(test_input_gauss_legendre, run_name, vperp_bc, expected_base, 1.0e-14, 1.0e-14 ) + end + @testset "Gauss Legendre no enforced regularity condition at vperp = 0" begin + run_name = "gausslegendre_pseudospectral_no_regularity" + vperp_bc = "zero-no-regularity" + run_test(test_input_gauss_legendre, run_name, vperp_bc, expected_no_regularity, + 1.0e-13, 1.0e-13) end end end From 6ec3cefcbdce49817a0cf8911a3b1a99ab0fc040 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Fri, 28 Jun 2024 14:31:28 +0100 Subject: [PATCH 03/34] Correct test expected values, refactor test inputs to make better use of existing functionality. --- .../fokker_planck_time_evolution_tests.jl | 75 +++++++++---------- 1 file changed, 36 insertions(+), 39 deletions(-) diff --git a/moment_kinetics/test/fokker_planck_time_evolution_tests.jl b/moment_kinetics/test/fokker_planck_time_evolution_tests.jl index e11563d24..669af32de 100644 --- a/moment_kinetics/test/fokker_planck_time_evolution_tests.jl +++ b/moment_kinetics/test/fokker_planck_time_evolution_tests.jl @@ -86,47 +86,47 @@ expected_data( [-3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0], [0.155051025721682, 0.644948974278318, 1.000000000000000, 1.500000000000000, 2.000000000000000, 2.500000000000000, 3.000000000000000], # Expected phi: - [-1.249777922692213, -1.249777922692476], + [-1.254259088243025, -1.254259088243286], # Expected n_ion: - [0.286568430139637, 0.286568430139562], + [0.285287142537587, 0.285287142537513], # Expected upar_ion: [0.0, 0.0], # Expected ppar_ion: - [0.182900778868619, 0.157103495982094], + [0.182220654804438, 0.156448883483764], # Expected pperp_ion - [0.143950337301367, 0.156848978744563], + [0.143306715174515, 0.156192600834786], # Expected qpar_ion [0.0, 0.0], # Expected v_t_ion - [1.046547864034549, 1.046547864034539], + [1.046701532502699, 1.046701532502689], # Expected dSdt - [0.000000000000000, -0.000000000409521], + [0.000000000000000, -0.000000000865997], # Expected f_ion: [0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 ; - 0.000709315517525 0.000479326544463 0.000267291885601 0.000076580407386 0.000013307679383 0.000001402619088 0.000000000000000 ; - 0.006729798254906 0.004547723633218 0.002535994801781 0.000726574675523 0.000126259746577 0.000013307679383 0.000000000000000 ; - 0.038727315040126 0.026170342585036 0.014593642470202 0.004181148571375 0.000726574675400 0.000076580407373 0.000000000000000 ; - 0.135165490857824 0.091339334982523 0.050934510844468 0.014592981682667 0.002535879973697 0.000267279782810 0.000000000000000 ; - 0.262731408806302 0.177543188036297 0.099005269067555 0.028365484502288 0.004929182099845 0.000519531971048 0.000000000000000 ; - 0.000976245858025 0.000659707199564 0.000367879441171 0.000105399224562 0.000018315638889 0.000001930454136 0.000000000000000 ; - 0.262731408806302 0.177543188036297 0.099005269067555 0.028365484502288 0.004929182099845 0.000519531971048 0.000000000000000 ; - 0.135165490857824 0.091339334982523 0.050934510844468 0.014592981682667 0.002535879973697 0.000267279782810 0.000000000000000 ; - 0.038727315040126 0.026170342585036 0.014593642470202 0.004181148571375 0.000726574675400 0.000076580407373 0.000000000000000 ; - 0.006729798254906 0.004547723633218 0.002535994801781 0.000726574675523 0.000126259746577 0.000013307679383 0.000000000000000 ; - 0.000709315517525 0.000479326544463 0.000267291885601 0.000076580407386 0.000013307679383 0.000001402619088 0.000000000000000 ; + 0.000706724195475 0.000477575434536 0.000266315395816 0.000076300638379 0.000013259062819 0.000001397494940 0.000000000000000 ; + 0.006705212475828 0.004531109564814 0.002526730124657 0.000723920301085 0.000125798485463 0.000013259062819 0.000000000000000 ; + 0.038585833650058 0.026074735222541 0.014540327934436 0.004165873701136 0.000723920300963 0.000076300638366 0.000000000000000 ; + 0.134671678398729 0.091005636629981 0.050748427134097 0.014539667807029 0.002526615411287 0.000266303305116 0.000000000000000 ; + 0.261709398369233 0.176852554997681 0.098620144126568 0.028255144359305 0.004910007858078 0.000517511020834 0.000000000000000 ; + 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 ; + 0.261709398369233 0.176852554997682 0.098620144126568 0.028255144359305 0.004910007858078 0.000517511020834 0.000000000000000 ; + 0.134671678398729 0.091005636629981 0.050748427134097 0.014539667807029 0.002526615411287 0.000266303305116 0.000000000000000 ; + 0.038585833650058 0.026074735222541 0.014540327934436 0.004165873701136 0.000723920300963 0.000076300638366 0.000000000000000 ; + 0.006705212475828 0.004531109564814 0.002526730124657 0.000723920301085 0.000125798485463 0.000013259062819 0.000000000000000 ; + 0.000706724195475 0.000477575434536 0.000266315395816 0.000076300638379 0.000013259062819 0.000001397494940 0.000000000000000 ; 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 ;;; 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 ; - 0.000474751734681 0.000190566474175 -0.000068453615894 -0.000220474853312 -0.000166467621731 -0.000070269966866 0.000000000000000 ; - 0.008012848857986 0.005724808854515 0.003340228975343 0.000928890127606 -0.000038607624651 -0.000175709447738 0.000000000000000 ; - 0.034159249923911 0.024646971077049 0.014803095427209 0.004994733815267 0.000871813908241 -0.000279350037847 0.000000000000000 ; - 0.095764438891876 0.069477778980504 0.042139603659514 0.014773029429542 0.003247498909427 -0.000201036843458 0.000000000000000 ; - 0.177293203999877 0.129346984538368 0.079203821037167 0.028272781052345 0.006554590511338 0.000118708200311 0.000000000000000 ; - 0.216600640826244 0.158225686082234 0.097107273632873 0.034812857798807 0.008154333928791 0.000303091030245 0.000000000000000 ; - 0.177293203999877 0.129346984538368 0.079203821037167 0.028272781052345 0.006554590511338 0.000118708200311 0.000000000000000 ; - 0.095764438891876 0.069477778980504 0.042139603659514 0.014773029429542 0.003247498909427 -0.000201036843458 0.000000000000000 ; - 0.034159249923911 0.024646971077049 0.014803095427209 0.004994733815267 0.000871813908241 -0.000279350037847 0.000000000000000 ; - 0.008012848857986 0.005724808854515 0.003340228975343 0.000928890127606 -0.000038607624651 -0.000175709447738 0.000000000000000 ; - 0.000474751734681 0.000190566474175 -0.000068453615894 -0.000220474853312 -0.000166467621731 -0.000070269966866 0.000000000000000 ; + 0.000473874379555 0.000190726186343 -0.000067413540342 -0.000219129923463 -0.000165609965457 -0.000069940808382 0.000000000000000 ; + 0.007980077110978 0.005701898246888 0.003327438535529 0.000925955494251 -0.000037963012435 -0.000174809935222 0.000000000000000 ; + 0.034010103035913 0.024541441396039 0.014741892515397 0.004975785949797 0.000869171347288 -0.000277723510529 0.000000000000000 ; + 0.095320265225675 0.069160937049138 0.041953067849510 0.014711859721633 0.003235165444330 -0.000199356792174 0.000000000000000 ; + 0.176443375955397 0.128736406466194 0.078839749105573 0.028150317244036 0.006528121645593 0.000119350436884 0.000000000000000 ; + 0.215552595150568 0.157471466895046 0.096656199821596 0.034660213647059 0.008120924795998 0.000303073103709 0.000000000000000 ; + 0.176443375955397 0.128736406466194 0.078839749105573 0.028150317244036 0.006528121645593 0.000119350436884 0.000000000000000 ; + 0.095320265225675 0.069160937049138 0.041953067849510 0.014711859721633 0.003235165444330 -0.000199356792174 0.000000000000000 ; + 0.034010103035913 0.024541441396039 0.014741892515397 0.004975785949797 0.000869171347288 -0.000277723510529 0.000000000000000 ; + 0.007980077110978 0.005701898246888 0.003327438535529 0.000925955494251 -0.000037963012435 -0.000174809935222 0.000000000000000 ; + 0.000473874379555 0.000190726186343 -0.000067413540342 -0.000219129923463 -0.000165609965457 -0.000069940808382 0.000000000000000 ; 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000]) ########################################################################################### # to modify the test, with a new expected f, print the new f using the following commands @@ -191,7 +191,7 @@ test_input_gauss_legendre = Dict("run_name" => "gausslegendre_pseudospectral", "vperp_nelement" => 3, "vperp_L" => 3.0, "vperp_discretization" => "gausslegendre_pseudospectral", - "split_operators" => false, + #"split_operators" => false, "ionization_frequency" => 0.0, "charge_exchange_frequency" => 0.0, "constant_ionization_rate" => false, @@ -238,21 +238,19 @@ test_input_gauss_legendre = Dict("run_name" => "gausslegendre_pseudospectral", Run a sound-wave test for a single set of parameters """ # Note 'name' should not be shared by any two tests in this file -function run_test(test_input, run_name, vperp_bc, expected, rtol, atol, upar_rtol=nothing; args...) +function run_test(test_input, expected, rtol, atol, upar_rtol=nothing; args...) # by passing keyword arguments to run_test, args becomes a Dict which can be used to # update the default inputs # Make a copy to make sure nothing modifies the input Dicts defined in this test # script. test_input = deepcopy(test_input) - test_input["vperp_bc"] = vperp_bc if upar_rtol === nothing upar_rtol = rtol end # Convert keyword arguments to a unique name - test_input["run_name"] = run_name name = test_input["run_name"] if length(args) > 0 name = string(name, "_", (string(k, "-", v, "_") for (k, v) in args)...) @@ -269,10 +267,8 @@ function run_test(test_input, run_name, vperp_bc, expected, rtol, atol, upar_rto # Update default inputs with values to be changed input = merge(test_input, modified_inputs) - #input = test_input input["run_name"] = name - print(input) # Suppress console output while running quietoutput() do # run simulation @@ -361,8 +357,6 @@ function run_test(test_input, run_name, vperp_bc, expected, rtol, atol, upar_rto ###################################### @test isapprox(expected.n_ion[tind], n_ion[tind], atol=atol) - println(expected.n_ion[tind]) - println(n_ion[tind]) @test isapprox(expected.upar_ion[tind], upar_ion[tind], atol=atol) @test isapprox(expected.ppar_ion[tind], ppar_ion[tind], atol=atol) @test isapprox(expected.pperp_ion[tind], pperp_ion[tind], atol=atol) @@ -394,13 +388,16 @@ function runtests() @testset "Gauss Legendre base" begin run_name = "gausslegendre_pseudospectral" vperp_bc = "zero" - run_test(test_input_gauss_legendre, run_name, vperp_bc, expected_base, 1.0e-14, 1.0e-14 ) + run_test(test_input_gauss_legendre, + expected_base, 1.0e-14, 1.0e-14; + vperp_bc="zero") end @testset "Gauss Legendre no enforced regularity condition at vperp = 0" begin run_name = "gausslegendre_pseudospectral_no_regularity" vperp_bc = "zero-no-regularity" - run_test(test_input_gauss_legendre, run_name, vperp_bc, expected_no_regularity, - 1.0e-13, 1.0e-13) + run_test(test_input_gauss_legendre, + expected_no_regularity, + 1.0e-14, 1.0e-14; vperp_bc="zero-no-regularity") end end end From adc1facafdacf62d5e4256767bea498dccc0c665 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Fri, 28 Jun 2024 14:34:15 +0100 Subject: [PATCH 04/34] Change input to match check-in test. --- examples/fokker-planck/fokker-planck-relaxation.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/fokker-planck/fokker-planck-relaxation.toml b/examples/fokker-planck/fokker-planck-relaxation.toml index f2d5651a1..0f633c12b 100644 --- a/examples/fokker-planck/fokker-planck-relaxation.toml +++ b/examples/fokker-planck/fokker-planck-relaxation.toml @@ -13,7 +13,7 @@ initial_temperature1 = 1.0 initial_density2 = 0.5 initial_temperature2 = 1.0 z_IC_option1 = "sinusoid" -z_IC_density_amplitude1 = 0.001 +z_IC_density_amplitude1 = 0.0 z_IC_density_phase1 = 0.0 z_IC_upar_amplitude1 = 0.0 z_IC_upar_phase1 = 0.0 From f2def768aa9d2237229a6ba61325dfae92e528a7 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Fri, 28 Jun 2024 14:55:41 +0100 Subject: [PATCH 05/34] Make sure both vperp_bc options are precompiled. --- util/precompile_run.jl | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/util/precompile_run.jl b/util/precompile_run.jl index f5a7fb512..e927edd36 100644 --- a/util/precompile_run.jl +++ b/util/precompile_run.jl @@ -66,17 +66,24 @@ for input ∈ [base_input, cheb_input, wall_bc_input, wall_bc_cheb_input] push!(inputs_list, x) end -collisions_input = merge(wall_bc_cheb_input, Dict("n_neutral_species" => 0, +collisions_input1 = merge(wall_bc_cheb_input, Dict("n_neutral_species" => 0, "krook_collisions" => Dict{String,Any}("use_krook" => true), "fokker_planck_collisions" => Dict{String,Any}("use_fokker_planck" => true, "self_collisions" => true, "slowing_down_test" => true), "vperp_discretization" => "gausslegendre_pseudospectral", "vpa_discretization" => "gausslegendre_pseudospectral", )) +collisions_input2 = merge(wall_bc_cheb_input, Dict("n_neutral_species" => 0, + "krook_collisions" => Dict{String,Any}("use_krook" => true), + "fokker_planck_collisions" => Dict{String,Any}("use_fokker_planck" => true, "self_collisions" => true, "slowing_down_test" => true), + "vperp_discretization" => "gausslegendre_pseudospectral", + "vpa_discretization" => "gausslegendre_pseudospectral", + "vperp_bc" => "zero-no-regularity", + )) # add an additional input for every geometry option available in addition to the default geo_input1 = merge(wall_bc_cheb_input, Dict("n_neutral_species" => 0, "geometry" => Dict{String,Any}("option" => "1D-mirror", "DeltaB" => 0.5, "pitch" => 0.5, "rhostar" => 1.0))) -push!(inputs_list, collisions_input, geo_input1) +push!(inputs_list, collisions_input1, collisions_input2, geo_input1) for input in inputs_list run_moment_kinetics(input) From 2afbbae96627056960779bc06cad40aa3187fd00 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Fri, 28 Jun 2024 14:58:41 +0100 Subject: [PATCH 06/34] Make conserving correction terms an optional feature through "use_conserving_corrections" in the Fokker Planck input namelist. --- moment_kinetics/src/fokker_planck.jl | 12 +++++++++--- moment_kinetics/src/input_structs.jl | 2 ++ 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/moment_kinetics/src/fokker_planck.jl b/moment_kinetics/src/fokker_planck.jl index 1f520e9f6..fb9c6ef19 100644 --- a/moment_kinetics/src/fokker_planck.jl +++ b/moment_kinetics/src/fokker_planck.jl @@ -90,6 +90,7 @@ function setup_fkpl_collisions_input(toml_input::Dict, reference_params) nuii = -1.0, frequency_option = "reference_parameters", self_collisions = true, + use_conserving_corrections = true, slowing_down_test = false, sd_density = 1.0, sd_temp = 0.01, @@ -257,6 +258,7 @@ function explicit_fp_collisions_weak_form_Maxwellian_cross_species!(pdf_out,pdf_ @boundscheck r.n == size(dSdt,2) || throw(BoundsError(dSdt)) @boundscheck n_ion_species == size(dSdt,3) || throw(BoundsError(dSdt)) + use_conserving_corrections = collisions.fkpl.use_conserving_corrections fkin = collisions.fkpl # masses charge numbers and collision frequencies mref = 1.0 # generalise if multiple ions evolved @@ -282,8 +284,10 @@ function explicit_fp_collisions_weak_form_Maxwellian_cross_species!(pdf_out,pdf_ # enforce the boundary conditions on CC before it is used for timestepping enforce_vpavperp_BCs!(fkpl_arrays.CC,vpa,vperp,vpa_spectral,vperp_spectral) # make sure that the cross-species terms conserve density - density_conserving_correction!(fkpl_arrays.CC, pdf_in[:,:,iz,ir,is], vpa, vperp, + if use_conserving_corrections + density_conserving_correction!(fkpl_arrays.CC, pdf_in[:,:,iz,ir,is], vpa, vperp, fkpl_arrays.S_dummy) + end # advance this part of s,r,z with the resulting sum_s' C[Fs,Fs'] begin_anyv_vperp_vpa_region() CC = fkpl_arrays.CC @@ -326,6 +330,7 @@ function explicit_fokker_planck_collisions_weak_form!(pdf_out,pdf_in,dSdt,compos nuref = collisions.fkpl.nuii # generalise! Zi = collisions.fkpl.Zi # generalise! nussp = nuref*(Zi^4) # include charge number factor for self collisions + use_conserving_corrections = collisions.fkpl.use_conserving_corrections # N.B. parallelisation using special 'anyv' region begin_s_r_z_anyv_region() @loop_s_r_z is ir iz begin @@ -336,9 +341,10 @@ function explicit_fokker_planck_collisions_weak_form!(pdf_out,pdf_in,dSdt,compos # enforce the boundary conditions on CC before it is used for timestepping enforce_vpavperp_BCs!(fkpl_arrays.CC,vpa,vperp,vpa_spectral,vperp_spectral) # make ad-hoc conserving corrections - conserving_corrections!(fkpl_arrays.CC, pdf_in[:,:,iz,ir,is], vpa, vperp, + if use_conserving_corrections + conserving_corrections!(fkpl_arrays.CC, pdf_in[:,:,iz,ir,is], vpa, vperp, fkpl_arrays.S_dummy) - + end # advance this part of s,r,z with the resulting C[Fs,Fs] begin_anyv_vperp_vpa_region() CC = fkpl_arrays.CC diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl index fe7661e9a..1c3f18f3f 100644 --- a/moment_kinetics/src/input_structs.jl +++ b/moment_kinetics/src/input_structs.jl @@ -365,6 +365,8 @@ Base.@kwdef struct fkpl_collisions_input nuii::mk_float # option to determine if self collisions are used (for physics test) self_collisions::Bool + # option to determine if ad-hoc moment_kinetics-style conserving corrections are used + use_conserving_corrections::Bool # option to determine if cross-collisions against fixed Maxwellians are used slowing_down_test::Bool # Setting to switch between different options for Fokker-Planck collision frequency input From 6253195929822c3129bd00e8f27aa4615ae4b318 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 1 Jul 2024 15:48:25 +0000 Subject: [PATCH 07/34] Input file for slowing-down Fokker-Planck relaxation, without enforcing dF/dvperp = 0. Simulation shows reasonable behaviour with density preserved and thermal speed hitting a steady-state value near 0.15, which is held for multiple collision times. Ti/Talpha = 1/40 was chosen to reduce strain on resolution. --- ...planck-relaxation-slowing-down-alphas.toml | 92 +++++++++++++++++++ 1 file changed, 92 insertions(+) create mode 100644 examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas.toml diff --git a/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas.toml b/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas.toml new file mode 100644 index 000000000..24aaf6e13 --- /dev/null +++ b/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas.toml @@ -0,0 +1,92 @@ +# cheap input file for a 0D2V relaxation to a collisional Maxwellian distribution with self-ion collisions and collisions with fixed Maxwellian background of cold ions and electrons. +n_ion_species = 1 +n_neutral_species = 0 +electron_physics = "boltzmann_electron_response" +evolve_moments_density = false +evolve_moments_parallel_flow = false +evolve_moments_parallel_pressure = false +evolve_moments_conservation = false +T_e = 1.0 +T_wall = 1.0 +initial_density1 = 1.0 +initial_temperature1 = 1.0 +z_IC_option1 = "sinusoid" +z_IC_density_amplitude1 = 0.0 +z_IC_density_phase1 = 0.0 +z_IC_upar_amplitude1 = 0.0 +z_IC_upar_phase1 = 0.0 +z_IC_temperature_amplitude1 = 0.0 +z_IC_temperature_phase1 = 0.0 +vpa_IC_option1 = "isotropic-beam" +vpa_IC_v01 = 1.0 +vpa_IC_vth01 = 0.1 +#vpa_IC_option1 = "directed-beam" +#vpa_IC_vpa01 = -1.5 +#vpa_IC_vperp01 = 0.0 +charge_exchange_frequency = 0.0 +ionization_frequency = 0.0 +constant_ionization_rate = false + +z_ngrid = 1 +z_nelement = 1 +z_nelement_local = 1 +z_bc = "wall" +z_discretization = "chebyshev_pseudospectral" +r_ngrid = 1 +r_nelement = 1 +r_nelement_local = 1 +r_bc = "periodic" +r_discretization = "chebyshev_pseudospectral" +vpa_ngrid = 5 +vpa_nelement = 32 +vpa_L = 3.0 +vpa_bc = "zero" +vpa_discretization = "gausslegendre_pseudospectral" +vperp_ngrid = 5 +vperp_nelement = 16 +vperp_L = 1.5 +vperp_discretization = "gausslegendre_pseudospectral" +vperp_bc = "zero-no-regularity" +# Fokker-Planck operator requires the "gausslegendre_pseudospectral +# options for the vpa and vperp grids + +[fokker_planck_collisions] +# nuii sets the normalised input C[F,F] Fokker-Planck collision frequency +nuii = 1.876334222e-3 #(1/nu_alphae, as computed from input diagnostic) +Zi = 2.0 +self_collisions = true +slowing_down_test = true +frequency_option = "manual" +use_fokker_planck = true +use_conserving_corrections = true +sd_density = 1.0 +sd_temp = 0.025 # TD/Ealpha +sd_mi = 0.5 # mD/malpha +sd_me = 0.000013616 # 0.25/1836.0 me/malpha +sd_q = 1.0 + +[ion_source] +active=false +source_strength=0.0 +source_T=0.005 +source_n=1.0 +r_profile="constant" +r_width=1.0 +r_relative_minimum=0.0 +z_profile="gaussian" +z_width=0.1 +z_relative_minimum=0.0 +source_v0 = 1.0 +#source_type="alphas" +source_type="alphas-with-losses" +#source_type="beam-with-losses" +#source_vpa0 = 1.0 +#source_vperp0 = 1.0 +sink_strength = 0.0 +sink_vth=0.07071067811865475 + +[timestepping] +nstep = 50000 +dt = 5.0e-4 +nwrite = 100 +nwrite_dfns = 100 From 5715036b45796d85e95d984ed7c6dfa461ab38b9 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Wed, 3 Jul 2024 12:57:05 +0000 Subject: [PATCH 08/34] Correct typo in input for vpa_IC_v0 parameter, no underscore assumed for other parameters and no underscore being used in the examples files, so remove it in the source for consistency. --- moment_kinetics/src/moment_kinetics_input.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/moment_kinetics/src/moment_kinetics_input.jl b/moment_kinetics/src/moment_kinetics_input.jl index 8e0772442..cddadfec0 100644 --- a/moment_kinetics/src/moment_kinetics_input.jl +++ b/moment_kinetics/src/moment_kinetics_input.jl @@ -548,7 +548,7 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) species.ion[is].vpa_IC.upar_amplitude, species.ion[is].vpa_IC.upar_phase, species.ion[is].vpa_IC.temperature_amplitude, species.ion[is].vpa_IC.temperature_phase, species.ion[is].vpa_IC.monomial_degree, - get(scan_input, "vpa_IC_v0_$is", 0.5*sqrt(vperp.L^2 + (0.5*vpa.L)^2)), + get(scan_input, "vpa_IC_v0$is", 0.5*sqrt(vperp.L^2 + (0.5*vpa.L)^2)), get(scan_input, "vpa_IC_vth0$is", 0.1*sqrt(vperp.L^2 + (0.5*vpa.L)^2)), get(scan_input, "vpa_IC_vpa0$is", 0.25*0.5*abs(vpa.L)), get(scan_input, "vpa_IC_vperp0$is", 0.5*abs(vperp.L))) From 7efa515554d09cfbacccfc9178e984c8b7053434 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Wed, 3 Jul 2024 14:14:45 +0000 Subject: [PATCH 09/34] Prevent animations being made if one of the dimensions has size 1. --- .../plots_post_processing/src/plots_post_processing.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/plots_post_processing/plots_post_processing/src/plots_post_processing.jl b/plots_post_processing/plots_post_processing/src/plots_post_processing.jl index 972b92620..353df4961 100644 --- a/plots_post_processing/plots_post_processing/src/plots_post_processing.jl +++ b/plots_post_processing/plots_post_processing/src/plots_post_processing.jl @@ -2858,7 +2858,7 @@ function plot_ion_pdf(run_name, run_name_label, vpa, vperp, z, r, z_local, r_loc spec_string = [string("_", spec_type)] end # make a gif animation of f(vpa,z,t) at a given (vperp,r) location - if pp.animate_f_vs_vpa_z + if pp.animate_f_vs_vpa_z && (vpa.n > 1 && z_local.n > 1) pdf = load_distributed_ion_pdf_slice(run_name, nblocks, itime_min:iskip:itime_max, n_species, r_local, z_local, vperp, vpa; @@ -2876,7 +2876,7 @@ function plot_ion_pdf(run_name, run_name_label, vpa, vperp, z, r, z_local, r_loc end end # make a gif animation of f(vpa,r,t) at a given (vperp,z) location - if pp.animate_f_vs_vpa_r + if pp.animate_f_vs_vpa_r && (vpa.n > 1 && r_local.n > 1) pdf = load_distributed_ion_pdf_slice(run_name, nblocks, itime_min:iskip:itime_max, n_species, r_local, z_local, vperp, vpa; @@ -2894,7 +2894,7 @@ function plot_ion_pdf(run_name, run_name_label, vpa, vperp, z, r, z_local, r_loc end end # make a gif animation of f(vperp,z,t) at a given (vpa,r) location - if pp.animate_f_vs_vperp_z + if pp.animate_f_vs_vperp_z && (vperp.n > 1 && z_local.n > 1) pdf = load_distributed_ion_pdf_slice(run_name, nblocks, itime_min:iskip:itime_max, n_species, r_local, z_local, vperp, vpa; ivpa=ivpa0, @@ -2908,7 +2908,7 @@ function plot_ion_pdf(run_name, run_name_label, vpa, vperp, z, r, z_local, r_loc end end # make a gif animation of f(vperp,r,t) at a given (vpa,z) location - if pp.animate_f_vs_vperp_r + if pp.animate_f_vs_vperp_r && (vperp.n > 1 && r_local.n > 1) pdf = load_distributed_ion_pdf_slice(run_name, nblocks, itime_min:iskip:itime_max, n_species, r_local, z_local, vperp, vpa; ivpa=ivpa0, @@ -2970,7 +2970,7 @@ function plot_ion_pdf(run_name, run_name_label, vpa, vperp, z, r, z_local, r_loc end end # make a gif animation of f(z,r,t) at a given (vpa,vperp) location - if pp.animate_f_vs_r_z + if pp.animate_f_vs_r_z && (z_local.n > 1 && r_local.n > 1) pdf = load_distributed_ion_pdf_slice(run_name, nblocks, itime_min:iskip:itime_max, n_species, r_local, z_local, vperp, vpa; ivpa=ivpa0, From 41e04fdedb8c6da537c8801b2cca5ef99f545417 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Wed, 3 Jul 2024 14:26:37 +0000 Subject: [PATCH 10/34] Reduce time step size of example input file to avoid numerical instability. --- .../fokker-planck-relaxation-slowing-down-alphas.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas.toml b/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas.toml index 24aaf6e13..dfab272e3 100644 --- a/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas.toml +++ b/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas.toml @@ -87,6 +87,6 @@ sink_vth=0.07071067811865475 [timestepping] nstep = 50000 -dt = 5.0e-4 +dt = 2.5e-4 nwrite = 100 nwrite_dfns = 100 From 797952d5eb632183ebf60e3e0f0f38c10d7127ba Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Thu, 4 Jul 2024 14:53:16 +0000 Subject: [PATCH 11/34] Prevent pointless plots being made when r and z dimension size is 1. --- .../src/plots_post_processing.jl | 63 ++++++++++--------- 1 file changed, 33 insertions(+), 30 deletions(-) diff --git a/plots_post_processing/plots_post_processing/src/plots_post_processing.jl b/plots_post_processing/plots_post_processing/src/plots_post_processing.jl index 353df4961..a29b57dd7 100644 --- a/plots_post_processing/plots_post_processing/src/plots_post_processing.jl +++ b/plots_post_processing/plots_post_processing/src/plots_post_processing.jl @@ -2806,13 +2806,13 @@ function plot_fields_rt(phi, delta_phi, time, itime_min, itime_max, nwrite_movie outfile = string(run_name, "_delta_phi(r0,z0)_vs_t.pdf") trysavefig(outfile) end - if pp.plot_phi_vs_z_t + if pp.plot_phi_vs_z_t && r.n > 1 # make a heatmap plot of ϕ(r,t) heatmap(time, r.grid, phi, xlabel="time", ylabel="r", title="ϕ", c = :deep) outfile = string(run_name, "_phi_vs_r_t.pdf") trysavefig(outfile) end - if pp.animate_phi_vs_z + if pp.animate_phi_vs_z && r.n > 1 # make a gif animation of ϕ(r) at different times anim = @animate for i ∈ itime_min:nwrite_movie:itime_max @views plot(r.grid, phi[:,i], xlabel="r", ylabel="ϕ", ylims = (phimin,phimax)) @@ -2826,10 +2826,11 @@ function plot_fields_rt(phi, delta_phi, time, itime_min, itime_max, nwrite_movie # plot!(exp.(-(phi[cld(nz,2),end] .- phi[izmid:end,end])) .* erfi.(sqrt.(abs.(phi[cld(nz,2),end] .- phi[izmid:end,end])))/sqrt(pi)/0.688, phi[izmid:end,end] .- phi[izmid,end], label = "analytical", linewidth=2) # outfile = string(run_name, "_harrison_comparison.pdf") # trysavefig(outfile) - plot(r.grid, phi[:,end], xlabel="r/Lr", ylabel="eϕ/Te", label="", linewidth=2) - outfile = string(run_name, "_phi(r)_final.pdf") - trysavefig(outfile) - + if r.n > 1 + plot(r.grid, phi[:,end], xlabel="r/Lr", ylabel="eϕ/Te", label="", linewidth=2) + outfile = string(run_name, "_phi(r)_final.pdf") + trysavefig(outfile) + end println("done.") end @@ -3067,15 +3068,16 @@ end function plot_fields_2D(phi, Ez, Er, time, z, r, iz0, ir0, itime_min, itime_max, nwrite_movie, run_name, pp, description) nr = size(r,1) + nz = size(z,1) print("Plotting fields data...") phimin = minimum(phi) phimax = maximum(phi) - if pp.plot_phi_vs_r0_z # plot last timestep phi[z,ir0] + if pp.plot_phi_vs_r0_z && nz > 1 # plot last timestep phi[z,ir0] @views plot(z, phi[:,ir0,end], xlabel=L"z/L_z", ylabel=L"\phi") outfile = string(run_name, "_phi"*description*"(r0,z)_vs_z.pdf") trysavefig(outfile) end - if pp.animate_phi_vs_r_z && nr > 1 + if pp.animate_phi_vs_r_z && nr > 1 && nz > 1 # make a gif animation of ϕ(z) at different times anim = @animate for i ∈ itime_min:nwrite_movie:itime_max @views heatmap(r, z, phi[:,:,i], xlabel="r", ylabel="z", c = :deep, interpolation = :cubic) @@ -3087,7 +3089,7 @@ function plot_fields_2D(phi, Ez, Er, time, z, r, iz0, ir0, end outfile = string(run_name, "_phi(zwall-)_vs_r.gif") trygif(anim, outfile, fps=5) - elseif pp.animate_phi_vs_r_z && nr == 1 # make a gif animation of ϕ(z) at different times + elseif pp.animate_phi_vs_r_z && nr == 1 && nz > 1 # make a gif animation of ϕ(z) at different times anim = @animate for i ∈ itime_min:nwrite_movie:itime_max @views plot(z, phi[:,1,i], xlabel="z", ylabel=L"\widetilde{\phi}", ylims = (phimin,phimax)) end @@ -3096,17 +3098,17 @@ function plot_fields_2D(phi, Ez, Er, time, z, r, iz0, ir0, end Ezmin = minimum(Ez) Ezmax = maximum(Ez) - if pp.plot_Ez_vs_r0_z # plot last timestep Ez[z,ir0] + if pp.plot_Ez_vs_r0_z && nz > 1 # plot last timestep Ez[z,ir0] @views plot(z, Ez[:,ir0,end], xlabel=L"z/L_z", ylabel=L"E_z") outfile = string(run_name, "_Ez"*description*"(r0,z)_vs_z.pdf") trysavefig(outfile) end - if pp.plot_wall_Ez_vs_r && nr > 1 # plot last timestep Ez[z_wall,r] + if pp.plot_wall_Ez_vs_r && nr > 1 && nz > 1 # plot last timestep Ez[z_wall,r] @views plot(r, Ez[1,:,end], xlabel=L"r/L_r", ylabel=L"E_z") outfile = string(run_name, "_Ez"*description*"(r,z_wall-)_vs_r.pdf") trysavefig(outfile) end - if pp.animate_Ez_vs_r_z && nr > 1 + if pp.animate_Ez_vs_r_z && nr > 1 && nz > 1 # make a gif animation of ϕ(z) at different times anim = @animate for i ∈ itime_min:nwrite_movie:itime_max @views heatmap(r, z, Ez[:,:,i], xlabel="r", ylabel="z", c = :deep, interpolation = :cubic) @@ -3118,7 +3120,7 @@ function plot_fields_2D(phi, Ez, Er, time, z, r, iz0, ir0, end outfile = string(run_name, "_Ez(zwall-)_vs_r.gif") trygif(anim, outfile, fps=5) - elseif pp.animate_Ez_vs_r_z && nr == 1 + elseif pp.animate_Ez_vs_r_z && nr == 1 && nz > 1 anim = @animate for i ∈ itime_min:nwrite_movie:itime_max @views plot(z, Ez[:,1,i], xlabel="z", ylabel=L"\widetilde{E}_z", ylims = (Ezmin,Ezmax)) end @@ -3127,12 +3129,12 @@ function plot_fields_2D(phi, Ez, Er, time, z, r, iz0, ir0, end Ermin = minimum(Er) Ermax = maximum(Er) - if pp.plot_Er_vs_r0_z # plot last timestep Er[z,ir0] + if pp.plot_Er_vs_r0_z && nz > 1 # plot last timestep Er[z,ir0] @views plot(z, Er[:,ir0,end], xlabel=L"z/L_z", ylabel=L"E_r") outfile = string(run_name, "_Er"*description*"(r0,z)_vs_z.pdf") trysavefig(outfile) end - if pp.plot_wall_Er_vs_r && nr > 1 # plot last timestep Er[z_wall,r] + if pp.plot_wall_Er_vs_r && nr > 1 && nz > 1 # plot last timestep Er[z_wall,r] @views plot(r, Er[1,:,end], xlabel=L"r/L_r", ylabel=L"E_r") outfile = string(run_name, "_Er"*description*"(r,z_wall-)_vs_r.pdf") trysavefig(outfile) @@ -3142,7 +3144,7 @@ function plot_fields_2D(phi, Ez, Er, time, z, r, iz0, ir0, outfile = string(run_name, "_Er(zwall-)_vs_r.gif") trygif(anim, outfile, fps=5) end - if pp.animate_Er_vs_r_z && nr > 1 + if pp.animate_Er_vs_r_z && nr > 1 && nz > 1 # make a gif animation of ϕ(z) at different times anim = @animate for i ∈ itime_min:nwrite_movie:itime_max @views heatmap(r, z, Er[:,:,i], xlabel="r", ylabel="z", c = :deep, interpolation = :cubic) @@ -3159,6 +3161,7 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure, time, z, r, iz0, ir0, n_ion_species, itime_min, itime_max, nwrite_movie, run_name, pp) nr = size(r,1) + nz = size(z,1) ntime = size(time,1) print("Plotting ion moments data...") for is in 1:n_ion_species @@ -3166,7 +3169,7 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure, # the density densitymin = minimum(density[:,:,is,:]) densitymax = maximum(density) - if pp.plot_density_vs_r0_z # plot last timestep density[z,ir0] + if pp.plot_density_vs_r0_z && nz > 1 # plot last timestep density[z,ir0] @views plot(z, density[:,ir0,is,end], xlabel=L"z/L_z", ylabel=L"n_i") outfile = string(run_name, "_density"*description*"(r0,z)_vs_z.pdf") trysavefig(outfile) @@ -3176,7 +3179,7 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure, outfile = string(run_name, "_density"*description*"(r,z_wall)_vs_r.pdf") trysavefig(outfile) end - if pp.animate_density_vs_r_z && nr > 1 + if pp.animate_density_vs_r_z && nr > 1 && nz > 1 # make a gif animation of ϕ(z) at different times anim = @animate for i ∈ itime_min:nwrite_movie:itime_max @views heatmap(r, z, density[:,:,is,i], xlabel="r", ylabel="z", c = :deep, interpolation = :cubic) @@ -3184,7 +3187,7 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure, outfile = string(run_name, "_density"*description*"_vs_r_z.gif") trygif(anim, outfile, fps=5) end - if pp.plot_density_vs_r_z && nr > 1 + if pp.plot_density_vs_r_z && nr > 1 && nz > 1 @views heatmap(r, z, density[:,:,is,end], xlabel=L"r", ylabel=L"z", c = :deep, interpolation = :cubic, windowsize = (360,240), margin = 15pt) outfile = string(run_name, "_density"*description*"_vs_r_z.pdf") @@ -3201,17 +3204,17 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure, # the parallel flow parallel_flowmin = minimum(parallel_flow[:,:,is,:]) parallel_flowmax = maximum(parallel_flow) - if pp.plot_parallel_flow_vs_r0_z # plot last timestep parallel_flow[z,ir0] + if pp.plot_parallel_flow_vs_r0_z && nz > 1 # plot last timestep parallel_flow[z,ir0] @views plot(z, parallel_flow[:,ir0,is,end], xlabel=L"z/L_z", ylabel=L"u_{i\|\|}") outfile = string(run_name, "_parallel_flow"*description*"(r0,z)_vs_z.pdf") trysavefig(outfile) end - if pp.plot_wall_parallel_flow_vs_r && nr > 1 # plot last timestep parallel_flow[z_wall,r] + if pp.plot_wall_parallel_flow_vs_r && nr > 1 && nz > 1# plot last timestep parallel_flow[z_wall,r] @views plot(r, parallel_flow[end,:,is,end], xlabel=L"r/L_r", ylabel=L"u_{i\|\|}") outfile = string(run_name, "_parallel_flow"*description*"(r,z_wall)_vs_r.pdf") trysavefig(outfile) end - if pp.animate_parallel_flow_vs_r_z && nr > 1 + if pp.animate_parallel_flow_vs_r_z && nr > 1 && nz > 1 # make a gif animation of ϕ(z) at different times anim = @animate for i ∈ itime_min:nwrite_movie:itime_max @views heatmap(r, z, parallel_flow[:,:,is,i], xlabel="r", ylabel="z", c = :deep, interpolation = :cubic) @@ -3219,7 +3222,7 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure, outfile = string(run_name, "_parallel_flow"*description*"_vs_r_z.gif") trygif(anim, outfile, fps=5) end - if pp.plot_parallel_flow_vs_r_z && nr > 1 + if pp.plot_parallel_flow_vs_r_z && nr > 1 && nz > 1 @views heatmap(r, z, parallel_flow[:,:,is,end], xlabel=L"r", ylabel=L"z", c = :deep, interpolation = :cubic, windowsize = (360,240), margin = 15pt) outfile = string(run_name, "_parallel_flow"*description*"_vs_r_z.pdf") @@ -3236,7 +3239,7 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure, # the parallel pressure parallel_pressuremin = minimum(parallel_pressure[:,:,is,:]) parallel_pressuremax = maximum(parallel_pressure) - if pp.plot_parallel_pressure_vs_r0_z # plot last timestep parallel_pressure[z,ir0] + if pp.plot_parallel_pressure_vs_r0_z && nz > 1 # plot last timestep parallel_pressure[z,ir0] @views plot(z, parallel_pressure[:,ir0,is,end], xlabel=L"z/L_z", ylabel=L"p_{i\|\|}") outfile = string(run_name, "_parallel_pressure"*description*"(r0,z)_vs_z.pdf") trysavefig(outfile) @@ -3246,7 +3249,7 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure, outfile = string(run_name, "_parallel_pressure"*description*"(r,z_wall)_vs_r.pdf") trysavefig(outfile) end - if pp.animate_parallel_pressure_vs_r_z && nr > 1 + if pp.animate_parallel_pressure_vs_r_z && nr > 1 && nz > 1 # make a gif animation of ϕ(z) at different times anim = @animate for i ∈ itime_min:nwrite_movie:itime_max @views heatmap(r, z, parallel_pressure[:,:,is,i], xlabel="r", ylabel="z", c = :deep, interpolation = :cubic) @@ -3254,7 +3257,7 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure, outfile = string(run_name, "_parallel_pressure"*description*"_vs_r_z.gif") trygif(anim, outfile, fps=5) end - if pp.plot_parallel_pressure_vs_r_z && nr > 1 + if pp.plot_parallel_pressure_vs_r_z && nr > 1 && nz > 1 @views heatmap(r, z, parallel_pressure[:,:,is,end], xlabel=L"r", ylabel=L"z", c = :deep, interpolation = :cubic, windowsize = (360,240), margin = 15pt) outfile = string(run_name, "_parallel_pressure"*description*"_vs_r_z.pdf") @@ -3264,17 +3267,17 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure, # Note factor of 2 here because currently temperatures are normalised by # Tref, while pressures are normalised by m*nref*c_ref^2=2*nref*Tref temperature = 2.0 * parallel_pressure ./ density - if pp.plot_parallel_temperature_vs_r0_z # plot last timestep parallel_temperature[z,ir0] + if pp.plot_parallel_temperature_vs_r0_z && nz > 1# plot last timestep parallel_temperature[z,ir0] @views plot(z, temperature[:,ir0,is,end], xlabel=L"z/L_z", ylabel=L"T_i") outfile = string(run_name, "_temperature"*description*"(r0,z)_vs_z.pdf") trysavefig(outfile) end - if pp.plot_wall_parallel_temperature_vs_r && nr > 1 # plot last timestep parallel_temperature[z_wall,r] + if pp.plot_wall_parallel_temperature_vs_r && nr > 1 && nz > 1 # plot last timestep parallel_temperature[z_wall,r] @views plot(r, temperature[end,:,is,end], xlabel=L"r/L_r", ylabel=L"T_i") outfile = string(run_name, "_temperature"*description*"(r,z_wall)_vs_r.pdf") trysavefig(outfile) end - if pp.animate_parallel_temperature_vs_r_z && nr > 1 + if pp.animate_parallel_temperature_vs_r_z && nr > 1 && nz > 1 # make a gif animation of T_||(z) at different times anim = @animate for i ∈ itime_min:nwrite_movie:itime_max @views heatmap(r, z, temperature[:,:,is,i], xlabel="r", ylabel="z", c = :deep, interpolation = :cubic) @@ -3282,7 +3285,7 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure, outfile = string(run_name, "_temperature"*description*"_vs_r_z.gif") trygif(anim, outfile, fps=5) end - if pp.plot_parallel_temperature_vs_r_z && nr > 1 + if pp.plot_parallel_temperature_vs_r_z && nr > 1 && nz > 1 @views heatmap(r, z, temperature[:,:,is,end], xlabel=L"r", ylabel=L"z", c = :deep, interpolation = :cubic, windowsize = (360,240), margin = 15pt) outfile = string(run_name, "_temperature"*description*"_vs_r_z.pdf") From e6bdb90023ea44ed8c5537d8a635448837c79f89 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Thu, 4 Jul 2024 16:07:30 +0000 Subject: [PATCH 12/34] Rename and add example file for slowing-down distribution test. Testing without self collisions is much faster. --- ...lowing-down-alphas-no-self-collisions.toml | 91 +++++++++++++++++++ ...lanck-source-sink-slowing-down-alphas.toml | 91 +++++++++++++++++++ 2 files changed, 182 insertions(+) create mode 100644 examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas-no-self-collisions.toml create mode 100644 examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas.toml diff --git a/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas-no-self-collisions.toml b/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas-no-self-collisions.toml new file mode 100644 index 000000000..e3070a6f7 --- /dev/null +++ b/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas-no-self-collisions.toml @@ -0,0 +1,91 @@ +# cheap input file for a 0D2V relaxation to a collisional Maxwellian distribution with self-ion collisions and collisions with fixed Maxwellian background of cold ions and electrons. +n_ion_species = 1 +n_neutral_species = 0 +electron_physics = "boltzmann_electron_response" +evolve_moments_density = false +evolve_moments_parallel_flow = false +evolve_moments_parallel_pressure = false +evolve_moments_conservation = false +T_e = 1.0 +T_wall = 1.0 +initial_density1 = 1.0 +initial_temperature1 = 1.0 +z_IC_option1 = "sinusoid" +z_IC_density_amplitude1 = 0.0 +z_IC_density_phase1 = 0.0 +z_IC_upar_amplitude1 = 0.0 +z_IC_upar_phase1 = 0.0 +z_IC_temperature_amplitude1 = 0.0 +z_IC_temperature_phase1 = 0.0 +vpa_IC_option1 = "isotropic-beam" +vpa_IC_v01 = 1.0 +vpa_IC_vth01 = 0.1 +#vpa_IC_option1 = "directed-beam" +#vpa_IC_vpa01 = -1.5 +#vpa_IC_vperp01 = 0.0 +charge_exchange_frequency = 0.0 +ionization_frequency = 0.0 +constant_ionization_rate = false + +z_ngrid = 1 +z_nelement = 1 +z_nelement_local = 1 +z_bc = "wall" +z_discretization = "chebyshev_pseudospectral" +r_ngrid = 1 +r_nelement = 1 +r_nelement_local = 1 +r_bc = "periodic" +r_discretization = "chebyshev_pseudospectral" +vpa_ngrid = 5 +vpa_nelement = 32 +vpa_L = 3.0 +vpa_bc = "zero" +vpa_discretization = "gausslegendre_pseudospectral" +vperp_ngrid = 5 +vperp_nelement = 16 +vperp_L = 1.5 +vperp_discretization = "gausslegendre_pseudospectral" +vperp_bc = "zero-no-regularity" +# Fokker-Planck operator requires the "gausslegendre_pseudospectral +# options for the vpa and vperp grids + +[fokker_planck_collisions] +# nuii sets the normalised input C[F,F] Fokker-Planck collision frequency +nuii = 1.876334222e-3 #(1/nu_alphae, as computed from input diagnostic) +Zi = 2.0 +self_collisions = false +slowing_down_test = true +frequency_option = "manual" +use_fokker_planck = true +sd_density = 1.0 +sd_temp = 0.0025 # TD/Ealpha +sd_mi = 0.5 # mD/malpha +sd_me = 0.000013616 # 0.25/1836.0 me/malpha +sd_q = 1.0 + +[ion_source] +active=true +source_strength=1.0 +source_T=0.005 +source_n=1.0 +r_profile="constant" +r_width=1.0 +r_relative_minimum=0.0 +z_profile="gaussian" +z_width=0.1 +z_relative_minimum=0.0 +source_v0 = 1.0 +#source_type="alphas" +source_type="alphas-with-losses" +#source_type="beam-with-losses" +#source_vpa0 = 1.0 +#source_vperp0 = 1.0 +sink_strength = 1.0 +sink_vth=0.07071067811865475 + +[timestepping] +nstep = 250000 +dt = 1.0e-4 +nwrite = 500 +nwrite_dfns = 500 diff --git a/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas.toml b/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas.toml new file mode 100644 index 000000000..841af3a1e --- /dev/null +++ b/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas.toml @@ -0,0 +1,91 @@ +# cheap input file for a 0D2V relaxation to a collisional Maxwellian distribution with self-ion collisions and collisions with fixed Maxwellian background of cold ions and electrons. +n_ion_species = 1 +n_neutral_species = 0 +electron_physics = "boltzmann_electron_response" +evolve_moments_density = false +evolve_moments_parallel_flow = false +evolve_moments_parallel_pressure = false +evolve_moments_conservation = false +T_e = 1.0 +T_wall = 1.0 +initial_density1 = 1.0 +initial_temperature1 = 1.0 +z_IC_option1 = "sinusoid" +z_IC_density_amplitude1 = 0.0 +z_IC_density_phase1 = 0.0 +z_IC_upar_amplitude1 = 0.0 +z_IC_upar_phase1 = 0.0 +z_IC_temperature_amplitude1 = 0.0 +z_IC_temperature_phase1 = 0.0 +vpa_IC_option1 = "isotropic-beam" +vpa_IC_v01 = 1.0 +vpa_IC_vth01 = 0.1 +#vpa_IC_option1 = "directed-beam" +#vpa_IC_vpa01 = -1.5 +#vpa_IC_vperp01 = 0.0 +charge_exchange_frequency = 0.0 +ionization_frequency = 0.0 +constant_ionization_rate = false + +z_ngrid = 1 +z_nelement = 1 +z_nelement_local = 1 +z_bc = "wall" +z_discretization = "chebyshev_pseudospectral" +r_ngrid = 1 +r_nelement = 1 +r_nelement_local = 1 +r_bc = "periodic" +r_discretization = "chebyshev_pseudospectral" +vpa_ngrid = 5 +vpa_nelement = 32 +vpa_L = 3.0 +vpa_bc = "zero" +vpa_discretization = "gausslegendre_pseudospectral" +vperp_ngrid = 5 +vperp_nelement = 16 +vperp_L = 1.5 +vperp_discretization = "gausslegendre_pseudospectral" +vperp_bc = "zero-no-regularity" +# Fokker-Planck operator requires the "gausslegendre_pseudospectral +# options for the vpa and vperp grids + +[fokker_planck_collisions] +# nuii sets the normalised input C[F,F] Fokker-Planck collision frequency +nuii = 1.876334222e-3 #(1/nu_alphae, as computed from input diagnostic) +Zi = 2.0 +self_collisions = true +slowing_down_test = true +frequency_option = "manual" +use_fokker_planck = true +sd_density = 1.0 +sd_temp = 0.0025 # TD/Ealpha +sd_mi = 0.5 # mD/malpha +sd_me = 0.000013616 # 0.25/1836.0 me/malpha +sd_q = 1.0 + +[ion_source] +active=true +source_strength=1.0 +source_T=0.005 +source_n=1.0 +r_profile="constant" +r_width=1.0 +r_relative_minimum=0.0 +z_profile="gaussian" +z_width=0.1 +z_relative_minimum=0.0 +source_v0 = 1.0 +#source_type="alphas" +source_type="alphas-with-losses" +#source_type="beam-with-losses" +#source_vpa0 = 1.0 +#source_vperp0 = 1.0 +sink_strength = 1.0 +sink_vth=0.07071067811865475 + +[timestepping] +nstep = 250000 +dt = 1.0e-4 +nwrite = 500 +nwrite_dfns = 500 From 8c84a1e8af5436678cba99843614e76a5965822d Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Fri, 5 Jul 2024 18:51:04 +0100 Subject: [PATCH 13/34] Change vperp numerical dissipation to use laplacian_derivative!(), support laplacian_derivative!() for chebyshev points with direct differentiation calculation. --- moment_kinetics/src/calculus.jl | 70 +++++++++++++++++++- moment_kinetics/src/numerical_dissipation.jl | 4 +- 2 files changed, 70 insertions(+), 4 deletions(-) diff --git a/moment_kinetics/src/calculus.jl b/moment_kinetics/src/calculus.jl index fbff04a2f..f4a73fb83 100644 --- a/moment_kinetics/src/calculus.jl +++ b/moment_kinetics/src/calculus.jl @@ -131,6 +131,70 @@ function second_derivative!(d2f, f, coord, spectral; handle_periodic=true) return nothing end +function laplacian_derivative!(d2f, f, coord, spectral; handle_periodic=true) + # computes (1/coord) d / coord ( coord d f / d(coord)) for vperp coordinate + # For spectral element methods, calculate second derivative by applying first + # derivative twice, with special treatment for element boundaries + + # First derivative + elementwise_derivative!(coord, f, spectral) + derivative_elements_to_full_grid!(coord.scratch3, coord.scratch_2d, coord) + if coord.name == "vperp" + # include the Jacobian + @. coord.scratch3 *= coord.grid + end + # MPI reconcile code here if used with z or r coords + + # Save elementwise first derivative result + coord.scratch2_2d .= coord.scratch_2d + + # Second derivative for element interiors + elementwise_derivative!(coord, coord.scratch3, spectral) + derivative_elements_to_full_grid!(d2f, coord.scratch_2d, coord) + if coord.name == "vperp" + # include the Jacobian + @. d2f /= coord.grid + end + # MPI reconcile code here if used with z or r coords + + # Add contribution to penalise discontinuous first derivatives at element + # boundaries. For smooth functions this would do nothing so should not affect + # convergence of the second derivative. Aims to stabilise numerical instability when + # spike develops at an element boundary. The coefficient is an arbitrary choice, it + # should probably be large enough for stability but as small as possible. + # + # Arbitrary numerical coefficient + C = 1.0 + function penalise_discontinuous_first_derivative!(d2f, imin, imax, df) + # Left element boundary + d2f[imin] += C * df[1] + + # Right element boundary + d2f[imax] -= C * df[end] + + return nothing + end + @views penalise_discontinuous_first_derivative!(d2f, 1, coord.imax[1], + coord.scratch2_2d[:,1]) + for ielement ∈ 2:coord.nelement_local + @views penalise_discontinuous_first_derivative!(d2f, coord.imin[ielement]-1, + coord.imax[ielement], + coord.scratch2_2d[:,ielement]) + end + + if coord.bc ∈ ("zero", "both_zero", "zero-no-regularity") + # For stability don't contribute to evolution at boundaries, in case these + # points are not set by a boundary condition. + # Full grid may be across processes and bc only applied to extreme ends of the + # domain. + if coord.irank == coord.nrank - 1 + d2f[end] = 0.0 + end + else + error("Unsupported bc '$(coord.bc)'") + end + return nothing +end """ mass_matrix_solve!(f, b, spectral::weak_discretization_info) @@ -176,9 +240,11 @@ function laplacian_derivative!(d2f, f, coord, spectral::weak_discretization_info # for all other coord.name, do exactly the same as second_derivative! above. mul!(coord.scratch, spectral.L_matrix, f) - if handle_periodic && coord.bc == "periodic" + if coord.bc == "periodic" && coord.name == "vperp" + error("laplacian_derivative!() cannot handle periodic boundaries for vperp") + elseif coord.bc == "periodic" if coord.nrank > 1 - error("second_derivative!() cannot handle periodic boundaries for a " + error("laplacian_derivative!() cannot handle periodic boundaries for a " * "distributed coordinate") end diff --git a/moment_kinetics/src/numerical_dissipation.jl b/moment_kinetics/src/numerical_dissipation.jl index 86e517229..f48a158c6 100644 --- a/moment_kinetics/src/numerical_dissipation.jl +++ b/moment_kinetics/src/numerical_dissipation.jl @@ -12,7 +12,7 @@ export setup_numerical_dissipation, vpa_boundary_buffer_decay!, using Base.Iterators: flatten using ..looping -using ..calculus: derivative!, second_derivative! +using ..calculus: derivative!, second_derivative!, laplacian_derivative! using ..derivatives: derivative_r!, derivative_z!, second_derivative_r!, second_derivative_z! using ..type_definitions: mk_float @@ -369,7 +369,7 @@ function vperp_dissipation!(f_out, f_in, vperp, spectral::T_spectral, dt, end @loop_s_r_z_vpa is ir iz ivpa begin - @views second_derivative!(vperp.scratch, f_in[ivpa,:,iz,ir,is], vperp, spectral) + @views laplacian_derivative!(vperp.scratch, f_in[ivpa,:,iz,ir,is], vperp, spectral) @views @. f_out[ivpa,:,iz,ir,is] += dt * diffusion_coefficient * vperp.scratch end From 4269c9ed3de61d5003458e2c84555d8311d28c91 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Mon, 8 Jul 2024 11:28:34 +0100 Subject: [PATCH 14/34] Addition of tests of laplacian_derivative!(). Note the large errors, especially associated with the Chebyshev pseudospectral method, which does not show good convergence even at large ngrid. The finite-element GaussLegendre method shows much better behaviour. --- moment_kinetics/test/calculus_tests.jl | 201 +++++++++++++++++++++++++ 1 file changed, 201 insertions(+) diff --git a/moment_kinetics/test/calculus_tests.jl b/moment_kinetics/test/calculus_tests.jl index e439ac4cf..ef13b75bb 100644 --- a/moment_kinetics/test/calculus_tests.jl +++ b/moment_kinetics/test/calculus_tests.jl @@ -5,6 +5,7 @@ include("setup.jl") using moment_kinetics.input_structs: grid_input, advection_input using moment_kinetics.coordinates: define_coordinate using moment_kinetics.calculus: derivative!, second_derivative!, integral +using moment_kinetics.calculus: laplacian_derivative! using MPI using Random @@ -1355,6 +1356,101 @@ function runtests() end end + @testset "Chebyshev pseudospectral cylindrical laplacian derivatives (4 argument), zero" verbose=false begin + @testset "$nelement $ngrid" for (nelement, ngrid, rtol) ∈ + ( + (4, 7, 2.e-1), + (4, 8, 2.e-1), + (4, 9, 4.e-2), + (4, 10, 4.e-2), + (4, 11, 5.e-3), + (4, 12, 5.e-3), + (4, 13, 5.e-3), + (4, 14, 5.e-3), + (4, 15, 5.e-3), + (4, 16, 5.e-3), + (4, 17, 5.e-3), + (4, 18, 5.e-3), + (4, 19, 5.e-3), + (4, 20, 5.e-3), + (4, 21, 5.e-3), + (4, 22, 5.e-3), + (4, 23, 5.e-3), + (4, 24, 4.e-3), + (4, 25, 4.e-3), + (4, 26, 4.e-3), + (4, 27, 4.e-3), + (4, 28, 4.e-3), + (4, 29, 4.e-3), + (4, 30, 4.e-3), + (4, 31, 4.e-3), + (4, 32, 4.e-3), + (4, 33, 4.e-3), + + (5, 7, 2.e-1), + (5, 8, 8.e-2), + (5, 9, 5.e-2), + (5, 10, 8.e-3), + (5, 11, 8.e-3), + (5, 12, 8.e-3), + (5, 13, 8.e-3), + (5, 14, 8.e-3), + (5, 15, 8.e-3), + (5, 16, 2.e-3), + (5, 17, 2.e-3), + (5, 18, 2.e-3), + (5, 19, 2.e-3), + (5, 20, 2.e-3), + (5, 21, 2.e-3), + (5, 22, 2.e-3), + (5, 23, 2.e-3), + (5, 24, 2.e-3), + (5, 25, 2.e-3), + (5, 26, 2.e-3), + (5, 27, 2.e-3), + (5, 28, 2.e-3), + (5, 29, 2.e-3), + (5, 30, 2.e-3), + (5, 31, 2.e-3), + (5, 32, 2.e-3), + (5, 33, 2.e-3), + ), cheb_option in ("FFT","matrix") + + # define inputs needed for the test + L = 6.0 + bc = "zero" + # fd_option and adv_input not actually used so given values unimportant + fd_option = "" + adv_input = advection_input("default", 1.0, 0.0, 0.0) + # create the 'input' struct containing input info needed to create a + # coordinate + nelement_local = nelement + nrank_per_block = 0 # dummy value + irank = 0 # dummy value + comm = MPI.COMM_NULL # dummy value + element_spacing_option = "uniform" + input = grid_input("vperp", ngrid, nelement, + nelement_local, nrank_per_block, irank, L, + "chebyshev_pseudospectral", fd_option, cheb_option, bc, adv_input, comm, + element_spacing_option) + # create the coordinate struct 'x' and info for derivatives, etc. + # This test runs effectively in serial, so use `ignore_MPI=true` to avoid + # errors due to communicators not being fully set up. + x, spectral = define_coordinate(input; ignore_MPI=true) + + f = @. exp(-x.grid^2) + expected_d2f = @. 4.0*(x.grid^2 - 1.0)*exp(-x.grid^2) + # create array for the derivative d2f/dx2 + d2f = similar(f) + + # differentiate f + laplacian_derivative!(d2f, f, x, spectral) + + @test isapprox(d2f, expected_d2f, rtol=rtol, atol=1.e-10, + norm=maxabs_norm) + end + end + @testset "GaussLegendre pseudospectral second derivatives (4 argument), periodic" verbose=false begin @testset "$nelement $ngrid" for (nelement, ngrid, rtol) ∈ ( @@ -1460,6 +1556,111 @@ function runtests() # differentiate f second_derivative!(d2f, f, x, spectral) + @test isapprox(d2f, expected_d2f, rtol=rtol, atol=1.e-10, + norm=maxabs_norm) + end + end + + @testset "GaussLegendre pseudospectral cylindrical laplacian derivatives (4 argument), zero" verbose=false begin + @testset "$nelement $ngrid" for (nelement, ngrid, rtol) ∈ + ( + (1, 8, 5.e-2), + (1, 9, 1.e-1), + (1, 10, 2.e-1), + (1, 11, 5.e-2), + (1, 12, 5.e-2), + (1, 13, 5.e-2), + (1, 14, 5.e-2), + (1, 15, 5.e-3), + (1, 16, 5.e-2), + (1, 17, 5.e-3), + + (2, 6, 8.e-2), + (2, 7, 8.e-2), + (2, 8, 5.e-2), + (2, 9, 5.e-2), + (2, 10, 5.e-2), + (2, 11, 5.e-3), + (2, 12, 5.e-3), + (2, 13, 5.e-4), + (2, 14, 5.e-4), + (2, 15, 5.e-4), + (2, 16, 5.e-4), + (2, 17, 5.e-4), + + (3, 6, 5.e-2), + (3, 7, 5.e-3), + (3, 8, 5.e-2), + (3, 9, 5.e-4), + (3, 10, 5.e-3), + (3, 11, 5.e-4), + (3, 12, 5.e-5), + (3, 13, 5.e-5), + (3, 14, 5.e-6), + (3, 15, 5.e-6), + (3, 16, 5.e-6), + (3, 17, 5.e-8), + + (4, 5, 5.e-2), + (4, 6, 2.e-2), + (4, 7, 2.e-2), + (4, 8, 2.e-3), + (4, 9, 1.e-3), + (4, 10, 1.e-4), + (4, 11, 8.e-5), + (4, 12, 8.e-6), + (4, 13, 8.e-6), + (4, 14, 8.e-7), + (4, 15, 8.e-7), + (4, 16, 8.e-8), + (4, 17, 8.e-8), + + (5, 5, 4.e-2), + (5, 6, 8.e-3), + (5, 7, 5.e-3), + (5, 8, 5.e-4), + (5, 9, 8.e-5), + (5, 10, 5.e-6), + (5, 11, 8.e-6), + (5, 12, 4.e-7), + (5, 13, 2.e-7), + (5, 14, 2.e-7), + (5, 15, 8.e-7), + (5, 16, 8.e-10), + (5, 17, 8.e-10), + ) + + # define inputs needed for the test + L = 6.0 + bc = "zero" + # fd_option and adv_input not actually used so given values unimportant + fd_option = "" + adv_input = advection_input("default", 1.0, 0.0, 0.0) + cheb_option = "" + # create the 'input' struct containing input info needed to create a + # coordinate + nelement_local = nelement + nrank_per_block = 0 # dummy value + irank = 0 # dummy value + comm = MPI.COMM_NULL # dummy value + element_spacing_option = "uniform" + input = grid_input("vperp", ngrid, nelement, + nelement_local, nrank_per_block, irank, L, + "gausslegendre_pseudospectral", fd_option, cheb_option, bc, adv_input, comm, + element_spacing_option) + # create the coordinate struct 'x' and info for derivatives, etc. + # This test runs effectively in serial, so use `ignore_MPI=true` to avoid + # errors due to communicators not being fully set up. + x, spectral = define_coordinate(input; ignore_MPI=true, collision_operator_dim=false) + + f = @. exp(-x.grid^2) + expected_d2f = @. 4.0*(x.grid^2 - 1.0)*exp(-x.grid^2) + # create array for the derivative d2f/dx2 + d2f = similar(f) + + # differentiate f + laplacian_derivative!(d2f, f, x, spectral) + @test isapprox(d2f, expected_d2f, rtol=rtol, atol=1.e-10, norm=maxabs_norm) end From 99e9525a5d388b1d06056ea6bd2ee03eaa03bc42 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Mon, 8 Jul 2024 12:31:25 +0100 Subject: [PATCH 15/34] Initial refactor moving Er_constant to geo struct. --- moment_kinetics/ext/manufactured_solns_ext.jl | 19 +++++++++++------ moment_kinetics/src/em_fields.jl | 9 ++++---- moment_kinetics/src/geo.jl | 21 ++++++++++++++----- moment_kinetics/src/input_structs.jl | 6 +++++- moment_kinetics/src/moment_kinetics_input.jl | 6 +++--- moment_kinetics/src/time_advance.jl | 8 +++---- 6 files changed, 46 insertions(+), 23 deletions(-) diff --git a/moment_kinetics/ext/manufactured_solns_ext.jl b/moment_kinetics/ext/manufactured_solns_ext.jl index de6d4e44b..b2bb21290 100644 --- a/moment_kinetics/ext/manufactured_solns_ext.jl +++ b/moment_kinetics/ext/manufactured_solns_ext.jl @@ -43,6 +43,8 @@ using IfElse # does not appear in a particular geometric coefficient, because # that coefficient is a constant. struct geometric_coefficients_sym{} + Er_constant::mk_float + Ez_constant::mk_float rhostar::mk_float Bzed::Union{mk_float,Num} Bzeta::Union{mk_float,Num} @@ -62,6 +64,8 @@ using IfElse option = geometry_input_data.option rhostar = geometry_input_data.rhostar pitch = geometry_input_data.pitch + Er_constant = geometry_input_data.Er_constant + Ez_constant = geometry_input_data.Ez_constant if option == "constant-helical" || option == "default" bzed = pitch bzeta = sqrt(1 - bzed^2) @@ -93,7 +97,8 @@ using IfElse else input_option_error("$option", option) end - geo_sym = geometric_coefficients_sym(rhostar,Bzed,Bzeta,Bmag,bzed,bzeta,dBdz,dBdr,jacobian) + geo_sym = geometric_coefficients_sym(Er_constant,Ez_constant, + rhostar,Bzed,Bzeta,Bmag,bzed,bzeta,dBdz,dBdr,jacobian) return geo_sym end @@ -272,7 +277,7 @@ using IfElse upari = 0.0 elseif z_bc == "wall" densi = densi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) - Er, Ez, phi = electric_fields(Lr,Lz,r_bc,z_bc,composition,nr,manufactured_solns_input,species) + Er, Ez, phi = electric_fields(Lr,Lz,r_bc,z_bc,composition,geometry,nr,manufactured_solns_input,species) Bzeta = geometry.Bzeta Bmag = geometry.Bmag rhostar = geometry.rhostar @@ -402,7 +407,7 @@ using IfElse return dfni end - function electric_fields(Lr, Lz, r_bc, z_bc, composition, nr, + function electric_fields(Lr, Lz, r_bc, z_bc, composition, geometry, nr, manufactured_solns_input, species) # define derivative operators @@ -500,11 +505,13 @@ using IfElse return manufactured_solns_list end - function manufactured_electric_fields(Lr, Lz, r_bc, z_bc, composition, nr, + function manufactured_electric_fields(Lr, Lz, r_bc, z_bc, composition, geometry_input_data::geometry_input, nr, manufactured_solns_input, species) + # calculate the geometry symbolically + geometry = geometry_sym(geometry_input_data,Lz,Lr,nr) # calculate the electric fields and the potential - Er, Ez, phi = electric_fields(Lr, Lz, r_bc, z_bc, composition, nr, + Er, Ez, phi = electric_fields(Lr, Lz, r_bc, z_bc, composition, geometry, nr, manufactured_solns_input, species) Er_func = build_function(Er, z, r, t, expression=Val{false}) @@ -603,7 +610,7 @@ using IfElse # calculate the electric fields and the potential Er, Ez, phi = electric_fields(r_coord.L, z_coord.L, r_coord.bc, z_coord.bc, - composition, r_coord.n, manufactured_solns_input, + composition, geometry, r_coord.n, manufactured_solns_input, ion_species) # the adiabatic invariant (for compactness) diff --git a/moment_kinetics/src/em_fields.jl b/moment_kinetics/src/em_fields.jl index 6591e9062..5c021099b 100644 --- a/moment_kinetics/src/em_fields.jl +++ b/moment_kinetics/src/em_fields.jl @@ -31,7 +31,7 @@ end """ update_phi updates the electrostatic potential, phi """ -function update_phi!(fields, fvec, vperp, z, r, composition, z_spectral, r_spectral, scratch_dummy, gyroavs::gyro_operators) +function update_phi!(fields, fvec, vperp, z, r, composition, geometry, z_spectral, r_spectral, scratch_dummy, gyroavs::gyro_operators) n_ion_species = composition.n_ion_species @boundscheck size(fields.phi,1) == z.n || throw(BoundsError(fields.phi)) @boundscheck size(fields.phi,2) == r.n || throw(BoundsError(fields.phi)) @@ -125,8 +125,8 @@ function update_phi!(fields, fvec, vperp, z, r, composition, z_spectral, r_spect end else @loop_r_z ir iz begin - fields.Er[iz,ir] = composition.Er_constant - # Er_constant defaults to 0.0 in moment_kinetics_input.jl + fields.Er[iz,ir] = geometry.input.Er_constant + # Er_constant defaults to 0.0 in geo.jl end end #Ez = - d phi / dz @@ -137,7 +137,8 @@ function update_phi!(fields, fvec, vperp, z, r, composition, z_spectral, r_spect z_spectral,z) else @serial_region begin - fields.Ez[:,:] .= 0.0 + fields.Ez[:,:] .= = geometry.input.Ez_constant + # Ez_constant defaults to 0.0 in geo.jl end end diff --git a/moment_kinetics/src/geo.jl b/moment_kinetics/src/geo.jl index 9ad6d4885..e9c188bb8 100644 --- a/moment_kinetics/src/geo.jl +++ b/moment_kinetics/src/geo.jl @@ -73,11 +73,22 @@ option = "" """ function setup_geometry_input(toml_input::Dict, reference_rhostar) - input_section = get(toml_input, "geometry", Dict{String,Any}()) - if !("rhostar" ∈ keys(input_section)) - # Set default rhostar with reference value - input_section["rhostar"] = get(input_section, "rhostar", reference_rhostar) - end + # read the input toml and specify a sensible default + input_section = set_defaults_and_check_section!(toml_input, "geometry", + # begin default inputs (as kwargs) + # rhostar ion (ref) + rhostar = reference_rhostar #used to premultiply ExB drift terms + # magnetic geometry option + option = "constant-helical" # "1D-mirror" + # pitch ( = Bzed/Bmag if geometry_option == "constant-helical") + pitch = 1.0, + # DeltaB ( = (Bzed(z=L/2) - Bzed(0))/Bref if geometry_option == "1D-mirror") + DeltaB = 0.0, + # constant for testing nonzero Er when nr = 1 + Er_constant = 0.0, + # constant for testing nonzero Ez when nz = 1 + Ez_constant = 0.0) + input = Dict(Symbol(k)=>v for (k,v) in input_section) #println(input) return geometry_input(; input...) diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl index 1c3f18f3f..634b2fadd 100644 --- a/moment_kinetics/src/input_structs.jl +++ b/moment_kinetics/src/input_structs.jl @@ -305,7 +305,7 @@ mutable struct species_composition # wall potential used if electron_physics=boltzmann_electron_response_with_simple_sheath phi_wall::mk_float # constant for testing nonzero Er - Er_constant::mk_float + #Er_constant::mk_float # ratio of the neutral particle mass to the ion mass mn_over_mi::mk_float # ratio of the electron particle mass to the ion mass @@ -414,6 +414,10 @@ Base.@kwdef struct geometry_input pitch::mk_float = 1.0 # DeltaB ( = (Bzed(z=L/2) - Bzed(0))/Bref if geometry_option == "1D-mirror") DeltaB::mk_float = 0.0 + # constant for testing nonzero Er when nr = 1 + Er_constant::mk_float + # constant for testing nonzero Ez when nz = 1 + Ez_constant::mk_float end @enum binary_format_type hdf5 netcdf diff --git a/moment_kinetics/src/moment_kinetics_input.jl b/moment_kinetics/src/moment_kinetics_input.jl index cddadfec0..393d5e55d 100644 --- a/moment_kinetics/src/moment_kinetics_input.jl +++ b/moment_kinetics/src/moment_kinetics_input.jl @@ -103,7 +103,7 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) # if false use true Knudsen cosine for neutral wall bc composition.use_test_neutral_wall_pdf = get(scan_input, "use_test_neutral_wall_pdf", false) # constant to be used to test nonzero Er in wall boundary condition - composition.Er_constant = get(scan_input, "Er_constant", 0.0) + #composition.Er_constant = get(scan_input, "Er_constant", 0.0) # The ion flux reaching the wall that is recycled as neutrals is reduced by # `recycling_fraction` to account for ions absorbed by the wall. composition.recycling_fraction = get(scan_input, "recycling_fraction", 1.0) @@ -1015,7 +1015,7 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) # wall potential at z = 0 phi_wall = 0.0 # constant to test nonzero Er - Er_constant = 0.0 + #Er_constant = 0.0 # ratio of the neutral particle mass to the ion particle mass mn_over_mi = 1.0 # ratio of the electron particle mass to the ion particle mass @@ -1025,7 +1025,7 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) recycling_fraction = 1.0 gyrokinetic_ions = false composition = species_composition(n_species, n_ion_species, n_neutral_species, - electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, Er_constant, + electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, #Er_constant, mn_over_mi, me_over_mi, recycling_fraction, gyrokinetic_ions, allocate_float(n_species)) species_ion = Array{species_parameters_mutable,1}(undef,n_ion_species) diff --git a/moment_kinetics/src/time_advance.jl b/moment_kinetics/src/time_advance.jl index 27eec4a0f..0fbc1d1ac 100644 --- a/moment_kinetics/src/time_advance.jl +++ b/moment_kinetics/src/time_advance.jl @@ -394,7 +394,7 @@ function setup_time_advance!(pdf, fields, vz, vr, vzeta, vpa, vperp, z, r, gyrop # initialize the electrostatic potential begin_serial_region() - update_phi!(fields, scratch[1], vperp, z, r, composition, z_spectral, r_spectral, scratch_dummy, gyroavs) + update_phi!(fields, scratch[1], vperp, z, r, composition, geometry, z_spectral, r_spectral, scratch_dummy, gyroavs) @serial_region begin # save the initial phi(z) for possible use later (e.g., if forcing phi) fields.phi0 .= fields.phi @@ -592,7 +592,7 @@ function setup_time_advance!(pdf, fields, vz, vr, vzeta, vpa, vperp, z, r, gyrop end end - update_phi!(fields, scratch[1], vperp, z, r, composition, z_spectral, r_spectral, + update_phi!(fields, scratch[1], vperp, z, r, composition, geometry, z_spectral, r_spectral, scratch_dummy, gyroavs) calculate_ion_moment_derivatives!(moments, scratch[1], scratch_dummy, z, z_spectral, ion_mom_diss_coeff) @@ -1572,7 +1572,7 @@ function rk_update!(scratch, pdf, moments, fields, boundary_distributions, vz, v end # update the electrostatic potential phi - update_phi!(fields, scratch[istage+1], vperp, z, r, composition, z_spectral, + update_phi!(fields, scratch[istage+1], vperp, z, r, composition, geometry, z_spectral, r_spectral, scratch_dummy, gyroavs) # _block_synchronize() here because phi needs to be read on different ranks than # it was written on, even though the loop-type does not change here. However, @@ -1603,7 +1603,7 @@ function rk_update!(scratch, pdf, moments, fields, boundary_distributions, vz, v end # update the electrostatic potential phi - update_phi!(fields, scratch[istage+1], vperp, z, r, composition, z_spectral, + update_phi!(fields, scratch[istage+1], vperp, z, r, composition, geometry, z_spectral, r_spectral, scratch_dummy, gyroavs) if !(( moments.evolve_upar || moments.evolve_ppar) && istage == length(scratch)-1) From 4cde46c91d85bedf618e2fe47d6aba823b789009 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Mon, 8 Jul 2024 12:35:21 +0100 Subject: [PATCH 16/34] Remove Er_constant from comments. --- .../fokker-planck-1D2V-even_nz-shorttest-nstep200.toml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/fokker-planck-1D2V/fokker-planck-1D2V-even_nz-shorttest-nstep200.toml b/examples/fokker-planck-1D2V/fokker-planck-1D2V-even_nz-shorttest-nstep200.toml index eb7d60b18..7422ef0b0 100644 --- a/examples/fokker-planck-1D2V/fokker-planck-1D2V-even_nz-shorttest-nstep200.toml +++ b/examples/fokker-planck-1D2V/fokker-planck-1D2V-even_nz-shorttest-nstep200.toml @@ -9,7 +9,6 @@ evolve_moments_parallel_flow = false evolve_moments_parallel_pressure = false evolve_moments_conservation = false #force_Er_zero_at_wall = false #true -#Er_constant = 0.0 #epsilon_offset = 0.1 #use_vpabar_in_mms_dfni = true T_e = 1.0 @@ -95,4 +94,4 @@ vperp_discretization = "gausslegendre_pseudospectral" [fokker_planck_collisions] use_fokker_planck = true nuii = 1.0 -frequency_option = "manual" \ No newline at end of file +frequency_option = "manual" From 705f4a3a0c66a646c62e86558188815722faf7eb Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Mon, 8 Jul 2024 12:37:03 +0100 Subject: [PATCH 17/34] Delete out of date input files. --- ..._new_nel_r_1_z_16_vpa_16_vperp_8_diss.toml | 98 ------------------ ...MS_new_nel_r_1_z_6_vpa_6_vperp_3_diss.toml | 99 ------------------- 2 files changed, 197 deletions(-) delete mode 100644 runs/1D-mirror_MMS_new_nel_r_1_z_16_vpa_16_vperp_8_diss.toml delete mode 100644 runs/1D-mirror_MMS_new_nel_r_1_z_6_vpa_6_vperp_3_diss.toml diff --git a/runs/1D-mirror_MMS_new_nel_r_1_z_16_vpa_16_vperp_8_diss.toml b/runs/1D-mirror_MMS_new_nel_r_1_z_16_vpa_16_vperp_8_diss.toml deleted file mode 100644 index 9a475346b..000000000 --- a/runs/1D-mirror_MMS_new_nel_r_1_z_16_vpa_16_vperp_8_diss.toml +++ /dev/null @@ -1,98 +0,0 @@ -n_ion_species = 1 -n_neutral_species = 0 -electron_physics = "boltzmann_electron_response" -#electron_physics = "boltzmann_electron_response_with_simple_sheath" -evolve_moments_density = false -evolve_moments_parallel_flow = false -evolve_moments_parallel_pressure = false -evolve_moments_conservation = false -force_Er_zero_at_wall = false #true -geometry_option="1D-mirror" -DeltaB=10.0 -#geometry_option="constant-helical" -#pitch=1.0 -Er_constant = 0.0 -T_e = 1.0 -T_wall = 1.0 -rhostar = 1.0 -initial_density1 = 0.5 -initial_temperature1 = 1.0 -initial_density2 = 0.5 -initial_temperature2 = 1.0 -z_IC_option1 = "sinusoid" -z_IC_density_amplitude1 = 0.001 -z_IC_density_phase1 = 0.0 -z_IC_upar_amplitude1 = 0.0 -z_IC_upar_phase1 = 0.0 -z_IC_temperature_amplitude1 = 0.0 -z_IC_temperature_phase1 = 0.0 -z_IC_option2 = "sinusoid" -z_IC_density_amplitude2 = 0.001 -z_IC_density_phase2 = 0.0 -z_IC_upar_amplitude2 = 0.0 -z_IC_upar_phase2 = 0.0 -z_IC_temperature_amplitude2 = 0.0 -z_IC_temperature_phase2 = 0.0 -charge_exchange_frequency = 0.0 -ionization_frequency = 0.0 -nstep = 2000 -dt = 0.0005 -nwrite = 200 -nwrite_dfns = 200 -use_semi_lagrange = false -n_rk_stages = 4 -split_operators = false -z_ngrid = 17 -z_nelement = 16 -z_nelement_local = 16 -z_bc = "wall" -z_discretization = "chebyshev_pseudospectral" -r_ngrid = 1 -r_nelement = 1 -r_nelement_local = 1 -r_bc = "periodic" -r_discretization = "chebyshev_pseudospectral" -vpa_ngrid = 17 -vpa_nelement = 16 -vpa_L = 12.0 -vpa_bc = "zero" -vpa_discretization = "chebyshev_pseudospectral" -vperp_ngrid = 17 -vperp_nelement = 8 -vperp_L = 6.0 -#vperp_discretization = "finite_difference" -vperp_discretization = "chebyshev_pseudospectral" -vperp_bc = "zero" - -vz_ngrid = 17 -vz_nelement = 4 -vz_L = 12.0 -vz_bc = "periodic" -vz_discretization = "chebyshev_pseudospectral" - -vr_ngrid = 17 -vr_nelement = 4 -vr_L = 12.0 -vr_bc = "periodic" -vr_discretization = "chebyshev_pseudospectral" - -vzeta_ngrid = 17 -vzeta_nelement = 4 -vzeta_L = 12.0 -vzeta_bc = "periodic" -vzeta_discretization = "chebyshev_pseudospectral" - -[manufactured_solns] - use_for_advance=true - use_for_init=true - # constant to be used to control Ez divergence in MMS tests - epsilon_offset=0.1 - # bool to control if dfni is a function of vpa or vpabar in MMS test - use_vpabar_in_mms_dfni=true - alpha_switch=1.0 - type="default" -[ion_numerical_dissipation] -vpa_dissipation_coefficient = 0.001 -vperp_dissipation_coefficient = 0.001 -#z_dissipation_coefficient = 0.1 -r_dissipation_coefficient = 0.0 diff --git a/runs/1D-mirror_MMS_new_nel_r_1_z_6_vpa_6_vperp_3_diss.toml b/runs/1D-mirror_MMS_new_nel_r_1_z_6_vpa_6_vperp_3_diss.toml deleted file mode 100644 index 347dd2b04..000000000 --- a/runs/1D-mirror_MMS_new_nel_r_1_z_6_vpa_6_vperp_3_diss.toml +++ /dev/null @@ -1,99 +0,0 @@ -n_ion_species = 1 -n_neutral_species = 0 -electron_physics = "boltzmann_electron_response" -#electron_physics = "boltzmann_electron_response_with_simple_sheath" -evolve_moments_density = false -evolve_moments_parallel_flow = false -evolve_moments_parallel_pressure = false -evolve_moments_conservation = false -force_Er_zero_at_wall = false #true -geometry_option="1D-mirror" -DeltaB=0.5 -#geometry_option="constant-helical" -#pitch=1.0 -Er_constant = 0.0 -T_e = 1.0 -T_wall = 1.0 -rhostar = 1.0 -initial_density1 = 0.5 -initial_temperature1 = 1.0 -initial_density2 = 0.5 -initial_temperature2 = 1.0 -z_IC_option1 = "sinusoid" -z_IC_density_amplitude1 = 0.001 -z_IC_density_phase1 = 0.0 -z_IC_upar_amplitude1 = 0.0 -z_IC_upar_phase1 = 0.0 -z_IC_temperature_amplitude1 = 0.0 -z_IC_temperature_phase1 = 0.0 -z_IC_option2 = "sinusoid" -z_IC_density_amplitude2 = 0.001 -z_IC_density_phase2 = 0.0 -z_IC_upar_amplitude2 = 0.0 -z_IC_upar_phase2 = 0.0 -z_IC_temperature_amplitude2 = 0.0 -z_IC_temperature_phase2 = 0.0 -charge_exchange_frequency = 0.0 -ionization_frequency = 0.0 -nstep = 2000 -dt = 0.0005 -nwrite = 200 -nwrite_dfns = 200 -use_semi_lagrange = false -n_rk_stages = 4 -split_operators = false -z_ngrid = 17 -z_nelement = 6 -z_nelement_local = 6 -z_bc = "wall" -z_discretization = "chebyshev_pseudospectral" -r_ngrid = 1 -r_nelement = 1 -r_nelement_local = 1 -r_bc = "periodic" -r_discretization = "chebyshev_pseudospectral" -vpa_ngrid = 17 -vpa_nelement = 6 -vpa_L = 12.0 -vpa_bc = "zero" -vpa_discretization = "chebyshev_pseudospectral" -vperp_ngrid = 17 -vperp_nelement = 3 -vperp_L = 6.0 -#vperp_discretization = "finite_difference" -#vperp_discretization = "chebyshev_pseudospectral" -vperp_discretization = "gausslegendre_pseudospectral" -vperp_bc = "zero" - -vz_ngrid = 17 -vz_nelement = 4 -vz_L = 12.0 -vz_bc = "periodic" -vz_discretization = "chebyshev_pseudospectral" - -vr_ngrid = 17 -vr_nelement = 4 -vr_L = 12.0 -vr_bc = "periodic" -vr_discretization = "chebyshev_pseudospectral" - -vzeta_ngrid = 17 -vzeta_nelement = 4 -vzeta_L = 12.0 -vzeta_bc = "periodic" -vzeta_discretization = "chebyshev_pseudospectral" - -[manufactured_solns] - use_for_advance=true - use_for_init=true - # constant to be used to control Ez divergence in MMS tests - epsilon_offset=0.1 - # bool to control if dfni is a function of vpa or vpabar in MMS test - use_vpabar_in_mms_dfni=true - alpha_switch=1.0 - type="default" -[ion_numerical_dissipation] -vpa_dissipation_coefficient = 0.001 -vperp_dissipation_coefficient = 0.001 -#z_dissipation_coefficient = 0.1 -r_dissipation_coefficient = 0.0 From d1d091f266e3abc98ac589b6398b9e4cd8a891a8 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Mon, 8 Jul 2024 13:06:21 +0100 Subject: [PATCH 18/34] Move Er_constant to geometry namelist. --- ...1D-mirror_MMS_ngrid_9_nel_r_1_z_16_vpa_16_vperp_8_diss.toml | 2 +- ...D-mirror_MMS_ngrid_9_nel_r_1_z_32_vpa_32_vperp_16_diss.toml | 2 +- runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_4_vpa_4_vperp_2_diss.toml | 2 +- ...D-mirror_MMS_ngrid_9_nel_r_1_z_64_vpa_64_vperp_32_diss.toml | 2 +- runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_8_vpa_8_vperp_4_diss.toml | 2 +- runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_diss.toml | 3 ++- runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml | 3 ++- runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml | 3 ++- ...D-mirror_MMS_ngrid_5_nel_r_16_z_16_vpa_16_vperp_8_diss.toml | 2 +- ...-mirror_MMS_ngrid_5_nel_r_32_z_32_vpa_32_vperp_16_diss.toml | 2 +- runs/2D-mirror_MMS_ngrid_5_nel_r_4_z_4_vpa_4_vperp_2_diss.toml | 2 +- runs/2D-mirror_MMS_ngrid_5_nel_r_8_z_8_vpa_8_vperp_4_diss.toml | 2 +- runs/2D-wall_MMS_nel_r_32_z_32_vpa_16_vperp_1_diss.toml | 2 +- 13 files changed, 16 insertions(+), 13 deletions(-) diff --git a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_16_vpa_16_vperp_8_diss.toml b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_16_vpa_16_vperp_8_diss.toml index 883a73275..f22e52894 100644 --- a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_16_vpa_16_vperp_8_diss.toml +++ b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_16_vpa_16_vperp_8_diss.toml @@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false evolve_moments_parallel_pressure = false evolve_moments_conservation = false force_Er_zero_at_wall = false #true -Er_constant = 0.0 T_e = 1.0 T_wall = 1.0 initial_density1 = 0.5 @@ -99,3 +98,4 @@ DeltaB=0.5 #option="constant-helical" pitch=1.0 rhostar = 1.0 +Er_constant = 0.0 diff --git a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_32_vpa_32_vperp_16_diss.toml b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_32_vpa_32_vperp_16_diss.toml index 50e878cb3..fefea7811 100644 --- a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_32_vpa_32_vperp_16_diss.toml +++ b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_32_vpa_32_vperp_16_diss.toml @@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false evolve_moments_parallel_pressure = false evolve_moments_conservation = false force_Er_zero_at_wall = false #true -Er_constant = 0.0 T_e = 1.0 T_wall = 1.0 initial_density1 = 0.5 @@ -99,3 +98,4 @@ DeltaB=0.5 #option="constant-helical" pitch=1.0 rhostar = 1.0 +Er_constant = 0.0 diff --git a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_4_vpa_4_vperp_2_diss.toml b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_4_vpa_4_vperp_2_diss.toml index 1624ca129..58051541c 100644 --- a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_4_vpa_4_vperp_2_diss.toml +++ b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_4_vpa_4_vperp_2_diss.toml @@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false evolve_moments_parallel_pressure = false evolve_moments_conservation = false force_Er_zero_at_wall = false #true -Er_constant = 0.0 T_e = 1.0 T_wall = 1.0 initial_density1 = 0.5 @@ -99,3 +98,4 @@ DeltaB=0.5 #option="constant-helical" pitch=1.0 rhostar = 1.0 +Er_constant = 0.0 diff --git a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_64_vpa_64_vperp_32_diss.toml b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_64_vpa_64_vperp_32_diss.toml index 74e3af388..0cc6b6ef1 100644 --- a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_64_vpa_64_vperp_32_diss.toml +++ b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_64_vpa_64_vperp_32_diss.toml @@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false evolve_moments_parallel_pressure = false evolve_moments_conservation = false force_Er_zero_at_wall = false #true -Er_constant = 0.0 T_e = 1.0 T_wall = 1.0 initial_density1 = 0.5 @@ -99,3 +98,4 @@ DeltaB=0.5 #option="constant-helical" pitch=1.0 rhostar = 1.0 +Er_constant = 0.0 diff --git a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_8_vpa_8_vperp_4_diss.toml b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_8_vpa_8_vperp_4_diss.toml index f4d4956ab..0f3a24c74 100644 --- a/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_8_vpa_8_vperp_4_diss.toml +++ b/runs/1D-mirror_MMS_ngrid_9_nel_r_1_z_8_vpa_8_vperp_4_diss.toml @@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false evolve_moments_parallel_pressure = false evolve_moments_conservation = false force_Er_zero_at_wall = false #true -Er_constant = 0.0 T_e = 1.0 T_wall = 1.0 initial_density1 = 0.5 @@ -99,3 +98,4 @@ DeltaB=0.5 #option="constant-helical" pitch=1.0 rhostar = 1.0 +Er_constant = 0.0 diff --git a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_diss.toml b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_diss.toml index 8bf508ee1..96611183c 100644 --- a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_diss.toml +++ b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_diss.toml @@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false evolve_moments_parallel_pressure = false evolve_moments_conservation = false force_Er_zero_at_wall = false #true -Er_constant = 0.0 T_e = 1.0 T_wall = 1.0 initial_density1 = 0.5 @@ -89,3 +88,5 @@ split_operators = false vpa_dissipation_coefficient = 0.01 #z_dissipation_coefficient = 0.1 #r_dissipation_coefficient = 0.0 +[geometry] +Er_constant = 0.0 diff --git a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml index 1abe96001..346996f69 100644 --- a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml +++ b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml @@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false evolve_moments_parallel_pressure = false evolve_moments_conservation = false force_Er_zero_at_wall = false #true -Er_constant = 0.0 T_e = 1.0 T_wall = 1.0 initial_density1 = 0.5 @@ -93,3 +92,5 @@ split_operators = false vpa_dissipation_coefficient = -1.0 z_dissipation_coefficient = -1.0 r_dissipation_coefficient = -1.0 +[geometry] +Er_constant = 0.0 diff --git a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml index d9ddc1aeb..3378a93b1 100644 --- a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml +++ b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml @@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false evolve_moments_parallel_pressure = false evolve_moments_conservation = false force_Er_zero_at_wall = false #true -Er_constant = 0.0 T_e = 1.0 T_wall = 1.0 initial_density1 = 0.5 @@ -96,3 +95,5 @@ split_operators = false vpa_dissipation_coefficient = -1.0 z_dissipation_coefficient = -1.0 r_dissipation_coefficient = -1.0 +[geometry] +Er_constant = 0.0 diff --git a/runs/2D-mirror_MMS_ngrid_5_nel_r_16_z_16_vpa_16_vperp_8_diss.toml b/runs/2D-mirror_MMS_ngrid_5_nel_r_16_z_16_vpa_16_vperp_8_diss.toml index a4a31ccbb..29542c43f 100644 --- a/runs/2D-mirror_MMS_ngrid_5_nel_r_16_z_16_vpa_16_vperp_8_diss.toml +++ b/runs/2D-mirror_MMS_ngrid_5_nel_r_16_z_16_vpa_16_vperp_8_diss.toml @@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false evolve_moments_parallel_pressure = false evolve_moments_conservation = false force_Er_zero_at_wall = false #true -Er_constant = 0.0 T_e = 1.0 T_wall = 1.0 initial_density1 = 0.5 @@ -99,3 +98,4 @@ DeltaB=0.5 #option="constant-helical" pitch=0.5 rhostar = 1.0 +Er_constant = 0.0 diff --git a/runs/2D-mirror_MMS_ngrid_5_nel_r_32_z_32_vpa_32_vperp_16_diss.toml b/runs/2D-mirror_MMS_ngrid_5_nel_r_32_z_32_vpa_32_vperp_16_diss.toml index 01dd1b487..3561f1d11 100644 --- a/runs/2D-mirror_MMS_ngrid_5_nel_r_32_z_32_vpa_32_vperp_16_diss.toml +++ b/runs/2D-mirror_MMS_ngrid_5_nel_r_32_z_32_vpa_32_vperp_16_diss.toml @@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false evolve_moments_parallel_pressure = false evolve_moments_conservation = false force_Er_zero_at_wall = false #true -Er_constant = 0.0 T_e = 1.0 T_wall = 1.0 initial_density1 = 0.5 @@ -99,3 +98,4 @@ DeltaB=0.5 #option="constant-helical" pitch=0.5 rhostar = 1.0 +Er_constant = 0.0 diff --git a/runs/2D-mirror_MMS_ngrid_5_nel_r_4_z_4_vpa_4_vperp_2_diss.toml b/runs/2D-mirror_MMS_ngrid_5_nel_r_4_z_4_vpa_4_vperp_2_diss.toml index b13205d70..a8add1f3e 100644 --- a/runs/2D-mirror_MMS_ngrid_5_nel_r_4_z_4_vpa_4_vperp_2_diss.toml +++ b/runs/2D-mirror_MMS_ngrid_5_nel_r_4_z_4_vpa_4_vperp_2_diss.toml @@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false evolve_moments_parallel_pressure = false evolve_moments_conservation = false force_Er_zero_at_wall = false #true -Er_constant = 0.0 T_e = 1.0 T_wall = 1.0 initial_density1 = 0.5 @@ -99,3 +98,4 @@ DeltaB=0.5 #option="constant-helical" pitch=0.5 rhostar = 1.0 +Er_constant = 0.0 diff --git a/runs/2D-mirror_MMS_ngrid_5_nel_r_8_z_8_vpa_8_vperp_4_diss.toml b/runs/2D-mirror_MMS_ngrid_5_nel_r_8_z_8_vpa_8_vperp_4_diss.toml index f19b91e39..2eb57f6ee 100644 --- a/runs/2D-mirror_MMS_ngrid_5_nel_r_8_z_8_vpa_8_vperp_4_diss.toml +++ b/runs/2D-mirror_MMS_ngrid_5_nel_r_8_z_8_vpa_8_vperp_4_diss.toml @@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false evolve_moments_parallel_pressure = false evolve_moments_conservation = false force_Er_zero_at_wall = false #true -Er_constant = 0.0 T_e = 1.0 T_wall = 1.0 initial_density1 = 0.5 @@ -99,3 +98,4 @@ DeltaB=0.5 #option="constant-helical" pitch=0.5 rhostar = 1.0 +Er_constant = 0.0 diff --git a/runs/2D-wall_MMS_nel_r_32_z_32_vpa_16_vperp_1_diss.toml b/runs/2D-wall_MMS_nel_r_32_z_32_vpa_16_vperp_1_diss.toml index 274fc59e6..fa1f6c220 100644 --- a/runs/2D-wall_MMS_nel_r_32_z_32_vpa_16_vperp_1_diss.toml +++ b/runs/2D-wall_MMS_nel_r_32_z_32_vpa_16_vperp_1_diss.toml @@ -7,7 +7,6 @@ evolve_moments_parallel_flow = false evolve_moments_parallel_pressure = false evolve_moments_conservation = false force_Er_zero_at_wall = false #true -Er_constant = 0.0 T_e = 1.0 T_wall = 1.0 initial_density1 = 0.5 @@ -95,3 +94,4 @@ r_dissipation_coefficient = 0.1 option="constant-helical" pitch=0.5 rhostar = 1.0 +Er_constant = 0.0 From 1ad159cadb03364a4bbabaa9bc7d9aec29a8db6c Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Mon, 8 Jul 2024 13:50:36 +0100 Subject: [PATCH 19/34] Fix various bugs in last commit so that refactored code runs and MMS tests can be executed. --- .../makie_post_processing/src/shared_utils.jl | 19 +++++-------------- moment_kinetics/ext/manufactured_solns_ext.jl | 8 ++++---- moment_kinetics/src/em_fields.jl | 4 ++-- moment_kinetics/src/geo.jl | 6 +++--- moment_kinetics/test/gyroaverage_tests.jl | 12 +++++------- .../src/plot_MMS_sequence.jl | 2 +- .../src/plots_post_processing.jl | 2 +- 7 files changed, 21 insertions(+), 32 deletions(-) diff --git a/makie_post_processing/makie_post_processing/src/shared_utils.jl b/makie_post_processing/makie_post_processing/src/shared_utils.jl index 031bcda75..f168bd457 100644 --- a/makie_post_processing/makie_post_processing/src/shared_utils.jl +++ b/makie_post_processing/makie_post_processing/src/shared_utils.jl @@ -5,7 +5,6 @@ export calculate_and_write_frequencies, get_geometry, get_composition using moment_kinetics.analysis: fit_delta_phi_mode using moment_kinetics.array_allocation: allocate_float using moment_kinetics.coordinates: define_coordinate -using moment_kinetics.geo: init_magnetic_geometry using moment_kinetics.input_structs: boltzmann_electron_response, boltzmann_electron_response_with_simple_sheath, grid_input, geometry_input, species_composition @@ -13,7 +12,7 @@ using moment_kinetics.moment_kinetics_input: get_default_rhostar, setup_referenc using moment_kinetics.type_definitions: mk_float, mk_int using moment_kinetics.reference_parameters: setup_reference_parameters using moment_kinetics.moment_kinetics_input: get_default_rhostar -using moment_kinetics.geo: init_magnetic_geometry +using moment_kinetics.geo: init_magnetic_geometry, setup_geometry_input using MPI """ @@ -90,7 +89,7 @@ function get_composition(scan_input) use_test_neutral_wall_pdf = get(scan_input, "use_test_neutral_wall_pdf", false) gyrokinetic_ions = get(scan_input, "gyrokinetic_ions", false) # constant to be used to test nonzero Er in wall boundary condition - Er_constant = get(scan_input, "Er_constant", 0.0) + #Er_constant = get(scan_input, "Er_constant", 0.0) recycling_fraction = get(scan_input, "recycling_fraction", 1.0) # constant to be used to control Ez divergences epsilon_offset = get(scan_input, "epsilon_offset", 0.001) @@ -106,7 +105,7 @@ function get_composition(scan_input) # ratio of the electron particle mass to the ion particle mass me_over_mi = 1.0/1836.0 composition = species_composition(n_species, n_ion_species, n_neutral_species, - electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, Er_constant, + electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, #Er_constant, mn_over_mi, me_over_mi, recycling_fraction, gyrokinetic_ions, allocate_float(n_species)) return composition @@ -114,18 +113,10 @@ end function get_geometry(scan_input,z,r) reference_params = setup_reference_parameters(scan_input) - # set geometry_input - # MRH need to get this in way that does not duplicate code - # MRH from moment_kinetics_input.jl - option = get(scan_input, "geometry_option", "constant-helical") #"1D-mirror" - pitch = get(scan_input, "pitch", 1.0) - rhostar = get(scan_input, "rhostar", get_default_rhostar(reference_params)) - DeltaB = get(scan_input, "DeltaB", 1.0) - geo_in = geometry_input(rhostar,option,pitch,DeltaB) + reference_rhostar = get_default_rhostar(reference_params) + geo_in = setup_geometry_input(scan_input, reference_rhostar) geometry = init_magnetic_geometry(geo_in,z,r) - return geometry - end end # shared_utils.jl diff --git a/moment_kinetics/ext/manufactured_solns_ext.jl b/moment_kinetics/ext/manufactured_solns_ext.jl index 8f1f42650..314853725 100644 --- a/moment_kinetics/ext/manufactured_solns_ext.jl +++ b/moment_kinetics/ext/manufactured_solns_ext.jl @@ -362,7 +362,7 @@ using IfElse if manufactured_solns_input.type == "default" # calculate the electric fields and the potential - Er, Ez, phi = electric_fields(Lr, Lz, r_bc, z_bc, composition, nr, + Er, Ez, phi = electric_fields(Lr, Lz, r_bc, z_bc, composition, geometry, nr, manufactured_solns_input, species) # get geometric/composition data @@ -417,7 +417,7 @@ using IfElse # get N_e factor for boltzmann response if composition.electron_physics == boltzmann_electron_response_with_simple_sheath && nr == 1 # so 1D MMS test with 3V neutrals where ion current can be calculated prior to knowing Er - jpari_into_LHS_wall = jpari_into_LHS_wall_sym(Lr, Lz, r_bc, z_bc, composition, + jpari_into_LHS_wall = jpari_into_LHS_wall_sym(Lr, Lz, r_bc, z_bc, manufactured_solns_input) N_e = -2.0*sqrt(pi*composition.me_over_mi)*exp(-composition.phi_wall/composition.T_e)*jpari_into_LHS_wall elseif composition.electron_physics == boltzmann_electron_response_with_simple_sheath && nr > 1 @@ -442,8 +442,8 @@ using IfElse # calculate the electric fields dense = densi # get the electron density via quasineutrality with Zi = 1 phi = composition.T_e*log(dense/N_e) # use the adiabatic response of electrons for me/mi -> 0 - Er = -Dr(phi)*rfac + composition.Er_constant - Ez = -Dz(phi) + Er = -Dr(phi)*rfac + geometry.Er_constant + Ez = -Dz(phi) + geometry.Ez_constant Er_expanded = expand_derivatives(Er) Ez_expanded = expand_derivatives(Ez) diff --git a/moment_kinetics/src/em_fields.jl b/moment_kinetics/src/em_fields.jl index 5c021099b..67534d3c3 100644 --- a/moment_kinetics/src/em_fields.jl +++ b/moment_kinetics/src/em_fields.jl @@ -136,8 +136,8 @@ function update_phi!(fields, fvec, vperp, z, r, composition, geometry, z_spectra scratch_dummy.buffer_r_3, scratch_dummy.buffer_r_4, z_spectral,z) else - @serial_region begin - fields.Ez[:,:] .= = geometry.input.Ez_constant + @loop_r_z ir iz begin + fields.Ez[iz,ir] = geometry.input.Ez_constant # Ez_constant defaults to 0.0 in geo.jl end end diff --git a/moment_kinetics/src/geo.jl b/moment_kinetics/src/geo.jl index e9c188bb8..75637425a 100644 --- a/moment_kinetics/src/geo.jl +++ b/moment_kinetics/src/geo.jl @@ -8,7 +8,7 @@ module geo export init_magnetic_geometry export setup_geometry_input -using ..input_structs: geometry_input +using ..input_structs: geometry_input, set_defaults_and_check_section! using ..file_io: input_option_error using ..array_allocation: allocate_float using ..type_definitions: mk_float, mk_int @@ -77,9 +77,9 @@ function setup_geometry_input(toml_input::Dict, reference_rhostar) input_section = set_defaults_and_check_section!(toml_input, "geometry", # begin default inputs (as kwargs) # rhostar ion (ref) - rhostar = reference_rhostar #used to premultiply ExB drift terms + rhostar = reference_rhostar, #used to premultiply ExB drift terms # magnetic geometry option - option = "constant-helical" # "1D-mirror" + option = "constant-helical",# "1D-mirror" # pitch ( = Bzed/Bmag if geometry_option == "constant-helical") pitch = 1.0, # DeltaB ( = (Bzed(z=L/2) - Bzed(0))/Bref if geometry_option == "1D-mirror") diff --git a/moment_kinetics/test/gyroaverage_tests.jl b/moment_kinetics/test/gyroaverage_tests.jl index 38a6c7085..21c43f0b7 100644 --- a/moment_kinetics/test/gyroaverage_tests.jl +++ b/moment_kinetics/test/gyroaverage_tests.jl @@ -10,7 +10,7 @@ using SpecialFunctions: besselj0 import moment_kinetics using moment_kinetics.input_structs using moment_kinetics.coordinates: define_coordinate -using moment_kinetics.geo: init_magnetic_geometry +using moment_kinetics.geo: init_magnetic_geometry, setup_geometry_input using moment_kinetics.communication using moment_kinetics.looping using moment_kinetics.array_allocation: allocate_float, allocate_shared_float @@ -142,11 +142,9 @@ function gyroaverage_test(absolute_error; rhostar=0.1, pitch=0.5, ngrid=5, kr=2, gyrophase, gyrophase_spectral = define_coordinate(gyrophase_input; collision_operator_dim=false, run_directory=test_output_directory) # create test geometry - #rhostar = 0.1 #rhostar of ions for ExB drift option = "constant-helical" - #pitch = 1.0 - DeltaB = 1.0 - geometry_in = geometry_input(rhostar,option,pitch,DeltaB) + inputdict = Dict("geometry" => Dict("option" => option, "rhostar" => rhostar, "pitch" => pitch)) + geometry_in = setup_geometry_input(inputdict, 0.1) geometry = init_magnetic_geometry(geometry_in,z,r) # create test composition @@ -243,7 +241,7 @@ function create_test_composition() # wall potential at z = 0 phi_wall = 0.0 # constant to test nonzero Er - Er_constant = 0.0 + #Er_constant = 0.0 # ratio of the neutral particle mass to the ion particle mass mn_over_mi = 1.0 # ratio of the electron particle mass to the ion particle mass @@ -253,7 +251,7 @@ function create_test_composition() recycling_fraction = 1.0 gyrokinetic_ions = true return composition = species_composition(n_species, n_ion_species, n_neutral_species, - electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, Er_constant, + electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, #Er_constant, mn_over_mi, me_over_mi, recycling_fraction, gyrokinetic_ions, allocate_float(n_species)) end diff --git a/plots_post_processing/plots_post_processing/src/plot_MMS_sequence.jl b/plots_post_processing/plots_post_processing/src/plot_MMS_sequence.jl index 84bfcf21b..92d780b27 100644 --- a/plots_post_processing/plots_post_processing/src/plot_MMS_sequence.jl +++ b/plots_post_processing/plots_post_processing/src/plot_MMS_sequence.jl @@ -255,7 +255,7 @@ function get_MMS_error_data(path_list,scan_type,scan_name) dfnn_func = manufactured_solns_list.dfnn_func densn_func = manufactured_solns_list.densn_func - manufactured_E_fields = manufactured_electric_fields(Lr_in,z.L,r_bc,z_bc,composition,r.n,manufactured_solns_input,species) + manufactured_E_fields = manufactured_electric_fields(Lr_in,z.L,r_bc,z_bc,composition,geo_in,r.n,manufactured_solns_input,species) Er_func = manufactured_E_fields.Er_func Ez_func = manufactured_E_fields.Ez_func phi_func = manufactured_E_fields.phi_func diff --git a/plots_post_processing/plots_post_processing/src/plots_post_processing.jl b/plots_post_processing/plots_post_processing/src/plots_post_processing.jl index a29b57dd7..81d1e3828 100644 --- a/plots_post_processing/plots_post_processing/src/plots_post_processing.jl +++ b/plots_post_processing/plots_post_processing/src/plots_post_processing.jl @@ -1097,7 +1097,7 @@ function analyze_and_plot_data(prefix...; run_index=nothing) densn_func = manufactured_solns_list.densn_func manufactured_E_fields = manufactured_electric_fields(Lr_in, z_global.L, r_global.bc, z_global.bc, - composition, r_global.n, + composition, geometry.input, r_global.n, manufactured_solns_input, species) Er_func = manufactured_E_fields.Er_func Ez_func = manufactured_E_fields.Ez_func 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 20/34] 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 From 8bff6a5658e8e3012fb8c0bef81c01c02760fb22 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Wed, 10 Jul 2024 11:56:02 +0100 Subject: [PATCH 21/34] Changes to permit testing of "none" boundary condition option for Fokker-Planck time evolution. In this mode of operation, boundary conditions on the fluxes are imposed through the assembly of the finite-element method, and no additional zero-ing of the distribution function is carried out. Initial tests suggest that this mode of operation is at least as good as other tested options. --- moment_kinetics/src/boundary_conditions.jl | 2 + moment_kinetics/src/fokker_planck_calculus.jl | 16 ++++-- .../fokker_planck_time_evolution_tests.jl | 56 +++++++++++++++++++ 3 files changed, 68 insertions(+), 6 deletions(-) diff --git a/moment_kinetics/src/boundary_conditions.jl b/moment_kinetics/src/boundary_conditions.jl index 7082c1383..27282d65e 100644 --- a/moment_kinetics/src/boundary_conditions.jl +++ b/moment_kinetics/src/boundary_conditions.jl @@ -954,6 +954,8 @@ function enforce_v_boundary_condition_local!(f, bc, speed, v_diffusion, v, v_spe elseif bc == "periodic" f[1] = 0.5*(f[1]+f[end]) f[end] = f[1] + elseif bc == "none" + # Do nothing else error("Unsupported boundary condition option '$bc' for $(v.name)") end diff --git a/moment_kinetics/src/fokker_planck_calculus.jl b/moment_kinetics/src/fokker_planck_calculus.jl index 7838608fe..4007e6356 100644 --- a/moment_kinetics/src/fokker_planck_calculus.jl +++ b/moment_kinetics/src/fokker_planck_calculus.jl @@ -2287,18 +2287,22 @@ function enforce_vpavperp_BCs!(pdf,vpa,vperp,vpa_spectral,vperp_spectral) D0 = vperp_spectral.radau.D0 # vpa boundary conditions # zero at infinity - begin_anyv_vperp_region() - @loop_vperp ivperp begin - pdf[1,ivperp] = 0.0 - pdf[nvpa,ivperp] = 0.0 + if vpa.bc == "zero" + begin_anyv_vperp_region() + @loop_vperp ivperp begin + pdf[1,ivperp] = 0.0 + pdf[nvpa,ivperp] = 0.0 + end end # vperp boundary conditions # zero boundary condition at infinity # set regularity condition d F / d vperp = 0 at vperp = 0 # adjust F(vperp = 0) so that d F / d vperp = 0 at vperp = 0 begin_anyv_vpa_region() - @loop_vpa ivpa begin - pdf[ivpa,nvperp] = 0.0 + if vperp.bc in ("zero", "zero-no-regularity") + @loop_vpa ivpa begin + pdf[ivpa,nvperp] = 0.0 + end end if vperp.bc == "zero" buffer = @view vperp.scratch[1:ngrid_vperp-1] diff --git a/moment_kinetics/test/fokker_planck_time_evolution_tests.jl b/moment_kinetics/test/fokker_planck_time_evolution_tests.jl index 669af32de..b5d5086d6 100644 --- a/moment_kinetics/test/fokker_planck_time_evolution_tests.jl +++ b/moment_kinetics/test/fokker_planck_time_evolution_tests.jl @@ -128,6 +128,54 @@ expected_data( 0.007980077110978 0.005701898246888 0.003327438535529 0.000925955494251 -0.000037963012435 -0.000174809935222 0.000000000000000 ; 0.000473874379555 0.000190726186343 -0.000067413540342 -0.000219129923463 -0.000165609965457 -0.000069940808382 0.000000000000000 ; 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000]) + +const expected_none_bc = +expected_data( + [-3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0], + [0.155051025721682, 0.644948974278318, 1.000000000000000, 1.500000000000000, 2.000000000000000, 2.500000000000000, 3.000000000000000], + # Expected phi: + [-1.254104813718096, -1.254104813718360], + # Expected n_ion: + [0.285331158471152, 0.285331158471077], + # Expected upar_ion: + [0.0, 0.0], + # Expected ppar_ion: + [0.182321263834348, 0.154595996906667], + # Expected pperp_ion + [0.143470130069393, 0.157332763533171], + # Expected qpar_ion + [0.0, 0.0], + # Expected v_t_ion + [1.047097792428007, 1.047097792428005], + # Expected dSdt + [0.000000000000000, 0.000000019115425], + # Expected f_ion: + [0.000045179366280 0.000030530376095 0.000017024973661 0.000004877736620 0.000000847623528 0.000000089338863 0.000000005711242 ; + 0.000706724195475 0.000477575434536 0.000266315395816 0.000076300638379 0.000013259062819 0.000001397494940 0.000000089338863 ; + 0.006705212475828 0.004531109564814 0.002526730124657 0.000723920301085 0.000125798485463 0.000013259062819 0.000000847623528 ; + 0.038585833650058 0.026074735222541 0.014540327934436 0.004165873701136 0.000723920300963 0.000076300638366 0.000004877736619 ; + 0.134671678398729 0.091005636629981 0.050748427134097 0.014539667807029 0.002526615411287 0.000266303305116 0.000017024200728 ; + 0.261709398369233 0.176852554997681 0.098620144126568 0.028255144359305 0.004910007858078 0.000517511020834 0.000033083372713 ; + 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 ; + 0.261709398369233 0.176852554997682 0.098620144126568 0.028255144359305 0.004910007858078 0.000517511020834 0.000033083372713 ; + 0.134671678398729 0.091005636629981 0.050748427134097 0.014539667807029 0.002526615411287 0.000266303305116 0.000017024200728 ; + 0.038585833650058 0.026074735222541 0.014540327934436 0.004165873701136 0.000723920300963 0.000076300638366 0.000004877736619 ; + 0.006705212475828 0.004531109564814 0.002526730124657 0.000723920301085 0.000125798485463 0.000013259062819 0.000000847623528 ; + 0.000706724195475 0.000477575434536 0.000266315395816 0.000076300638379 0.000013259062819 0.000001397494940 0.000000089338863 ; + 0.000045179366280 0.000030530376095 0.000017024973661 0.000004877736620 0.000000847623528 0.000000089338863 0.000000005711242 ;;; + 0.000447615535468 0.000364801746561 0.000270093752229 0.000138454732861 0.000052917434226 0.000014331973588 -0.000000022222325 ; + 0.000530676740860 0.000337644325418 0.000161654162384 0.000044642657051 0.000030848860156 0.000021605350900 0.000014496785601 ; + 0.006652401806692 0.004725322877613 0.002773765135978 0.000908975538425 0.000210924069594 0.000037798409136 0.000057517483301 ; + 0.030674916797523 0.021423584760802 0.012211671337997 0.003838328661082 0.000895825759878 0.000045643481270 0.000159762122747 ; + 0.097751022951371 0.068620792366922 0.039365625553724 0.012397445282748 0.002818297719268 0.000149127033057 0.000311439525732 ; + 0.194564361178625 0.137570544970062 0.079900546066037 0.025549350260855 0.005673491663727 0.000408759661293 0.000445141805523 ; + 0.242684098382075 0.171867841736931 0.100113917461427 0.032148650920148 0.007091277425658 0.000554833974109 0.000500018986081 ; + 0.194564361178625 0.137570544970062 0.079900546066037 0.025549350260855 0.005673491663727 0.000408759661293 0.000445141805523 ; + 0.097751022951371 0.068620792366922 0.039365625553724 0.012397445282748 0.002818297719268 0.000149127033057 0.000311439525732 ; + 0.030674916797523 0.021423584760802 0.012211671337997 0.003838328661082 0.000895825759878 0.000045643481270 0.000159762122747 ; + 0.006652401806692 0.004725322877613 0.002773765135978 0.000908975538425 0.000210924069594 0.000037798409136 0.000057517483301 ; + 0.000530676740860 0.000337644325418 0.000161654162384 0.000044642657051 0.000030848860156 0.000021605350900 0.000014496785601 ; + 0.000447615535468 0.000364801746561 0.000270093752229 0.000138454732861 0.000052917434226 0.000014331973588 -0.000000022222325]) ########################################################################################### # to modify the test, with a new expected f, print the new f using the following commands # in an interative Julia REPL. The path is the path to the .dfns file. @@ -399,6 +447,14 @@ function runtests() expected_no_regularity, 1.0e-14, 1.0e-14; vperp_bc="zero-no-regularity") end + @testset "Gauss Legendre no (explicitly) enforced boundary conditions" begin + run_name = "gausslegendre_pseudospectral_none_bc" + vperp_bc = "none" + vpa_bc = "none" + run_test(test_input_gauss_legendre, + expected_none_bc, + 1.0e-14, 1.0e-14; vperp_bc=vperp_bc, vpa_bc=vpa_bc) + end end end From 37777b0cd8382843fc95be4df3ee7bb02059b1f4 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Tue, 30 Jul 2024 16:56:01 +0000 Subject: [PATCH 22/34] Remove old example file. --- ...ation-slowing-down-alphas-source-sink.toml | 91 ------------------- 1 file changed, 91 deletions(-) delete mode 100644 examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas-source-sink.toml diff --git a/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas-source-sink.toml b/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas-source-sink.toml deleted file mode 100644 index 098dc787a..000000000 --- a/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas-source-sink.toml +++ /dev/null @@ -1,91 +0,0 @@ -# cheap input file for a 0D2V relaxation to a collisional Maxwellian distribution with self-ion collisions and collisions with fixed Maxwellian background of cold ions and electrons. -n_ion_species = 1 -n_neutral_species = 0 -electron_physics = "boltzmann_electron_response" -evolve_moments_density = false -evolve_moments_parallel_flow = false -evolve_moments_parallel_pressure = false -evolve_moments_conservation = false -T_e = 1.0 -T_wall = 1.0 -initial_density1 = 1.0 -initial_temperature1 = 1.0 -z_IC_option1 = "sinusoid" -z_IC_density_amplitude1 = 0.0 -z_IC_density_phase1 = 0.0 -z_IC_upar_amplitude1 = 0.0 -z_IC_upar_phase1 = 0.0 -z_IC_temperature_amplitude1 = 0.0 -z_IC_temperature_phase1 = 0.0 -vpa_IC_option1 = "isotropic-beam" -vpa_IC_v01 = 1.0 -vpa_IC_vth01 = 0.1 -#vpa_IC_option1 = "directed-beam" -#vpa_IC_vpa01 = -1.5 -#vpa_IC_vperp01 = 0.0 -charge_exchange_frequency = 0.0 -ionization_frequency = 0.0 -constant_ionization_rate = false - -z_ngrid = 1 -z_nelement = 1 -z_nelement_local = 1 -z_bc = "wall" -z_discretization = "chebyshev_pseudospectral" -r_ngrid = 1 -r_nelement = 1 -r_nelement_local = 1 -r_bc = "periodic" -r_discretization = "chebyshev_pseudospectral" -vpa_ngrid = 5 -vpa_nelement = 32 -vpa_L = 3.0 -vpa_bc = "zero" -vpa_discretization = "gausslegendre_pseudospectral" -vperp_ngrid = 5 -vperp_nelement = 16 -vperp_L = 1.5 -vperp_discretization = "gausslegendre_pseudospectral" -vperp_bc = "zero" -# Fokker-Planck operator requires the "gausslegendre_pseudospectral -# options for the vpa and vperp grids - -[fokker_planck_collisions] -# nuii sets the normalised input C[F,F] Fokker-Planck collision frequency -nuii = 1.876334222e-3 #(1/nu_alphae, as computed from input diagnostic) -Zi = 2.0 -self_collisions = true -slowing_down_test = true -frequency_option = "manual" -use_fokker_planck = true -sd_density = 1.0 -sd_temp = 0.0025 # TD/Ealpha -sd_mi = 0.5 # mD/malpha -sd_me = 0.000013616 # 0.25/1836.0 me/malpha -sd_q = 1.0 - -[ion_source] -active=true -source_strength=1.0 -source_T=0.005 -source_n=1.0 -r_profile="constant" -r_width=1.0 -r_relative_minimum=0.0 -z_profile="gaussian" -z_width=0.1 -z_relative_minimum=0.0 -source_v0 = 1.0 -#source_type="alphas" -source_type="alphas-with-losses" -#source_type="beam-with-losses" -#source_vpa0 = 1.0 -#source_vperp0 = 1.0 -sink_strength = 1.0 -sink_vth=0.07071067811865475 - -[timestepping] -nstep = 50000 -dt = 5.0e-4 -nwrite = 100 -nwrite_dfns = 100 From 368302276d375f2eb5ad5b4ce697f983921ae3ff Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Tue, 30 Jul 2024 16:57:25 +0000 Subject: [PATCH 23/34] Add plot of pperp(z) at final time. --- moment_kinetics/src/input_structs.jl | 2 ++ .../src/plots_post_processing.jl | 30 ++++++++++++++----- .../src/post_processing_input.jl | 4 ++- 3 files changed, 27 insertions(+), 9 deletions(-) diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl index cddb11226..f59fbb049 100644 --- a/moment_kinetics/src/input_structs.jl +++ b/moment_kinetics/src/input_structs.jl @@ -553,6 +553,8 @@ struct pp_input animate_parallel_flow_vs_r_z::Bool # if plot_parallel_pressure_vs_r0_z = true plot last timestep parallel_pressure[z,ir0] plot_parallel_pressure_vs_r0_z::Bool + # if plot_perpendicular_pressure_vs_r0_z = true plot last timestep perpendicular_pressure[z,ir0] + plot_perpendicular_pressure_vs_r0_z::Bool # if plot_wall_parallel_pressure_vs_r = true plot last timestep parallel_pressure[z_wall,r] plot_wall_parallel_pressure_vs_r::Bool # if plot_parallel_pressure_vs_r_z = true plot parallel_pressure vs r z at last timestep diff --git a/plots_post_processing/plots_post_processing/src/plots_post_processing.jl b/plots_post_processing/plots_post_processing/src/plots_post_processing.jl index 81d1e3828..8fdb8fd81 100644 --- a/plots_post_processing/plots_post_processing/src/plots_post_processing.jl +++ b/plots_post_processing/plots_post_processing/src/plots_post_processing.jl @@ -3308,33 +3308,47 @@ function plot_ion_moments_2D(density, parallel_flow, parallel_pressure, outfile = string(run_name, "_delta_perpendicular_pressure"*description*"(iz0,ir0)_vs_t.pdf") trysavefig(outfile) end + if pp.plot_perpendicular_pressure_vs_r0_z && nz > 1 # plot last timestep perpendicular_pressure[z,ir0] + @views plot(z, perpendicular_pressure[:,ir0,is,end], xlabel=L"z/L_z", ylabel=L"p_{i\perp}",label="") + outfile = string(run_name, "_perpendicular_pressure"*description*"(r0,z)_vs_z.pdf") + trysavefig(outfile) + end # the total pressure + total_pressure = copy(parallel_pressure) + @. total_pressure = (2.0/3.0)*perpendicular_pressure + (1.0/3.0)*parallel_pressure if pp.plot_ppar0_vs_t && pp.plot_pperp0_vs_t + @views plot([time, time, time] , [parallel_pressure[iz0,ir0,is,:], perpendicular_pressure[iz0,ir0,is,:], - (2.0/3.0).*perpendicular_pressure[iz0,ir0,is,:] .+ (1.0/3.0).*parallel_pressure[iz0,ir0,is,:]], - xlabel=L"t/ (L_{ref}/c_{ref})", ylabel="", label = [L"p_{i\|\|}(t)" L"p_{i\perp}(t)" L"p_{i}(t)"]) + total_pressure[iz0,ir0,is,:]], + xlabel=L"t/ (L_{ref}/c_{ref})", ylabel="", label = [L"p_{i\|\|}(t)" L"p_{i\perp}(t)" L"p_{i}(t)"], + foreground_color_legend = nothing, background_color_legend = nothing) outfile = string(run_name, "_pressures"*description*"(iz0,ir0)_vs_t.pdf") trysavefig(outfile) @views plot([time, time, time] , [parallel_pressure[iz0,ir0,is,:] .- parallel_pressure[iz0,ir0,is,1], perpendicular_pressure[iz0,ir0,is,:] .- perpendicular_pressure[iz0,ir0,is,1], - (2.0/3.0).*(perpendicular_pressure[iz0,ir0,is,:] .- perpendicular_pressure[iz0,ir0,is,1]).+ - (1.0/3.0).*(parallel_pressure[iz0,ir0,is,:] .- parallel_pressure[iz0,ir0,is,1])], - xlabel=L"t/ (L_{ref}/c_{ref})", ylabel="", label = [L"p_{i\|\|}(t) - p_{i\|\|}(0)" L"p_{i\perp}(t) - p_{i\perp}(0)" L"p_{i}(t) - p_{i}(0)"]) + total_pressure[iz0,ir0,is,:] .- total_pressure[iz0,ir0,is,1]], + xlabel=L"t/ (L_{ref}/c_{ref})", ylabel="", label = [L"p_{i\|\|}(t) - p_{i\|\|}(0)" L"p_{i\perp}(t) - p_{i\perp}(0)" L"p_{i}(t) - p_{i}(0)"], + foreground_color_legend = nothing, background_color_legend = nothing) outfile = string(run_name, "_delta_pressures"*description*"(iz0,ir0)_vs_t.pdf") trysavefig(outfile) @views plot([time] , - [(2.0/3.0).*perpendicular_pressure[iz0,ir0,is,:] .+ (1.0/3.0).*parallel_pressure[iz0,ir0,is,:]], + [total_pressure[iz0,ir0,is,:]], xlabel=L"t/ (L_{ref}/c_{ref})", ylabel=L"p_{i}(t)", label = "") outfile = string(run_name, "_pressure"*description*"(iz0,ir0)_vs_t.pdf") trysavefig(outfile) @views plot([time] , - [(2.0/3.0).*(perpendicular_pressure[iz0,ir0,is,:] .- perpendicular_pressure[iz0,ir0,is,1]).+ - (1.0/3.0).*(parallel_pressure[iz0,ir0,is,:] .- parallel_pressure[iz0,ir0,is,1])], + [total_pressure[iz0,ir0,is,:] .- total_pressure[iz0,ir0,is,1]], xlabel=L"t/ (L_{ref}/c_{ref})", ylabel=L"p_{i}(t) - p_{i}(0)", label = "") outfile = string(run_name, "_delta_pressure"*description*"(iz0,ir0)_vs_t.pdf") trysavefig(outfile) end + if pp.plot_perpendicular_pressure_vs_r0_z && pp.plot_parallel_pressure_vs_r0_z && nz > 1 # plot last timestep perpendicular_pressure[z,ir0] + @views plot([z,z,z], [parallel_pressure[:,ir0,is,end],perpendicular_pressure[:,ir0,is,end],total_pressure[:,ir0,is,end]], xlabel=L"z/L_z", ylabel="", + label = [L"p_{i\|\|}" L"p_{i\perp}" L"p_{i}"], foreground_color_legend = nothing, background_color_legend = nothing) + outfile = string(run_name, "_all_pressures"*description*"(r0,z)_vs_z.pdf") + trysavefig(outfile) + end # the thermal speed if pp.plot_vth0_vs_t @views plot(time, thermal_speed[iz0,ir0,is,:], xlabel=L"t / (L_{ref}/c_{ref})", ylabel=L"v_{i,th}(t)", label = "") diff --git a/plots_post_processing/plots_post_processing/src/post_processing_input.jl b/plots_post_processing/plots_post_processing/src/post_processing_input.jl index 5c3b6f58d..08cad55f8 100644 --- a/plots_post_processing/plots_post_processing/src/post_processing_input.jl +++ b/plots_post_processing/plots_post_processing/src/post_processing_input.jl @@ -102,6 +102,7 @@ const plot_wall_density_vs_r = true # plot last timestep density[z_wall,r] const plot_density_vs_r_z = true const animate_density_vs_r_z = true const plot_parallel_flow_vs_r0_z = true # plot last timestep parallel_flow[z,ir0] +const plot_perpendicular_flow_vs_r0_z = true # plot last timestep perpendicular_flow[z,ir0] const plot_wall_parallel_flow_vs_r = true # plot last timestep parallel_flow[z_wall,r] const plot_parallel_flow_vs_r_z = true const animate_parallel_flow_vs_r_z = true @@ -171,7 +172,8 @@ pp = pp_input(calculate_frequencies, plot_phi0_vs_t, plot_phi_vs_z_t, animate_ph animate_Er_vs_r_z, animate_Ez_vs_r_z, animate_phi_vs_r_z, plot_phi_vs_r0_z, plot_Ez_vs_r0_z, plot_wall_Ez_vs_r, plot_Er_vs_r0_z, plot_wall_Er_vs_r, plot_density_vs_r0_z, plot_wall_density_vs_r, plot_density_vs_r_z, - animate_density_vs_r_z, plot_parallel_flow_vs_r0_z, plot_wall_parallel_flow_vs_r, + animate_density_vs_r_z, plot_parallel_flow_vs_r0_z, plot_perpendicular_flow_vs_r0_z, + plot_wall_parallel_flow_vs_r, plot_parallel_flow_vs_r_z, animate_parallel_flow_vs_r_z, plot_parallel_pressure_vs_r0_z, plot_wall_parallel_pressure_vs_r, plot_parallel_pressure_vs_r_z, animate_parallel_pressure_vs_r_z, From 276d6827960fd19019ec195260e86a2aad9af63b Mon Sep 17 00:00:00 2001 From: mrhardman <29800382+mrhardman@users.noreply.github.com> Date: Wed, 31 Jul 2024 11:44:05 +0100 Subject: [PATCH 24/34] Update moment_kinetics/src/input_structs.jl Co-authored-by: John Omotani --- moment_kinetics/src/input_structs.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl index 1554dbd4c..b740615c9 100644 --- a/moment_kinetics/src/input_structs.jl +++ b/moment_kinetics/src/input_structs.jl @@ -318,7 +318,6 @@ mutable struct species_composition # wall potential used if electron_physics=boltzmann_electron_response_with_simple_sheath phi_wall::mk_float # constant for testing nonzero Er - #Er_constant::mk_float # ratio of the neutral particle mass to the ion mass mn_over_mi::mk_float # ratio of the electron particle mass to the ion mass From 9193a563425d9b766514663ac5fe1c6cf216e621 Mon Sep 17 00:00:00 2001 From: mrhardman <29800382+mrhardman@users.noreply.github.com> Date: Wed, 31 Jul 2024 11:44:18 +0100 Subject: [PATCH 25/34] Update makie_post_processing/makie_post_processing/src/shared_utils.jl Co-authored-by: John Omotani --- makie_post_processing/makie_post_processing/src/shared_utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/makie_post_processing/makie_post_processing/src/shared_utils.jl b/makie_post_processing/makie_post_processing/src/shared_utils.jl index f168bd457..9c7f6303d 100644 --- a/makie_post_processing/makie_post_processing/src/shared_utils.jl +++ b/makie_post_processing/makie_post_processing/src/shared_utils.jl @@ -105,7 +105,7 @@ function get_composition(scan_input) # ratio of the electron particle mass to the ion particle mass me_over_mi = 1.0/1836.0 composition = species_composition(n_species, n_ion_species, n_neutral_species, - electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, #Er_constant, + electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, mn_over_mi, me_over_mi, recycling_fraction, gyrokinetic_ions, allocate_float(n_species)) return composition From 185598a7b0cfb47e336386f096dc76871698ee95 Mon Sep 17 00:00:00 2001 From: mrhardman <29800382+mrhardman@users.noreply.github.com> Date: Wed, 31 Jul 2024 11:44:25 +0100 Subject: [PATCH 26/34] Update makie_post_processing/makie_post_processing/src/shared_utils.jl Co-authored-by: John Omotani --- makie_post_processing/makie_post_processing/src/shared_utils.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/makie_post_processing/makie_post_processing/src/shared_utils.jl b/makie_post_processing/makie_post_processing/src/shared_utils.jl index 9c7f6303d..32acaf1ec 100644 --- a/makie_post_processing/makie_post_processing/src/shared_utils.jl +++ b/makie_post_processing/makie_post_processing/src/shared_utils.jl @@ -89,7 +89,6 @@ function get_composition(scan_input) use_test_neutral_wall_pdf = get(scan_input, "use_test_neutral_wall_pdf", false) gyrokinetic_ions = get(scan_input, "gyrokinetic_ions", false) # constant to be used to test nonzero Er in wall boundary condition - #Er_constant = get(scan_input, "Er_constant", 0.0) recycling_fraction = get(scan_input, "recycling_fraction", 1.0) # constant to be used to control Ez divergences epsilon_offset = get(scan_input, "epsilon_offset", 0.001) From 126cea21c1e0e55923589a9934f052d43edb6e34 Mon Sep 17 00:00:00 2001 From: mrhardman <29800382+mrhardman@users.noreply.github.com> Date: Wed, 31 Jul 2024 11:44:32 +0100 Subject: [PATCH 27/34] Update moment_kinetics/src/moment_kinetics_input.jl Co-authored-by: John Omotani --- moment_kinetics/src/moment_kinetics_input.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/moment_kinetics/src/moment_kinetics_input.jl b/moment_kinetics/src/moment_kinetics_input.jl index 9eb9d3b25..68b2366d3 100644 --- a/moment_kinetics/src/moment_kinetics_input.jl +++ b/moment_kinetics/src/moment_kinetics_input.jl @@ -1027,7 +1027,7 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) recycling_fraction = 1.0 gyrokinetic_ions = false composition = species_composition(n_species, n_ion_species, n_neutral_species, - electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, #Er_constant, + electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, mn_over_mi, me_over_mi, recycling_fraction, gyrokinetic_ions, allocate_float(n_species)) species_ion = Array{species_parameters_mutable,1}(undef,n_ion_species) From 8e041b7331c16503ccaadd4a4194bba07616d275 Mon Sep 17 00:00:00 2001 From: mrhardman <29800382+mrhardman@users.noreply.github.com> Date: Wed, 31 Jul 2024 11:44:40 +0100 Subject: [PATCH 28/34] Update moment_kinetics/src/moment_kinetics_input.jl Co-authored-by: John Omotani --- moment_kinetics/src/moment_kinetics_input.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/moment_kinetics/src/moment_kinetics_input.jl b/moment_kinetics/src/moment_kinetics_input.jl index 68b2366d3..7d9f79a12 100644 --- a/moment_kinetics/src/moment_kinetics_input.jl +++ b/moment_kinetics/src/moment_kinetics_input.jl @@ -1016,7 +1016,6 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) T_wall = 1.0 # wall potential at z = 0 phi_wall = 0.0 - # constant to test nonzero Er #Er_constant = 0.0 # ratio of the neutral particle mass to the ion particle mass mn_over_mi = 1.0 From 1a89877da4fe825ac2cdc0589c819d18162c7ad4 Mon Sep 17 00:00:00 2001 From: mrhardman <29800382+mrhardman@users.noreply.github.com> Date: Wed, 31 Jul 2024 11:44:47 +0100 Subject: [PATCH 29/34] Update moment_kinetics/src/moment_kinetics_input.jl Co-authored-by: John Omotani --- moment_kinetics/src/moment_kinetics_input.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/moment_kinetics/src/moment_kinetics_input.jl b/moment_kinetics/src/moment_kinetics_input.jl index 7d9f79a12..a80c03751 100644 --- a/moment_kinetics/src/moment_kinetics_input.jl +++ b/moment_kinetics/src/moment_kinetics_input.jl @@ -1016,7 +1016,6 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) T_wall = 1.0 # wall potential at z = 0 phi_wall = 0.0 - #Er_constant = 0.0 # ratio of the neutral particle mass to the ion particle mass mn_over_mi = 1.0 # ratio of the electron particle mass to the ion particle mass From a6e13a981c895db011bb6ce2fe6ba54c5f8a6f5e Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Wed, 31 Jul 2024 11:05:23 +0000 Subject: [PATCH 30/34] Suggested change to dirichlet_bc option for Gauss Legendre, do not impose bc on lower endpoint of a cylindrical coordinate. --- moment_kinetics/src/gauss_legendre.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/moment_kinetics/src/gauss_legendre.jl b/moment_kinetics/src/gauss_legendre.jl index 539a5fd86..02eaf8b80 100644 --- a/moment_kinetics/src/gauss_legendre.jl +++ b/moment_kinetics/src/gauss_legendre.jl @@ -887,9 +887,12 @@ function setup_global_weak_form_matrix!(QQ_global::Array{mk_float,2}, if dirichlet_bc # Make matrix diagonal for first/last grid points so it does not change the values # there - if coord.irank == 0 - QQ_global[1,:] .= 0.0 - QQ_global[1,1] = 1.0 + if !(coord.name == "vperp") + # modify lower endpoint if not a radial/cylindrical coordinate + if coord.irank == 0 + QQ_global[1,:] .= 0.0 + QQ_global[1,1] = 1.0 + end end if coord.irank == coord.nrank - 1 QQ_global[end,:] .= 0.0 From 85c2ee746669a4953d791b6dd961afbaccf1b4b2 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Wed, 31 Jul 2024 11:34:57 +0000 Subject: [PATCH 31/34] Change names of vperp boundary conditions for sensible behaviour. "zero" now implies that zeros are imposed at the ends of the domain, "zero-impose-regularity" imposes dFdvperp = 0 at vperp = 0. --- ...er-planck-relaxation-slowing-down-alphas.toml | 2 +- ...k-slowing-down-alphas-no-self-collisions.toml | 2 +- ...r-planck-source-sink-slowing-down-alphas.toml | 2 +- moment_kinetics/src/boundary_conditions.jl | 6 +++--- moment_kinetics/src/fokker_planck_calculus.jl | 4 ++-- moment_kinetics/src/moment_kinetics_input.jl | 10 +++++----- .../test/fokker_planck_time_evolution_tests.jl | 16 ++++++++-------- 7 files changed, 21 insertions(+), 21 deletions(-) diff --git a/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas.toml b/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas.toml index dfab272e3..e6b75e29c 100644 --- a/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas.toml +++ b/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas.toml @@ -46,7 +46,7 @@ vperp_ngrid = 5 vperp_nelement = 16 vperp_L = 1.5 vperp_discretization = "gausslegendre_pseudospectral" -vperp_bc = "zero-no-regularity" +vperp_bc = "zero" # Fokker-Planck operator requires the "gausslegendre_pseudospectral # options for the vpa and vperp grids diff --git a/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas-no-self-collisions.toml b/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas-no-self-collisions.toml index e3070a6f7..1a95fc210 100644 --- a/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas-no-self-collisions.toml +++ b/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas-no-self-collisions.toml @@ -46,7 +46,7 @@ vperp_ngrid = 5 vperp_nelement = 16 vperp_L = 1.5 vperp_discretization = "gausslegendre_pseudospectral" -vperp_bc = "zero-no-regularity" +vperp_bc = "zero" # Fokker-Planck operator requires the "gausslegendre_pseudospectral # options for the vpa and vperp grids diff --git a/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas.toml b/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas.toml index 841af3a1e..cd49c82e7 100644 --- a/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas.toml +++ b/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas.toml @@ -46,7 +46,7 @@ vperp_ngrid = 5 vperp_nelement = 16 vperp_L = 1.5 vperp_discretization = "gausslegendre_pseudospectral" -vperp_bc = "zero-no-regularity" +vperp_bc = "zero" # Fokker-Planck operator requires the "gausslegendre_pseudospectral # options for the vpa and vperp grids diff --git a/moment_kinetics/src/boundary_conditions.jl b/moment_kinetics/src/boundary_conditions.jl index 27282d65e..b014abe99 100644 --- a/moment_kinetics/src/boundary_conditions.jl +++ b/moment_kinetics/src/boundary_conditions.jl @@ -974,7 +974,7 @@ function enforce_vperp_boundary_condition!(f::AbstractArray{mk_float,5}, bc, vpe end function enforce_vperp_boundary_condition!(f::AbstractArray{mk_float,4}, bc, vperp, vperp_spectral, vperp_advect, diffusion) - if bc == "zero" || bc == "zero-no-regularity" + if bc == "zero" || bc == "zero-impose-regularity" nvperp = vperp.n ngrid = vperp.ngrid # set zero boundary condition @@ -984,7 +984,7 @@ function enforce_vperp_boundary_condition!(f::AbstractArray{mk_float,4}, bc, vpe end end # set regularity condition d F / d vperp = 0 at vperp = 0 - if bc == "zero" && (vperp.discretization == "gausslegendre_pseudospectral" || vperp.discretization == "chebyshev_pseudospectral") + if bc == "zero-impose-regularity" && (vperp.discretization == "gausslegendre_pseudospectral" || vperp.discretization == "chebyshev_pseudospectral") D0 = vperp_spectral.radau.D0 buffer = @view vperp.scratch[1:ngrid-1] @loop_r_z_vpa ir iz ivpa begin @@ -994,7 +994,7 @@ function enforce_vperp_boundary_condition!(f::AbstractArray{mk_float,4}, bc, vpe f[ivpa,1,iz,ir] = -sum(buffer)/D0[1] end end - elseif bc == "zero-no-regularity" + elseif bc == "zero" # do nothing else println("vperp.bc=\"$bc\" not supported by discretization " diff --git a/moment_kinetics/src/fokker_planck_calculus.jl b/moment_kinetics/src/fokker_planck_calculus.jl index 4007e6356..f11dcf43c 100644 --- a/moment_kinetics/src/fokker_planck_calculus.jl +++ b/moment_kinetics/src/fokker_planck_calculus.jl @@ -2299,12 +2299,12 @@ function enforce_vpavperp_BCs!(pdf,vpa,vperp,vpa_spectral,vperp_spectral) # set regularity condition d F / d vperp = 0 at vperp = 0 # adjust F(vperp = 0) so that d F / d vperp = 0 at vperp = 0 begin_anyv_vpa_region() - if vperp.bc in ("zero", "zero-no-regularity") + if vperp.bc in ("zero", "zero-impose-regularity") @loop_vpa ivpa begin pdf[ivpa,nvperp] = 0.0 end end - if vperp.bc == "zero" + if vperp.bc == "zero-impose-regularity" buffer = @view vperp.scratch[1:ngrid_vperp-1] @loop_vpa ivpa begin @views @. buffer = D0[2:ngrid_vperp] * pdf[ivpa,2:ngrid_vperp] diff --git a/moment_kinetics/src/moment_kinetics_input.jl b/moment_kinetics/src/moment_kinetics_input.jl index a80c03751..3bdba79f8 100644 --- a/moment_kinetics/src/moment_kinetics_input.jl +++ b/moment_kinetics/src/moment_kinetics_input.jl @@ -1196,8 +1196,8 @@ function check_coordinate_input(coord, coord_name, io) println(io,">$coord_name.bc = 'constant'. enforcing constant incoming BC in $coord_name.") elseif coord.bc == "zero" println(io,">$coord_name.bc = 'zero'. enforcing zero incoming BC in $coord_name. Enforcing zero at both boundaries if diffusion operator is present.") - elseif coord.bc == "zero-no-regularity" - println(io,">$coord_name.bc = 'zero'. enforcing zero incoming BC in $coord_name. Enforcing zero at both boundaries if diffusion operator is present. Do not enforce dF/dcoord = 0 at origin if coord = vperp.") + elseif coord.bc == "zero-impose-regularity" + println(io,">$coord_name.bc = 'zero'. enforcing zero incoming BC in $coord_name. Enforcing zero at both boundaries if diffusion operator is present. Enforce dF/dcoord = 0 at origin if coord = vperp.") elseif coord.bc == "zero_gradient" println(io,">$coord_name.bc = 'zero_gradient'. enforcing zero gradients at both limits of $coord_name domain.") elseif coord.bc == "both_zero" @@ -1216,9 +1216,9 @@ function check_coordinate_input(coord, coord_name, io) coord.nelement_global, " elements across the $coord_name domain [", 0.0, ",", coord.L, "].") - if coord.bc != "zero" && coord.n_global > 1 && global_rank[] == 0 - println("WARNING: regularity condition (df/dvperp=0 at vperp=0) not being " - * "imposed. Collisions or vperp-diffusion may be unstable.") + if coord.bc == "zero-impose-regularity" && coord.n_global > 1 && global_rank[] == 0 + println("WARNING: regularity condition (df/dvperp=0 at vperp=0) being " + * "imposed explicitly.") end else println(io,">using ", coord.ngrid, " grid points per $coord_name element on ", diff --git a/moment_kinetics/test/fokker_planck_time_evolution_tests.jl b/moment_kinetics/test/fokker_planck_time_evolution_tests.jl index b5d5086d6..a02fddf70 100644 --- a/moment_kinetics/test/fokker_planck_time_evolution_tests.jl +++ b/moment_kinetics/test/fokker_planck_time_evolution_tests.jl @@ -34,7 +34,7 @@ struct expected_data f_ion::Array{mk_float, 3} # vpa, vperp, time end -const expected_base = +const expected_zero_impose_regularity = expected_data( [-3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0], [0.155051025721682, 0.644948974278318, 1.000000000000000, 1.500000000000000, 2.000000000000000, 2.500000000000000, 3.000000000000000], @@ -81,7 +81,7 @@ const expected_base = 0.005877384425022965 0.00466291282312835 0.0028530945920350282 0.0008130106752860425 2.4399432485774198e-5 -9.743480904212476e-5 0.0; 0.00017108987342944414 7.097261590252536e-5 -7.822316408657793e-5 -0.0001534897546375058 -9.087984332447822e-5 -3.3079379573126077e-5 0.0; 0.0 0.0 0.0 0.0 0.0 0.0 0.0]) -const expected_no_regularity = +const expected_zero = expected_data( [-3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0], [0.155051025721682, 0.644948974278318, 1.000000000000000, 1.500000000000000, 2.000000000000000, 2.500000000000000, 3.000000000000000], @@ -435,17 +435,17 @@ function runtests() # Benchmark data is taken from this run (GaussLegendre) @testset "Gauss Legendre base" begin run_name = "gausslegendre_pseudospectral" - vperp_bc = "zero" + vperp_bc = "zero-impose-regularity" run_test(test_input_gauss_legendre, - expected_base, 1.0e-14, 1.0e-14; - vperp_bc="zero") + expected_zero_impose_regularity, 1.0e-14, 1.0e-14; + vperp_bc=vperp_bc) end @testset "Gauss Legendre no enforced regularity condition at vperp = 0" begin run_name = "gausslegendre_pseudospectral_no_regularity" - vperp_bc = "zero-no-regularity" + vperp_bc = "zero" run_test(test_input_gauss_legendre, - expected_no_regularity, - 1.0e-14, 1.0e-14; vperp_bc="zero-no-regularity") + expected_zero, + 1.0e-14, 1.0e-14; vperp_bc=vperp_bc) end @testset "Gauss Legendre no (explicitly) enforced boundary conditions" begin run_name = "gausslegendre_pseudospectral_none_bc" From 6c35c7999976f3611ef3cb15ebc11219ef0505a6 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Wed, 31 Jul 2024 11:47:32 +0000 Subject: [PATCH 32/34] Change input option in util/precompile_run.jl --- util/precompile_run.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/util/precompile_run.jl b/util/precompile_run.jl index e927edd36..43f293ee5 100644 --- a/util/precompile_run.jl +++ b/util/precompile_run.jl @@ -77,7 +77,7 @@ collisions_input2 = merge(wall_bc_cheb_input, Dict("n_neutral_species" => 0, "fokker_planck_collisions" => Dict{String,Any}("use_fokker_planck" => true, "self_collisions" => true, "slowing_down_test" => true), "vperp_discretization" => "gausslegendre_pseudospectral", "vpa_discretization" => "gausslegendre_pseudospectral", - "vperp_bc" => "zero-no-regularity", + "vperp_bc" => "zero-impose-regularity", )) # add an additional input for every geometry option available in addition to the default geo_input1 = merge(wall_bc_cheb_input, Dict("n_neutral_species" => 0, From 29f630b6d67487616a882f6c2338ddd51c64319c Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Thu, 1 Aug 2024 07:59:15 +0000 Subject: [PATCH 33/34] Change default so that dirichlet_bc=false, and feature permitting ODE solve with K and L matrices is suppressed (but can be revived in a custom script with the dirichlet_bc=true optional argument in setup_gausslegendre_pseudospectral()). --- moment_kinetics/src/coordinates.jl | 3 +-- moment_kinetics/src/gauss_legendre.jl | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/moment_kinetics/src/coordinates.jl b/moment_kinetics/src/coordinates.jl index d034a2dd0..b761ea1fa 100644 --- a/moment_kinetics/src/coordinates.jl +++ b/moment_kinetics/src/coordinates.jl @@ -252,8 +252,7 @@ function define_coordinate(input, parallel_io::Bool=false; run_directory=nothing elseif input.discretization == "gausslegendre_pseudospectral" # create arrays needed for explicit GaussLegendre pseudospectral treatment in this # coordinate and create the matrices for differentiation - spectral = setup_gausslegendre_pseudospectral(coord, collision_operator_dim=collision_operator_dim, - dirichlet_bc=occursin("zero", coord.bc)) + spectral = setup_gausslegendre_pseudospectral(coord, collision_operator_dim=collision_operator_dim) # obtain the local derivatives of the uniform grid with respect to the used grid derivative!(coord.duniform_dgrid, coord.uniform_grid, coord, spectral) else diff --git a/moment_kinetics/src/gauss_legendre.jl b/moment_kinetics/src/gauss_legendre.jl index 02eaf8b80..9a0afb6b3 100644 --- a/moment_kinetics/src/gauss_legendre.jl +++ b/moment_kinetics/src/gauss_legendre.jl @@ -100,7 +100,7 @@ struct gausslegendre_info{TSparse, TLU} <: weak_discretization_info Qmat::Array{mk_float,2} end -function setup_gausslegendre_pseudospectral(coord; collision_operator_dim=true, dirichlet_bc=true) +function setup_gausslegendre_pseudospectral(coord; collision_operator_dim=true, dirichlet_bc=false) lobatto = setup_gausslegendre_pseudospectral_lobatto(coord,collision_operator_dim=collision_operator_dim) radau = setup_gausslegendre_pseudospectral_radau(coord,collision_operator_dim=collision_operator_dim) From 0ed436918e7a99263f7deda014f53d0b0e7a9c54 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Tue, 6 Aug 2024 12:42:59 +0100 Subject: [PATCH 34/34] Comment out mms_tests.jl from debug_test/runtests.jl due to missing Symbolics package. --- moment_kinetics/debug_test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/moment_kinetics/debug_test/runtests.jl b/moment_kinetics/debug_test/runtests.jl index 94b5dbbc7..7317ee974 100644 --- a/moment_kinetics/debug_test/runtests.jl +++ b/moment_kinetics/debug_test/runtests.jl @@ -8,7 +8,7 @@ function runtests() include(joinpath(@__DIR__, "fokker_planck_collisions_tests.jl")) include(joinpath(@__DIR__, "wall_bc_tests.jl")) include(joinpath(@__DIR__, "harrisonthompson.jl")) - include(joinpath(@__DIR__, "mms_tests.jl")) + #include(joinpath(@__DIR__, "mms_tests.jl")) include(joinpath(@__DIR__, "restart_interpolation_tests.jl")) include(joinpath(@__DIR__, "recycling_fraction_tests.jl")) include(joinpath(@__DIR__, "gyroaverage_tests.jl"))