Snapshot from this documentation build. SHA-256: bc7acb9e295ebeabafcbc99d28cf5326952c8b40374b7b3371d1d5acd68dc5e7
1 #!/usr/bin/env julia
2 # Warm MHD-kernel throughput. --mhd-solver=glm (default) keeps the historical behaviour;
3 # ct measures face-staggered constrained transport (always the :ref integrator); both runs
4 # GLM (default integrator), GLM on :ref (if that differs: the like-for-like integrator for
5 # CT) and CT, and prints one table. Single-solver and combined output default to results/.
6 using PMFGLM
7 using KernelAbstractions
8 using TOML
9 using Dates
10 using FFTW
11 using Printf
12
13 be_name = :cpu
14 n = 32
15 iterations = 10
16 solver = "glm"
17 output = nothing
18 for arg in ARGS
19 key,value = split(arg,'=';limit=2)
20 key == "--backend" ? (global be_name=Symbol(value)) :
21 key == "--n" ? (global n=parse(Int,value)) :
22 key == "--iterations" ? (global iterations=parse(Int,value)) :
23 key == "--mhd-solver" ? (global solver=lowercase(value)) :
24 key == "--output" ? (global output=abspath(value)) : error("unknown option $key")
25 end
26 solver in ("glm","ct","both") || error("--mhd-solver must be glm, ct, or both")
27 if output === nothing
28 name = solver == "glm" ? "benchmark.toml" : solver == "ct" ? "benchmark-ct.toml" : "benchmark-compare.toml"
29 output = joinpath(@__DIR__,"..","results",name)
30 end
31 n >= 8 && n % 4 == 0 || error("benchmark n must be a multiple of four >=8")
32 iterations > 0 || error("iterations must be positive")
33 if be_name == :metal
34 @eval using Metal
35 Metal.functional() || error("Metal hardware is unavailable")
36 elseif be_name == :cuda
37 @eval using CUDA
38 CUDA.functional() || error("CUDA hardware is unavailable")
39 elseif be_name != :cpu
40 error("backend must be cpu, metal, or cuda")
41 end
42 be = backend(be_name)
43 device = be_name == :metal ? string(Metal.device()) :
44 be_name == :cuda ? CUDA.name(CUDA.device()) : first(Sys.cpu_info()).model
45 default_integrator = be_name == :cpu ? :ref : :cube
46
47 function measure(mhd_solver::Symbol,integrator::Symbol)
48 s = allocate_state(be,Float32,(n,n,n);dx=1/n,recon=:ppm,riemann=:hlld,mhd_solver)
49 init!() = (MHDKernels.init_turb_field!(s); mhd_solver === :ct && MHDKernels.ct_faces_from_cells!(s))
50 init!()
51 work = MHDKernels.device_zeros(be,Float32,(n^3,))
52 dt,smax = MHDKernels.compute_dt!(work,s;cfl=0.2)
53 advance() = step!(s,dt;ch=smax,integrator) # ch is the GLM wave speed; CT ignores it
54 advance(); KernelAbstractions.synchronize(be)
55 init!() # Warm compilation must not advance the measured IC.
56 GC.gc()
57 KernelAbstractions.synchronize(be)
58 start = time_ns()
59 for _ in 1:iterations
60 advance()
61 end
62 KernelAbstractions.synchronize(be)
63 elapsed = (time_ns()-start)/1e9
64 host = fields_to_host(s)
65 finite = all(all(isfinite,a) for a in host)
66 positive_density = minimum(host[1]) > 0
67 brms = sqrt(sum(v->sum(abs2,Float64.(host[v])),6:8)/n^3)
68 divb = Float64(max_divb(s))*Float64(s.dx)/brms # centred-difference diagnostic, any solver
69 face = mhd_solver === :ct ? Float64(MHDKernels.max_divb_face(s))*Float64(s.dx)/brms : NaN
70 payload = Dict("schema_version"=>1,"backend"=>string(be_name),"n"=>n,
71 "iterations"=>iterations,"integrator"=>string(integrator),"mhd_solver"=>string(mhd_solver),
72 "wall_seconds"=>elapsed,"mcells_per_second"=>n^3*iterations/elapsed/1e6,
73 "julia_threads"=>Threads.nthreads(),"julia_version"=>string(VERSION),
74 "created_utc"=>string(now(UTC)),"source_fingerprint_sha256"=>source_fingerprint(),
75 "machine"=>gethostname(),"device"=>device,"platform"=>string(Sys.KERNEL,"/",Sys.ARCH),
76 "finite"=>finite,"positive_density"=>positive_density,
77 "divb_centred_dx_over_brms"=>divb,
78 "status"=>finite && positive_density ? "PASS" : "FAIL",
79 "scope"=>"warm MHD kernel only; excludes chemistry, FFTs, gravity, diagnostics, startup")
80 isnan(face) || (payload["divb_face_dx_over_brms"] = face)
81 payload
82 end
83
84 function write_payload(payload,path)
85 mkpath(dirname(path))
86 open(io->TOML.print(io,payload;sorted=true),path,"w")
87 end
88
89 if solver != "both"
90 payload = measure(Symbol(solver),solver == "ct" ? :ref : default_integrator)
91 write_payload(payload,output)
92 TOML.print(stdout,payload;sorted=true)
93 payload["status"] == "PASS" || error("benchmark ended with an invalid MHD state")
94 else
95 cases = [(:glm,default_integrator,"glm")]
96 default_integrator == :cube && push!(cases,(:glm,:ref,"glm-ref"))
97 push!(cases,(:ct,:ref,"ct"))
98 stem,ext = splitext(output)
99 payloads = Dict{String,Any}[]
100 skipped = Tuple[]
101 for (mhd_solver,integrator,label) in cases
102 result = try
103 measure(mhd_solver,integrator)
104 catch err
105 # e.g. a kernel that does not compile on this backend; say so, do not hide it.
106 reason = first(split(sprint(showerror,err),'\n'))
107 println(stderr,"case ",label," NOT_RUN: ",reason)
108 push!(skipped,(label,mhd_solver,integrator,reason))
109 continue
110 end
111 result["case_label"] = label
112 write_payload(result,string(stem,"-",label,ext))
113 push!(payloads,result)
114 end
115 println("| solver | integrator | backend | N | iterations | Mcell/s | divB centred dx/Brms | divB face dx/Brms | status |")
116 println("|---|---|---|---|---|---|---|---|---|")
117 for p in payloads
118 @printf("| %s | %s | %s | %d | %d | %.3f | %.2e | %s | %s |\n",p["mhd_solver"],p["integrator"],
119 p["backend"],p["n"],p["iterations"],p["mcells_per_second"],p["divb_centred_dx_over_brms"],
120 haskey(p,"divb_face_dx_over_brms") ? @sprintf("%.2e",p["divb_face_dx_over_brms"]) : "n/a",p["status"])
121 end
122 for (label,mhd_solver,integrator,reason) in skipped
123 println("| ",mhd_solver," | ",integrator," | ",be_name," | ",n," | ",iterations," | NOT_RUN | n/a | n/a | ",reason," |")
124 end
125 println("Scope: ",first(payloads)["scope"],". GLM rows pass ch=smax; CT is not GLM-limited. Files: ",stem,"-<label>",ext)
126 all(p["status"]=="PASS" for p in payloads) || error("benchmark ended with an invalid MHD state")
127 end