Snapshot from this documentation build. SHA-256: 57e6c18e19a23820fe55fe22bdd38d2f2e2c1e460d90c9dc349df39f246df227
1 #!/usr/bin/env julia
2 using PMFGLM
3 using Serialization
4 using TOML
5 using Dates
6
7 name = "cpu"
8 solver = "glm"
9 output = nothing
10 usage = "usage: scripts/validate_production.jl --backend=cpu|metal|cuda [--mhd-solver=glm|ct] [--output=DIR]"
11 for arg in ARGS
12 if startswith(arg,"--backend=")
13 global name = split(arg,'=';limit=2)[2]
14 elseif startswith(arg,"--mhd-solver=")
15 global solver = lowercase(split(arg,'=';limit=2)[2])
16 elseif startswith(arg,"--output=")
17 global output = abspath(split(arg,'=';limit=2)[2])
18 else
19 error(usage)
20 end
21 end
22 name in ("cpu","metal","cuda") || error("unsupported backend")
23 solver in ("glm","ct") || error("--mhd-solver must be glm or ct")
24 # Without --output this refreshes the tracked docs/data/production-BACKEND.toml (GLM only).
25 # CT evidence never goes to tracked files implicitly: --output=DIR is required.
26 solver == "ct" && output === nothing && error("--mhd-solver=ct requires --output=DIR (results/... keeps tracked evidence untouched)")
27 root = dirname(@__DIR__)
28 mkpath(joinpath(root,"results"))
29 if output === nothing
30 out = mktempdir(joinpath(root,"results"); prefix="production-$name-",cleanup=false)
31 else
32 ispath(output) && !isempty(readdir(output)) && error("--output must be new or empty: $output")
33 mkpath(output)
34 out = output
35 end
36
37 function checkpoint_arrays(path)
38 open(path) do io
39 magic = deserialize(io)
40 magic == "MHDKernels PMF checkpoint v1" || error("unexpected checkpoint schema")
41 config = deserialize(io)
42 runtime = deserialize(io)
43 arrays = Dict{String,Any}()
44 for _ in 1:deserialize(io)
45 key = deserialize(io)
46 arrays[key] = deserialize(io)
47 end
48 (; config, runtime, arrays)
49 end
50 end
51
52 function launch(config,path;expected="COMPLETED")
53 get!(config.settings,"run",Dict{String,Any}())["mhd_solver"] = solver
54 result = run_simulation(config;output_dir=path)
55 result.status == expected || error("production case has unexpected endpoint: $result")
56 result
57 end
58
59 function streaming_initial_check(path)
60 lines = readlines(joinpath(path,"streaming.tsv"))
61 header = Base.split(lines[1],'\t'); row = Base.split(lines[2],'\t')
62 values = Dict(header .=> row)
63 initial = parse.(Float64,[values["vbc_x_cms"],values["vbc_y_cms"],values["vbc_z_cms"]])
64 maximum(abs,initial./[1e5,0.5e5,0.25e5].-1) <= 2e-6 ||
65 error("streaming IC units or sign mismatch")
66 modes = readlines(joinpath(path,"modes.tsv"))
67 columns = Base.split(modes[1],'\t'); mode = Base.split(modes[2],'\t')
68 k = parse(Float64,mode[findfirst(==("k_comoving_mpc_inv"),columns)])
69 isapprox(k,2pi*1000*0.674/2;rtol=1e-12) || error("physical mode box mismatch")
70 end
71
72 function last_history_value(path,column)
73 lines = readlines(joinpath(path,"history.tsv"))
74 header = Base.split(lines[1],'\t')
75 index = findfirst(==(column),header)
76 index === nothing ? NaN : parse(Float64,Base.split(lines[end],'\t')[index])
77 end
78
79 started = time()
80 smoke = load_config(joinpath(root,"examples","cosmology-smoke.toml"))
81 smoke.settings["run"]["backend"] = name
82 full = launch(smoke,joinpath(out,"cosmology"))
83
84 # Deliberately stop at a cycle budget so the three trajectories have identical
85 # timesteps and final scale factor. Include chemistry and persistent photons.
86 base = load_config(joinpath(root,"examples","streaming-planar-smoke.toml"))
87 base.settings["run"]["backend"] = name
88 # This also exercises physical box != diagnostic coordinate box, coherent
89 # streaming, planar symmetry, and their checkpoint preservation.
90 base.settings["cosmology"]["zend"] = 3190.0
91 base.settings["cosmology"]["max_dln_a"] = 1e-4
92 base.settings["run"]["max_cycles"] = 4
93 base.settings["run"]["checkpoint_final"] = true
94 continuous = launch(base,joinpath(out,"continuous");expected="STOPPED_AT_LIMIT")
95 base.settings["run"]["max_cycles"] = 2
96 split = launch(base,joinpath(out,"resumed");expected="STOPPED_AT_LIMIT")
97 base.settings["run"]["max_cycles"] = 4
98 base.settings["run"]["resume"] = true
99 resumed = launch(base,joinpath(out,"resumed");expected="STOPPED_AT_LIMIT")
100 a = checkpoint_arrays(joinpath(continuous.output_dir,"checkpoint.jls"))
101 b = checkpoint_arrays(joinpath(resumed.output_dir,"checkpoint.jls"))
102 streaming_initial_check(continuous.output_dir)
103 planar_error = maximum(1:9) do i
104 field = reshape(a.arrays["mhd_U_$i"],16,16,16)
105 maximum(abs,field.-circshift(field,(0,0,1)))/max(1,maximum(abs,field))
106 end
107 planar_error <= 5e-6 || error("planar symmetry lost in coupled evolution: $planar_error")
108 keys(a.arrays) == keys(b.arrays) || error("checkpoint field sets differ")
109 a.runtime.cycle == b.runtime.cycle == 4 || error("restart did not reach cycle four")
110 a.runtime.a == b.runtime.a || error("restart changed the cosmological clock")
111 errors = Dict{String,Float64}()
112 for (key,x) in a.arrays
113 y = b.arrays[key]
114 all(isfinite,x) && all(isfinite,y) || error("nonfinite state: $key")
115 size(x) == size(y) || error("checkpoint field sizes differ: $key")
116 wide = eltype(x) <: Complex ? ComplexF64 : Float64
117 errors[key] = maximum(abs, wide.(x) .- wide.(y)) /
118 max(maximum(abs,wide.(x)),1.0)
119 end
120 # A plotting/diagnostic request must not alter the physical trajectory.
121 sparse_config = RunConfig(deepcopy(base.settings))
122 sparse_config.settings["run"]["resume"] = false
123 sparse_config.settings["diagnostics"]["every"] = 0
124 sparse_run = launch(sparse_config,joinpath(out,"sparse-diagnostics");expected="STOPPED_AT_LIMIT")
125 sparse_state = checkpoint_arrays(joinpath(sparse_run.output_dir,"checkpoint.jls"))
126 a.runtime.a == sparse_state.runtime.a || error("diagnostic cadence changed the clock")
127 cadence_error = maximum(keys(a.arrays)) do key
128 x=a.arrays[key]; y=sparse_state.arrays[key]
129 wide=eltype(x) <: Complex ? ComplexF64 : Float64
130 maximum(abs,wide.(x).-wide.(y))/max(maximum(abs,wide.(x)),1.0)
131 end
132 cadence_error <= 5e-6 || error("diagnostic cadence changed evolution: $cadence_error")
133 worst = maximum(values(errors))
134 passed = worst <= 5e-6
135 record = Dict("schema_version"=>1,"backend"=>name,"mhd_solver"=>solver,"status"=>passed ? "PASS" : "FAIL",
136 "created_utc"=>string(now(UTC)),"machine"=>gethostname(),
137 "source_fingerprint_sha256"=>source_fingerprint(),"julia_version"=>string(VERSION),
138 "julia_threads"=>Threads.nthreads(),"wall_seconds"=>time()-started,
139 "scope"=>"full-physics endpoint smoke plus four-cycle planar streaming full-MHD/photon/H-He-chemistry/gravity restart and diagnostic-cadence parity",
140 "cosmology_status"=>full.status,
141 "cosmology_final_divb_centred_dx_over_brms"=>last_history_value(full.output_dir,"divBdx_over_brms"),
142 "cosmology_final_divb_face_dx_over_brms"=>last_history_value(full.output_dir,"divb_face_dx_over_brms"),"planar_max_scaled_error"=>planar_error,
143 "streaming_initial_relative_tolerance"=>2e-6,"diagnostic_cadence_max_scaled_error"=>cadence_error,"restart_fields"=>length(errors),
144 "restart_max_scaled_error"=>worst,"restart_error_tolerance"=>5e-6,
145 "run_directory"=>relpath(out,root),"divb_note"=>solver == "ct" ? "CT: exact face divergence is in each run's history.tsv column divb_face_dx_over_brms" : "GLM: centred divB column divBdx_over_brms","restart_field_errors"=>errors)
146 destination = output === nothing ? joinpath(root,"docs","data","production-$name.toml") :
147 joinpath(output,solver == "glm" ? "production-$name.toml" : "production-$name-$solver.toml")
148 open(io->TOML.print(io,record;sorted=true),destination,"w")
149 TOML.print(stdout,record;sorted=true)
150 passed || error("production restart parity failed")