From ec573db2404a8799fcd91e62555045bcd61bc1c9 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Sat, 23 Nov 2024 17:18:40 +0000 Subject: [PATCH] Fix adaptive-timestep IMEX for kinetic electrons Can (?) skip the problematic solve, that was run to update the electron shape function without updating the electron pressure, on the 'explicit first stage' of ESDIRK schemes. This might have the effect of reducing the order of accuracy of the scheme somehow, as the qpar_e used for the 'explicit' calculation of the time derivative of electron_ppar is taken from the most recent implicit solve rather than a solve updated with the new ion/electron profiles. However, the change is probably small (?) and at least the solver does run now - it is useful to have an adaptive ion timestep as it may let the code recover by reducing the ion timestep when an electron implicit solve fails to converge. --- moment_kinetics/src/time_advance.jl | 4 +- .../test/kinetic_electron_tests.jl | 392 +++++++++--------- 2 files changed, 202 insertions(+), 194 deletions(-) diff --git a/moment_kinetics/src/time_advance.jl b/moment_kinetics/src/time_advance.jl index f3592dc80..1ea997da5 100644 --- a/moment_kinetics/src/time_advance.jl +++ b/moment_kinetics/src/time_advance.jl @@ -3104,9 +3104,7 @@ end || (istage == n_rk_stages && t_params.implicit_coefficient_is_zero[1]) || t_params.implicit_coefficient_is_zero[istage+1]) update_electrons = (t_params.rk_coefs_implicit === nothing - || !(t_params.implicit_electron_advance || t_params.implicit_electron_ppar) - || (istage < n_rk_stages && t_params.implicit_coefficient_is_zero[istage+1]) - || (istage == n_rk_stages && t_params.implicit_coefficient_is_zero[1])) + || !(t_params.implicit_electron_advance || t_params.implicit_electron_ppar)) diagnostic_moments = diagnostic_checks && istage == n_rk_stages success = apply_all_bcs_constraints_update_moments!( scratch[istage+1], pdf, moments, fields, boundary_distributions, diff --git a/moment_kinetics/test/kinetic_electron_tests.jl b/moment_kinetics/test/kinetic_electron_tests.jl index 33738748d..d20b8dccd 100644 --- a/moment_kinetics/test/kinetic_electron_tests.jl +++ b/moment_kinetics/test/kinetic_electron_tests.jl @@ -123,9 +123,14 @@ kinetic_input["nonlinear_solver"] = OptionsDict("nonlinear_max_iterations" => 10 "rtol" => 1.0e-8, "atol" => 1.0e-14, "linear_restart" => 5, - "preconditioner_update_interval" => 1000, + "preconditioner_update_interval" => 100, ) +kinetic_input_adaptive_timestep = deepcopy(kinetic_input) +kinetic_input_adaptive_timestep["output"]["run_name"] = "kinetic_electron_adaptive_timestep_test" +kinetic_input_adaptive_timestep["timestepping"]["type"] = "KennedyCarpenterARK324" +kinetic_input_adaptive_timestep["timestepping"]["maximum_dt"] = 1.0e-5 + """ Run a test for a single set of parameters @@ -136,204 +141,209 @@ function run_test() this_boltzmann_input = deepcopy(boltzmann_input) this_boltzmann_input["output"]["base_directory"] = test_output_directory - this_kinetic_input = deepcopy(kinetic_input) - this_kinetic_input["output"]["base_directory"] = test_output_directory - - # Provide some progress info - println(" - testing kinetic electrons") - - # Suppress console output while running? Test is pretty long, so maybe better to leave - # intermediate output visible. Leaving `quietoutput()` commented out for now... + # Suppress console output while running. quietoutput() do run_moment_kinetics(this_boltzmann_input) + end - restart_from_directory = joinpath(this_boltzmann_input["output"]["base_directory"], this_boltzmann_input["output"]["run_name"]) - restart_from_file_pattern = this_boltzmann_input["output"]["run_name"] * ".dfns*.h5" - restart_from_file = glob(restart_from_file_pattern, restart_from_directory)[1] + for (this_kinetic_input, label, tol) ∈ ((deepcopy(kinetic_input), "", 1.0e-6), + (deepcopy(kinetic_input_adaptive_timestep), "adaptive timestep", 1.0e-4)) + # Provide some progress info + println(" - testing kinetic electrons $label") - # run kinetic electron simulation - run_moment_kinetics(this_kinetic_input; restart=restart_from_file) - end + this_kinetic_input["output"]["base_directory"] = test_output_directory - if global_rank[] == 0 - # Load and analyse output - ######################### - - path = joinpath(realpath(this_kinetic_input["output"]["base_directory"]), this_kinetic_input["output"]["run_name"]) - - # open the output file(s) - run_info = get_run_info_no_setup(path, dfns=true) - - # load fields data - Ez = postproc_load_variable(run_info, "Ez")[:,1,:] - vthe = postproc_load_variable(run_info, "electron_thermal_speed")[:,1,:] - electron_advance_linear_iterations = postproc_load_variable(run_info, "electron_advance_linear_iterations")[end] - - close_run_info(run_info) - - # Regression test - # Benchmark data generated in serial on Linux - if global_size[] == 1 - # Serial solves use LU preconditioner - expected_Ez = [-0.5990683230706185 -1.136483186157602; - -0.4944296396481284 -0.9873296990705788; - -0.30889032954504736 -0.6694380824928302; - -0.2064830747303776 -0.4471331690708596; - -0.21232457328748663 -0.423069171542538; - -0.18233875912042674 -0.3586467595624931; - -0.16711429522309232 -0.3018272987758344; - -0.16920776495088916 -0.27814384649305496; - -0.1629417555658927 -0.26124630661090814; - -0.16619150334079993 -0.2572789330163811; - -0.15918194883360942 -0.23720078037362732; - -0.14034706409006803 -0.20520396656341475; - -0.12602184032280567 -0.1827016549071128; - -0.10928716440800472 -0.15808919669899502; - -0.07053969674257217 -0.10137753767917096; - -0.0249577746169536 -0.0358411459260082; - -2.8327303308330514e-15 -2.0803303361189427e-5; - 0.024957774616960776 0.03584490974053962; - 0.07053969674257636 0.1013692898656727; - 0.10928716440799909 0.15807862358546687; - 0.1260218403227975 0.18263049748179466; - 0.1403470640900294 0.20516566362571026; - 0.1591819488336015 0.23711236692241613; - 0.16619150334082114 0.257126146434857; - 0.16294175556587748 0.2609881259705107; - 0.16920776495090983 0.2778978154805798; - 0.1671142952230893 0.3015349192528757; - 0.1823387591204167 0.3585291689672981; - 0.21232457328753865 0.4231179549656996; - 0.20648307473037922 0.44816400221269476; - 0.3088903295450278 0.6716787105435247; - 0.4944296396481271 0.9861165590258743; - 0.5990683230705801 1.1300034111861956] - expected_vthe = [22.64555285302391 22.485481713141688; - 23.763411647653097 23.63281883616836; - 25.26907160117684 25.181703459470448; - 26.17920352818247 26.12461016686916; - 26.514772631426933 26.476018852279974; - 26.798783188585713 26.774387562937218; - 27.202255545479264 27.203662204308202; - 27.50424749120107 27.527732850637264; - 27.630498656270504 27.6642323848215; - 27.748483758260697 27.79134809261204; - 27.933760382468346 27.990808336620802; - 28.08611508251559 28.153978618442775; - 28.14959662643782 28.221734439130564; - 28.207730844115044 28.283677711828023; - 28.28567669896009 28.36634261525836; - 28.32728392065335 28.410489883644782; - 28.331064506972027 28.41437629072209; - 28.32729968986601 28.41050992096321; - 28.285678151542136 28.366352683865195; - 28.207765527709956 28.28373408727703; - 28.149604559462947 28.221771261090687; - 28.086248527111163 28.154158507899695; - 27.933979289064936 27.991103719847732; - 27.74906125092813 27.792046191405188; - 27.631210333523736 27.66508092926101; - 27.505479130159543 27.529115937508752; - 27.20422756527604 27.20578114592589; - 26.801712351383053 26.77740066591359; - 26.517644511297203 26.478915386575462; - 26.18176436913143 26.127099000267552; - 25.26635932097994 25.178676836919877; - 23.756593489029708 23.625697695979085; - 22.64390166090378 22.48400980852866] - else - # Parallel solves, which here use only shared-memory parallelism, use the ADI - # preconditioner, which should be as accurate, but may give different results - # within Newton-Krylov tolerances. - expected_Ez = [-0.5990683230706185 -1.136484793603861; - -0.4944296396481284 -0.9873300031440772; - -0.30889032954504736 -0.6694378168618197; - -0.2064830747303776 -0.447133132132065; - -0.21232457328748663 -0.42306913446372424; - -0.18233875912042674 -0.3586467771727455; - -0.16711429522309232 -0.30182728110160495; - -0.16920776495088916 -0.27814382747995164; - -0.1629417555658927 -0.2612463784138094; - -0.16619150334079993 -0.25727894258000966; - -0.15918194883360942 -0.23720078814350573; - -0.14034706409006803 -0.20520397188041256; - -0.12602184032280567 -0.18270162474892546; - -0.10928716440800472 -0.1580892035790512; - -0.07053969674257217 -0.10137753682381391; - -0.0249577746169536 -0.03584114725793184; - -2.8327303308330514e-15 -2.0802378395589373e-5; - 0.024957774616960776 0.0358449101669449; - 0.07053969674257636 0.10136928934666747; - 0.10928716440799909 0.15807862867071673; - 0.1260218403227975 0.18263047522175488; - 0.1403470640900294 0.20516566756031385; 0.1591819488336015 0.2371123741024713; - 0.16619150334082114 0.2571261543920033; - 0.16294175556587748 0.2609882062708652; - 0.16920776495090983 0.27789779494370415; - 0.1671142952230893 0.30153489797658445; - 0.1823387591204167 0.35852918516786003; - 0.21232457328753865 0.42311789840457864; - 0.20648307473037922 0.44816400062147066; - 0.3088903295450278 0.6716785459169026; - 0.4944296396481271 0.9861167610959626; - 0.5990683230705801 1.1300045383907789] - expected_vthe = [22.64555338227396 22.48548119549829; - 23.76341164436594 23.632819782771243; - 25.26907163394297 25.18170391887767; - 26.179203467285365 26.12461016927763; - 26.514772629327332 26.47601877788725; - 26.79878318858447 26.774387534342114; - 27.20225551034186 27.20366217166485; - 27.504247525601926 27.527732760234755; - 27.630498605068166 27.66423228184859; - 27.748483763235846 27.791348082529804; - 27.933760371994826 27.990808308571204; - 28.08611509938479 28.153978648601132; - 28.149596610550738 28.221734405417436; - 28.207730848524463 28.28367771694209; - 28.28567670146647 28.366342613061416; - 28.32728392764203 28.410489892675102; - 28.331064498175866 28.414376282256146; - 28.327299695349158 28.41050992979778; - 28.285678155424083 28.366352683054103; - 28.207765532359442 28.28373409338897; - 28.149604554344048 28.22177123547944; - 28.086248537316628 28.154158532699547; - 27.933979285563435 27.991103698041254; - 27.749061255285646 27.79204618050744; - 27.63121031067771 27.665080846653012; - 27.505479148983177 27.529115838548574; - 27.204227550854288 27.205781129997607; - 26.801712356957204 26.777400644678224; - 26.517644516966772 26.478915353716097; - 26.181764354679014 26.12709901369174; - 25.266359355820907 25.178677080491074; - 23.756593465755735 23.625698257711747; - 22.64390180335094 22.48400934735562] - end + # Suppress console output while running. + quietoutput() do + restart_from_directory = joinpath(this_boltzmann_input["output"]["base_directory"], this_boltzmann_input["output"]["run_name"]) + restart_from_file_pattern = this_boltzmann_input["output"]["run_name"] * ".dfns*.h5" + restart_from_file = glob(restart_from_file_pattern, restart_from_directory)[1] - if expected_Ez == nothing - # Error: no expected input provided - println("data tested would be: Ez=", Ez) - @test false - else - @test elementwise_isapprox(Ez, expected_Ez, rtol=0.0, atol=2.0e-6) - end - if expected_vthe == nothing - # Error: no expected input provided - println("data tested would be: vthe=", vthe) - @test false - else - @test elementwise_isapprox(vthe, expected_vthe, rtol=1.0e-6, atol=0.0) + # run kinetic electron simulation + run_moment_kinetics(this_kinetic_input; restart=restart_from_file) end - # Iteration counts are fairly inconsistent, but it's good to check that they at - # least don't unexpectedly increase by an order of magnitude. - # Expected iteration count is from a serial run on Linux. - expected_electron_advance_linear_iterations = 49307 - @test electron_advance_linear_iterations < 2 * expected_electron_advance_linear_iterations - if !(electron_advance_linear_iterations < 2 * expected_electron_advance_linear_iterations) - println("electron_advance_linear_iterations=$electron_advance_linear_iterations was greater than twice the expected $expected_electron_advance_linear_iterations.") + if global_rank[] == 0 + # Load and analyse output + ######################### + + path = joinpath(realpath(this_kinetic_input["output"]["base_directory"]), this_kinetic_input["output"]["run_name"]) + + # open the output file(s) + run_info = get_run_info_no_setup(path, dfns=true) + + # load fields data + Ez = postproc_load_variable(run_info, "Ez")[:,1,:] + vthe = postproc_load_variable(run_info, "electron_thermal_speed")[:,1,:] + electron_advance_linear_iterations = postproc_load_variable(run_info, "electron_advance_linear_iterations")[end] + + close_run_info(run_info) + + # Regression test + # Benchmark data generated in serial on Linux + if global_size[] == 1 + # Serial solves use LU preconditioner + expected_Ez = [-0.5990683230706185 -1.3888270203716615; + -0.4944296396481284 -1.0551455974484043; + -0.30889032954504736 -0.6767520362404918; + -0.2064830747303776 -0.485090580118428; + -0.21232457328748663 -0.41493279708476516; + -0.18233875912042674 -0.37176844952064725; + -0.16711429522309232 -0.31570138238223894; + -0.16920776495088916 -0.29068053599780513; + -0.1629417555658927 -0.27188002998776034; + -0.16619150334079993 -0.2687203611465835; + -0.15918194883360942 -0.2458365285082124; + -0.14034706409006803 -0.212578741420612; + -0.12602184032280567 -0.1884105208926891; + -0.10928716440800472 -0.16329672769037926; + -0.07053969674257217 -0.10430785352778167; + -0.0249577746169536 -0.037002796148943216; + -2.8327303308330514e-15 1.057063009693631e-13; + 0.024957774616960776 0.0370027961489565; + 0.07053969674257636 0.10430785352779906; + 0.10928716440799909 0.1632967276903652; + 0.1260218403227975 0.18841052089259683; + 0.1403470640900294 0.2125787414205574; + 0.1591819488336015 0.24583652850820256; + 0.16619150334082114 0.268720361146576; + 0.16294175556587748 0.27188002998764893; + 0.16920776495090983 0.2906805359978753; + 0.1671142952230893 0.3157013823822701; + 0.1823387591204167 0.3717684495205846; + 0.21232457328753865 0.4149327970849557; + 0.20648307473037922 0.48509058011840156; + 0.3088903295450278 0.6767520362404763; + 0.4944296396481271 1.055145597448456; + 0.5990683230705801 1.388827020371195] + expected_vthe = [22.060784019895635 21.9342060834967; + 23.494436916773658 23.38112321216884; + 25.089111111668707 25.013074485041443; + 26.093509589070717 26.045983400236945; + 26.471857106550107 26.436430919559548; + 26.75791751426904 26.737415688470954; + 27.212567708836758 27.21475541254089; + 27.542095657620997 27.56546581494082; + 27.70118448108486 27.73292369826694; + 27.81848812706628 27.85943826644967; + 28.036861843934524 28.09121450556879; + 28.19767723444817 28.262744860955333; + 28.28507878259119 28.353950302304714; + 28.33330668235518 28.405797687881684; + 28.428130279157052 28.505369272106943; + 28.464772163084138 28.544309158365877; + 28.48502646639033 28.564850718833164; + 28.464772163084135 28.544309158365746; + 28.428130279157056 28.50536927210684; + 28.333306682355172 28.40579768788156; + 28.285078782591192 28.353950302304696; + 28.19767723444816 28.262744860955298; + 28.036861843934517 28.091214505568818; + 27.81848812706627 27.859438266449644; + 27.701184481084855 27.73292369826705; + 27.542095657620987 27.565465814940936; + 27.212567708836758 27.214755412540853; + 26.75791751426904 26.737415688470907; + 26.471857106550132 26.436430919559665; + 26.093509589070734 26.04598340023685; + 25.08911111166873 25.013074485041553; + 23.494436916773665 23.3811232121686; + 22.060784019895582 21.934206083496772] + else + # Parallel solves, which here use only shared-memory parallelism, use the ADI + # preconditioner, which should be as accurate, but may give different results + # within Newton-Krylov tolerances. + expected_Ez = [-0.5990683230706185 -1.3888273991885798; + -0.4944296396481284 -1.055145755247293; + -0.30889032954504736 -0.6767521121103302; + -0.2064830747303776 -0.4850905804127858; + -0.21232457328748663 -0.4149327380676377; + -0.18233875912042674 -0.3717684678982668; + -0.16711429522309232 -0.3157013804917729; + -0.16920776495088916 -0.29068053867861765; + -0.1629417555658927 -0.2718800247898976; + -0.16619150334079993 -0.26872036244271424; + -0.15918194883360942 -0.24583652878817755; + -0.14034706409006803 -0.212578741555737; + -0.12602184032280567 -0.18841052477585926; + -0.10928716440800472 -0.1632967261364965; + -0.07053969674257217 -0.10430785526375647; + -0.0249577746169536 -0.03700279305551095; + -2.8327303308330514e-15 2.1437357110269042e-10; + 0.024957774616960776 0.03700279301018719; + 0.07053969674257636 0.10430785521602967; + 0.10928716440799909 0.1632967263726813; + 0.1260218403227975 0.18841052495380603; + 0.1403470640900294 0.21257874155389292; + 0.1591819488336015 0.2458365287482765; + 0.16619150334082114 0.2687203623073682; + 0.16294175556587748 0.2718800250829084; + 0.16920776495090983 0.2906805385932972; + 0.1671142952230893 0.315701380421582; + 0.1823387591204167 0.3717684678638724; + 0.21232457328753865 0.41493273816921206; + 0.20648307473037922 0.48509058034936786; + 0.3088903295450278 0.6767521120678559; + 0.4944296396481271 1.0551457553034789; + 0.5990683230705801 1.3888273974565115] + expected_vthe = [22.060784026809287 21.93420509723251; + 23.494436910748018 23.38112275273686; + 25.089111112557333 25.013074367833635; + 26.09350958918034 26.04598344924481; + 26.471857102756474 26.436430851855782; + 26.757917514391167 26.737415676638154; + 27.2125677086776 27.214755411271334; + 27.542095657834963 27.565465815602288; + 27.70118448110617 27.732923699730097; + 27.81848812739529 27.859438265048865; + 28.036861844197087 28.09121450892992; + 28.197677235241777 28.26274486181854; + 28.285078783322888 28.35395030826896; + 28.333306682806334 28.405797691637765; + 28.42813028005704 28.505369276012765; + 28.46477216719669 28.544309163991375; + 28.48502647392552 28.564850713212525; + 28.46477216632498 28.544309163777076; + 28.428130280444595 28.505369276099227; + 28.333306682325304 28.405797691420656; + 28.285078783465817 28.35395030775343; + 28.19767723446788 28.262744861162105; + 28.036861843447365 28.09121450837069; + 27.818488127101517 27.85943826461028; + 27.701184480925914 27.73292369967636; + 27.54209565717233 27.565465814986407; + 27.212567708728393 27.21475541125295; + 26.757917514100008 26.737415676347563; + 26.471857103367785 26.4364308523701; + 26.093509588815845 26.045983448888492; + 25.089111113033802 25.013074368227194; + 23.494436910466618 23.38112275249835; + 22.060784027938254 21.93420509818273] + end + + if expected_Ez == nothing + # Error: no expected input provided + println("data tested would be: Ez=", Ez) + @test false + else + @test elementwise_isapprox(Ez, expected_Ez, rtol=0.0, atol=2.0*tol) + end + if expected_vthe == nothing + # Error: no expected input provided + println("data tested would be: vthe=", vthe) + @test false + else + @test elementwise_isapprox(vthe, expected_vthe, rtol=tol, atol=0.0) + end + + # Iteration counts are fairly inconsistent, but it's good to check that they at + # least don't unexpectedly increase by an order of magnitude. + # Expected iteration count is from a serial run on Linux. + expected_electron_advance_linear_iterations = 59514 + @test electron_advance_linear_iterations < 2 * expected_electron_advance_linear_iterations + if !(electron_advance_linear_iterations < 2 * expected_electron_advance_linear_iterations) + println("electron_advance_linear_iterations=$electron_advance_linear_iterations was greater than twice the expected $expected_electron_advance_linear_iterations.") + end end end