Snapshot from this documentation build. SHA-256: f50ba02d4c21c0c8772f427a90911a9ec2340d47d0f4949bab22ad95b8257e42
1 #!/usr/bin/env julia
2 using PMFGLM
3 using FFTW
4 using TOML
5 using SHA
6
7 length(ARGS) == 2 || error("usage: julia --project=. -t 4 scripts/generate_ic.jl CONFIG.toml OUTPUT_DIR")
8 config = load_config(ARGS[1])
9 out = abspath(ARGS[2])
10 ispath(out) && !isempty(readdir(out)) && error("IC output must be new or empty")
11 mkpath(out)
12 env = run_environment(config;output_dir=out)
13 n = parse(Int,env["MHD_PMF_N"])
14 get(env,"MHD_PMF_INIT","batchelor") == "batchelor" || error("this generator writes Batchelor fields; use tools/reference for matched linear ICs")
15 FFTW.set_num_threads(parse(Int,get(ENV,"PMF_FFTW_THREADS",string(min(4,Threads.nthreads())))))
16 brms = parse(Float64,get(env,"MHD_PMF_BRMS","0.001"))
17 if haskey(env,"MHD_PMF_B0_NG")
18 z = parse(Float64,env["MHD_ZSTART"])
19 h = parse(Float64,get(env,"MHD_H","0.674"))
20 om = parse(Float64,get(env,"MHD_OM","0.315"))
21 ob = parse(Float64,get(env,"MHD_OB",string(0.17037*om)))
22 vunit = parse(Float64,get(env,"MHD_CHEM_VUNIT","8e5"))
23 brms = magnetic_code_amplitude(parse(Float64,env["MHD_PMF_B0_NG"]),z;
24 h,omega_b=ob,velocity_unit_cms=vunit)
25 end
26 kwargs = Dict{Symbol,Any}(:seed=>parse(Int,get(env,"MHD_PMF_SEED","42")),
27 :kcut=>parse(Float64,get(env,"MHD_PMF_KCUT",string(0.625pi*n))),
28 :brms=>brms,
29 :box=>1.0)
30 for (name,key,convert) in (
31 ("MHD_PMF_KMAX",:kmax,x->parse(Float64,x)),
32 ("MHD_PMF_KMIN",:kmin,x->parse(Float64,x)),
33 ("MHD_PMF_MODE_LOCK_N",:mode_lock_n,x->parse(Int,x)),
34 ("MHD_PMF_CURL_REFERENCE_N",:curl_reference_n,x->parse(Int,x)),
35 ("MHD_PMF_CURL_SYMBOL",:curl_symbol,Symbol),
36 ("MHD_PMF_NORMALIZATION_KIND",:normalization_kind,Symbol),
37 ("MHD_PMF_PRESERVE_LOW_KMAX",:preserve_low_kmax,x->parse(Float64,x)),
38 ("MHD_PMF_PRESERVE_LOW_KCUT",:preserve_low_kcut,x->parse(Float64,x)),
39 ("MHD_PMF_MIN_CELLS_PER_WAVELENGTH",:minimum_cells_per_wavelength,x->parse(Float64,x)),
40 ("MHD_PMF_FIXED_MODE_POWER",:fixed_mode_power,x->x in ("1","true")),
41 ("MHD_PMF_MODE_REALIZATION",:mode_realization,Symbol),
42 ("MHD_PMF_REALIZATION_BOX",:realization_box,x->parse(Float64,x)),
43 )
44 haskey(env,name) && (kwargs[key]=convert(env[name]))
45 end
46 Bx,By,Bz,stats = pmf_batchelor_field((n,n,n);kwargs...)
47 hashes = Dict{String,String}()
48 for (name,field) in zip(("bx.f32","by.f32","bz.f32"),(Bx,By,Bz))
49 path = joinpath(out,name)
50 open(io->write(io,Float32.(field)),path,"w")
51 hashes[name] = bytes2hex(sha256(read(path)))
52 end
53 unit = "MHD code magnetic units at the configured initial redshift"
54 metadata = Dict("n"=>n,"unit"=>unit,"order"=>"Julia column-major (x fastest)",
55 "dtype"=>"Float32","endianness"=>ENDIAN_BOM==0x04030201 ? "little" : "big",
56 "source_fingerprint_sha256"=>source_fingerprint(),"sha256"=>hashes,
57 "total_brms"=>Float64(stats.total_brms_measured),
58 "reference_band_brms"=>Float64(stats.reference_band_brms_measured),
59 "normalization_kind"=>stats.normalization_kind,
60 "mhd_solver"=>get(env,"MHD_SOLVER","glm"),
61 "mhd_solver_note"=>"The field is solver-independent. A CT run builds its face field from this cell field on import/start (see CT_INITIAL_FACES); GLM uses the cell field directly.","configuration"=>config.settings)
62 open(io->TOML.print(io,metadata;sorted=true),joinpath(out,"ic.toml"),"w")
63 println("Wrote ",out,"; B units: ",unit,"; total Brms=",stats.total_brms_measured)
64 println("Import at the same start redshift/units with MHD_INITIAL_B_IMPORT_DIR; raw import preserves each Fourier amplitude. Regenerate for a different normalization epoch.")