|
| 1 | +# # Alternative Kappa Calibration Example |
| 2 | +# ## Overview |
| 3 | +#= |
| 4 | +In this example, just like in kappa_calibration.jl, we use the inverse problem to calibrate the von-karman constant, κ in |
| 5 | +the equation: u(z) = u^* / κ log (z / z0), |
| 6 | +which represents the wind profile in Monin-Obukhov |
| 7 | +Similarity Theory (MOST) formulations. We use the same dataset: https://turbulence.pha.jhu.edu/Channel_Flow.aspx |
| 8 | +
|
| 9 | +Instead of using u^* as an observable, we use the dataset's u, and each ensemble member will estimate u |
| 10 | +through the profile equation u(z) = u^* / κ log (z / z0). |
| 11 | +=# |
| 12 | + |
| 13 | +# ## Prerequisites |
| 14 | +#= |
| 15 | +[EnsembleKalmanProcess.jl](https://github.com/CliMA/EnsembleKalmanProcesses.jl), |
| 16 | +=# |
| 17 | + |
| 18 | +# ## Example |
| 19 | + |
| 20 | +# First, we import relevant modules. |
| 21 | +using LinearAlgebra, Random |
| 22 | +using Distributions, Plots |
| 23 | +using EnsembleKalmanProcesses |
| 24 | +using EnsembleKalmanProcesses.ParameterDistributions |
| 25 | +const EKP = EnsembleKalmanProcesses |
| 26 | + |
| 27 | +using Downloads |
| 28 | +using DelimitedFiles |
| 29 | + |
| 30 | +FT = Float64 |
| 31 | + |
| 32 | +mkpath(joinpath(@__DIR__, "data")) # create data folder if not exists |
| 33 | +web_datafile_path = "https://turbulence.oden.utexas.edu/channel2015/data/LM_Channel_5200_mean_prof.dat" |
| 34 | +localfile = "data/profiles.dat" |
| 35 | +Downloads.download(web_datafile_path, localfile) |
| 36 | +data_mean_velocity = readdlm("data/profiles.dat", skipstart = 112) ## We skip 72 lines (header) and 40(laminar layer) |
| 37 | + |
| 38 | +web_datafile_path = "https://turbulence.oden.utexas.edu/channel2015/data/LM_Channel_5200_mean_stdev.dat" |
| 39 | +localfile = "data/vel_stdev.dat" |
| 40 | +Downloads.download(web_datafile_path, localfile) |
| 41 | +# We skip 72 lines (header) and 40(laminar layer) |
| 42 | +data_stdev_velocity = readdlm("data/vel_stdev.dat", skipstart = 112) |
| 43 | + |
| 44 | +# We extract the required info for this problem |
| 45 | +u_star_obs = 4.14872e-02 # add noise later |
| 46 | +z0 = FT(0.0001) |
| 47 | +κ = 0.4 |
| 48 | + |
| 49 | +# turn u into distributions |
| 50 | +u = data_mean_velocity[:, 3] * u_star_obs |
| 51 | +z = data_mean_velocity[:, 1] |
| 52 | +u = u[1:(length(u) - 1)] # filter out last element because σᵤ is only of length 727, not 728 |
| 53 | +z = z[1:(length(z) - 1)] |
| 54 | + |
| 55 | +σᵤ = data_stdev_velocity[:, 4] * u_star_obs |
| 56 | +dist_u = Array{Normal{Float64}}(undef, length(u)) |
| 57 | +for i in 1:length(u) |
| 58 | + dist_u[i] = Normal(u[i], σᵤ[i]) |
| 59 | +end |
| 60 | + |
| 61 | +# u(z) = u^* / κ log (z / z0) |
| 62 | +function physical_model(parameters, inputs) |
| 63 | + κ = parameters[1] # this is being updated by the EKP iterator |
| 64 | + (; u_star_obs, z, z0) = inputs |
| 65 | + u_profile = u_star_obs ./ κ .* log.(z ./ z0) |
| 66 | + return u_profile |
| 67 | +end |
| 68 | + |
| 69 | +function G(parameters, inputs, u_profile = nothing) |
| 70 | + if (isnothing(u_profile)) |
| 71 | + u_profile = physical_model(parameters, inputs) |
| 72 | + end |
| 73 | + return [maximum(u_profile) - minimum(u_profile), mean(u_profile)] |
| 74 | +end |
| 75 | + |
| 76 | +Γ = 0.0001 * I |
| 77 | +η_dist = MvNormal(zeros(2), Γ) |
| 78 | +noisy_u_profile = [rand(dist_u[i]) for i in 1:length(u)] |
| 79 | +y = G(nothing, nothing, noisy_u_profile) |
| 80 | + |
| 81 | +parameters = (; κ) |
| 82 | +inputs = (; u_star_obs, z, z0) |
| 83 | +# y = G(parameters, inputs) .+ rand(η_dist) |
| 84 | + |
| 85 | +# Assume that users have prior knowledge of approximate truth. |
| 86 | +# (e.g. via physical models / subset of obs / physical laws.) |
| 87 | +prior_u1 = constrained_gaussian("κ", 0.35, 0.25, 0, Inf); |
| 88 | +prior = combine_distributions([prior_u1]) |
| 89 | + |
| 90 | +# Set up the initial ensembles |
| 91 | +N_ensemble = 5; |
| 92 | +N_iterations = 10; |
| 93 | + |
| 94 | +rng_seed = 41 |
| 95 | +rng = Random.MersenneTwister(rng_seed) |
| 96 | +initial_ensemble = EKP.construct_initial_ensemble(rng, prior, N_ensemble); |
| 97 | + |
| 98 | +# Define EKP and run iterative solver for defined number of iterations |
| 99 | +ensemble_kalman_process = EKP.EnsembleKalmanProcess(initial_ensemble, y, Γ, Inversion(); rng = rng) |
| 100 | + |
| 101 | +for n in 1:N_iterations |
| 102 | + params_i = get_ϕ_final(prior, ensemble_kalman_process) |
| 103 | + G_ens = hcat([G(params_i[:, m], inputs) for m in 1:N_ensemble]...) |
| 104 | + EKP.update_ensemble!(ensemble_kalman_process, G_ens) |
| 105 | +end |
| 106 | + |
| 107 | +# Mean values in final ensemble for the two parameters of interest reflect the "truth" within some degree of |
| 108 | +# uncertainty that we can quantify from the elements of `final_ensemble`. |
| 109 | +final_ensemble = get_ϕ_final(prior, ensemble_kalman_process) |
| 110 | +mean(final_ensemble[1, :]) # [param, ens_no] |
| 111 | + |
| 112 | + |
| 113 | +ENV["GKSwstype"] = "nul" |
| 114 | +zrange = z |
| 115 | +plot(zrange, noisy_u_profile, c = :black, label = "Truth", linewidth = 2, legend = :bottomright) |
| 116 | +plot!(zrange, physical_model(parameters, inputs), c = :green, label = "Model truth", linewidth = 2)# reshape to convert from vector to matrix) |
| 117 | +plot!( |
| 118 | + zrange, |
| 119 | + [physical_model(get_ϕ(prior, ensemble_kalman_process, 1)[:, i]..., inputs) for i in 1:N_ensemble], |
| 120 | + c = :red, |
| 121 | + label = reshape(vcat(["Initial ensemble"], ["" for i in 1:(N_ensemble - 1)]), 1, N_ensemble), # reshape to convert from vector to matrix |
| 122 | +) |
| 123 | +plot!( |
| 124 | + zrange, |
| 125 | + [physical_model(final_ensemble[:, i]..., inputs) for i in 1:N_ensemble], |
| 126 | + c = :blue, |
| 127 | + label = reshape(vcat(["Final ensemble"], ["" for i in 1:(N_ensemble - 1)]), 1, N_ensemble), |
| 128 | +) |
| 129 | +xlabel!("Z") |
| 130 | +ylabel!("U") |
| 131 | +png("profile_plot") |
0 commit comments