From ba1bdcdafb4277eaaaf1a9cc0ee9b1332cd431a2 Mon Sep 17 00:00:00 2001 From: Arnau Quera-Bofarull Date: Thu, 10 Mar 2022 12:06:10 +0000 Subject: [PATCH 1/4] removed old scripts --- README.md | 2 +- scripts/config.yaml | 52 ----- scripts/debug.jl | 25 -- scripts/debug_hpc.jl | 101 -------- scripts/example_run.jl | 122 ---------- scripts/interpolation_tests.jl | 162 ------------- scripts/intro.jl | 3 - scripts/make_models.jl | 19 -- scripts/plots.jl | 54 ----- scripts/radiation_transfer_trials.jl | 334 --------------------------- scripts/re.jl | 83 ------- scripts/refine_quadtree.jl | 25 -- scripts/run.jl | 80 ------- scripts/run_parallel.jl | 31 --- scripts/ryota_meeting.jl | 160 ------------- scripts/scattering.jl | 15 -- scripts/scratch.jl | 0 scripts/test_interpolation.jl | 49 ---- scripts/tests.jl | 58 ----- 19 files changed, 1 insertion(+), 1374 deletions(-) delete mode 100644 scripts/config.yaml delete mode 100644 scripts/debug.jl delete mode 100644 scripts/debug_hpc.jl delete mode 100644 scripts/example_run.jl delete mode 100644 scripts/interpolation_tests.jl delete mode 100644 scripts/intro.jl delete mode 100644 scripts/make_models.jl delete mode 100644 scripts/plots.jl delete mode 100644 scripts/radiation_transfer_trials.jl delete mode 100644 scripts/re.jl delete mode 100644 scripts/refine_quadtree.jl delete mode 100644 scripts/run.jl delete mode 100644 scripts/run_parallel.jl delete mode 100644 scripts/ryota_meeting.jl delete mode 100644 scripts/scattering.jl delete mode 100644 scripts/scratch.jl delete mode 100644 scripts/test_interpolation.jl delete mode 100644 scripts/tests.jl diff --git a/README.md b/README.md index da0399e..73bc6e9 100644 --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@ ![CI](https://github.com/arnauqb/Qwind.jl/workflows/CI/badge.svg) [![codecov](https://codecov.io/gh/arnauqb/Qwind.jl/branch/main/graph/badge.svg?token=KQPtxMDMAm)](https://codecov.io/gh/arnauqb/Qwind.jl) -[![Documentation](https://github.com/arnauqb/Qwind.jl/actions/workflows/docs.yml/badge.svg)](https://arnauqb.github.io/Qwind.jl/) +[![Documentation](https://github.com/arnauqb/Qwind.jl/actions/workflows/docs.yml/badge.svg)](https://arnauqb.github.io/Qwind.jl/dev/) [![TagBot](https://github.com/arnauqb/Qwind.jl/actions/workflows/tagbot.yml/badge.svg)](https://github.com/arnauqb/Qwind.jl/actions/workflows/tagbot.yml) [![DOI](https://zenodo.org/badge/271830057.svg)](https://zenodo.org/badge/latestdoi/271830057) diff --git a/scripts/config.yaml b/scripts/config.yaml deleted file mode 100644 index 3da12d4..0000000 --- a/scripts/config.yaml +++ /dev/null @@ -1,52 +0,0 @@ -system: - qwind_path: /cosma7/data/dp004/dc-quer1/Qwind.jl - n_cpus: 112 - partition: cosma7-shm - account: dp004 - max_time: 72:00:00 - job_name: qwind - -black_hole: - M: 1e8 # solar masses - mdot: 0.5 # Eddington units - spin: 0.0 # spin parameter (-1,1) - -radiation: - mode: QsosedRadiation # RERadiation: Risaliti & Elvis, QsosedRadiation: Done - relativistic: true - f_x: 0.15 - n_r: 1000 - z_xray: 0.0 - -radiative_transfer: - mode: RegularGrid - density_interpolator: - type: WindInterpolator - vacuum_density: 1e2 - n_timesteps: 10000 - nz: 100 #"@vary grid 50,100,500,1000" - nr: "auto" - -grid: - # units of Rg - r_min: 6.0 - r_max: 10000.0 - z_min: 0.0 - z_max: 10000.0 - -initial_conditions: - mode: CAKIC #UniformIC - r_in: 50.0 - r_fi: 1300.0 - n_lines: "auto" #"@vary grid 50,100,500,1000,5000,10000" #2000 - log_spaced: true - z_0: 0.0 # Rg - K: 0.03 - alpha: 0.6 - mu: 0.5 - -integrator: - n_iterations: 5 - atol: 1e-8 - rtol: 1e-3 - save_path: "papers/paper2/data" diff --git a/scripts/debug.jl b/scripts/debug.jl deleted file mode 100644 index 9257cc0..0000000 --- a/scripts/debug.jl +++ /dev/null @@ -1,25 +0,0 @@ -using Distributed -@everywhere using DrWatson -@everywhere @quickactivate "Qwind" -@everywhere using Qwind - -using PyPlot -#using YAML, HDF5, CSV, DataFrames -include("scripts/plotting.jl") -#using Profile, PProf, TimerOutputs, BenchmarkTools, ProgressMeter -LogNorm = matplotlib.colors.LogNorm -Normalize = matplotlib.colors.Normalize - -function get_model(config) - model = Model(config) - try - #mv(model.parameters.save_path, "/cosma7/data/dp004/dc-quer1/backup", force = true) - mv(model.parameters.save_path, "./backup", force = true) - catch - end - model = Model(config) - iterations_dict = Dict() - return model, iterations_dict -end -model, iterations_dict = get_model("./to_debug.yaml"); -run!(model, iterations_dict) diff --git a/scripts/debug_hpc.jl b/scripts/debug_hpc.jl deleted file mode 100644 index 552402b..0000000 --- a/scripts/debug_hpc.jl +++ /dev/null @@ -1,101 +0,0 @@ -ENV["JULIA_WORKER_TIMEOUT"] = 250 -using Distributed, ClusterManagers -pids = addprocs_slurm( - 100, - topology = :master_worker, - p = "cosma7-shm", - A = "dp004", - t = "04:00:00", - job_file_loc = "cpu_logs", -) - -@everywhere using Pkg -@everywhere Pkg.activate("/cosma/home/dp004/dc-quer1/Qwind.jl") -#@everywhere pushfirst!(Base.DEPOT_PATH, "/tmp/julia.cache") -println("Running on $(nprocs()) cores.") -@everywhere using LinearAlgebra, YAML, CSV, DataFrames -@everywhere BLAS.set_num_threads(1) -@info "Compiling Qwind..." -flush(stdout) -flush(stderr) -@everywhere using Qwind, YAML -@info "Done" -flush(stdout) -flush(stderr) - -using ProgressMeter, PyPlot -LogNorm = matplotlib.colors.LogNorm - -function get_model(config) - model = Model(config) - try - mv(model.config[:integrator][:save_path], "backup", force = true) - catch - end - model = Model(config) - iterations_dict = Dict() - return model, iterations_dict -end -model, iterations_dict = get_model("./configs/debug.yaml"); -run!(model, iterations_dict, parallel = true); - -@everywhere function compute_xi(rad, r, z; include_scattering = true) - n = get_density(rad.wi.density_grid, r, z) - tau_x = compute_tau_xray(rad, ri = 0, phii = 0, zi = 0, rf = r, zf = z, phif = 0) - ret = compute_ionization_parameter( - r = r, - z = z, - vr = 0, - vz = 0, - number_density = n, - tau_x = tau_x, - xray_luminosity = rad.xray_luminosity, - Rg = rad.bh.Rg, - include_scattering = include_scattering, - density_grid = rad.wi.density_grid, - absorption_opacity = rad.xray_opacity, - zh = rad.z_xray, - mu_electron = rad.mu_electron, - mu_nucleon = rad.mu_nucleon, - scattered_luminosity_grid = rad.wi.scattered_lumin_grid, - ) - return ret -end - -@everywhere rr = collect(10 .^ range(log10(6), log10(1000), length = 200)) -@everywhere zz = 10 .^ range(log10(1e-3), log10(1000), length = 201) -@everywhere rets = zeros(length(rr), length(zz)) -@showprogress for (i, r) in enumerate(rr) - rets[i, :] = pmap( - z -> compute_xi(model.rad, r, z, include_scattering = true), - zz, - batch_size = 10, - ) -end - - -fig, ax = plt.subplots() -cm = ax.pcolormesh(rr, zz, rets', norm = LogNorm(vmin=1e-2, vmax=1e7), shading="auto") -cbar = plt.colorbar(cm, ax = ax) -ax.set_xlabel("R [Rg]") -ax.set_ylabel("R [Rg]") -cbar.set_label("Ionization parameter", rotation=270) -#fig.savefig("./tests_noscatt.pdf") -fig.savefig("./tests.pdf") - -fig, ax = plt.subplots() -cm = ax.pcolormesh(model.rad.wi.density_grid.r_range, model.rad.wi.density_grid.z_range, model.rad.wi.density_grid.grid', norm=LogNorm()) -fig.savefig("./density_tests.pdf") - -@everywhere function f(i, model) - ret = zeros(length(zz)) - r = rr[i] - for (j, z) in enumerate(zz) - end - return ret -end -xi_scatt = @showprogress pmap(i -> f(i, model), 1:length(rr)) -xi_scatt = hcat(xi_scatt)'; - -#iterations_dict[1] = Dict() -#integrators, streamlines = Qwind.run_integrators!(model, iterations_dict, it_num = 1, parallel = true); diff --git a/scripts/example_run.jl b/scripts/example_run.jl deleted file mode 100644 index f211fe7..0000000 --- a/scripts/example_run.jl +++ /dev/null @@ -1,122 +0,0 @@ -## -using DrWatson -@quickactivate "Qwind" -using Qwind - -using PyPlot -LogNorm = matplotlib.colors.LogNorm -using RegionTrees - - -# create black hole -## -M = 1e8 * M_SUN -mdot = 0.5 -spin = 0.0 -black_hole = BlackHole(M, mdot, spin) - -# create radiation model -## -#fuv = 0.85 -fx = 0.15 -#radiation = RERadiation(black_hole, fuv, fx) -radiation = QsosedRadiation(black_hole, 300, fx) - -# radiative transfer model -shielding_density = 2e8 -rin = 50.0 -radiative_transfer = RERadiativeTransfer(radiation, shielding_density, rin) - -# create simulation grid -r_min = 0.0 -r_max = 10000 -z_min = 1.0 -z_max = 10000 -sim_grid = Grid(r_min, r_max, z_min, z_max) - -# initial conditions -z0 = 1.0 -n0 = 2e8 -v0 = 5e7 / C -rfi = 1600 -nlines = 50 -logspaced = true -initial_conditions = UniformIC(rin, rfi, nlines, z0, n0, v0, logspaced) - - -## -integrators = initialize_integrators( - radiative_transfer, - sim_grid, - initial_conditions, - atol = 1e-8, - rtol = 1e-3, -) - -iterations = Dict() -iterations[0] = integrators - -iterations_dict = Dict(0 => integrators) -## -run_integrators!(integrators) - -# plot lines -#fig, ax = plot_streamlines(integrators, xh = 10000, yh = 2000) -#gcf() - -# iteration 2 -radiation = QsosedRadiation(black_hole, 300, fx) - -tau_max_std = 0.01 - -adaptive_mesh = AdaptiveMesh(radiation, integrators, tau_max_std = tau_max_std) - -## - -initial_conditions = - CAKIC(radiation, black_hole, rin, rfi, nlines, z0, 0.03, 0.6, 0.5, logspaced) - -integrators = initialize_integrators( - adaptive_mesh, - sim_grid, - initial_conditions, - atol = 1e-8, - rtol = 1e-3, -) - -iterations[1] = integrators - - -run_integrators!(integrators) - -adaptive_mesh = AdaptiveMesh(radiation, integrators, tau_max_std = tau_max_std) - -integrators = initialize_integrators( - adaptive_mesh, - sim_grid, - initial_conditions, - atol = 1e-8, - rtol = 1e-3, -) - -iterations[2] = integrators - -run_integrators!(integrators) - -adaptive_mesh = AdaptiveMesh(radiation, integrators, tau_max_std = tau_max_std) - -# iteration 3 - -#end -# start iterations -#n_iterations = 3 -# -#for i in 1:n_iterations -# global adaptive_mesh -# global iterations -# global tau_max_std -# integrators = initialize_integrators(adaptive_mesh, sim_grid, initial_conditions, atol=1e-8, rtol=1e-3) -# iterations[i] = integrators -# run_integrators!(integrators) -# adaptive_mesh = AdaptiveMesh(radiation, integrators, tau_max_std=tau_max_std) -#end diff --git a/scripts/interpolation_tests.jl b/scripts/interpolation_tests.jl deleted file mode 100644 index af9e131..0000000 --- a/scripts/interpolation_tests.jl +++ /dev/null @@ -1,162 +0,0 @@ -using Qwind -using Test -using PyPlot -LogNorm = matplotlib.colors.LogNorm - -# create black hole -M = 1e8 * M_SUN -mdot = 0.5 -spin = 0.0 -black_hole = BlackHole(M, mdot, spin) - -# create radiation -f_uv = 0.85 -f_x = 0.15 -shielding_density = 2e8 -r_in = 20.0 -radiation = SimpleRadiation(black_hole, f_uv, f_x, shielding_density, r_in) - -# create simulation grid -r_min = 0.0 -r_max = 10000 -z_min = 0.0 -z_max = 10000 -vacuum_density = 1e2 -grid = Grid(r_min, r_max, z_min, z_max) - -# initialize streamlines -r_fi = 1600 -n_lines = 80 -velocity_generator(r_0) = 5e7 -density_generator(r_0) = 2e8 -streamlines = Streamlines( - r_in, - r_fi, - n_lines, - velocity_generator, - density_generator, - black_hole.M; - z_0 = 1e-2, - log_spaced = true, -) - -# initialize quadtree - -quadtree = initialize_quadtree(grid) - -#run simulation - -solvers = [] -for line in streamlines.lines - parameters = Parameters(line, grid, radiation, quadtree) - solver = initialize_solver(line, parameters) - push!(solvers, solver) - run_solver!(solver) -end - -function gaussian_smoothing(x, mu, sigma) - return 1 / (sigma * 2π) * exp(-0.5 * ((x - mu) / sigma)^2) -end -function my_distance(dp, d1, d2) - ret = d2^2 - ((d1^2 - dp^2 - d2^2) / (2dp))^2 - return sqrt(ret) -end - -function initialize_kdtree(solvers) - solver = solvers[1] - t_range = 10 .^ range(-3, log10(solver.t), length = 10000) - dense_solution = solver.sol(t_range) - r_dense_array = dense_solution[1, :] - z_dense_array = dense_solution[2, :] - v_r_dense_array = dense_solution[3, :] - v_z_dense_array = dense_solution[4, :] - z_max_array = maximum(z_dense_array) .* ones((length(z_dense_array))) - widths_array = r_dense_array .* solver.p.line.width_norm - if length(solvers) >= 2 - for solver in solvers[2:end] - t_range = - 10 .^ range(-3, log10(solver.sol.t[end-2]), length = 10000) - dense_solution = solver.sol(t_range) - r_dense_array = vcat(r_dense_array, dense_solution[1, :]) - z_dense_array = vcat(z_dense_array, dense_solution[2, :]) - v_r_dense_array = vcat(v_r_dense_array, dense_solution[3, :]) - v_z_dense_array = vcat(v_z_dense_array, dense_solution[4, :]) - widths_array = vcat( - widths_array, - dense_solution[1, :] .* solver.p.line.width_norm, - ) - z_max_array = vcat( - z_max_array, - maximum(dense_solution[2, :]) .* - ones((length(dense_solution[2, :]))), - ) - end - end - v_t = @. sqrt(v_r_dense_array^2 + v_z_dense_array^2) - points = hcat(r_dense_array, z_dense_array)' - density_dense_array = - compute_density.(r_dense_array, z_dense_array, v_t, Ref(solver.p.line)) - kdtree = KDTree(points) - return kdtree, - density_dense_array, - widths_array, - z_max_array, - r_dense_array, - z_dense_array -end - -function get_density_from_grid( - kdtree, - r, - z, - densities, - widths, - r_dense_array, - z_dense_array, -) - idcs, distances = knn(kdtree, [r, z], 2) - p1 = [r_dense_array[idcs[1]], z_dense_array[idcs[1]]] - d1 = distances[1] - p2 = [r_dense_array[idcs[2]], z_dense_array[idcs[2]]] - d2 = distances[2] - dp = evaluate(Euclidean(), p1, p2) - distance = my_distance(dp, d1, d2) - width = sum(widths[idcs]) / 4 - density = sum(densities[idcs]) / 2 - if distance <= width - ret = density #max(density, density_grid[i,j]) - else - ret = max(density * exp(-(distance - width)^2 / 2), 1e2) - end - return ret -end - -kdtree, densities, widths, z_array, r_dense_array, z_dense_array = - initialize_kdtree(solvers) -MIN_R = 0 -MIN_Z = 0 -MAX_R = 1000 -MAX_Z = 1000#600#0.21 -r_range = range(MIN_R, MAX_R, length = convert(Int, 500)) -z_range = range(MIN_Z, MAX_Z, length = convert(Int, 501)) -density_grid = 1e2 .* ones((length(r_range), length(z_range))) -for (i, r) in enumerate(r_range) - for (j, z) in enumerate(z_range) - density_grid[i,j] = get_density_from_grid(kdtree, r, z, densities, widths, r_dense_array, z_dense_array) - end -end -fig, ax = plt.subplots() -cm = ax.pcolormesh(r_range, z_range, density_grid', norm = LogNorm()) -plt.colorbar(cm, ax = ax) -#ax.plot(r_dense_array, z_dense_array, "o", markersize=0.1) -plot_streamlines(streamlines.lines, fig, ax) -#for i in 1:length(streamlines) -# idx = (i-1)*1000+1:i*1000 -# ax.plot(r_dense_array[idx], z_dense_array[idx], "o", markersize=1) -#end -ax.set_xlim(MIN_R, MAX_R) -ax.set_ylim(MIN_Z, MAX_Z) -#ax.set_yscale("log") -gcf() -#xlims!(0,1000) -#ylims!(0,100) diff --git a/scripts/intro.jl b/scripts/intro.jl deleted file mode 100644 index b4e5bd1..0000000 --- a/scripts/intro.jl +++ /dev/null @@ -1,3 +0,0 @@ -using DrWatson -@quickactivate "Qwind" -DrWatson.greet() diff --git a/scripts/make_models.jl b/scripts/make_models.jl deleted file mode 100644 index 67cd7aa..0000000 --- a/scripts/make_models.jl +++ /dev/null @@ -1,19 +0,0 @@ -using DrWatson -#@quickactivate "/cosma7/data/dp004/dc-quer1/Qwind.jl" -quickactivate("/cosma7/data/dp004/dc-quer1/Qwind.jl") -using Qwind -using ArgParse - -s= ArgParseSettings() -@add_arg_table s begin - "--config", "-c" - help = "Config path" - arg_type = String - required = true -end -parsed_args = parse_args(ARGS, s) - -config = parsed_args["config"] - -create_models_folders(config) - diff --git a/scripts/plots.jl b/scripts/plots.jl deleted file mode 100644 index f553c1d..0000000 --- a/scripts/plots.jl +++ /dev/null @@ -1,54 +0,0 @@ -using Qwind -using PyPlot -LogNorm = matplotlib.colors.LogNorm - -it_num = 2 -rt = iterations_dict[it_num]["radiative_transfer"] -sl = iterations_dict[it_num]["integrators"] - -fig, ax = plt.subplots() -plot_streamlines(sl, fig, ax) - -r_range = 10 .^ range(-1, 3, length=50) -z_range = 10 .^ range(-6, 3, length=50) - -rfield = zeros((length(r_range), length(z_range), 2)) -for (i, r) in enumerate(r_range) - for (j, z) in enumerate(z_range) - rfield[i, j, :] .= compute_disc_radiation_field(rt, r, z) - end -end - -force = compute_disc_radiation_field(rt, 6.250551925273973, 0.49417133613238345, atol=0, rtol=1e-3) - -tau_uv = compute_uv_tau(rt, 28.96510171628526, 6.250551925273973, 0.49417133613238345) - -fig, ax = plt.subplots() -cm = ax.pcolormesh(r_range, z_range, rfield[:,:,1]', norm=LogNorm()) -plt.colorbar(cm, ax=ax) - -fig, ax = plt.subplots() -cm = ax.pcolormesh(r_range, z_range, rfield[:,:,2]', norm=LogNorm(1e-10, 1e-5)) -plt.colorbar(cm, ax=ax) - -r_range = range(0, 200, length=250) -z_range = range(0, 100, length=250) -density_field = zeros((length(r_range), length(z_range))) -for (i, r) in enumerate(r_range) - for (j, z) in enumerate(z_range) - density_field[i, j] = get_density(rt.density_interpolator, r, z) - end -end -fig, ax = plt.subplots() -cm = ax.pcolormesh(r_range, z_range, density_field', norm=LogNorm()) -#plot_streamlines(sl, fig, ax) -plt.colorbar(cm, ax=ax) -ax.set_xlim(r_range[1], r_range[end]) -ax.set_ylim(z_range[1], z_range[end]) - -fig, ax = plt.subplots() -cm = ax.pcolormesh(rt.density_interpolator.r_range, rt.density_interpolator.z_range, rt.density_interpolator.grid', norm=LogNorm()) -#plot_streamlines(sl, fig, ax) -plt.colorbar(cm, ax=ax) -ax.set_xlim(0, 200) -ax.set_ylim(0, 1) diff --git a/scripts/radiation_transfer_trials.jl b/scripts/radiation_transfer_trials.jl deleted file mode 100644 index 0f184b9..0000000 --- a/scripts/radiation_transfer_trials.jl +++ /dev/null @@ -1,334 +0,0 @@ -using Qwind -using Test -using PyPlot -using RegionTrees -using StaticArrays -using Roots -import RegionTrees: AbstractRefinery, needs_refinement, refine_data -LogNorm = matplotlib.colors.LogNorm - -# create black hole -M = 1e8 * M_SUN -mdot = 0.5 -spin = 0.0 -black_hole = BlackHole(M, mdot, spin) - -# create radiation -f_uv = 0.85 -f_x = 0.15 -shielding_density = 2e8 -r_in = 20.0 -radiation = SimpleRadiation(black_hole, f_uv, f_x, shielding_density, r_in) - -# create simulation grid -r_min = 0.0 -r_max = 10000 -z_min = 0.0 -z_max = 10000 -vacuum_density = 1e2 -grid = Grid(r_min, r_max, z_min, z_max) - -# initialize streamlines -r_fi = 1600 -n_lines = 80 -velocity_generator(r_0) = 5e7 -density_generator(r_0) = 2e8 -streamlines = Streamlines( - r_in, - r_fi, - n_lines, - velocity_generator, - density_generator, - black_hole.M; - z_0 = 1e-2, - log_spaced = true, -) - -# initialize quadtree - -quadtree = initialize_quadtree(grid) - -#run simulation - -solvers = [] -for line in streamlines.lines - parameters = Parameters(line, grid, radiation, quadtree) - solver = initialize_solver(line, parameters) - push!(solvers, solver) - run_solver!(solver) -end - -windtree = initialize_wind_tree(solvers) - -MIN_R = 0 -MIN_Z = 0 -MAX_R = 1000 -MAX_Z = 1000#600#0.21 -r_range = range(MIN_R, MAX_R, length = convert(Int, 500)) -z_range = range(MIN_Z, MAX_Z, length = convert(Int, 501)) -density_grid = 1e2 .* ones((length(r_range), length(z_range))) -xi_grid = 1e2 .* ones((length(r_range), length(z_range))) -for (i, r) in enumerate(r_range) - for (j, z) in enumerate(z_range) - #xi_grid[i, j] = compute_ionization_parameter( - # windtree, - # r, - # z, - # xray_luminosity(radiation), - # black_hole.Rg, - # n_points=500, - #) - density_grid[i, j] = get_density_from_tree(windtree, r, z) - end -end -fig, ax = plt.subplots() -cm = ax.pcolormesh(r_range, z_range, density_grid', norm = LogNorm()) -#cm = ax.pcolormesh(r_range, z_range, xi_grid', norm = LogNorm()) -plt.colorbar(cm, ax = ax) -plot_streamlines(streamlines.lines, fig, ax) -ax.set_xlim(MIN_R, MAX_R) -ax.set_ylim(MIN_Z, MAX_Z) -#ax.set_yscale("log") -gcf() - - - - -function getdata(cell, child_indices, Rg, windtree) - r, z = quadtree.boundary.origin + quadtree.boundary.widths / 2 - density = get_density_from_tree(windtree, r, z) - return density -end - -getdata(cell, child_indices) = - getdata(cell, child_indices, black_hole.Rg, windtree) - -struct MyRefinery <: AbstractRefinery - delta_tau_max::Float64 - windtree::Any - Rg::Any -end - -function needs_refinement(refinery::MyRefinery, cell) - r, z = cell.boundary.origin + cell.boundary.widths / 2 - rs = range( - cell.boundary.origin[1], - cell.boundary.origin[1] + cell.boundary.widths[1], - length = 10, - ) - zs = range( - cell.boundary.origin[2], - cell.boundary.origin[2] + cell.boundary.widths[2], - length = 10, - ) - densities = [] - for rp in rs - for zp in zs - push!(densities, get_density_from_tree(refinery.windtree, rp, zp)) - end - end - delta_d = sqrt(sum(cell.boundary.widths .^ 2)) - taus = densities .* delta_d * refinery.Rg * SIGMA_T - #densities = get_density_from_tree.(Ref(refinery.windtree), rs, zs) - println("r : $r z : $z") - println("std: $(std(taus))") - refine_condition = std(taus) > 0.01 || delta_d > 100 - if !refine_condition - cell.data = mean(densities) - end - refine_condition - - #density = maximum(densities) - - #density = get_density_from_tree(refinery.windtree, r, z) - #tau = delta_d * refinery.Rg * SIGMA_T * density - #tau > refinery.delta_tau_max# || delta_d > 50 -end - -function refine_data(refinery::MyRefinery, cell::Cell, indices) - rs = range( - cell.boundary.origin[1], - cell.boundary.origin[1] + cell.boundary.widths[1], - length = 10, - ) - zs = range( - cell.boundary.origin[2], - cell.boundary.origin[2] + cell.boundary.widths[2], - length = 10, - ) - densities = [] - for rp in rs - for zp in zs - push!(densities, get_density_from_tree(refinery.windtree, rp, zp)) - end - end - mean(densities) -end - -quadtree = Cell(SVector(0.0, 0.0), SVector(1000.0, 1000.0), 1e2) - -refinery = MyRefinery(1, windtree, black_hole.Rg) - -adaptivesampling!(quadtree, refinery) - - -function compute_xray_tau_leaf2( - quadtree::Cell, - leaf::Cell, - point, - intersection, - taux0, - xray_luminosity, - Rg, -) - deltad = evaluate(Euclidean(), point, intersection) * Rg # cell size - d = evaluate(Euclidean(), [0.0, 0.0], intersection) * Rg # distance from the center - #density = interpolate_density(leaf.data.line_id, point, leaf, wind) #leaf.data[1] / leaf.data[2] - #deltatau = compute_tau_cell(quadtree, leaf, point, intersection, 1e2) * Rg - #n_eff = quadtree_effective_density(quadtree, point, intersection) - density = leaf.data - deltatau = density * deltad - xi0 = xray_luminosity / (density * d^2) - f(t) = - t - log(xi0) + taux0 + min(40, deltatau * compute_xray_opacity(exp(t))) - if f(20) < 0 - xi = xi0 - elseif f(-20) > 0 - xi = 1e-20 - else - t = find_zero(f, (-20, 20), Bisection(), atol = 0, rtol = 0.1) - xi = exp(t) - end - #xi = 1e6 - taux = compute_xray_opacity(xi) * deltatau - return taux -end - -function compute_xray_tau2(quadtree::Cell, r, z, z_0, xray_luminosity, Rg) - #return 40 - z = max(z, get_zmin(quadtree)) - point1 = [0.0, z_0] - #coords_list = [point1] - point1leaf = findleaf(quadtree, point1) - point2 = [r, z] - point2leaf = findleaf(quadtree, point2) - taux = 0.0 - if point1leaf == point2leaf - #push!(coords_list, point2) - taux = compute_xray_tau_leaf2( - quadtree, - point1leaf, - point1, - point2, - 0.0, - xray_luminosity, - Rg, - ) - return taux - end - intersection = compute_cell_intersection2(point1leaf, point1, point1, point2) - taux = compute_xray_tau_leaf2( - quadtree, - point1leaf, - point1, - intersection, - taux, - xray_luminosity, - Rg, - ) - currentpoint = intersection - #push!(coords_list, intersection) - currentleaf = findleaf(quadtree, currentpoint) - previousleaf = copy(currentleaf) - while currentleaf != point2leaf - intersection = - compute_cell_intersection2(currentleaf, currentpoint, point1, point2) - #push!(coords_list, intersection) - taux += compute_xray_tau_leaf2( - quadtree, - currentleaf, - currentpoint, - intersection, - taux, - xray_luminosity, - Rg, - ) - #if taux > 40 - # return 40.0 - #end - currentpoint = intersection - currentleaf = findleaf(quadtree, currentpoint) - if currentleaf == previousleaf - break - end - previousleaf = copy(currentleaf) - end - if currentpoint[2] < point2[2] - taux += compute_xray_tau_leaf2( - quadtree, - currentleaf, - currentpoint, - point2, - taux, - xray_luminosity, - Rg, - ) - end - #push!(coords_list, point2) - return taux - #if return_coords - # return taux, coords_list - #else - # return taux - #end -end - -function compute_cell_intersection2( - cell::Cell, - currentpoint, - initialpoint, - finalpoint, -) - r, z = currentpoint - ri, zi = initialpoint - rf, zf = finalpoint - direction_r = sign(rf - ri) - direction_z = sign(zf - zi) - cell_vertices = vertices(cell) - direction_r == 1 ? cell_r_boundary = cell_vertices[2, 2][1] : - cell_r_boundary = cell_vertices[1, 1][1] - direction_z == 1 ? cell_z_boundary = cell_vertices[2, 2][2] : - cell_z_boundary = cell_vertices[1, 1][2] - ri == rf ? lambda_r = Inf : lambda_r = (cell_r_boundary - r) / (rf - ri) - zi == zf ? lambda_z = Inf : lambda_z = (cell_z_boundary - z) / (zf - zi) - lambda = min(lambda_r, lambda_z) - intersection = [r, z] + lambda .* [rf - ri, zf - zi] - try - @assert lambda >= 0.0 - catch - println("lambda: $lambda") - println("current point : $currentpoint") - println("cell: $cell") - println("initial point: $initialpoint") - println("final point: $finalpoint") - throw(DomainError) - end - return intersection -end - - -r_range = range(0, 1000, length = 100) -z_range = range(0, 1000, length = 100) - -tau_x_grid = zeros(length(r_range), length(z_range)) - -for (i, r) in enumerate(r_range) - for (j, z) in enumerate(z_range) - tau_x_grid[i,j] = compute_xray_tau2(quadtree, r, z, 0, xlumin, rg) - end -end - -fig, ax = plt.subplots() -cm = ax.pcolormesh(r_range, z_range, tau_x_grid', norm=LogNorm(1e2, 5e3)) -plt.colorbar(cm ,ax=ax) -gcf() diff --git a/scripts/re.jl b/scripts/re.jl deleted file mode 100644 index e6c5d5a..0000000 --- a/scripts/re.jl +++ /dev/null @@ -1,83 +0,0 @@ -using DrWatson -using Cubature -using PyPlot -@quickactivate "Qwind" -using Qwind - -function compute_fr(integrator) - data = integrator.p.data - fr_r = data[:disc_radiation_field_r] .* exp.(-data[:tauuv]) .* (1 .+ data[:fm]) - fr_z = data[:disc_radiation_field_z] .* exp.(-data[:tauuv]) .* (1 .+ data[:fm]) - return hcat(fr_r, fr_z) -end - -model = Model("scripts/config_re.yaml") -#iterations_dict = Dict() -#run!(model, iterations_dict) - - -function compute_vertical_flux(r) - disk_flux_norel(model.bh, r) -end - -function estimate_from_force(r, z=1) - acc_units = model.bh.M *G /model.bh.Rg/model.bh.Rg; #Acceleration scaling; - res, err = integrate_radiation_force_integrand( - model.rt, - r, - z, - 6, - 1400, - atol = 0, - rtol = 1e-3, - norm = Cubature.INDIVIDUAL, - maxevals = 0, - ) - radiation_constant = compute_radiation_constant(model.rt.radiation) / model.rt.radiation.fuv; - force = z * radiation_constant .* res[2]; - a_z = force * acc_units; - F = a_z * C / (SIGMA_E) - return F -end - -function find_closest_fraction(number) - distance = Inf - num = 0 - den = 0 - for i in 1:9 - for j in 1:9 - dist = abs(number - i/j) - if dist < distance - distance = dist - num = i - den = j - end - end - end - return "$num/$den" -end - -r_range = range(200, 1000, length=50); -flux = compute_vertical_flux.(r_range); -estimate = estimate_from_force.(r_range, 1); -fig, ax = plt.subplots() -ax.plot(r_range, flux, label="SS Flux", linewidth=2) -ax.plot(r_range, estimate, label="Qwind", linewidth=4, alpha=0.5) -#ratios = estimate ./ flux -#ax.plot(r_range, ratios, label="Qwind") -#ax.set_yscale("log") -ax.legend() -#ax.set_ylim(0.0, 2) - -# lets check we recover total luminosity. - -function S(r) - return 3 * 0.5 / (8 * π * model.bh.efficiency) / r^2 # * ( 1- sqrt(6/r)) / r^2 -end -using QuadGK -integral, err = quadgk(S, 6, 2000, rtol=1e-8) -println(integral * 2 * π * 0.85) - - -println(compute_bolometric_luminosity(model.bh)) - diff --git a/scripts/refine_quadtree.jl b/scripts/refine_quadtree.jl deleted file mode 100644 index 1aff60b..0000000 --- a/scripts/refine_quadtree.jl +++ /dev/null @@ -1,25 +0,0 @@ -using DrWatson -@quickactivate "Qwind" -using RegionTrees, DataFrames, CSV, YAML, Printf, DataFrames -using Qwind - -kdtree_data = CSV.read("kdtree_data.csv", DataFrame) - -bh = BlackHole(1e8, 0.5, 0.5) - -wind_kdtree = create_wind_kdtree( - kdtree_data.r, - kdtree_data.z, - kdtree_data.zmax, - kdtree_data.z_0, - kdtree_data.width, - kdtree_data.n, -) - -quadtree = create_and_refine_quadtree( - wind_kdtree, - Rg = bh.Rg, - atol = 1e-2, - rtol = 1e-2, - cell_min_size = 1e-6, -) diff --git a/scripts/run.jl b/scripts/run.jl deleted file mode 100644 index c75cd24..0000000 --- a/scripts/run.jl +++ /dev/null @@ -1,80 +0,0 @@ -using Distributed -@everywhere using DrWatson -@everywhere @quickactivate "Qwind" -using YAML, Profile, PProf, PyCall, ProgressMeter, PyPlot, JLD2 -using Qwind -LogNorm = matplotlib.colors.LogNorm -include("scripts/plotting.jl") - -config_path = "configs/rin_scan.yaml" -config = YAML.load_file(config_path, dicttype = Dict{Symbol,Any}) -try - mv(config[:integrator][:save_path], "backup", force = true) -catch -end -model = Model(config); -#iterations_dict = Dict(); -#run!(model, iterations_dict) -#lr, lw = Qwind.compute_lines_range(model) -#fig, ax = plt.subplots() -#for l in lr -# ax.axvline(l) -#end - - - -integrators = iterations_dict[1]["integrators"]; -max_times = get_intersection_times(integrators); - -fig, ax =plt.subplots() -for integrator in integrators - ax.plot(integrator.p.data[:r], integrator.p.data[:z]) -end - - -dense_integrators = DenseIntegrator.(integrators); - -max_indices = Dict() -for integrator in dense_integrators - max_indices[integrator.id] = searchsorted_nearest(integrator.t, max_times[integrator.id]) -end - - -sliced_integs = [] -for integ in dense_integrators - push!(sliced_integs, Qwind.slice_integrator(integ, fi = max_indices[integ.id])) -end - -fig, ax =plt.subplots() -for integ in sliced_integs - ax.plot(integ.r, integ.z) -end - -integrators_interpolated_linear = interpolate_integrators( - integrators, - max_times = max_times, - n_timesteps = 100, - log = true, -) - - -fig, ax = plt.subplots() -#for integ in integrators -# ax.plot(integ.p.data[:r], integ.p.data[:z]) -#end -for integ in integrators_interpolated_linear - ax.plot(integ.r, integ.z) -end - -max_times = get_intersection_times(integrators); - -integs2 = Qwind.filter_close_trajectories(integrators, 0.5); -length(integs2) - -points = Hull(integrators, max_times); - -QwindPlotting.plot_wind_hull(hull, zmax=25); - -fig, ax = plt.subplots() -ps = reduce(hcat, points) -ax.scatter(ps[1,:], ps[2,:]) diff --git a/scripts/run_parallel.jl b/scripts/run_parallel.jl deleted file mode 100644 index bf9ce9f..0000000 --- a/scripts/run_parallel.jl +++ /dev/null @@ -1,31 +0,0 @@ -using Distributed -@everywhere using DrWatson -@everywhere @quickactivate "Qwind" -using Qwind -using YAML, Profile, PProf, PyCall -include("scripts/plotting.jl") - -config_path = "configs/debug.yaml" -config = YAML.load_file(config_path, dicttype = Dict{Symbol,Any}) -try - mv(config[:integrator][:save_path], "backup", force = true) -catch -end -model = Model(config_path); -iterations_dict = Dict(); -#run_iteration!(model, iterations_dict, it_num=1, parallel=true); -#run!(model, iterations_dict, parallel=true) - - -fig, ax = QwindPlotting.plot_streamlines(iterations_dict[1]["integrators"], linestyle="-") - - -using PyPlot -LogNorm = matplotlib.colors.LogNorm; - - -interp = iterations_dict[2]["radiative_transfer"].interpolator; -vgrid = interp.velocity_grid -fig, ax = plt.subplots() -cm = ax.pcolormesh(vgrid.r_range, vgrid.z_range, vgrid.vr_grid') -plt.colorbar(cm, ax=ax) diff --git a/scripts/ryota_meeting.jl b/scripts/ryota_meeting.jl deleted file mode 100644 index 31aba58..0000000 --- a/scripts/ryota_meeting.jl +++ /dev/null @@ -1,160 +0,0 @@ -using DrWatson -@quickactivate "Qwind" -using RegionTrees, DataFrames, CSV, YAML, Printf -using Qwind -using Profile -using PProf -using PyPlot -LogNorm = matplotlib.colors.LogNorm -using ColorSchemes -using ColorSchemes: colorschemes - - -function plot_density_grid_nn( - windkdtree::KDTree; - xl = 0, - xh = 5000, - yl = 0, - yh = 5000, - nr = 500, - nz = 501, - vmin = nothing, - vmax = nothing, - fig = nothing, - ax = nothing, -) - if fig === nothing || ax === nothing - fig, ax = plt.subplots() - end - grid = zeros(Float64, nr, nz) - r_range = range(xl, stop = xh, length = nr) - z_range = range(yl, stop = yh, length = nz) - for (i, r) in enumerate(r_range[1:(end - 1)]) - for (j, z) in enumerate(z_range[1:(end - 1)]) - grid[i, j] = get_density(windkdtree, r, z) - end - end - cm = ax.pcolormesh(r_range, z_range, grid', norm = LogNorm(vmin = vmin, vmax = vmax)) - cbar = plt.colorbar(cm, ax = ax) - cbar.ax.set_ylabel("Density [cm^-3]", labelpad = 10, rotation = -90) - ax.set_xlabel("r [Rg]") - ax.set_ylabel("z [Rg]") - ax.set_xlim(xl, xh) - ax.set_ylim(yl, yh) - return fig, ax -end - - -function plot_density_grid( - quadtree::Cell; - xl = 0, - xh = 5000, - yl = 0, - yh = 5000, - nr = 500, - nz = 501, - vmin = nothing, - vmax = nothing, - fig = nothing, - ax = nothing, -) - if fig === nothing || ax === nothing - fig, ax = plt.subplots() - end - grid = zeros(Float64, nr, nz) - r_range = range(xl, stop = xh, length = nr) - z_range = range(yl, stop = yh, length = nz) - for (i, r) in enumerate(r_range[1:(end - 1)]) - for (j, z) in enumerate(z_range[1:(end - 1)]) - grid[i, j] = get_density(quadtree, r, z) - end - end - cm = ax.pcolormesh(r_range, z_range, grid', norm = LogNorm()) - plt.colorbar(cm, ax = ax) - ax.set_xlabel("r [Rg]") - ax.set_ylabel("z [Rg]") - ax.set_xlim(xl, xh) - ax.set_ylim(yl, yh) - return fig, ax -end - - -config = YAML.load_file("scripts/config_uniform_ic.yaml"); - -black_hole = BlackHole(config); - -radiation = @eval $(Symbol(config["radiation"]["mode"]))(black_hole, config); - -radiative_transfer = - @eval $(Symbol(config["radiative_transfer"]["mode"]))(radiation, config); - -grid = Grid(config); - -initial_conditions = - @eval $(Symbol(config["initial_conditions"]["mode"]))(radiation, black_hole, config); - -integrators = initialize_integrators( - radiative_transfer, - grid, - initial_conditions, - atol = config["integrator"]["atol"], - rtol = config["integrator"]["rtol"], -); -run_integrators!(integrators); - -radiative_transfer = update_radiative_transfer(radiative_transfer, integrators); -qt = radiative_transfer.quadtree; -nn = radiative_transfer.windkdtree; - -fig, ax = subplots() -Qwind.plot_streamlines(integrators, fig, ax, colorscheme = "matter"); -savefig("results/ryota_meeting/streamlines.png", dpi = 300, bbox_inches = "tight") -show() - -fig, ax = subplots() -plot_density_grid_nn(nn, fig = fig, ax = ax, vmin = 1e2, nr = 500, nz = 501); -savefig("results/ryota_meeting/nn_density.png", dpi = 300, bbox_inches = "tight") -show() - -fig, ax = subplots() -plot_density_grid_nn(nn, fig = fig, ax = ax, nr = 1000, nz = 1001); -Qwind.plot_streamlines(integrators, fig, ax, color = "white", alpha = 0.5); -savefig("results/ryota_meeting/nn_density_with_lines.png", dpi = 300, bbox_inches = "tight") -show() - -fig, ax = subplots() -plot_density_grid(qt, fig = fig, ax = ax, vmin = 1e2, nr = 500, nz = 501); -savefig("results/ryota_meeting/qt_density.png", dpi = 300, bbox_inches = "tight") -show() - -fig, ax = subplots() -plot_density_grid(qt, fig = fig, ax = ax, vmin = 1e2, nr = 500, nz = 501); -plot_grid(qt, 7, fig, ax) -savefig("results/ryota_meeting/qt_density_with_grid.png", dpi = 300, bbox_inches = "tight") -show() - -fig, ax = subplots() -plot_density_grid(qt, fig = fig, ax = ax, vmin = 1e2, nr = 500, nz = 501); -plot_grid(qt, 7, fig, ax) -Qwind.plot_streamlines(integrators, fig, ax, color = "white", alpha = 0.5); -ax.set_xlim(0, 5000) -ax.set_ylim(0, 5000) -savefig( - "results/ryota_meeting/qt_density_with_grid_and_lines.png", - dpi = 300, - bbox_inches = "tight", -) -show() - -fig, ax = subplots() -plot_density_grid(qt, fig = fig, ax = ax, vmin = 1e2, nr = 500, nz = 501); -plot_grid(qt, 10, fig, ax) -tau, coords = - compute_xray_tau(qt, 1000.0, 250.0, 0.0, 0.0, radiation.xray_luminosity, black_hole.Rg) -coords = hcat(coords...) -ax.plot(coords[1, :], coords[2, :], "o-", color = "blue", linewidth = 1, markersize = 3) -ax.set_xlim(0, 1200) -ax.set_ylim(0, 300) -savefig("results/ryota_meeting/xray_tau_example.png", dpi = 300, bbox_inches = "tight") -show() - diff --git a/scripts/scattering.jl b/scripts/scattering.jl deleted file mode 100644 index 65db1cb..0000000 --- a/scripts/scattering.jl +++ /dev/null @@ -1,15 +0,0 @@ -using Distributed -@everywhere using DrWatson -@everywhere @quickactivate "Qwind" -@everywhere using Qwind -using PyPlot -LogNorm = matplotlib.colors.LogNorm -Normalize = matplotlib.colors.Normalize -using Profile, PProf - -function read_table(file) - df = CSV.read(file, DataFrame, header=false); - return Matrix(df) -end - -ry_values = "/home/arnau/Downloads/GX13p1_density_ps_05Ledd/np_1000.txt" diff --git a/scripts/scratch.jl b/scripts/scratch.jl deleted file mode 100644 index e69de29..0000000 diff --git a/scripts/test_interpolation.jl b/scripts/test_interpolation.jl deleted file mode 100644 index 762996c..0000000 --- a/scripts/test_interpolation.jl +++ /dev/null @@ -1,49 +0,0 @@ -using DrWatson -@quickactivate "Qwind" -using PyPlot -LogNorm = matplotlib.colors.LogNorm -using Qwind -using Profile -using PProf - -model = Model("paper2/plot_nn.yaml"); -iterations_dict = Dict(); - -do_iteration!(model, iterations_dict, it_num=1); - -integrators = iterations_dict[1]["integrators"] -plot_streamlines(integrators) - -lines_kdtrees = Qwind.create_lines_kdtrees(iterations_dict[1]["integrators"], n_timesteps=10000); - -fig, ax = plt.subplots() -for line in lines_kdtrees - ax.plot(line.r, line.z) -end - -nr = 150; -nz = 151; -#r_range, z_range = Qwind.get_spatial_grid(lines_kdtrees, nr, nz); -r_range = range(10, 20, length=nr) -z_range = range(0.001, 0.2, length=nz) -#r_range = range(31.5, 33.5, length=nr) -#z_range = range(0.1, 0.2, length=nz) -function get_density_grid(lines_kdtrees, r_range, z_range) - density_grid = zeros(Float64, nr, nz); - for (i, r) in enumerate(r_range) - for (j, z) in enumerate(z_range) - density_grid[i, j] = get_density(lines_kdtrees, r, z); - end - end - return density_grid -end -@time density_grid = get_density_grid(lines_kdtrees, r_range, z_range); - -fig, ax= plt.subplots(); -cm = ax.pcolormesh(r_range, z_range, density_grid', norm=LogNorm()); -plt.colorbar(cm, ax=ax); -for line in lines_kdtrees - ax.scatter(line.r, line.z, color="white", s=0.1) -end -ax.set_xlim(r_range[1], r_range[end]) -ax.set_ylim(z_range[1], z_range[end]) diff --git a/scripts/tests.jl b/scripts/tests.jl deleted file mode 100644 index b68d317..0000000 --- a/scripts/tests.jl +++ /dev/null @@ -1,58 +0,0 @@ -using Distributed -@everywhere using DrWatson -@everywhere @quickactivate "Qwind" -using YAML, Profile, PProf, PyCall, ProgressMeter, JLD2, PyPlot -using Qwind -LogNorm = matplotlib.colors.LogNorm -include("scripts/plotting.jl") - -config_path = "configs/debug.yaml" -config = YAML.load_file(config_path, dicttype = Dict{Symbol,Any}) -try - mv(config[:integrator][:save_path], "backup", force = true) -catch -end -model = Model(config); -iterations_dict = Dict(); -run!(model, iterations_dict) - -integrator = initialize_integrator(model.rt, model.wind_grid, model.ic, 25.0, 1); -solve!(integrator) - -fig, ax = plt.subplots() -#ax.loglog(integrator.p.data[:z], integrator.p.data[:n]) -ax.plot(integrator.p.data[:r], integrator.p.data[:z]) - - - -#update_radiative_transfer(model.rt, iterations_dict[1]["integrators"]) -integrators = iterations_dict[6]["integrators"]; -fig, ax = plt.subplots() -for integ in integrators - #ax.loglog(integ.p.data[:z], integ.p.data[:n]) - ax.plot(integ.p.data[:r], integ.p.data[:z]) -end - - -r0 = [integ.p.r0 for integ in integrators]; -max_times = get_intersection_times(integrators); -hull = Hull(integrators, max_times); - -density_grid_old = DensityGrid(integrators, max_times, hull, nr=1000, nz=500); - -fig, ax = plt.subplots() -ax.pcolormesh(density_grid.r_range, density_grid.z_range, density_grid.grid', norm=LogNorm()) - -r_grid = density_grid.r_range' .* ones(length(density_grid.z_range)); -z_grid = density_grid.z_range .* ones(length(density_grid.r_range))'; - -#@load "iterations_dict.jld2" iterations_dict - -integrators = load_trajectories("./results.hdf5", 2); - -fig, ax = plt.subplots() -for integ in integrators[1400:1500] - ax.plot(integ.r, integ.z) -end - - From 5a162b318d69b4567fc7edbc73f59018809042af Mon Sep 17 00:00:00 2001 From: Arnau Quera-Bofarull Date: Thu, 10 Mar 2022 12:07:50 +0000 Subject: [PATCH 2/4] removed old configs --- configs/bh_scan.yaml | 60 --------------------------------------- configs/debug.yaml | 49 -------------------------------- configs/proga.yaml | 65 ------------------------------------------- configs/rin_scan.yaml | 59 --------------------------------------- 4 files changed, 233 deletions(-) delete mode 100644 configs/bh_scan.yaml delete mode 100644 configs/debug.yaml delete mode 100644 configs/proga.yaml delete mode 100644 configs/rin_scan.yaml diff --git a/configs/bh_scan.yaml b/configs/bh_scan.yaml deleted file mode 100644 index 6314066..0000000 --- a/configs/bh_scan.yaml +++ /dev/null @@ -1,60 +0,0 @@ -system: - qwind_path: /cosma/home/dp004/dc-quer1/Qwind.jl - sysimage_path: /cosma/home/dp004/dc-quer1/sys_everywhere.so - n_cpus: 16 - partition: cosma - account: durham - max_time: 72:00:00 - job_name: qwind - -black_hole: - M: "@vary grid 1e6,1e7,1e8,1e9,1e10" # solar masses - mdot: "@vary log 0.05 0.5 10" - spin: 0.0 # spin parameter (-1,1) - -radiation: - relativistic: true - f_uv: disk - f_x: 0.15 - n_r: 2000 - disk_r_in: 6.0 - z_xray: 0.0 - disk_height: 0.0 - uv_opacity: thomson - xray_opacity: boost - xray_scattering: false - tau_uv_calculation: disk - vacuum_density: 1e2 - update_grid_method: "average" - grid_nz: 500 - grid_nr: auto - mu_nucleon: 0.61 - -initial_conditions: - mode: cak - r_in: "@vary grid 10,25,50" - r_fi: 1500.0 - n_trajs: auto - z_0: 0.0 - v_0: thermal - K: 0.03 - alpha: 0.6 - use_precalculated: false - -integrator: - n_iterations: 40 - save_path: "/cosma5/data/durham/covid19/arnau/qwind/revised/bh_scan_tauuv" - integrator_r_min: 6.0 - integrator_r_max: 10000.0 - integrator_z_min: 0.0 - integrator_z_max: 10000.0 - -numerical_tolerances: - disk_integral_atol: 0 - disk_integral_rtol: 1e-3 - disk_integral_maxevals: 10000 - integrator_atol: 1e-8 - integrator_rtol: 1e-3 - scattering_atol: 0 - scattering_rtol: 1e-3 - diff --git a/configs/debug.yaml b/configs/debug.yaml deleted file mode 100644 index 5bf2ca7..0000000 --- a/configs/debug.yaml +++ /dev/null @@ -1,49 +0,0 @@ -black_hole: - M: 1e8 # solar masses - mdot: 0.5 # Eddington units - spin: 0.0 # spin parameter (-1,1) - -radiation: - relativistic: true - f_uv: qsosed - f_x: qsosed - n_r: 2000 - disk_r_in: corona_radius - z_xray: 0.0 - disk_height: 0.0 - uv_opacity: thomson - xray_opacity: boost - xray_scattering: false - tau_uv_calculation: no_tau_uv - vacuum_density: 1e2 - update_grid_method: "average" - grid_nz: 500 - grid_nr: auto - -initial_conditions: - mode: cak - r_in: warm_radius - r_fi: 1500.0 - n_trajs: 20 - trajs_spacing: log - z0: 0.0 - K: 0.03 - alpha: 0.6 - use_precalculated: true - -integrator: - n_iterations: 5 - save_path: "./tests" - integrator_r_min: 6.0 - integrator_r_max: 10000.0 - integrator_z_min: 0.0 - integrator_z_max: 10000.0 - -numerical_tolerances: - disk_integral_atol: 0 - disk_integral_rtol: 5e-4 - disk_integral_maxevals: 10000 - integrator_atol: 1e-8 - integrator_rtol: 1e-3 - scattering_atol: 0 - scattering_rtol: 1e-3 diff --git a/configs/proga.yaml b/configs/proga.yaml deleted file mode 100644 index 4ba7122..0000000 --- a/configs/proga.yaml +++ /dev/null @@ -1,65 +0,0 @@ -system: - qwind_path: /cosma/home/dp004/dc-quer1/Qwind.jl - sysimage_path: /cosma/home/dp004/dc-quer1/sys_everywhere.so - n_cpus: 16 - partition: cosma6 - account: dp004 - max_time: 72:00:00 - job_name: proga04 - -black_hole: - M: 1e8 # solar masses - mdot: 0.5 # Eddington units - spin: 0.0 # spin parameter (-1,1) - -radiation: - relativistic: false - f_uv: 0.9 - f_x: 0.1 - n_r: 2000 - disk_r_in: 6.0 - disk_r_out: 3000.0 - z_xray: 1.5 - disk_height: 1.5 - xray_opacity: thomson - uv_opacity: thomson - xray_scattering: false - tau_uv_calculation: no_tau_uv - fm_interp_method_flag: analytic - vacuum_density: 1e2 - update_grid_method: "average" - grid_nz: 250 - grid_nr: 250 - mu_nucleon: 1.0 - -initial_conditions: - mode: cak - r_in: 60.0 - r_fi: 1600.0 - n_trajs: auto - trajs_spacing: log - z_0: 1.5 - v_0: thermal - K: 0.03 - alpha: 0.6 - z_disk: 0.0 - use_precalculated: false - -integrator: - n_iterations: 50 - #save_path: "/cosma7/data/dp004/dc-quer1/qwind/proga_honest" - save_path: "./proga_honest2" - integrator_r_min: 6.0 - integrator_r_max: 3000.0 - integrator_z_min: 0.0 - integrator_z_max: 3000.0 - -numerical_tolerances: - disk_integral_atol: 0 - disk_integral_rtol: 5e-4 - disk_integral_maxevals: 10000 - integrator_atol: 1e-8 - integrator_rtol: 1e-3 - scattering_atol: 0 - scattering_rtol: 1e-3 - diff --git a/configs/rin_scan.yaml b/configs/rin_scan.yaml deleted file mode 100644 index 9872643..0000000 --- a/configs/rin_scan.yaml +++ /dev/null @@ -1,59 +0,0 @@ -system: - qwind_path: /cosma/home/dp004/dc-quer1/Qwind.jl - sysimage_path: /cosma/home/dp004/dc-quer1/sys_everywhere.so - n_cpus: 16 - partition: cosma - account: durham - max_time: 72:00:00 - job_name: qwind - -black_hole: - M: 1e8 # solar masses - mdot: 0.5 # Eddington units - spin: 0.0 # spin parameter (-1,1) - -radiation: - relativistic: true - f_uv: disk - f_x: 0.15 - n_r: 2000 - disk_r_in: 6.0 - z_xray: 0.0 - disk_height: 0.0 - uv_opacity: thomson - xray_opacity: boost - xray_scattering: false - tau_uv_calculation: disk - vacuum_density: 1e2 - update_grid_method: "average" - grid_nz: 250 - grid_nr: auto - mu_nucleon: 0.61 - -initial_conditions: - mode: cak - r_in: "@vary log 10 250 10" - r_fi: 1500.0 - n_trajs: auto - z_0: 0.0 - v_0: thermal - K: 0.03 - alpha: 0.6 - use_precalculated: false - -integrator: - n_iterations: 40 - save_path: "/cosma5/data/durham/covid19/arnau/qwind/revised/rin_scan_tauuv" - integrator_r_min: 6.0 - integrator_r_max: 10000.0 - integrator_z_min: 0.0 - integrator_z_max: 10000.0 - -numerical_tolerances: - disk_integral_atol: 0 - disk_integral_rtol: 1e-3 - disk_integral_maxevals: 10000 - integrator_atol: 1e-8 - integrator_rtol: 1e-3 - scattering_atol: 0 - scattering_rtol: 1e-3 From f4c579a2eb14753095031d6beb938cf87e29ea9f Mon Sep 17 00:00:00 2001 From: Arnau Quera-Bofarull Date: Thu, 10 Mar 2022 13:05:24 +0000 Subject: [PATCH 3/4] testing cluster docs --- configs/parameter_scan_example.yaml | 58 +++++++++++++++++++++++++++++ src/cluster_utils.jl | 11 ++++-- 2 files changed, 66 insertions(+), 3 deletions(-) create mode 100644 configs/parameter_scan_example.yaml diff --git a/configs/parameter_scan_example.yaml b/configs/parameter_scan_example.yaml new file mode 100644 index 0000000..eac55b5 --- /dev/null +++ b/configs/parameter_scan_example.yaml @@ -0,0 +1,58 @@ +system: + qwind_path: /cosma/home/dp004/dc-quer1/Qwind.jl + n_cpus: 16 + partition: cosma + account: durham + max_time: 72:00:00 + job_name: qwind + +black_hole: + M: "@vary grid 1e8,1e9,1e10" # take the values specified + mdot: "@vary log 0.1 0.5 5" # from 0.1 to 0.5, 5 values log spaced + spin: 0.0 + +radiation: + relativistic: true + f_uv: 0.85 + f_x: 0.15 + n_r: 2000 + disk_r_in: 6.0 + z_xray: 0.0 + disk_height: 0.0 + uv_opacity: thomson + xray_opacity: boost + xray_scattering: false + tau_uv_calculation: no_tau_uv + vacuum_density: 1e2 + update_grid_method: "average" + grid_nz: 250 + grid_nr: auto + mu_nucleon: 0.61 + +initial_conditions: + mode: cak + r_in: 20.0 + r_fi: 1500.0 + n_trajs: 50 + z0: 0.0 + K: 0.03 + alpha: 0.6 + use_precalculated: false + +integrator: + n_iterations: 5 + save_path: "./parameter_scan_example" + integrator_r_min: 6.0 + integrator_r_max: 10000.0 + integrator_z_min: 0.0 + integrator_z_max: 10000.0 + +numerical_tolerances: + disk_integral_atol: 0 + disk_integral_rtol: 1e-3 + disk_integral_maxevals: 10000 + integrator_atol: 1e-8 + integrator_rtol: 1e-3 + scattering_atol: 0 + scattering_rtol: 1e-3 + diff --git a/src/cluster_utils.jl b/src/cluster_utils.jl index c2f7f19..74056bb 100644 --- a/src/cluster_utils.jl +++ b/src/cluster_utils.jl @@ -14,7 +14,7 @@ end function cosma_script(; qwind_path, - sysimage_path, + sysimage_path=nothing, n_cpus = 1, job_name = "qwind", script_number = 1, @@ -25,6 +25,11 @@ function cosma_script(; save_path = nothing, run_script_path = nothing, ) + if sysimage_path === nothing + sysimage_str = "" + else + sysimage_str = "--sysimage \"$sysimage_path\"" + end script = """ #!/bin/bash -l #SBATCH --ntasks $n_cpus @@ -37,7 +42,7 @@ function cosma_script(; export JULIA_PROJECT=$qwind_path - julia --sysimage "$sysimage_path" $run_script_path -m $script_number + julia $sysimage_str $run_script_path -m $script_number """ open(save_path, "w") do io write(io, script) @@ -47,7 +52,7 @@ end function make_cosma_scripts( n_models; qwind_path, - sysimage_path, + sysimage_path=nothing, n_cpus = 1, path = nothing, job_name = "qwind", From c6d190e1a139993c331789f06f751418f3294325 Mon Sep 17 00:00:00 2001 From: Arnau Quera-Bofarull Date: Thu, 10 Mar 2022 13:05:47 +0000 Subject: [PATCH 4/4] added cluster docs --- docs/src/run_parameter_grid.md | 99 ++++++++++++++++++++++++++++++++++ 1 file changed, 99 insertions(+) create mode 100644 docs/src/run_parameter_grid.md diff --git a/docs/src/run_parameter_grid.md b/docs/src/run_parameter_grid.md new file mode 100644 index 0000000..deaf9dc --- /dev/null +++ b/docs/src/run_parameter_grid.md @@ -0,0 +1,99 @@ +# Running a parameter grid + +The main goal of Qwind is to be able to run fast parameter searches. To make it easy, the user can specify ranges to vary in the parameter file. + +Let's illustrate this with an example changing BH mass and mass accretion rate. The parameter file we use (located at `configs/parameter_scan_example.yaml` ) is + + + +```yaml +system: + qwind_path: /path/to/Qwind.jl + n_cpus: 16 + partition: cosma + account: durham + max_time: 72:00:00 + job_name: qwind + +black_hole: + M: "@vary grid 1e8,1e9,1e10" # take the values specified + mdot: "@vary log 0.1 0.5 5" # from 0.1 to 0.5, 5 values log spaced + spin: 0.0 + +radiation: + relativistic: true + f_uv: 0.85 + f_x: 0.15 + n_r: 2000 + disk_r_in: 6.0 + z_xray: 0.0 + disk_height: 0.0 + uv_opacity: thomson + xray_opacity: boost + xray_scattering: false + tau_uv_calculation: no_tau_uv + vacuum_density: 1e2 + update_grid_method: "average" + grid_nz: 250 + grid_nr: auto + mu_nucleon: 0.61 + +initial_conditions: + mode: cak + r_in: 20.0 + r_fi: 1500.0 + n_trajs: 50 + z0: 0.0 + K: 0.03 + alpha: 0.6 + use_precalculated: false + +integrator: + n_iterations: 5 + save_path: "./parameter_scan_example" + integrator_r_min: 6.0 + integrator_r_max: 10000.0 + integrator_z_min: 0.0 + integrator_z_max: 10000.0 + +numerical_tolerances: + disk_integral_atol: 0 + disk_integral_rtol: 1e-3 + disk_integral_maxevals: 10000 + integrator_atol: 1e-8 + integrator_rtol: 1e-3 + scattering_atol: 0 + scattering_rtol: 1e-3 +``` + +The first config blog is the `system` configuration block, which specifies the path to our Qwind installation, the number of cpus to use for each model, the partition and account of our HPC cluster, the maximum time of the job, and the job name. Note that this only works for a SLURM scheduler. + +We specify the BH mass and mass accretion rate to vary in a certain range. For the BH mass we have specified a grid of values, and for mdot we specify a log range. Another possible option could have been `@vary linear 0.1 0.5 5` which would have varied the value linearly rather than in log-space. + +We can now create a folder hierarchy structure with the variations like this + +```julia +using Qwind +create_models_folders("configs/parameter_scan_example.yaml") +``` + +This will create a folder called `parameter_scan_example` with the following structure + +``` +├── all_configs.yaml +| +├── model_001 +│ ├── config.yaml +│ ├── submit.sh +├── model_002 +│ ├── config.yaml +│ ├── submit.sh +... +├── run_model.jl +├── stdout +├── submit_all.sh +└── all_configs.yaml +``` + +In total we can see 15 models, as we would expect from 3 x 5. Each model has its own `config.yaml` plus a submit script for our system. We can then submit the scripts one by one, or submit them all at once with `bash submit_all.sh`. The logs of the jobs will be stored in the `stdout` folder. +