From 1ff204e15623dabc870b099940d498e96b2fb66c 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 0cd59185d..8f2d7095b 100644 --- a/moment_kinetics/src/time_advance.jl +++ b/moment_kinetics/src/time_advance.jl @@ -3116,9 +3116,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..63c6b6ee2 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.1053138725180998; + -0.4944296396481284 -0.9819332128466166; + -0.30889032954504736 -0.6745656961983237; + -0.2064830747303776 -0.4459531272930669; + -0.21232457328748663 -0.4253218487528007; + -0.18233875912042674 -0.3596054334022437; + -0.16711429522309232 -0.3021381799340685; + -0.16920776495088916 -0.2784335484692499; + -0.1629417555658927 -0.2612551389558109; + -0.16619150334079993 -0.2574841927015592; + -0.15918194883360942 -0.23740132549636406; + -0.14034706409006803 -0.20534503972256973; + -0.12602184032280567 -0.1827098539044343; + -0.10928716440800472 -0.1582133200686042; + -0.07053969674257217 -0.10145491369831482; + -0.0249577746169536 -0.03585934915825971; + -2.8327303308330514e-15 3.742211718942586e-14; + 0.024957774616960776 0.03585934915827381; + 0.07053969674257636 0.10145491369829167; + 0.10928716440799909 0.15821332006862954; + 0.1260218403227975 0.18270985390445083; + 0.1403470640900294 0.20534503972250218; + 0.1591819488336015 0.23740132549634094; + 0.16619150334082114 0.2574841927015898; + 0.16294175556587748 0.261255138955811; + 0.16920776495090983 0.2784335484692798; + 0.1671142952230893 0.3021381799340713; + 0.1823387591204167 0.3596054334022252; + 0.21232457328753865 0.4253218487528467; + 0.20648307473037922 0.44595312729305947; + 0.3088903295450278 0.6745656961983009; + 0.4944296396481271 0.9819332128466268; + 0.5990683230705801 1.1053138725180645] + expected_vthe = [22.654024448490784 22.494016350356883; + 23.744503682730446 23.61361063067715; + 25.26061134578617 25.173128418725682; + 26.177253875120066 26.122412383901523; + 26.510545637302872 26.47158368991228; + 26.798827552847246 26.77429043464489; + 27.202535498354287 27.2038739551587; + 27.506373594650846 27.529813468465488; + 27.631027625644876 27.664719606410365; + 27.750902611036295 27.793759280909274; + 27.935780521313532 27.992775960575692; + 28.089380398280714 28.157198480516957; + 28.15152314377127 28.223553488629253; + 28.211115085781678 28.2870195116558; + 28.28856778918977 28.369130039283018; + 28.330972960680672 28.41411592647979; + 28.33351348538364 28.416680586218863; + 28.330972960680675 28.41411592647976; + 28.288567789189763 28.369130039283064; + 28.211115085781678 28.287019511655785; + 28.15152314377127 28.223553488629236; + 28.089380398280724 28.157198480516957; + 27.93578052131354 27.992775960575713; + 27.750902611036295 27.79375928090935; + 27.63102762564488 27.664719606410383; + 27.506373594650853 27.529813468465495; + 27.202535498354287 27.2038739551587; + 26.79882755284725 26.774290434644872; + 26.510545637302886 26.471583689912283; + 26.177253875120083 26.122412383901523; + 25.26061134578619 25.173128418725696; + 23.744503682730446 23.613610630677236; + 22.65402444849082 22.494016350356937] + 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.1053137071260657; + -0.4944296396481284 -0.9819330928307715; + -0.30889032954504736 -0.6745656725019216; + -0.2064830747303776 -0.44595313784207047; + -0.21232457328748663 -0.425321828548; + -0.18233875912042674 -0.3596054340570364; + -0.16711429522309232 -0.30213818089568956; + -0.16920776495088916 -0.27843354821637; + -0.1629417555658927 -0.2612551385019989; + -0.16619150334079993 -0.2574841930766524; + -0.15918194883360942 -0.23740132557788143; + -0.14034706409006803 -0.20534504018275174; + -0.12602184032280567 -0.18270985430997166; + -0.10928716440800472 -0.1582133189704785; + -0.07053969674257217 -0.101454914566153; + -0.0249577746169536 -0.035859347929368034; + -2.8327303308330514e-15 -4.536628997349189e-9; + 0.024957774616960776 0.035859348624052545; + 0.07053969674257636 0.10145491474282464; + 0.10928716440799909 0.15821331955573922; + 0.1260218403227975 0.18270985667178208; + 0.1403470640900294 0.2053450392202274; + 0.1591819488336015 0.23740132578753803; + 0.16619150334082114 0.25748419283426127; + 0.16294175556587748 0.2612551396310432; + 0.16920776495090983 0.2784335479625835; + 0.1671142952230893 0.3021381809909585; + 0.1823387591204167 0.35960543399747075; + 0.21232457328753865 0.4253218286915096; + 0.20648307473037922 0.44595313782295487; + 0.3088903295450278 0.6745656725300222; + 0.4944296396481271 0.9819330927685747; + 0.5990683230705801 1.1053137082172033] + expected_vthe = [22.654024454479018 22.494016869931663; + 23.74450367962989 23.61361086266046; + 25.260611341892094 25.173128419566062; + 26.17725387357487 26.122412390676395; + 26.510545632956767 26.47158369227529; + 26.7988275507785 26.774290427357606; + 27.20253549703805 27.20387395613098; + 27.506373594719115 27.529813465559865; + 27.63102762567087 27.6647196112545; + 27.75090260968854 27.79375927764987; + 27.935780521822277 27.992775962652605; + 28.08938039775227 28.157198478502867; + 28.151523156278788 28.223553495610926; + 28.211115080270424 28.28701950947455; + 28.288567793141777 28.369130040934596; + 28.330972955353705 28.414115925374524; + 28.333513456094945 28.41668058720323; + 28.330972961606466 28.414115929999316; + 28.288567792143006 28.369130041232697; + 28.211115083430062 28.287019512466056; + 28.15152314952673 28.223553491119628; + 28.089380398299795 28.157198479157458; + 27.93578052229754 27.99277596224337; + 27.750902609816293 27.79375927871885; + 27.631027625671482 27.664719609967122; + 27.50637359506551 27.52981346582775; + 27.20253549697429 27.203873955958308; + 26.798827550864885 26.77429042759387; + 26.510545632587316 26.471583691722795; + 26.177253873758893 26.122412390844207; + 25.26061134158348 25.17312841929966; + 23.7445036798294 23.613610862832093; + 22.654024453873603 22.494016869407307] + 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 = 48716 + @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