diff --git a/.gitignore b/.gitignore index d189a7a..d16ba09 100644 --- a/.gitignore +++ b/.gitignore @@ -27,3 +27,9 @@ config/config.json plots/* .DS_Store + +# Tom's files +--temp/* +.DS_Store + +.ipynb* diff --git a/Project.toml b/Project.toml index 48e2171..5d5a3b0 100644 --- a/Project.toml +++ b/Project.toml @@ -7,13 +7,17 @@ version = "0.1.0-DEV" ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63" Atom = "c52e3926-4ff0-5f6e-af25-54175e0327b1" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Callbacks = "db1e321a-d383-57b4-a664-0144fd54e973" ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" +DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94" DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" GEMPIC = "b6d65c3a-4a4e-11e9-25d0-d309dc85ddeb" @@ -26,6 +30,7 @@ JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" Juno = "e5e0dc1b-0480-54bc-9374-aad01c23163d" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" @@ -35,6 +40,7 @@ Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" PkgTemplates = "14b8a8f1-9102-5b29-a752-f990bacb7fe1" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" ProfileView = "c46f51b8-102a-5cf2-8d2c-8597cb0e0da7" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" diff --git a/analysis/SIR_analysis.jl b/analysis/SIR_analysis.jl new file mode 100644 index 0000000..59c77e4 --- /dev/null +++ b/analysis/SIR_analysis.jl @@ -0,0 +1,1316 @@ +using DataFrames, DelimitedFiles +using LinearAlgebra +import Statistics as Stc +using GLM +using LaTeXStrings +# import Plots as plt +using Colors +using Random, Distributions +using PlotlyJS +using SpecialFunctions +import PlotlyBase: make_subplots +using Pkg +Pkg.add(url="https://github.com/cossio/TruncatedNormal.jl") +import TruncatedNormal as tn + +## Defining useful functions + +function gaussD(mu, sigma, x1) + value = 1/((2*pi)^0.5 * sigma^0.5) * exp(-0.5*(x1-mu)*sigma^(-1)*(x1-mu)) + return value +end + +function print_mat(M) + a = size(M)[1] + b = size(M)[2] + print("[") + for i in 1:a + if i != 1 println() end + print(M[i,1]) + for j in 2:b + print(", ") + print(M[i,j]) + end + end + println("]") +end + +function gaussMultiD(mu, sigma, sigma1, x1) + dim = length(mu) + value = 1/((2*pi)^(dim/2) * det(sigma)^0.5) * exp(-0.5*transpose(x1-mu)*sigma1*(x1-mu)) + return value +end + +function spectrum_at_point(T, PartStates, xyPoint, sigma, kxrange, kyrange, nkX, nkY) + kxs = (0:(nkX-1))/nkX*(kxrange[2]-kxrange[1]).+kxrange[1] + kys = (0:(nkY-1))/nkY*(kyrange[2]-kyrange[1]).+kyrange[1] + + sigma1 = inv(sigma) + influence_radius_2 = sqrt(sigma[1,1])^2+sqrt(sigma[2,2])^2 + + particles = PartStates[T,"particleList"] + nPart = (size(particles))[1] + + result = zeros(nkX,nkY) + + for i in 1:nPart + if (particles[i,"x"]-xyPoint[1])^2+(particles[i,"y"]-xyPoint[2])^2 < influence_radius_2*1.5 + for x in 1:nkX + for y in 1:nkY + result[x,y] += gaussMultiD([xyPoint[1],xyPoint[2], kxs[x], kys[y]], sigma, sigma1, [particles[i,"x"],particles[i,"y"],particles[i,"cx"],particles[i,"cy"]]) + end + end + end + end + return result +end + +function compute_spectrum_parameters(T, PartStates, xyPoint, sigma, aggregation_type) + all_particles = PartStates[T,"particleList"][:,["cx","cy","x","y"]] + + influence_radius = (sqrt(sigma[1,1])^2+sqrt(sigma[2,2])^2) + sigma = [influence_radius 0; 0 influence_radius] + sigma1 = inv(sigma) + result = DataFrame(cx=[],cy=[],x=[],y=[]) + weights = [] + + for i in 1:size(all_particles)[1] + temp_pos = Vector(all_particles[i,["x","y"]]) + dist = (temp_pos-xyPoint)[1]^2+(temp_pos-xyPoint)[2]^2 + if dist < influence_radius* (10)^2 + append!(result,DataFrame(all_particles[i,:])) + if aggregation_type == "gauss" + append!(weights, gaussMultiD(temp_pos, sigma, sigma1, xyPoint)) + elseif aggregation_type == "const" + append!(weights, 1) + elseif aggregation_type == "lin" + if abs(temp_pos[1]-xyPoint[1]) <= sqrt(sigma[1,1]) && abs(temp_pos[2]-xyPoint[2]) <= sqrt(sigma[2,2]) + append!(weights,(abs(temp_pos[1]-xyPoint[1])/sqrt(sigma[1,1])*abs(temp_pos[2]-xyPoint[2])/sqrt(sigma[2,2]))) + else + append!(weights, 0) + end + end + end + end + + if length(weights) == 0 + return [0 0; 0 0], 0 + end + sumWeights=sum(weights) + weights=weights./sumWeights + cBar_cx = sum(Matrix(weights.*result[:,["cx"]])) + cBar_cy = sum(Matrix(weights.*result[:,["cy"]])) + + cov_cxcx = 1/(1-sum(weights.^2))*sum(weights.*(Matrix(result[:,["cx"]]).-cBar_cx).^2) + cov_cycy = 1/(1-sum(weights.^2))*sum(weights.*(Matrix(result[:,["cy"]]).-cBar_cy).^2) + cov_cxcy = 1/(1-sum(weights.^2))*sum(weights.*(Matrix(result[:,["cx"]]).-cBar_cx).*(Matrix(result[:,["cy"]]).-cBar_cy)) + result_cov = [cov_cxcx cov_cxcy; cov_cxcy cov_cycy] + + return result_cov, sumWeights +end + +function compute_derivative(y, t, step) + results = zeros(length(t)-2*step) + times = t[(step+1):(end-step)] + for i in 1:(length(t)-2*step) + results[i] = (y[i+2*step]-y[i])/(t[i+2*step]-t[i]) + end + return times, results +end + +function spectrum_parameters_at_point(T, PartStates, xyPoint, influenceZone) + peak_speeds = Matrix{} +end + +function peak_spectrum(T, PartStates, xyPoint, sigma, sigma1) + particles = PartStates[T,"particleList"] + nPart = (size(particles))[1] + influence_radius_2 = (sigma[1,1])+(sigma[2,2]) + + value = [0; 0] + weights = 0 + for i in 1:nPart + if (particles[i,"x"]-xyPoint[1])^2+(particles[i,"y"]-xyPoint[2])^2 < influence_radius_2*(25)^2 + local_weight = gaussMultiD(xyPoint, sigma, sigma1, [particles[i, "x"]; particles[i, "y"]]) + value += [particles[i,"cx"]; particles[i,"cy"]] * local_weight + weights += local_weight + end + end + if weights == 0 + weights = 1 + end + return value/weights +end + +function find_geometrical_center(T, PartStates) + particles = PartStates[T,"particleList"] + return [Stc.mean(particles[:,"x"]), Stc.mean(particles[:,"y"])] +end + +function unit_dec(a) + unit=string(Int64(floor((a)))) + dec=string(Int64(floor(10*(a-floor((a)))))) + return unit*","*dec +end + +function plot_state_and_error_points(State, xs, ys, cRange, widthPlot, heightPlot; showTitle = true) + custom_color_scale = [(0.0,"rgb(229, 229, 229)"),(0.3,"rgb(255, 230, 44)"), (1.0, "rgb(255, 0, 0)")] + + traceq1 = heatmap(x=xs, y=ys, z=Matrix(State[:, :]), colorscale=custom_color_scale + ,colorbar=attr(tickformat=".1e",cmin = cRange[1], cmax = cRange[2]) + ,zmin = cRange[1], zmax = cRange[2] + ) + + LatexFont = attr(color="black",family="Computer Modern",size=22) + + if showTitle + q1=plot(traceq1, + Layout( + title=attr(text="\$ \\Large{\\text{Energy distribution in space}}\$", + x = 0.5, + xanchor="center", + font=attr(color="black")), + font=LatexFont, + width=widthPlot, height=heightPlot, + xaxis=attr(title=attr(text="\$\\Large{\\text{x position in } km}\$",font=attr(color="black"))), + yaxis=attr(title=attr(text="\$\\Large{\\text{y position in } km}\$",font=attr(color="black")),scaleanchor="x")) + ) + else + q1=plot(traceq1, + Layout( + font=LatexFont, + width=widthPlot, height=heightPlot, + xaxis=attr(title=attr(text="\$\\Large{\\text{x position in } km}\$",font=attr(color="black"))), + yaxis=attr(title=attr(text="\$\\Large{\\text{y position in } km}\$",font=attr(color="black")),scaleanchor="x")) + ) + end + return q1 +end + +function compute_translated_gaussian(mu_init, sigma_init, M, time) + translated_mu = M(time)*mu_init + translated_sigma = M(time)*transpose(sigma_init)*transpose(M(time)) + return translated_mu, translated_sigma +end + +function pearson(data, predicted) + meanValue = Stc.mean(data) + return 1-sum((data-predicted).^2)/sum((data.-meanValue).^2) +end + +## Reading the data +localpath = pwd()*"/.." + +# parentPath = localpath * "/plots/tests/paper/test_case_1/v2_pi" +# parentPath = localpath * "/plots/tests/paper/test_case_2_local/v2_pi_4" +# parentPath = localpath * "/plots/tests/paper/test_case_2/v2_pi_4" +# parentPath = localpath * "/plots/tests/paper/test_case_3/v2_half_cell" +# parentPath = localpath * "/plots/tests/paper/test_case_3/v2_half_cell" +# parentPath = localpath * "/plots/tests/paper/test_case_2/v2_pi_4_operationnel" +# parentPath = localpath * "/plots/tests/paper/test_debug/v2_half_cell" + +println("------------ Reading the data from : "* parentPath * " ------------") +path = parentPath*"/data/simu_infos.csv" +path3 = parentPath*"/data/sigma.csv" +data, header = readdlm(path, ',', header=true) +data3, header3 = readdlm(path3, ',', header=true) + +simu_infos = DataFrame(data, vec(header)) +init_sigma_df = DataFrame(data3, vec(header3)) + +init_sigma = Matrix(init_sigma_df) +Nx = simu_infos.Nx[1] +Ny = simu_infos.Ny[1] +xmin = simu_infos.xmin[1] +xmax = simu_infos.xmax[1] +ymin = simu_infos.ymin[1] +ymax = simu_infos.ymax[1] +lne_source = simu_infos.lne_source[1] +c_x_source = simu_infos.c_x_source[1] +c_y_source = simu_infos.c_y_source[1] +init_mu_C = [c_x_source, c_y_source] +x_source = simu_infos.x_source[1] +y_source = simu_infos.y_source[1] +init_mu_X = [x_source,y_source] +angular_spread_source = simu_infos.angular_spread_source[1] +Δt = simu_infos.Δt[1] +stop_time = simu_infos.stop_time[1] + +nTimes = Int64(ceil(stop_time/Δt))+1 + +# times = (1:(nTimes))*Δt +times = (0:(nTimes-1))*Δt +timesMinutes = times/60 + +timesMinutesString = unit_dec.(timesMinutes) + +ParticleStates = DataFrame([[],[]],["times", "particleList"]) +MeshState = DataFrame([[],[]],["times", "mesh_state"]) + +for i in 1:nTimes + path1 = parentPath*"/data/particles_"*timesMinutesString[i]*".csv" + path2 = parentPath*"/data/mesh_values_"*timesMinutesString[i]*".csv" + data1, header1 = readdlm(path1, ',', header=true) + data2, header2 = readdlm(path2, ',', header=true) + + read_df1 = DataFrame(data1, vec(header1)) + read_df2 = DataFrame(data2, vec(header2)) + + temp_df1 = DataFrame([[timesMinutes[i]], [read_df1]], ["times", "particleList"]) + temp_df2 = DataFrame([[timesMinutes[i]], [read_df2]], ["times", "mesh_state"]) + append!(ParticleStates, temp_df1) + append!(MeshState, temp_df2) +end + +## Building the spectrum at different points and different times + +nParticles = size(ParticleStates[1,"particleList"])[1] + +nkX = 15 +nkY = 15 + +kxscale = [-10, 30] +kyscale = [-10, 30] + +alpha = 5 +alpha = alpha /180 * pi +rotMatPos = [cos(alpha) -sin(alpha); sin(alpha) cos(alpha)] +rotMatNeg = [cos(-alpha) -sin(-alpha); sin(-alpha) cos(-alpha)] + +kxs = (0:(nkX-1))/nkX*(kxscale[2]-kxscale[1]).+kxscale[1] +kys = (0:(nkY-1))/nkY*(kyscale[2]-kyscale[1]).+kyscale[1] + +sigma_x = (xmax-xmin)/(Nx-1) +sigma_y = (ymax-ymin)/(Ny-1) +sigma_kx = (kxscale[2]-kxscale[1])/nkX +sigma_ky = (kyscale[2]-kyscale[1])/nkY + +kernel = [sigma_x^2 0 0 0; + 0 sigma_y^2 0 0; + 0 0 sigma_kx^2 0; + 0 0 0 sigma_ky^2] + +sdx_end = sqrt(Stc.var(ParticleStates[nTimes, "particleList"][:,"x"])) +sdy_end = sqrt(Stc.var(ParticleStates[nTimes, "particleList"][:,"y"])) + +xyCenter_1 = find_geometrical_center(1, ParticleStates) +xyCenter_2 = find_geometrical_center(Int64(floor(nTimes/2)), ParticleStates) +xyCenter_end = find_geometrical_center(nTimes, ParticleStates) +xy_right = rotMatNeg * xyCenter_end +xy_left = rotMatPos * xyCenter_end +xy_forw = xyCenter_end + 1 * [sdx_end,sdy_end] +xy_back = xyCenter_end - 1 * [sdx_end,sdy_end] +println() +println("------------------------------------------------") +println(" POSITIONS :") +println() + +println("Center at t=0 : (x,y) = ("*string(xyCenter_1[1])*","*string(xyCenter_1[2])*")") +println("Center at t="*timesMinutesString[nTimes]*" : (x,y) = ("*string(xyCenter_end[1])*","*string(xyCenter_end[2])*")") +println("Right at t="*timesMinutesString[nTimes]*" : (x,y) = ("*string(xy_right[1])*","*string(xy_right[2])*")") +println("Left at t="*timesMinutesString[nTimes]*" : (x,y) = ("*string(xy_left[1])*","*string(xy_left[2])*")") +println("Forward at t="*timesMinutesString[nTimes]*" : (x,y) = ("*string(xy_forw[1])*","*string(xy_forw[2])*")") +println("Backward at t="*timesMinutesString[nTimes]*" : (x,y) = ("*string(xy_back[1])*","*string(xy_back[2])*")") + +compute_spectra_bool = false + +if compute_spectra_bool + + println() + println("------------------------------------------------") + println(" COMPUTING SPECTRA :") + println() + + spectrum_1 = spectrum_at_point(1, ParticleStates, xyCenter_1, kernel, kxscale, kyscale, nkX, nkY) + println("Spectrum 1 over") + spectrum_2 = spectrum_at_point(Int64(floor(nTimes/2)), ParticleStates, xyCenter_2, kernel, kxscale, kyscale, nkX, nkY) + println("Spectrum 1 over") + spectrum_3 = spectrum_at_point(nTimes, ParticleStates, xyCenter_end, kernel, kxscale, kyscale, nkX, nkY) + println("Spectrum 2 over") + #spectrum_right = spectrum_at_point(nTimes, ParticleStates, xy_right, kernel, kxscale, kyscale, nkX, nkY) + spectrum_end = zeros(nkX, nkY) + spectrum_right = zeros(nkX, nkY) + println("Spectrum 3 over") + # spectrum_left = spectrum_at_point(nTimes, ParticleStates, xy_left, kernel, kxscale, kyscale, nkX, nkY) + spectrum_left = zeros(nkX, nkY) + println("Spectrum 4 over") + spectrum_forw = spectrum_at_point(nTimes, ParticleStates, xy_forw, kernel, kxscale, kyscale, nkX, nkY) + println("Spectrum 5 over") + spectrum_back = spectrum_at_point(nTimes, ParticleStates, xy_back, kernel, kxscale, kyscale, nkX, nkY) + println("Spectrum 6 over") + println() + + # FIRST PLOT + axis_template = attr(range = [-10,30], autorange = false, + showgrid = false, zeroline = false, + linecolor = "black", showticklabels = true, + ticks = "" ) + + sp1_trace = heatmap(x=kxs, y=kys, z=spectrum_1) + sp1_trace2 = scatter(x=[0, 0],y=[-20, 40], line=attr(color="black", width=3), showticklabels = false, ticks = "", name="") + sp1_trace3 = scatter(y=[0, 0],x=[-20, 40], line=attr(color="black", width=3), showticklabels = false, ticks = "", name="") + sp1 = plot([sp1_trace, sp1_trace2, sp1_trace3], Layout( + legend=:none, + title = "Energy Spectrum (in velocity) at t=0", + xaxis=axis_template, + yaxis=axis_template, + yaxis_title="c_y", + xaxis_title="c_x",width=1000, height=1000, + font=attr(size=28))) + + # SECOND PLOT + axis_template = attr(range = [-10,30], autorange = false, + showgrid = false, zeroline = false, + linecolor = "black", showticklabels = true, + ticks = "" ) + + sp2_trace = heatmap(x=kxs, y=kys, z=spectrum_2) + sp2_trace2 = scatter(x=[0, 0],y=[-20, 40], line=attr(color="black", width=3), showticklabels = false, ticks = "", name="") + sp2_trace3 = scatter(y=[0, 0],x=[-20, 40], line=attr(color="black", width=3), showticklabels = false, ticks = "", name="") + sp2 = plot([sp2_trace, sp2_trace2, sp2_trace3], Layout( + legend=:none, + title = "Energy Spectrum (in velocity) at t=0", + xaxis=axis_template, + yaxis=axis_template, + yaxis_title="c_y", + xaxis_title="c_x",width=1000, height=1000, + font=attr(size=28))) + + # THIRD PLOT + axis_template = attr(range = [-10,30], autorange = false, + showgrid = false, zeroline = false, + linecolor = "black", showticklabels = true, + ticks = "" ) + + sp3_trace = heatmap(x=kxs, y=kys, z=spectrum_3) + sp3_trace2 = scatter(x=[0, 0],y=[-20, 40], line=attr(color="black", width=3), showticklabels = false, ticks = "", name="") + sp3_trace3 = scatter(y=[0, 0],x=[-20, 40], line=attr(color="black", width=3), showticklabels = false, ticks = "", name="") + sp3 = plot([sp3_trace, sp3_trace2, sp3_trace3], Layout( + legend=:none, + title = "Energy Spectrum (in velocity) at t=0", + xaxis=axis_template, + yaxis=axis_template, + yaxis_title="c_y", + xaxis_title="c_x",width=1000, height=1000, + font=attr(size=28))) + + # NEXT PLOTS + p2 = plt.heatmap(kxs, kys, spectrum_end, aspect_ratio=:equal, size=(2160, 2500)) + plt.plot!(legend=:none, + title = "Energy Spectrum in velocity at t="*timesMinutesString[nTimes], + ylabel="c_y", + xlabel="c_x") + + p3 = plt.heatmap(kxs, kys, spectrum_right, aspect_ratio=:equal, size=(2160, 2500)) + plt.plot!(legend=:none, + title = "Energy Spectrum in velocity at t="*timesMinutesString[nTimes], + ylabel="c_y", + xlabel="c_x") + + p4 = plt.heatmap(kxs, kys, spectrum_left, aspect_ratio=:equal, size=(2160, 2500)) + plt.plot!(legend=:none, + title = "Energy Spectrum in velocity at t="*timesMinutesString[nTimes], + ylabel="c_y", + xlabel="c_x") + + p5 = plt.heatmap(kxs, kys, spectrum_forw, aspect_ratio=:equal, size=(2160, 2500)) + plt.plot!(legend=:none, + title = "Energy Spectrum in velocity at t="*timesMinutesString[nTimes], + ylabel="c_y", + xlabel="c_x") + + p6 = plt.heatmap(kxs, kys, spectrum_back, aspect_ratio=:equal, size=(2160, 2500)) + plt.plot!(legend=:none, + title = "Energy Spectrum in velocity at t="*timesMinutesString[nTimes], + ylabel="c_y", + xlabel="c_x") + + + max_speed_id = argmax(spectrum_1) + max_speed = [kxs[max_speed_id[1]],kys[max_speed_id[2]]] + + proba_spectrum_1 = spectrum_1 / sum(spectrum_1) + proba_spectrum_1[max_speed_id[1],max_speed_id[2]] += 1-sum(proba_spectrum_1) +end + +println("DONE") + + +xs = (0:(Nx-1))./(Nx-1)*(xmax-xmin).+xmin +ys = (0:(Ny-1))./(Ny-1)*(ymax-ymin).+ymin + +# Plotting test case 1 plots + +plot_test_case_1 = false + +if plot_test_case_1 + plotting_times = [2, Int64(floor(size(MeshState)[1]/2))+1, size(MeshState)[1]] + + halfMeshStates = DataFrame([[],[]], ["times", "mesh_state"]) + + + yid_to_plot = 1:(size(MeshState[1, "mesh_state"])[1]) + xid_to_plot = 1:Int64(floor(size(MeshState[1, "mesh_state"])[1]/2)) + for i in 1:length(plotting_times) + values = MeshState[plotting_times[i], "mesh_state"] + values = values[xid_to_plot,yid_to_plot] + df = DataFrame([[MeshState[plotting_times[i],"times"]], [values]], ["times", "mesh_state"]) + append!(halfMeshStates, df) + end + + cRange = [0, maximum(Matrix(halfMeshStates[2, "mesh_state"]))] + + tc1_1 = plot_state_and_error_points(halfMeshStates[1, "mesh_state"], xs[xid_to_plot]./1000, ys[yid_to_plot]./1000, cRange, 1130, 573) + tc1_2 = plot_state_and_error_points(halfMeshStates[2, "mesh_state"], xs[xid_to_plot]./1000, ys[yid_to_plot]./1000, cRange, 1130, 573) + tc1_3 = plot_state_and_error_points(halfMeshStates[3, "mesh_state"], xs[xid_to_plot]./1000, ys[yid_to_plot]./1000, cRange, 1130, 573) + + save_plots = true + if save_plots + savefig(tc1_1,"test_case_1_t1.html") + savefig(tc1_2,"test_case_1_t2.html") + savefig(tc1_3,"test_case_1_t3.html") + println("Plotting times for the test case 1 energy distribution : "*string(MeshState[plotting_times[1], "times"])*", "*string(MeshState[plotting_times[2], "times"])*", "*string(MeshState[plotting_times[3], "times"])) + end + + # Plotting the energy decay rate + + axis_template = attr(autorange = true, + showgrid = true, zeroline = false, + linecolor = "black", showticklabels = true, + ticks = "" ) + + LatexFont = attr(color="black",family="Computer Modern",size=22) + + ed_layout = Layout(yaxis_type="log",xaxis_type="log", + title=attr(text="Energy density evolution with distance",font=LatexFont) + ,width=800, height=800 + ,xaxis_title=L"\Large{\text{distance in } km}" + ,yaxis_title=L"\Large{\text{Energy density (arbitrary units)}}" + ,legend=attr(borderwidth=1,x=0.1, y=0.4,font=LatexFont) + ,font=LatexFont + ,xaxis=axis_template + ,yaxis=axis_template + ,margin=attr(l=90,r=80,t=80,b=80) + ,template="simple_white") + + initPos = [x_source; y_source] + + argmaxValues = [argmax(transpose(Matrix(MeshState[i,"mesh_state"][:,:]))) for i in 2:(size(MeshState)[1])] + max_energy = [MeshState[i,"mesh_state"][argmaxValues[i-1][2],argmaxValues[i-1][1]] for i in 2:(size(MeshState)[1])] + maxPos = [[xmin.+(argmaxValues[i][1]-1) .* sigma_x; ymin.+(argmaxValues[i][2]-1) .* sigma_y] for i in 1:length(argmaxValues)] + maxDist = [sqrt((maxPos[i][1]-x_source)^2+(maxPos[i][2]-y_source)^2) for i in 1:length(maxPos)] + + coeff = 300 + perfect_1_r = [coeff/distance for distance in maxDist[12:end]] + + trace_energy_decay = scatter(x=maxDist./1000,y=max_energy, + name="\$\\text{Energy density}\$", + mode="markers" + ) + trace_perfect_1_r = scatter(x=maxDist[12:end]./1000,y=perfect_1_r, + name="\$\\frac{1}{r} \\text{ curve for reference}\$" + ) + tc1_ed = plot([trace_energy_decay, trace_perfect_1_r], ed_layout) + + if save_plots + savefig(tc1_ed,"test_case_1_energy_decay.html") + end + +end + +# Plotting test case 2 plots + +plot_test_case_2 = true + +if plot_test_case_2 + # Plotting the energy distribution + + plotting_times = [2, Int64(floor(size(MeshState)[1]/2))+2, size(MeshState)[1]] + + plotLog = false + colorbarDelay = 0 + if plotLog + cRange = [-10, maximum(Matrix(log.(MeshState[plotting_times[1]+colorbarDelay, "mesh_state"])))] + + tc2_1 = plot_state_and_error_points(log.(MeshState[plotting_times[1], "mesh_state"]), xs./1000, ys./1000, cRange, 1000, 885) + tc2_2 = plot_state_and_error_points(log.(MeshState[plotting_times[2], "mesh_state"]), xs./1000, ys./1000, cRange, 1000, 885) + tc2_3 = plot_state_and_error_points(log.(MeshState[plotting_times[3], "mesh_state"]), xs./1000, ys./1000, cRange, 1000, 885) + else + cRange = [0, maximum(Matrix(MeshState[plotting_times[2]+colorbarDelay, "mesh_state"]))] + + tc2_1 = plot_state_and_error_points(MeshState[plotting_times[1], "mesh_state"], xs./1000, ys./1000, cRange, 1000, 885; showTitle=false) + tc2_2 = plot_state_and_error_points(MeshState[plotting_times[2], "mesh_state"], xs./1000, ys./1000, cRange, 1000, 885; showTitle=false) + tc2_3 = plot_state_and_error_points(MeshState[plotting_times[3], "mesh_state"], xs./1000, ys./1000, cRange, 1000, 885; showTitle=false) + end + + save_plots = true + if save_plots + savefig(tc2_1,"test_case_2_t1.html") + savefig(tc2_2,"test_case_2_t2.html") + savefig(tc2_3,"test_case_2_t3.html") + println("Plotting times for the test case 2 energy distribution : "*string(MeshState[plotting_times[1], "times"])*", "*string(MeshState[plotting_times[2], "times"])*", "*string(MeshState[plotting_times[3], "times"])) + end + + # Plotting the two analysis plots + + xaxis_template = attr(#autorange = true, + showgrid = true, zeroline = false, + linecolor = "black", showticklabels = true, + ticks = "", + range=[xmin,xmax]./1000 + ) + + yaxis_template = attr(#autorange = true, + showgrid = true, zeroline = false, + linecolor = "black", showticklabels = true, + ticks = "", + range=[ymin,ymax]./1000 + ) + + LatexFont = attr(color="black",family="Computer Modern",size=22) + + ed_layout = Layout(#yaxis_type="log",xaxis_type="log", + title=attr(text="Particle positions at t=60sec",font=LatexFont) + ,width=1600, height=374 + ,xaxis_title=attr(text="\$\\Large{\\text{x position in } km}\$",font=attr(color="black")) + ,yaxis_title=attr(text="\$\\Large{\\text{y position in } km}\$",font=attr(color="black")) + ,legend=attr(borderwidth=1,x=0.1, y=0.4,font=LatexFont) + ,font=LatexFont + ,xaxis=xaxis_template + ,yaxis=yaxis_template + ,margin=attr(l=90,r=80,t=80,b=80) + ,template="simple_white") + + + cov_end = Stc.cov(Matrix(ParticleStates[1,"particleList"][:,["cx","cy","x","y"]])) + + # q1 = plot_state_and_error_points(MeshState[nTimes, "mesh_state"], xs, ys, xmin, xmax, ymin, ymax) + println("ici 1") + # particleSpeeds = sqrt.([sum((Matrix(ParticleStates[nTimes,"particleList"][:,["cx","cy"]]).^2)[i,:]) for i in 1:nParticles]) + particleSpeeds=zeros(nParticles) + plotting_time = Int64(floor(nTimes/2)) + for i in 1:nParticles + particleSpeeds[i] = sqrt(ParticleStates[plotting_time,"particleList"][:,"cx"][i]^2+ParticleStates[plotting_time,"particleList"][:,"cy"][i]^2) + end + println("ici 2") + particleXs = ParticleStates[plotting_time,"particleList"][:,"x"] + particleYs = ParticleStates[plotting_time,"particleList"][:,"y"] + M(t) = [Diagonal(ones(2)) Diagonal(zeros(2)); + t*Diagonal(ones(2)) Diagonal(ones(2))] + println("ici 3") + result_mu, result_sigma = compute_translated_gaussian([c_x_source, c_y_source, x_source, y_source], init_sigma, M, 60*timesMinutes[plotting_time]) + println("ici 4") + timeOfInterest = Int64(floor(nTimes/2)) + translated_mu_C = [c_x_source, c_y_source] + translated_mu_XY = [x_source, y_source] + timesMinutes[timeOfInterest]*60 * [c_x_source, c_y_source] + translated_sigma_XY = init_sigma[3:4,3:4] + (timesMinutes[timeOfInterest]*60)^2 * init_sigma[1:2,1:2] + translated_cross_cov = timesMinutes[timeOfInterest]*60 * init_sigma[1:2,1:2] + # translated_sigma_C = init_sigma[1:2,1:2] + # translated_cov_total=[translated_sigma_C translated_cross_cov; + # transpose(translated_cross_cov) translated_sigma_XY] + # println("ici 5") + spectrum(x,y) = gaussMultiD(translated_mu_XY,translated_sigma_XY,inv(translated_sigma_XY),[x,y]) + println("ici 6") + z = @. spectrum(xs', ys) + nPartSubset = 20000 + subsetPartID = Int64.(floor.(rand(nPartSubset)*length(particleXs))).+1 + color1=[255 0 0] + color2=[255 230 44] + minSpeed=minimum(particleSpeeds) + maxSpeed=maximum(particleSpeeds) + colorPart=color1.*((particleSpeeds[subsetPartID].-minSpeed)./maxSpeed).+color2.*(1 .-(particleSpeeds[subsetPartID].-minSpeed)./maxSpeed) + transparencyPoints = 1 + trace_tc_2_1 = scatter(x=particleXs[subsetPartID]./1000,y=particleYs[subsetPartID]./1000, + mode="markers", + marker=attr(size=3, + color=particleSpeeds[subsetPartID], + # color="rgba(".*string.(Int64.(floor.(colorPart[:,1]))).*", ".*string.(Int64.(floor.(colorPart[:,2]))).*", ".*string.(Int64.(floor.(colorPart[:,3]))).*", ".*string(transparencyPoints).*")", + showscale=true, + colorscale=[[0,"rgb(255, 230, 44)"],[1,"rgb(255, 0, 0)"]] + ), + zorder=0) + trace_tc_2_2 = contour(x=xs./1000,y=ys./1000,z=z,contours_coloring="lines",line_width=4,showscale=false,colorscale=[[0,"black"],[0.5,"black"],[1,"black"]],zorder=1) + + q1 = plot([trace_tc_2_1,trace_tc_2_2], ed_layout) + + if save_plots + savefig(q1,"test_case_2_particle_contour.html") + end + + line_plot_nPoints = 350 + line_plot_xy = zeros(line_plot_nPoints,2) + line_plot_xy[:,1] = (0:(line_plot_nPoints-1))./(line_plot_nPoints-1)*(xmax-xmin).+xmin + line_plot_xy[:,2] = (y_source) * ones(line_plot_nPoints) + line_plot_length = [sqrt((line_plot_xy[i,1]-line_plot_xy[1,1])^2+(line_plot_xy[i,2]-line_plot_xy[1,2])^2) for i in 1:line_plot_nPoints] + kernel1 = inv(kernel) + println("Start measuring peak_spectrum") + + # Compute model output curve + line_plot_vector = [peak_spectrum(timeOfInterest, ParticleStates, line_plot_xy[i,:],kernel[1:2,1:2],kernel1[1:2,1:2]) for i in 1:line_plot_nPoints] + # line_plot_values = [sqrt(sum(line_plot_vector[i].^2)) for i in 1:line_plot_nPoints] + line_plot_values = [line_plot_vector[i][1] for i in 1:line_plot_nPoints] + println("Finished measuring peak_spectrum") + + # Compute theoretical curve + init_mu_C_measured = [Stc.mean(Matrix(ParticleStates[1,"particleList"][:,["cx"]])), Stc.mean(Matrix(ParticleStates[1,"particleList"][:,["cy"]]))] + # init_mu_C_measured = init_mu_C + init_mu_X_measured = [Stc.mean(Matrix(ParticleStates[3,"particleList"][:,["x"]])), Stc.mean(Matrix(ParticleStates[3,"particleList"][:,["y"]]))] - Δt*2*init_mu_C_measured + init_sigma_X_measured = Stc.cov(Matrix(ParticleStates[3,"particleList"][:,["x","y"]])) + # init_sigma_X_measured = init_sigma[3:4,3:4] + # init_mu_X_measured = init_mu_X + translated_sigma_XY1 = inv(translated_sigma_XY) + translated_mu_XY = init_mu_X_measured+ timesMinutes[timeOfInterest]*60 * [c_x_source, c_y_source] + translated_cross_cov = timesMinutes[timeOfInterest]*60 * init_sigma[1:2,1:2] + line_plot_theory_vector = [translated_mu_C+translated_cross_cov*translated_sigma_XY1*(line_plot_xy[i,:]-translated_mu_XY) for i in 1:line_plot_nPoints] + # line_plot_theory_values = [sqrt(sum(line_plot_theory_vector[i].^2)) for i in 1:line_plot_nPoints] + line_plot_theory_values = [line_plot_theory_vector[i][1] for i in 1:line_plot_nPoints] + + # Compute diffused curve + coef = 0.5 + init_mu_C_measured = [Stc.mean(Matrix(ParticleStates[1,"particleList"][:,["cx"]])), Stc.mean(Matrix(ParticleStates[1,"particleList"][:,["cy"]]))] + # init_mu_C_measured = init_mu_C + init_mu_X_measured = [Stc.mean(Matrix(ParticleStates[3,"particleList"][:,["x"]])), Stc.mean(Matrix(ParticleStates[3,"particleList"][:,["y"]]))] - Δt*2*init_mu_C_measured + init_sigma_X_measured = Stc.cov(Matrix(ParticleStates[3,"particleList"][:,["x","y"]])) + init_sigma_C_measured = Stc.cov(Matrix(ParticleStates[1,"particleList"][:,["cx","cy"]])) + # init_sigma_X_measured = init_sigma[3:4,3:4] + # init_mu_X_measured = init_mu_X + D = coef*sigma_x * coef*sigma_x / Δt * I[1:2,1:2] + translated_diff_sigma_XY = init_sigma_X_measured + (timesMinutes[timeOfInterest]*60)^2 * init_sigma_C_measured + timesMinutes[timeOfInterest]*60 * D + translated_diff_sigma_XY1 = inv(translated_diff_sigma_XY) + translated_cross_cov = timesMinutes[timeOfInterest]*60 * init_sigma_C_measured# + (D) * timesMinutes[timeOfInterest]*60/Δt + translated_mu_XY = init_mu_X_measured + timesMinutes[timeOfInterest]*60 * init_mu_C_measured + line_plot_diff_theory_vector = [translated_mu_C+translated_cross_cov*translated_diff_sigma_XY1*(line_plot_xy[i,:]-translated_mu_XY) for i in 1:line_plot_nPoints] + # line_plot_diff_theory_vector = [translated_mu_C+translated_cross_cov*translated_diff_sigma_XY1*1.09*(line_plot_xy[i,:]-translated_mu_XY) for i in 1:line_plot_nPoints] + #line_plot_diff_theory_values = [sqrt(sum(line_plot_diff_theory_vector[i].^2)) for i in 1:line_plot_nPoints] + line_plot_diff_theory_values = [line_plot_diff_theory_vector[i][1] for i in 1:line_plot_nPoints] + + # Compute linereg through the data and plot it + startPlot = 1 + peakSpeedDataFrame = DataFrame(x=line_plot_length[startPlot:end]./1000, y=line_plot_values[startPlot:end]) + peak_speed_linreg = lm(@formula(y~x), peakSpeedDataFrame) + predict_peak_speed_lm = predict(peak_speed_linreg) + + axis_template = attr(autorange = true, + showgrid = true, zeroline = false, + linecolor = "black", showticklabels = true, + ticks = "" + ) + + LatexFont = attr(color="black",family="Computer Modern",size=28) + LatexFont2 = attr(color="black",family="Computer Modern",size=22) + + psd_layout = Layout(title=attr(text=L"\Large{\text{Peak speed along the axis of propagation}}",font=LatexFont) + ,width=900, height=900 + ,xaxis_title=attr(text=L"\Large{\text{Position along the axis of propagation in } km}",font=LatexFont) + ,yaxis_title=attr(text=L"\Large{\text{Peak speed of particles in } m/s}",font=LatexFont) + ,legend=attr(borderwidth=1,x=0.1, y=0.9,font=LatexFont2) + ,font=LatexFont + ,xaxis=axis_template + ,yaxis=axis_template + ,margin=attr(l=90,r=80,t=80,b=100) + ,template="simple_white" + ) + + + q2 = plot([scatter(x=line_plot_length./1000,y=line_plot_values,name="Model output"), + # scatter(x=line_plot_length./1000,y=line_plot_theory_values, name="Theory without diffusion"), + scatter(x=line_plot_length./1000,y=line_plot_diff_theory_values, name="Theory with diffusion"), + scatter(x=line_plot_length[startPlot:end]./1000,y=predict_peak_speed_lm, name="Linear regression through the data")],psd_layout + ) + + if save_plots + savefig(q2,"test_case_2_peak_speed_distribution.html") + end + +end + +## Computing statistical parameters + +compute_spectrum_bool = false + +if compute_spectrum_bool + counts_1 = Int64.(round.(proba_spectrum_1*10000)) + nSample = sum(counts_1) + sample_1=Array{Any,1}() + + for i in 1:nkX + for j in 1:nkY + for k in 1:counts_1[i,j] + push!(sample_1, [kxs[i],kys[j]]) + end + end + end + + proba_spectrum_end = spectrum_end / sum(spectrum_end) + proba_spectrum_end[max_speed_id[1],max_speed_id[2]] += 1-sum(proba_spectrum_end) + + counts_end = Int64.(round.(proba_spectrum_end*10000)) + nSample = sum(counts_end) + sample_end=Array{Any,1}() + + for i in 1:nkX + for j in 1:nkY + for k in 1:counts_end[i,j] + push!(sample_end, [kxs[i],kys[j]]) + end + end + end + + # proba_spectrum_right = spectrum_right / sum(spectrum_right) + # proba_spectrum_right[max_speed_id[1],max_speed_id[2]] += 1-sum(proba_spectrum_right) + + # counts_right = Int64.(round.(proba_spectrum_right*10000)) + # nSample = sum(counts_right) + # sample_right=Array{Any,1}() + + # for i in 1:nkX + # for j in 1:nkY + # for k in 1:counts_right[i,j] + # push!(sample_right, [kxs[i],kys[j]]) + # end + # end + # end + + # proba_spectrum_left = spectrum_left / sum(spectrum_left) + # proba_spectrum_left[max_speed_id[1],max_speed_id[2]] += 1-sum(proba_spectrum_left) + + # counts_left = Int64.(round.(proba_spectrum_left*10000)) + # nSample = sum(counts_left) + # sample_left=Array{Any,1}() + + # for i in 1:nkX + # for j in 1:nkY + # for k in 1:counts_left[i,j] + # push!(sample_left, [kxs[i],kys[j]]) + # end + # end + # end + + proba_spectrum_forw = spectrum_forw / sum(spectrum_forw) + proba_spectrum_forw[max_speed_id[1],max_speed_id[2]] += 1-sum(proba_spectrum_forw) + + counts_forw = Int64.(round.(proba_spectrum_forw*10000)) + nSample = sum(counts_forw) + sample_forw=Array{Any,1}() + + for i in 1:nkX + for j in 1:nkY + for k in 1:counts_forw[i,j] + push!(sample_forw, [kxs[i],kys[j]]) + end + end + end + + proba_spectrum_back = spectrum_back / sum(spectrum_back) + proba_spectrum_back[max_speed_id[1],max_speed_id[2]] += 1-sum(proba_spectrum_back) + + counts_back = Int64.(round.(proba_spectrum_back*10000)) + nSample = sum(counts_back) + sample_back=Array{Any,1}() + + for i in 1:nkX + for j in 1:nkY + for k in 1:counts_back[i,j] + push!(sample_back, [kxs[i],kys[j]]) + end + end + end + + std_1 = Stc.cov(sample_1) + mean_1 = [Stc.mean([sample_1[i][j] for i in 1:length(sample_1)]) for j in 1:2] + std_end = Stc.cov(sample_end) + mean_end = [Stc.mean([sample_end[i][j] for i in 1:length(sample_end)]) for j in 1:2] + # std_right = Stc.cov(sample_right) + # mean_right = [Stc.mean([sample_right[i][j] for i in 1:length(sample_right)]) for j in 1:2] + # std_left = Stc.cov(sample_left) + # mean_left = [Stc.mean([sample_left[i][j] for i in 1:length(sample_left)]) for j in 1:2] + std_forw = Stc.cov(sample_forw) + mean_forw = [Stc.mean([sample_forw[i][j] for i in 1:length(sample_forw)]) for j in 1:2] + std_back = Stc.cov(sample_back) + mean_back = [Stc.mean([sample_back[i][j] for i in 1:length(sample_back)]) for j in 1:2] + + + + println() + println("------------------------------------------------") + println(" RESULTS :") + println() + + ## T=0 + print("At t=0, the center energy is at (") + print(xyCenter_1[1]) + print(", ") + print(xyCenter_1[2]) + println(")") + print("The mean of the spectrum at this point is (") + print(mean_1[1]) + print(", ") + print(mean_1[2]) + println(")") + println("And the standard deviation is :") + print_mat(std_1) + println() + println() + + ## T=end -> Center + print("At t=60, the center energy is at (") + print(xyCenter_end[1]) + print(", ") + print(xyCenter_end[2]) + println(")") + print("The mean of the spectrum at this point is (") + print(mean_end[1]) + print(", ") + print(mean_end[2]) + println(")") + println("And the standard deviation is :") + print_mat(std_end) + println() + println() + + ## T=end -> Right + # print("At t=60, the right point is at (") + # print(xy_right[1]) + # print(", ") + # print(xy_right[2]) + # println(")") + # print("The mean of the spectrum at this point is (") + # print(mean_right[1]) + # print(", ") + # print(mean_right[2]) + # println(")") + # println("And the standard deviation is :") + # print_mat(std_right) + # println() + # println() + + # ## T=end -> Left + # print("At t=60, the center energy is at (") + # print(xy_left[1]) + # print(", ") + # print(xy_left[2]) + # println(")") + # print("The mean of the spectrum at this point is (") + # print(mean_left[1]) + # print(", ") + # print(mean_left[2]) + # println(")") + # println("And the standard deviation is :") + # print_mat(std_left) + # println() + # println() + + ## T=end -> Forward + print("At t=60, the center energy is at (") + print(xy_forw[1]) + print(", ") + print(xy_forw[2]) + println(")") + print("The mean of the spectrum at this point is (") + print(mean_forw[1]) + print(", ") + print(mean_forw[2]) + println(")") + println("And the standard deviation is :") + print_mat(std_forw) + println() + println() + + ## T=end -> Backward + print("At t=60, the center energy is at (") + print(xy_back[1]) + print(", ") + print(xy_back[2]) + println(")") + print("The mean of the spectrum at this point is (") + print(mean_back[1]) + print(", ") + print(mean_back[2]) + println(")") + println("And the standard deviation is :") + print_mat(std_back) + println() + println() +end + +plot_test_case_3 = true + +if plot_test_case_3 + sigmasComputeType = "multiplePoints" # "singlePoint" or "multiplePoints" + recompute_sigmas_bool = true + timesSeconds = 60*timesMinutes + + if recompute_sigmas_bool + nSamplePoints = 30 + + sigmas_cx = zeros(nTimes) + sigmas_cy = zeros(nTimes) + sigmas_cxy = zeros(nTimes) + sigmas_cx_1 = zeros(nTimes) + sigmas_cy_1 = zeros(nTimes) + sigmas_cxy_1 = zeros(nTimes) + sigmas_C = [zeros(2,2) for _ in 1:nTimes] + sigmas_C_1 = [zeros(2,2) for _ in 1:nTimes] + + for i in 1:nTimes + timeOfInterest = i + translated_mu_C = [c_x_source, c_y_source] + translated_mu_XY = [x_source, y_source] + timesMinutes[timeOfInterest]*60 * [c_x_source, c_y_source] + translated_sigma_XY = init_sigma[3:4,3:4] + (timesMinutes[timeOfInterest]*60)^2 * init_sigma[1:2,1:2] + translated_sigma_C = init_sigma[1:2,1:2] + if sigmasComputeType=="multiplePoints" # Version 2 + xy_center = translated_mu_XY + xy_sigma = translated_sigma_XY + # xy_center = timesSeconds[i].*init_mu_C+init_mu_X + # xy_sigma = timesSeconds[i]^2 .* init_sigma[1:2,1:2] + init_sigma[3:4,3:4] + d = MvNormal(xy_center, xy_sigma) + sampleXYs = rand(d,nSamplePoints) + sampleValues = [] + sampleWeights = [] + # print("computing for sample point number ") + for j in 1:nSamplePoints + # print(", "*string(j)) + res, weight = compute_spectrum_parameters(i, ParticleStates, sampleXYs[:,j], kernel, "gauss") + if res != [0 0; 0 0] + append!(sampleValues, [res]) + append!(sampleWeights, [weight]) + else + print("FOUND ZERO; ") + end + end + sigmas_C[i] = sum(sampleWeights.*sampleValues)/sum(sampleWeights) + elseif sigmasComputeType=="singlePoint" # Version 1 + xy_center = translated_mu_XY + xy_sigma = translated_sigma_XY + xy_chosen = xy_center + sqrt.([translated_sigma_XY[1,1], translated_sigma_XY[2,2]]) + sigmas_C[i] = compute_spectrum_parameters(i, ParticleStates, xy_chosen, kernel, "gauss")[1] + end + + sigmas_C_1[i] = inv(sigmas_C[i]) + + sigmas_cx[i] = sigmas_C[i][1,1] + sigmas_cy[i] = sigmas_C[i][2,2] + sigmas_cxy[i] = sigmas_C[i][2,1] + + sigmas_cx_1[i] = sigmas_C_1[i][1,1] + sigmas_cy_1[i] = sigmas_C_1[i][2,2] + sigmas_cxy_1[i] = sigmas_C_1[i][2,1] + if i%10==0 || i==nTimes + println("Computing sigma for t = "*string(timesMinutes[i])) + end + end + end + + time_long_index = 32 + + data_sigma_cx = DataFrame(x=timesSeconds[time_long_index:end], y=sigmas_cx[time_long_index:end].^-1) + ab_sigma_cx = lm(@formula(y~x), data_sigma_cx) + data_sigma_cy = DataFrame(x=timesSeconds[time_long_index:end], y=sigmas_cy[time_long_index:end].^-1) + ab_sigma_cy = lm(@formula(y~x), data_sigma_cy) + + pearson_sigma_cx = pearson(sigmas_cx[time_long_index:end].^-1, predict(ab_sigma_cx)) + pearson_sigma_cy = pearson(sigmas_cy[time_long_index:end].^-1, predict(ab_sigma_cy)) + + trace1 = scatter(x=timesSeconds,y=sigmas_cx.^-1,name="sigma c_x") + trace2 = scatter(x=timesSeconds,y=sigmas_cy.^-1,name="sigma c_y") + trace3 = scatter(x=timesSeconds[time_long_index:end],y=predict(ab_sigma_cx),name="linreg through sigma c_x; R2="*string(round(pearson_sigma_cx,digits=3))) + trace4 = scatter(x=timesSeconds[time_long_index:end],y=predict(ab_sigma_cy),name="linreg through sigma c_y; R2="*string(round(pearson_sigma_cy,digits=3))) + + trace5 = scatter(x=timesSeconds,y=sigmas_cx.^1,name=L"\text{Model output }(\sigma_{c_x})^2", mode="markers", marker=attr(color="blue")) + trace6 = scatter(x=timesSeconds,y=sigmas_cy.^1,name=L"\text{Model output }(\sigma_{c_y})^2", mode="markers", marker=attr(color="orange")) + + init_sigma_1 = inv(init_sigma) + + init_sigma_C_1 = init_sigma_1[1:2,1:2] + # init_sigma_C_1 = 1 ./(init_sigma[1:2,1:2]) + init_sigma_X_1 = init_sigma_1[3:4,3:4] + # init_sigma_X_1 = 1 ./(init_sigma[3:4,3:4]) + + sigmas_C_t = [inv(t^2*init_sigma_X_1^1 .+ init_sigma_C_1) for t in timesSeconds] + sigmas_cx_t = [sigmas_C_t[t][1,1] for t in 1:length(timesSeconds)] + sigmas_cy_t = [sigmas_C_t[t][2,2] for t in 1:length(timesSeconds)] + + trace7 = scatter(x=timesSeconds,y=sigmas_cx_t, name=L"\text{Theoretical } (\sigma_{c_x})^2", mode="lines", line=attr(dash="dash", color="blue")) + trace8 = scatter(x=timesSeconds,y=sigmas_cy_t, name=L"\text{Theoretical } (\sigma_{c_y})^2", mode="lines", line=attr(dash="dash", color="orange")) + + log_sigma_cx = DataFrame(x=log.(timesSeconds[time_long_index:end]), y=log.(sigmas_cx[time_long_index:end])) + ab_log_sigma_cx = lm(@formula(y~x), log_sigma_cx) + r2_log_sigma_cx = pearson(log_sigma_cx.y,predict(ab_log_sigma_cx)) + log_sigma_cy = DataFrame(x=log.(timesSeconds[time_long_index:end]), y=log.(sigmas_cy[time_long_index:end])) + ab_log_sigma_cy = lm(@formula(y~x), log_sigma_cy) + r2_log_sigma_cy = pearson(log_sigma_cy.y,predict(ab_log_sigma_cy)) + + trace9 = scatter(x=timesSeconds[time_long_index:end], y=exp.(predict(ab_log_sigma_cx))) + trace10 = scatter(x=timesSeconds[time_long_index:end], y=exp.(predict(ab_log_sigma_cy))) + + trace11 = scatter(x=timesSeconds[time_long_index:end], y=20000000 .* timesSeconds[time_long_index:end].^-2) + + # TESTING OTHER CURVES + + # Version 1 + # alpha_mu_Cx = init_mu_C[1]/(sigma_x/Δt) - floor(init_mu_C[1]/(sigma_x/Δt)) + # beta_sigma_Cx = sqrt(init_sigma[1,1])/(sigma_x/Δt) - floor(sqrt(init_sigma[1,1])/(sigma_x/Δt)) + # alpha_mu_Cy = init_mu_C[2]/(sigma_y/Δt) - floor(init_mu_C[2]/(sigma_y/Δt)) + # beta_sigma_Cy = sqrt(init_sigma[2,2])/(sigma_y/Δt) - floor(sqrt(init_sigma[2,2])/(sigma_y/Δt)) + + # coefx = (alpha_mu_Cx-(alpha_mu_Cx^2+beta_sigma_Cx^2)) + # coefy = (alpha_mu_Cy*(1-alpha_mu_Cy/2)+beta_sigma_Cy*(sqrt(2/pi)-beta_sigma_Cy/2)) + # D = sigma_x*sigma_x / Δt * [coefx 0; 0 coefy] + + # Version 2 + # init_mu_cx_corrected = (init_mu_C[1]*Δt/sigma_x - floor.(init_mu_C[1]*Δt/sigma_x))/Δt*sigma_x + # init_mu_cy_corrected = (init_mu_C[2]*Δt/sigma_y - floor.(init_mu_C[2]*Δt/sigma_y))/Δt*sigma_y + # init_sigma_cx_corrected = sqrt(init_sigma[1,1]) + # init_sigma_cy_corrected = sqrt(init_sigma[2,2]) + # # init_sigma_cx_corrected = ((sqrt(init_sigma[1,1])*Δt/sigma_x - floor.(sqrt(init_sigma[1,1])*Δt/sigma_x))/Δt*sigma_x)^2 + # # init_sigma_cy_corrected = ((sqrt(init_sigma[2,2])*Δt/sigma_y - floor.(sqrt(init_sigma[2,2])*Δt/sigma_y))/Δt*sigma_y)^2 + # localx = -init_mu_cx_corrected/(init_sigma_cx_corrected) + # localy = -init_mu_cy_corrected/(init_sigma_cy_corrected) + # phi_x = 1/sqrt(2*pi)*exp(-0.5*(localx)^2) + # phi_y = 1/sqrt(2*pi)*exp(-0.5*(localy)^2) + # Phi_x = 0.5*(1+erf(localx/sqrt(2))) + # Phi_y = 0.5*(1+erf(localy/sqrt(2))) + + # mu_c_upper_x = init_mu_cx_corrected + ((init_sigma_cx_corrected) * phi_x)/(1-Phi_x) + # mu_c_lower_x = init_mu_cx_corrected - ((init_sigma_cx_corrected) * phi_x)/(Phi_x) + # sigma_c_upper_x = init_sigma_cx_corrected * sqrt(1 + (localx* phi_x)/(1 - Phi_x) - (phi_x/(1 - Phi_x))^2) + # sigma_c_lower_x = init_sigma_cx_corrected * sqrt(1 - (localx* phi_x)/(Phi_x) - (phi_x/(Phi_x))^2) + + # mu_c_upper_y = init_mu_cy_corrected + (sqrt(init_sigma_cy_corrected) * phi_y)/(1-Phi_y) + # mu_c_lower_y = init_mu_cy_corrected - (sqrt(init_sigma_cy_corrected) * phi_y)/(Phi_y) + # sigma_c_upper_y = sqrt(init_sigma_cy_corrected * (1 - (localy* phi_y)/(1 - Phi_y) - (phi_y/(1-Phi_y))^2)) + # sigma_c_lower_y = sqrt(init_sigma_cy_corrected * (1 + (localy* phi_y)/(Phi_y) - (phi_y/(Phi_y))^2)) + + # coefx = (1-Phi_x)*(sigma_x*mu_c_upper_x-Δt*(sigma_c_upper_x^2+mu_c_upper_x^2))+Phi_x*(-sigma_x*mu_c_lower_x-Δt*(sigma_c_lower_x^2+mu_c_lower_x^2)) + # coefy = (1-Phi_y)*(sigma_y*mu_c_upper_y-Δt*(sigma_c_upper_y^2+mu_c_upper_y^2))+Phi_y*(-sigma_y*mu_c_lower_y-Δt*(sigma_c_lower_y^2+mu_c_lower_y^2)) + # D = [coefx 0; 0 coefy] + + # Version 3 + # Phi = (x,mu,sigma) -> 0.5*(1+erf((x-mu)/(sigma*sqrt(2)))) + # phi = (x,mu,sigma) -> 1/(sqrt(2*pi))*exp(-0.5*(x-mu)^2/(sigma^2)) + # p=0.0025 + + + # cmin_x = init_mu_C[1] + sqrt(init_sigma[1,1]*2)*erfinv(2*p-1) + # cmax_x = init_mu_C[1] + sqrt(init_sigma[1,1]*2)*erfinv(2*(1-p)-1) + # k_x = Int64(floor(cmin_x*Δt/sigma_x)) + # l_x = Int64(floor(cmax_x*Δt/sigma_x)) + # pos_cx = [(i)*sigma_x/Δt for i in k_x:(l_x+1)] + + # indexes_x = (k_x:l_x) + # probs_x = [(Phi(pos_cx[i+1], init_mu_C[1], sqrt(init_sigma[1,1]))-Phi(pos_cx[i], init_mu_C[1], sqrt(init_sigma[1,1]))) for i in 1:(l_x-k_x+1)] + # means_x = [tn.tnmean(pos_cx[i],pos_cx[i+1],init_mu_C[1],sqrt(init_sigma[1,1])) for i in 1:(l_x-k_x+1)] + # vars_x = [tn.tnvar(pos_cx[i],pos_cx[i+1],init_mu_C[1],sqrt(init_sigma[1,1])) for i in 1:(l_x-k_x+1)] + # coefs_x = probs_x.*( + # sigma_x .* (means_x - indexes_x.*(sigma_x/Δt)) + # - Δt .* ( + # vars_x .^ 2 + means_x .^ 2 + # - 2 .* (indexes_x) .* (Δt / sigma_x) .* means_x + # + ((indexes_x) * sigma_x / Δt) .^ 2 + # ) + # ) + # coefx = sum(coefs_x) + + + # cmin_y = init_mu_C[2] + sqrt(init_sigma[2,2]*2)*erfinv(2*p-1) + # cmax_y = init_mu_C[2] + sqrt(init_sigma[2,2]*2)*erfinv(2*(1-p)-1) + # k_y = Int64(floor(cmin_y*Δt/sigma_y)) + # l_y = Int64(floor(cmax_y*Δt/sigma_y)) + # pos_cy = [(i)*sigma_y/Δt for i in k_y:(l_y+1)] + + # probs_y = [(Phi(pos_cy[i+1], init_mu_C[2], sqrt(init_sigma[2,2]))-Phi(pos_cy[i], init_mu_C[2], sqrt(init_sigma[2,2]))) for i in 1:(l_y-k_y+1)] + # means_y = [tn.tnmean(pos_cy[i],pos_cy[i+1],init_mu_C[2],sqrt(init_sigma[2,2])) for i in 1:(l_y-k_y+1)] + # # means_y = (means_y.*Δt./sigma_y - floor.(means_y.*Δt./sigma_y))./Δt.*sigma_y + # vars_y = [tn.tnvar(pos_cy[i],pos_cy[i+1],init_mu_C[2],sqrt(init_sigma[2,2])) for i in 1:(l_y-k_y+1)] + # coefs_y = 1 + # coefy = sum(coefs_y) + + # # coefx = 884.5 + # # coefy = 599.96 + # D = [coefx 0; 0 coefy] + + # Version 4 + p=1e-5 + nIntPoints = 100000 + f = (c) -> c - floor(c) + g = (c) -> c^2 - 2 * c * floor(c) + floor(c)^2 + compute_Dc = (c,DeltaX,DeltaT) -> DeltaX * f(c/DeltaX*DeltaT) * (DeltaX/DeltaT) - DeltaT * g(c*DeltaT/DeltaX) * (DeltaX/DeltaT)^2 + + cmin_x = init_mu_C[1] + sqrt(init_sigma[1,1]*2)*erfinv(2*p-1) + cmax_x = init_mu_C[1] + sqrt(init_sigma[1,1]*2)*erfinv(2*(1-p)-1) + pos_cx = (((0:nIntPoints)./nIntPoints).*(cmax_x-cmin_x)).+cmin_x + DeltaCx = pos_cx[2]-pos_cx[1] + probas_cx = gaussD.(init_mu_C[1],init_sigma[1,1],pos_cx) + D_cx_values = compute_Dc.(pos_cx, sigma_x, Δt) + D_cx_integrand = probas_cx .* D_cx_values + coefx = sum((D_cx_integrand[1:(end-1)]+D_cx_integrand[2:(end)])/2)*DeltaCx + + cmin_y = init_mu_C[2] + sqrt(init_sigma[2,2]*2)*erfinv(2*p-1) + cmax_y = init_mu_C[2] + sqrt(init_sigma[2,2]*2)*erfinv(2*(1-p)-1) + pos_cy = (((0:nIntPoints)./nIntPoints).*(cmax_y-cmin_y)).+cmin_y + DeltaCy = pos_cy[2]-pos_cy[1] + probas_cy = gaussD.(init_mu_C[2],init_sigma[2,2],pos_cy) + D_cy_values = compute_Dc.(pos_cy, sigma_y, Δt) + D_cy_integrand = probas_cy .* D_cy_values + coefy = sum((D_cy_integrand[1:(end-1)]+D_cy_integrand[2:(end)])/2)*DeltaCy + + D = [coefx 0 ; 0 coefy] + + + timesSeconds2 = (0:(nTimes-1))./(nTimes-1).*(timesSeconds[end]) + + sigmas_C_t_new = [inv(t^2 .* inv(t .* D+ init_sigma[3:4,3:4]) .+ init_sigma_C_1) for t in timesSeconds2] + sigmas_cx_t_new = [sigmas_C_t_new[t][1,1] for t in 1:length(timesSeconds2)] + sigmas_cy_t_new = [sigmas_C_t_new[t][2,2] for t in 1:length(timesSeconds2)] + + derivSigmas_D_C_t_th = [(log.(sigmas_C_t_new[t]).-log.(sigmas_C_t_new[t-1]))./(log(timesSeconds2[t]).-log(timesSeconds2[t-1])) for t in 2:length(timesSeconds2)] + derivSigmas_D_cx_t_th = [derivSigmas_D_C_t_th[t][1,1] for t in 1:(length(timesSeconds2)-1)] + derivSigmas_D_cy_t_th = [derivSigmas_D_C_t_th[t][2,2] for t in 1:(length(timesSeconds2)-1)] + + newTimes, derivSigmas_D_cx_t = compute_derivative(log.(sigmas_cx),log.(timesSeconds), 1) + newTimes, derivSigmas_D_cy_t = compute_derivative(log.(sigmas_cy),log.(timesSeconds), 1) + newTimes = exp.(newTimes) + + trace12 = scatter(x=timesSeconds2, y=sigmas_cx_t_new, name=L"\text{Diffused } (\sigma_{c_x})^2", marker=attr(color="rgb(0, 98, 255)")) + trace13 = scatter(x=timesSeconds2, y=sigmas_cy_t_new, name=L"\text{Diffused } (\sigma_{c_y})^2", marker=attr(color="rgb(226, 45, 0)")) + + trace14 = scatter(x=timesSeconds2[2:end], y=derivSigmas_D_cx_t_th) + trace15 = scatter(x=timesSeconds2[2:end], y=derivSigmas_D_cy_t_th) + + trace16 = scatter(x=newTimes, y=derivSigmas_D_cx_t) + trace17 = scatter(x=newTimes, y=derivSigmas_D_cy_t) + + α = 1000 + β = 30000000 + γ = 17500000 + perfect_1_r = [α/t for t in timesSeconds2[10:end]] + perfect_1_r_2 = [β/t^2 for t in timesSeconds2[10:end]] + perfect_1_r_2_2 = [γ/t^2 for t in timesSeconds2[30:end]] + + trace_perfect_1_r = scatter(x=timesSeconds2[10:end],y=perfect_1_r, + name="\$\\frac{1}{t} \\text{ curve for reference}\$" + , mode="lines", line=attr(dash="dot", color="rgb(112, 187, 112)",width=4) + ) + trace_perfect_1_r_2 = scatter(x=timesSeconds2[10:end],y=perfect_1_r_2, + name="\$\\frac{1}{t^2} \\text{ curve for reference}\$" + , mode="lines", line=attr(dash="dot", color="rgb(181, 105, 178)",width=4) + ) + trace_perfect_1_r_2_2 = scatter(x=timesSeconds2[30:end],y=perfect_1_r_2_2, + name="\$\\frac{1}{t^2} \\text{ curve for reference}\$" + , mode="lines", line=attr(dash="dot", color="rgb(181, 105, 178)",width=4) + ) + + + # END TESTING OTHER CURVES + + axis_template = attr(autorange = true, + showgrid = true, zeroline = false, + linecolor = "black", showticklabels = true, + ticks = "" + ) + + axis_template_x_2 = attr(showgrid = true, zeroline = false, + linecolor = "black", showticklabels = true, + ticks = "", range = log10.([2500, 7500]) + ) + + axis_template_y_2 = attr(showgrid = true, zeroline = false, + linecolor = "black", showticklabels = true, + ticks = "", range = log10.([0.12, 1.34]) + ) + + LatexFont = attr(color="black",family="Computer Modern",size=28) + LatexFontlegend = attr(color="black",family="Computer Modern",size=40) + + sn_layout = Layout(yaxis_type="log",xaxis_type="log", + title=attr(text=L"\LARGE{\text{Evolution of the local } C_x \text{ energy spectrum variance in time}}",font=LatexFont) + ,width=1600, height=1000 + ,xaxis_title=L"\LARGE{\text{time } (sec)}" + ,yaxis_title=L"\LARGE{(\sigma_{c})^2}" + #,legend=attr(borderwidth=1,x=0.5, y=0.5,font=LatexFont) + ,legend=attr(borderwidth=1,x=0.05, y=0.05,font=LatexFontlegend) + ,font=LatexFont + ,xaxis=axis_template + ,yaxis=axis_template + ,margin=attr(l=90,r=80,t=80,b=140) + ,template="simple_white" + ) + + sn_layout2 = Layout(yaxis_type="log",xaxis_type="log", + title=attr(text=L"\LARGE{\text{Evolution of the local } C_y \text{ energy spectrum variance in time}}",font=LatexFont) + ,width=1600, height=1000 + ,xaxis_title=L"\LARGE{\text{time } (sec)}" + ,yaxis_title=L"\LARGE{(\sigma_{c})^2}" + #,legend=attr(borderwidth=1,x=0.5, y=0.5,font=LatexFont) + ,legend=attr(borderwidth=1,x=0.05, y=0.05,font=LatexFontlegend) + ,font=LatexFont + ,xaxis=axis_template + ,yaxis=axis_template + ,margin=attr(l=90,r=80,t=80,b=140) + ,template="simple_white" + ) + + tc3_1 = plot([trace5, trace7, trace12, trace_perfect_1_r, trace_perfect_1_r_2 + ],sn_layout + ) + + tc3_2 = plot([trace6, trace8, trace13, trace_perfect_1_r, trace_perfect_1_r_2_2 + ],sn_layout2 + ) + + save_plots = true + if save_plots + savefig(tc3_1,"test_case_2_spectral_peak_narrowing_cx.html") + savefig(tc3_2,"test_case_2_spectral_peak_narrowing_cy.html") + # savefig(tc3_3,"test_case_2_t3.html") + end + + tc3_3 = plot([trace14, trace16, trace17] + ,Layout(xaxis_type="log" + ,title="LogSlope of the spectrum variance in time" + ) + ) + + tc3_4 = plot([trace1, trace2, trace3, trace4],Layout(#yaxis_type="log",xaxis_type="log", + title="Evolution of the local energy spectrum variance in time" + ,width=1000, height=1000 + ,xaxis_title="\\text{time} (sec)" + ,yaxis_title="sigma_c" + ,legend=attr(borderwidth=3,x=0.1, y=0.8,font=attr(size=18)) + ,font=attr(size=20) + ,template="simple_white" + ) + ) + + # Plotting the space distribution + + plotting_times = [2, Int64(floor(size(MeshState)[1]/2))+2, size(MeshState)[1]] + + plotLog = false + colorbarDelay = 20 + width = 930 + height = 300 + showTitle = false + if plotLog + cRange = [-10, maximum(Matrix(log.(MeshState[plotting_times[1]+colorbarDelay, "mesh_state"])))] + + tc3_t1 = plot_state_and_error_points(log.(MeshState[plotting_times[1], "mesh_state"]), xs./1000, ys./1000, cRange, width, height; showTitle) + tc3_t2 = plot_state_and_error_points(log.(MeshState[plotting_times[2], "mesh_state"]), xs./1000, ys./1000, cRange, width, height; showTitle) + tc3_t3 = plot_state_and_error_points(log.(MeshState[plotting_times[3], "mesh_state"]), xs./1000, ys./1000, cRange, width, height; showTitle) + else + cRange = [0, maximum(Matrix(MeshState[plotting_times[2]+colorbarDelay, "mesh_state"]))] + + tc3_t1 = plot_state_and_error_points(MeshState[plotting_times[1], "mesh_state"], xs./1000, ys./1000, cRange, width, height; showTitle) + tc3_t2 = plot_state_and_error_points(MeshState[plotting_times[2], "mesh_state"], xs./1000, ys./1000, cRange, width, height; showTitle) + tc3_t3 = plot_state_and_error_points(MeshState[plotting_times[3], "mesh_state"], xs./1000, ys./1000, cRange, width, height; showTitle) + end + + save_plots = true + if save_plots + savefig(tc3_t1,"test_case_3_t1.html") + savefig(tc3_t2,"test_case_3_t2.html") + savefig(tc3_t3,"test_case_3_t3.html") + println("Plotting times for the test case 3 energy distribution : "*string(MeshState[plotting_times[1], "times"])*", "*string(MeshState[plotting_times[2], "times"])*", "*string(MeshState[plotting_times[3], "times"])) + end + +end \ No newline at end of file diff --git a/src/Architectures.jl b/src/Architectures.jl index f1596a9..b177856 100644 --- a/src/Architectures.jl +++ b/src/Architectures.jl @@ -4,7 +4,7 @@ using SharedArrays using StaticArrays -export AbstractGrid, AbstractODESettings, AbstractParticleInstance, AbstractMarkedParticleInstance, Abstract1DModel, Abstract2DModel, AbstractModel, AbstractStore, AbstractParticleSystem, StateTypeL1, IDConstantsInstance, ScgConstantsInstance, CartesianGrid, CartesianGrid1D, CartesianGrid2D, TripolarGrid, Grid2D, MeshGrids, MeshGridStatistics +export AbstractGrid, AbstractODESettings, AbstractParticleInstance, AbstractMarkedParticleInstance, Abstract1DModel, Abstract2DModel, Abstract2DStochasticModel, AbstractModel, AbstractStore, AbstractParticleSystem, StateTypeL1, IDConstantsInstance, ScgConstantsInstance, CartesianGrid, CartesianGrid1D, CartesianGrid2D, TripolarGrid, Grid2D, MeshGrids, MeshGridStatistics, AbstractStochasticParticleInstance export StandardRegular1D_old, StandardRegular2D_old @@ -39,6 +39,7 @@ MeshGridStatistics = Union{CartesianGridStatistics,TripolarGridStatistics,Spheri abstract type AbstractODESettings end abstract type AbstractParticleInstance end +abstract type AbstractStochasticParticleInstance <: AbstractParticleInstance end abstract type AbstractMarkedParticleInstance end abstract type IDConstantsInstance end @@ -53,6 +54,7 @@ abstract type AbstractModel{TS} end abstract type Abstract1DModel <: AbstractModel{Nothing} end abstract type Abstract2DModel <: AbstractModel{Nothing} end +abstract type Abstract2DStochasticModel <: Abstract2DModel end #All posiible types of a single-layer StateVectors StateTypeL1 = Union{SharedArray{Float64,3},MArray} diff --git a/src/Models/GeometricalOpticsModels.jl b/src/Models/GeometricalOpticsModels.jl new file mode 100644 index 0000000..1cd791f --- /dev/null +++ b/src/Models/GeometricalOpticsModels.jl @@ -0,0 +1,318 @@ +module GeometricalOpticsModels + +export GeometricalOptics, init_particles! +export fields + +using ...Architectures + +using ModelingToolkit: get_states, ODESystem + +#using core_1D: MarkedParticleInstance +using ...ParticleMesh: OneDGrid, OneDGridNotes, TwoDGrid, TwoDGridNotes + +using ...Operators.core_2D_spread: ParticleDefaults as ParticleDefaults2D +#using core_2D: SeedParticle! as SeedParticle2D! +using ...Operators.mapping_2D + +using SharedArrays +# using DistributedArrays +using Printf + +import Oceananigans: fields +using Oceananigans.TimeSteppers: Clock +using ...FetchRelations + +using LinearAlgebra + +using ..WaveGrowthModels2D +# include("WaveGrowthModels2D.jl") + +# TO DO : check if all the modules imported above are necessary, for now I just copied the ones from WaveGrowthModels2D.jl + + + +""" + GeometricalOptics{Grid, Lay, Tim, Clo, stat, PC, Ovar, Osys, Oses, Odev, bl_flag, bl, wnds, cur} + + ADD HERE A DESCRIPTION OF THE MODEL + +""" +mutable struct GeometricalOptics{Grid<:AbstractGrid, + Lay, + Tim, + Clo, + Int, + stat, + Pan, + PCollection, + PPool, + FPC, + Ovar, + Osys, + Oses, + Odev, + MinPar, + MinStat, + bl_flag, + bl, + bl_type, + as_type, + as_thresh, + pb_cov_init, + plt_steps, + plt_path, + sve_part, + nPart, + wnds, + cur, + Mstat} <: Abstract2DStochasticModel where {Mstat<:Union{Nothing,stat}, PCollection<:Union{Vector,Array}} + #Union{Vector,DArray} + grid::Grid + layers::Lay # number of layers used in the model, 1 is eneough + timestepper::Tim # silly Oceananigans + clock::Clo + dims::Int # number of dimensions + + State::stat # state of of the model at each grid point, for each layer it contains, energy, positions, group speed and directional spreading + ParticlesAtNode::Pan # list of the particles to regrid at each node + ParticleCollection::PCollection # Collection (list) of Particles + ParticlePool::PPool # Pool of particles (used for non parametric mode) + FailedCollection::FPC # Collection (list) of Particles that failed to integrate + + ODEvars::Ovar # list of variables in ODE system, have type Num from OrdinaryDiffEq and ModelingToolkit + ODEsystem::Osys # the ODE system used at each particle + ODEsettings::Oses # All setting needed to solve the ODEsystem + ODEdefaults::Odev # Dict{NUm, Float64} ODE defaults + minimal_particle::MinPar + minimal_state::MinStat + + periodic_boundary::bl_flag # If true we use a period boundary + boundary::bl # List of boundary points + boundary_defaults::bl_type # Dict{NUm, Float64} ODE defaults + + angular_spreading_thresh::as_thresh # threshold for splitting + angular_spreading_type::as_type # type of angular spreading used + proba_covariance_init::pb_cov_init # initial probability density 4x4 covariance matrix of particles + + plot_steps::plt_steps # Whether or not to plot the steps + plot_savepath::plt_path # The path to use for saving plots + + save_particles::sve_part # Whether or not to save particle data + + n_particles_launch::nPart # Number of particles to launch from each node when using a nonparametric spreading type + + winds::wnds # u, v, if needed u_x, u_y + currents::cur # u, v, currents + + MovieState::Mstat # state of of the model. Only used for producing movieframes + +end + + + + +## 2D version + + +""" +GeometricalOptics(; grid, winds, ODEsys, ODEvars, layers, clock, ODEsets, ODEdefaults, currents, periodic_boundary, CBsets) + This is the constructor for the GeometricalOptics model. The inputs are: + grid : the grid used in the model, + winds : the wind interpolation function here only 1D, + ODEsys : the ODE system used in the model, + ODEvars : the variables in the ODE system, + layers : the number of layers used in the model (default 1), + clock : the clock used in the model, + ODEsets : the ODE settings (type: ODESettings), + ODEdefaults : the ODE defaults (type: ParticleDefaults), + currents : the currents (not implimented yet), + periodic_boundary: if true we use a periodic boundary (default true), + CBsets : the callback settings (not implimented yet). +""" +function GeometricalOptics(; grid::GG, + winds::NamedTuple{(:u, :v)}, + ODEsys, + ODEvars=nothing, #needed for MTK for ODEsystem. will be depriciated later + layers::Int=1, + clock=Clock{eltype(grid)}(;time=0, last_Δt=0, last_stage_Δt=1), + ODEsets::AbstractODESettings=nothing, # ODE_settings + ODEinit_type::PP= "wind_sea", # default_ODE_parameters + minimal_particle=nothing, # minimum particle the model falls back to if a particle fails to integrate + minimal_state=nothing, # minimum state threshold needed for the state to be advanced + currents=nothing, # + periodic_boundary=true, + boundary_type="same", # or "minimal", "same", default is same, only used if periodic_boundary is false + angular_spreading_thresh=π/8, + angular_spreading_type="stochast", # or "geometrical" or "nonparametric" + proba_covariance_init = Matrix(I,4,4)*10^(-50), + plot_steps=false, + plot_savepath="", + save_particles=false, + n_particles_launch=150, + CBsets=nothing, + movie=false) where {PP<:Union{ParticleDefaults2D,String},GG<:AbstractGrid} + + # initialize state {SharedArray} given grid and layers + # Number of state variables + Nstate = 4 + State = init_StateArray(grid, Nstate, layers) + # if layers > 1 + # State = SharedArray{Float64,4}(grid.stats.Nx.N, grid.stats.Ny.N, Nstate, layers) + # else + # State = SharedArray{Float64,3}(grid.stats.Nx.N, grid.stats.Ny.N, Nstate) + # end + + if layers > 1 + ParticlesAtNode = Array{Array{Array{Array{Any,1},1},1},1}() + for i in 1:grid.stats.Nx.N + push!(ParticlesAtNode, []) + for j in 1:grid.stats.Ny.N + push!(ParticlesAtNode[i], []) + for _ in 1:layers + push!(ParticlesAtNode[i][j], []) + end + end + end + ParticlePool = Array{Array{Any,1},1}() + else + ParticlesAtNode = Array{Array{Array{Any,1},1},1}() + for i in 1:grid.stats.Nx.N + push!(ParticlesAtNode, []) + for _ in 1:grid.stats.Ny.N + push!(ParticlesAtNode[i], []) + end + end + ParticlePool = Array{Any,1}() + end + + if ODEinit_type isa ParticleDefaults2D + ODEdefaults = ODEinit_type + elseif ODEinit_type == "wind_sea" + ODEdefaults = nothing + elseif ODEinit_type == "mininmal" + ODEdefaults = ParticleDefaults2D(-11.0, 1e-3, 0.0) + else + error("ODEinit_type must be either 'wind_sea','mininmal', or ParticleDefaults2D instance ") + end + + if isnothing(minimal_particle) + @info "initalize minimum particle" + minimal_particle = FetchRelations.MinimalParticle(2, 2, ODEsets.timestep) + else + @info "use minimal particle" + end + + if isnothing(minimal_state) + @info "initalize minimum state" + minimal_state = FetchRelations.MinimalState(2, 2, ODEsets.timestep) + else + @info "use minimal state" + end + + # initliaze boundary points (periodic_boundary) + if ~periodic_boundary # if false, define boundary points here: + boundary = boundary = mark_boundary(grid) + else + boundary = [] + end + + if boundary_type == "wind_sea" + boundary_defaults = nothing + @info "use wind_sea boundary" + elseif boundary_type == "mininmal" + @info "use 'mininmal' boundary (1min with 2m/s)" + #FetchRelations.get_minimal_windsea(u(0, 0), ODEsets.DT) + WindSeamin = FetchRelations.get_minimal_windsea(1, 1, 5*60) # 5 min with 2 m/s + #WindSeamin = FetchRelations.get_minimal_windsea(u(0, 0, 0), v(0, 0, 0), DT / 2) + #WindSeamin = FetchRelations.get_initial_windsea(u(0, 0, 0), v(0, 0, 0), DT/5) + lne_local = log(WindSeamin["E"]) + cg_u_local = WindSeamin["cg_bar_x"] + cg_v_local = WindSeamin["cg_bar_y"] + boundary_defaults = copy(ParticleDefaults2D(lne_local, cg_u_local, cg_v_local, 0.0, 0.0)) + + elseif boundary_type == "same" + @info "use default value boundary" + boundary_defaults = ODEdefaults + else + error("boundary_type must be either 'wind_sea','mininmal', or 'same' ") + end + + ParticleCollection = [] + FailedCollection = Vector{AbstractMarkedParticleInstance}([]) + # particle initialization is not done in the init_particle! method + # for i in range(1,length = grid.Nx) + # SeedParticle!(ParticleCollection, State, i, + # ODEsys, ODEdev, ODEsets, + # gridnotes, winds, ODEsets.timestep, grid.Nx, boundary, periodic_boundary) + # end + + # check dimensions of all components: + # check Nstate ODEvars + if movie + Mstat = State + else + Mstat = nothing + end + + # return WaveGrowth2D structure + return GeometricalOptics( + grid, + layers, + nothing, + clock, # ??? + 2, # This is a 2D model + State, + ParticlesAtNode, + ParticleCollection, + ParticlePool, + FailedCollection, + ODEvars, + ODEsys, + ODEsets, + ODEdefaults, + minimal_particle, + minimal_state, + periodic_boundary, + boundary, boundary_defaults, + angular_spreading_thresh, + angular_spreading_type, + proba_covariance_init, + plot_steps, + plot_savepath, + save_particles, + n_particles_launch, + winds, + currents, + Mstat) +end + + + + +function Base.show(io::IO, ow::GeometricalOptics) + + if ow.ODEsystem isa ODESystem + sys_print = get_states(ow.ODEsystem) + else + sys_print = ow.ODEsystem + end + print(io, "GeometricalOptics ", "\n", + "├── grid: ", ow.grid, "\n", + "├── layers: ", ow.layers, "\n", + "├── clock: ", ow.clock, + "├── State: ", size(ow.State), "\n", + "├── ParticleCollection size: ", length(ow.ParticleCollection), "\n", + "├── ODEs \n", + "| ├── System: ", sys_print, "\n", + "| ├── Defaults: ", ow.ODEdefaults, "\n", + "| └── Settings: \n", ow.ODEsettings, "\n", + "├── winds ", ow.winds, "\n", + "├── currents ", ow.currents, "\n", + "└── Perdiodic Boundary ", ow.periodic_boundary, "\n") +end + + + + +# end of module +end \ No newline at end of file diff --git a/src/Models/Models.jl b/src/Models/Models.jl index daa8a40..b77c47a 100644 --- a/src/Models/Models.jl +++ b/src/Models/Models.jl @@ -1,11 +1,13 @@ module Models -export WaveGrowthModels1D, WaveGrowthModels2D, reset_boundary! +export WaveGrowthModels1D, WaveGrowthModels2D, GeometricalOpticsModels, reset_boundary! include("WaveGrowthModels1D.jl") -include("WaveGrowthModels2D.jl") - using .WaveGrowthModels1D +include("WaveGrowthModels2D.jl") using .WaveGrowthModels2D +include("GeometricalOpticsModels.jl") +using .GeometricalOpticsModels + end \ No newline at end of file diff --git a/src/Models/WaveGrowthModels2D.jl b/src/Models/WaveGrowthModels2D.jl index bc06fc0..b1dd20f 100644 --- a/src/Models/WaveGrowthModels2D.jl +++ b/src/Models/WaveGrowthModels2D.jl @@ -1,6 +1,6 @@ module WaveGrowthModels2D -export WaveGrowth2D, init_particles! +export WaveGrowth2D, init_particles!, init_StateArray, mark_boundary, reset_boundary!, fields export fields using ...Architectures @@ -352,11 +352,11 @@ end Return a flattened `NamedTuple` of the State vector for a `WaveGrowth2D` model. """ -fields(model::WaveGrowth2D) = (State=model.State,) +fields(model::Abstract2DModel) = (State=model.State,) # # Oceananigans.Simulations interface # fields(m::ContinuumIceModel) = merge(m.velocities, m.stresses) -function reset_boundary!(model::WaveGrowth2D) +function reset_boundary!(model::Abstract2DModel) if model.periodic_boundary # if false, define boundary points here: boundary = [] diff --git a/src/Operators/Operators.jl b/src/Operators/Operators.jl index d3a36bc..403a29f 100644 --- a/src/Operators/Operators.jl +++ b/src/Operators/Operators.jl @@ -1,6 +1,6 @@ module Operators -export core_1D, core_2D, custom_structures, mapping_1D, mapping_2D, TimeSteppers +export core_1D, core_2D, core_2D_spread, custom_structures, mapping_1D, mapping_2D, TimeSteppers export init_z0_to_State! using SharedArrays using StaticArrays @@ -24,6 +24,7 @@ include("utils.jl") include("initialize.jl") include("core_1D.jl") include("core_2D.jl") +include("core_2D_spread.jl") include("mapping_1D.jl") include("mapping_2D.jl") @@ -34,6 +35,7 @@ include("TimeSteppers.jl") using .core_1D using .core_2D +using .core_2D_spread using .mapping_1D using .mapping_2D using .TimeSteppers diff --git a/src/Operators/TimeSteppers.jl b/src/Operators/TimeSteppers.jl index 9ccc9bb..10402cb 100644 --- a/src/Operators/TimeSteppers.jl +++ b/src/Operators/TimeSteppers.jl @@ -10,8 +10,14 @@ using Statistics using Base.Threads using Printf +using Random, Distributions +using ..core_2D_spread: ParticleDefaults, InitParticleInstance +using ...custom_structures: ParticleInstance2D + using Oceananigans.TimeSteppers: tick! +using DataFrames, CSV, Dates + function mean_of_state(model::Abstract2DModel) return mean(model.State[:, :, 1]) end @@ -36,6 +42,34 @@ end ################# 1D #################### +function write_particles_to_csv(wave_model, Δt) + sec=string(Int64(floor((wave_model.clock.time+Δt)/60))) + dec=string(Int64(floor(10*((wave_model.clock.time+Δt)/60-floor((wave_model.clock.time+Δt)/60))))) + save_path = wave_model.plot_savepath + + nParticles = wave_model.n_particles_launch + @info wave_model.clock.time/60 + + parts = wave_model.ParticleCollection[(end-nParticles+1):end] + + logE = zeros(nParticles) + cx = zeros(nParticles) + cy = zeros(nParticles) + x = zeros(nParticles) + y = zeros(nParticles) + for i in 1:nParticles + logE[i] = parts[i].ODEIntegrator[1] + cx[i] = parts[i].ODEIntegrator[2] + cy[i] = parts[i].ODEIntegrator[3] + x[i] = parts[i].ODEIntegrator[4] + y[i] = parts[i].ODEIntegrator[5] + end + + data = DataFrame(id=1:nParticles, logE = logE, cx = cx, cy = cy, x = x, y = y) + data2 = Tables.table(transpose(wave_model.State[:, :, 1])) + CSV.write(save_path*"/data/particles_"*sec*","*dec*".csv", data) + CSV.write(save_path*"/data/mesh_values_"*sec*","*dec*".csv", data2) +end """ time_step!(model, Δt; callbacks=nothing) @@ -106,7 +140,7 @@ clock is ticked by Δt callbacks are not implimented yet """ -function time_step!(model::Abstract2DModel, Δt::Float64; callbacks=nothing, debug=false) +function time_step!(model::Abstract2DStochasticModel, Δt::Float64; callbacks=nothing, debug=false) # temporary FailedCollection to store failed particles FailedCollection = Vector{AbstractMarkedParticleInstance}([]) @@ -117,6 +151,199 @@ function time_step!(model::Abstract2DModel, Δt::Float64; callbacks=nothing, deb @info maximum(model.State[:, :, 1]), maximum(model.State[:, :, 2]), maximum(model.State[:, :, 3]) model.FailedCollection = FailedCollection end + for (i,j) in [(i,j) for i in 1:model.grid.stats.Nx.N for j in 1:model.grid.stats.Ny.N] + for k in 1:length(model.ParticlesAtNode[i][j]) + pop!(model.ParticlesAtNode[i][j]) + end + end + + @threads for a_particle in model.ParticleCollection + #@info a_particle.position_ij + mapping_2D.advance!( a_particle, model.ParticlesAtNode, model.State, FailedCollection, + model.grid, model.winds, Δt, + model.ODEsettings.log_energy_maximum, + model.ODEsettings.wind_min_squared, + model.periodic_boundary, + model.ODEdefaults) + end + + if model.save_particles && length(model.ParticleCollection) > 0 + write_particles_to_csv(model, Δt) + end + + for (i,j) in [(i,j) for i in 1:model.grid.stats.Nx.N for j in 1:model.grid.stats.Ny.N] + weights = [model.ParticlesAtNode[i][j][k][1] for k in 1:length(model.ParticlesAtNode[i][j])] + values = [model.ParticlesAtNode[i][j][k][2] for k in 1:length(model.ParticlesAtNode[i][j])] + if length(weights) > 0 + model.State[i,j,4] = sum(weights .* values) / sum(weights) + end + end + + if debug + print("mean energy after advance ", mean_of_state(model), "\n") + + @info "advanced: " + @info maximum(model.State[:, :, 1]), maximum(model.State[:, :, 2]), maximum(model.State[:, :, 3]) + #@info model.State[8:12, 1], model.State[8:12, 2] + @info model.clock.time, model.ParticleCollection[10].ODEIntegrator.t + @info "winds:", model.winds.u(model.ParticleCollection[10].ODEIntegrator.u[4], model.ParticleCollection[10].ODEIntegrator.u[5], model.ParticleCollection[10].ODEIntegrator.t) + end + + #@printf "re-mesh" + if model.angular_spreading_type == "stochast" + @threads for a_particle in model.ParticleCollection + mapping_2D.remesh!(a_particle, model.State, + model.winds, model.clock.time, + model.ODEsettings, Δt, + model.minimal_particle, + model.minimal_state, + model.ODEdefaults) + end + elseif model.angular_spreading_type == "geometrical" + i = 1 + particlesToBeReset = [] + + for _ in 1:length(model.ParticleCollection) + pos_ij = model.ParticleCollection[i].position_ij + big_enough = model.State[pos_ij[1], pos_ij[2],2]^2+model.State[pos_ij[1], pos_ij[2],3]^2>model.minimal_state[2] + if model.ParticleCollection[i].boundary || ~big_enough + i=i+1 + elseif !(model.ParticleCollection[i].position_ij in particlesToBeReset) + push!(particlesToBeReset, model.ParticleCollection[i].position_ij) + deleteat!(model.ParticleCollection, i) + end + end + @threads for a_particle in model.ParticleCollection + mapping_2D.remesh!(a_particle, model.State, + model.winds, model.clock.time, + model.ODEsettings, Δt, + model.minimal_particle, + model.minimal_state, + model.ODEdefaults) + end + @threads for (i,j) in particlesToBeReset + mapping_2D.remesh!(i, j, model, Δt) + end + + elseif model.angular_spreading_type == "nonparametric" + + # This next part needs to change + + nPreviousParticles = 0 + model.ParticlePool = [] + + @threads for (i,j) in [(i,j) for i in 1:model.grid.stats.Nx.N for j in 1:model.grid.stats.Ny.N] + locallen = length(model.ParticlesAtNode[i][j]) + nPreviousParticles += locallen + @threads for k in 1:locallen + + Upart=model.ParticlesAtNode[i][j][k][3].ODEIntegrator.u + # newodeint=init(model.ParticlesAtNode[i][j][k][3].ODEIntegrator.t,[deepcopy(model.ParticlesAtNode[i][j][k][3].ODEIntegrator.u[1]), deepcopy(model.ParticlesAtNode[i][j][k][3].ODEIntegrator.u[2]), deepcopy(model.ParticlesAtNode[i][j][k][3].ODEIntegrator.u[3]), deepcopy(model.ParticlesAtNode[i][j][k][3].ODEIntegrator.u[4]), deepcopy(model.ParticlesAtNode[i][j][k][3].ODEIntegrator.u[5]), deepcopy(model.ParticlesAtNode[i][j][k][3].ODEIntegrator.u[6])]) + # println(newodeint) + # println() + + z_init = ParticleDefaults(Upart[1],Upart[2],Upart[3],Upart[4],Upart[5],Upart[6]) + newpart=InitParticleInstance(model.ODEsystem,z_init,model.ODEsettings,(i,j),false,true) + push!(model.ParticlePool, [newpart, i, j, exp(model.ParticlesAtNode[i][j][k][3].ODEIntegrator.u[1])*model.ParticlesAtNode[i][j][k][1]]) + + #newpart2=[deepcopy(model.ParticlesAtNode[i][j][k][3]), i, j, exp(model.ParticlesAtNode[i][j][k][3].ODEIntegrator.u[1])*model.ParticlesAtNode[i][j][k][1]] + #push!(model.ParticlePool, newpart2) + end + end + + ParticuleEnergiesNorm = sum([model.ParticlePool[k][4] for k in 1:length(model.ParticlePool)]) + for k in 1:length(model.ParticlePool) + model.ParticlePool[k][4] = model.ParticlePool[k][4] / ParticuleEnergiesNorm + end + + energiesNormed = [model.ParticlePool[k][4] for k in 1:length(model.ParticlePool)] + a = Categorical(energiesNormed) + nNewParticles = model.n_particles_launch + particlesDrawn = rand(a, nNewParticles) + #new_energy_tot = sum([energiesNormed[k] for k in particlesDrawn]) + #energy_factor = totEnergyDomain[1] / new_energy_tot + + i = 1 + j = 0 + counter = 0 + for _ in 1:length(model.ParticleCollection) + if model.ParticleCollection[i].on + if counter >= nNewParticles + deleteat!(model.ParticleCollection, i) + j=j+1 + else + i=i+1 + end + counter +=1 + else + i=i+1 + end + end + #@info "removed particles :", j + #@info "new energy :", ParticuleEnergiesNorm + + #@info "particlesDrawn", model.ParticlePool[particlesDrawn[k]][1].ODEIntegrator.u[2] + + for k in 1:nNewParticles + i = model.ParticlePool[particlesDrawn[k]][2] + j = model.ParticlePool[particlesDrawn[k]][3] + x = model.grid.stats.xmin + model.grid.stats.dx*(i-1) + y = model.grid.stats.ymin + model.grid.stats.dy*(j-1) + #log_energy = log(totEnergyDomain[1]/nPreviousParticles) + log_energy = log(ParticuleEnergiesNorm/nNewParticles) + c_x = model.ParticlePool[particlesDrawn[k]][1].ODEIntegrator.u[2] + c_y = model.ParticlePool[particlesDrawn[k]][1].ODEIntegrator.u[3] + spreading = model.ParticlePool[particlesDrawn[k]][1].ODEIntegrator.u[6] + + #z_init = ParticleDefaults(log_energy, c_x, c_y, x, y, spreading) + model.ParticleCollection[end-k+1].position_ij = (i,j) + model.ParticleCollection[end-k+1].position_xy = (x,y) + model.ParticleCollection[end-k+1].ODEIntegrator.u[1] = log_energy + model.ParticleCollection[end-k+1].ODEIntegrator.u[2] = c_x + model.ParticleCollection[end-k+1].ODEIntegrator.u[3] = c_y + model.ParticleCollection[end-k+1].ODEIntegrator.u[4] = x + model.ParticleCollection[end-k+1].ODEIntegrator.u[5] = y + model.ParticleCollection[end-k+1].ODEIntegrator.u[6] = spreading + #push!(model.ParticleCollection,InitParticleInstance(model.ODEsystem, z_init, model.ODEsettings, (i,j), false, true)) + end + + #@info nPreviousParticles + #@info length(model.ParticleCollection) + + """ + @threads for k in eachindex(particlesToBeResetIndex) + i, j = particlesToBeResetIndex[k] + mapping_2D.remesh!(i, j, model, Δt) + end + """ + end + + if debug + @info "remeshed: " + #@info model.State[8:12, 1], model.State[8:12, 2] + @info maximum(model.State[:, :, 1]), maximum(model.State[:, :, 2]), maximum(model.State[:, :, 3]) + @info model.clock.time, model.ParticleCollection[10].ODEIntegrator.t + + end + #print("mean energy after remesh ", mean_of_state(model), "\n") + + # @printf("------- max state E=%.4e cgx=%.4e cgy=%.4e \n", max_energy(model), max_cgx(model), max_cgy(model)) + tick!(model.clock, Δt) + + +end + +function time_step!(model::Abstract2DModel, Δt::Float64; callbacks=nothing, debug=false) + + # temporary FailedCollection to store failed particles + FailedCollection = Vector{AbstractMarkedParticleInstance}([]) + + #print("mean energy before advance ", mean_of_state(model), "\n") + if debug + @info "before advance" + @info maximum(model.State[:, :, 1]), maximum(model.State[:, :, 2]), maximum(model.State[:, :, 3]) + model.FailedCollection = FailedCollection + end time_step!_advance(model, Δt, FailedCollection) # @threads for a_particle in model.ParticleCollection[model.ocean_points] @@ -128,7 +355,7 @@ function time_step!(model::Abstract2DModel, Δt::Float64; callbacks=nothing, deb # model.periodic_boundary, # model.ODEdefaults) # end - + if debug print("mean energy after advance ", mean_of_state(model), "\n") diff --git a/src/Operators/core_2D.jl b/src/Operators/core_2D.jl index ecd32b4..6154486 100644 --- a/src/Operators/core_2D.jl +++ b/src/Operators/core_2D.jl @@ -273,7 +273,7 @@ function InitParticleValues( #c̄_x = u_min[2] #c̄_y = u_min[3] - particle_on = false + particle_on = true end # initialize particle instance based on above devfined values diff --git a/src/Operators/core_2D_spread.jl b/src/Operators/core_2D_spread.jl new file mode 100644 index 0000000..84120fc --- /dev/null +++ b/src/Operators/core_2D_spread.jl @@ -0,0 +1,561 @@ +module core_2D_spread + +using DifferentialEquations: OrdinaryDiffEq.ODEIntegrator, OrdinaryDiffEq.ODEProblem, init +using ModelingToolkit: ODESystem ## depriciate when MTK is removed + +using SharedArrays +using DocStringExtensions + +# Particle-Node interaction +export GetParticleEnergyMomentum, GetVariablesAtVertex, Get_u_FromShared, ParticleDefaults, InitParticleInstance +export InitParticleValues + +#include("../Utils/FetchRelations.jl") +using ...FetchRelations + +using ...Architectures: AbstractParticleInstance, AbstractMarkedParticleInstance, StateTypeL1 +using ...Architectures: AbstractGridStatistics + +# using ..particle_waves_v3: init_vars +# t, x, y, c̄_x, c̄_y, lne, Δn, Δφ_p, r_g, C_α, C_φ, g, C_e = init_vars() + +using ...custom_structures: ParticleInstance1D, StochasticParticleInstance2D, MarkedParticleInstance + +export init_z0_to_State! +include("initialize.jl") + +# Callbacks +export wrap_pos!, periodic_BD_single_PI!, show_pos!, periodic_condition_x +include("utils.jl") + +### particle defaults ### +""" +ParticleDefaults +Structure holds default particles +# Fields +$(DocStringExtensions.FIELDS) +""" +mutable struct ParticleDefaults + "log energy" + lne::Float64 + "zonal velocity" + c̄_x::Float64 + "meridional velocity" + c̄_y::Float64 + "x Position" + x::Float64 + "y Position" + y::Float64 + "angular spreading" + angular_σ::Float64 +end + +Base.copy(s::ParticleDefaults) = ParticleDefaults(s.lne, s.c̄_x, s.c̄_y, s.x, s.y, s.angular_σ) +initParticleDefaults(s::ParticleDefaults) = [s.lne, s.c̄_x, s.c̄_y, s.x, s.y, s.angular_σ] + +mutable struct PointSource + "particle to launch" + particleLaunch::ParticleDefaults + "first launch time" + firstTimeLaunched::Float64 +end + +Base.copy(s::PointSource) = PointSource(s.particleLaunch, s.lastTimeLaunched) +PointSource(s::PointSource) = [s.particleLaunch, s.lastTimeLaunched] + +""" + + +""" + +######### Particle --> Node | 2D ########## + + +""" +GetParticleEnergyMomentum(PI) + +""" +function GetParticleEnergyMomentum(PI) + ui_lne, ui_c̄_x, ui_c̄_y, ui_x, ui_y, ui_angular_σ = PI.ODEIntegrator.u + + ui_e = exp(ui_lne) + c_speed = speed(ui_c̄_x, ui_c̄_y) + m_x = ui_c̄_x * ui_e / c_speed^2 / 2 + m_y = ui_c̄_y * ui_e / c_speed^2 / 2 + + return ui_e, m_x, m_y, ui_angular_σ +end + +function GetParticleEnergyMomentum(PI::AbstractParticleInstance) + return GetParticleEnergyMomentum(PI.ODEIntegrator.u) +end + +function GetParticleEnergyMomentum(zi::Dict) + return GetParticleEnergyMomentum([zi[lne], zi[c̄_x], zi[c̄_y], zi[x], zi[y], zi[angular_σ]]) +end + +function GetParticleEnergyMomentum(zi::ParticleDefaults) + return GetParticleEnergyMomentum([zi.lne, zi.c̄_x, zi.c̄_y, zi.x, zi.y, zi.angular_σ]) +end + +""" +GetParticleEnergyMomentum(PI) + +""" +function GetParticleEnergyMomentum(z0::Vector{Float64}) + + ui_lne, ui_c̄_x, ui_c̄_y, ui_x, ui_y, ui_angular_σ =z0 + ui_e = exp(ui_lne) + c_speed = speed(ui_c̄_x, ui_c̄_y) + m_x = ui_c̄_x * ui_e / c_speed^2 / 2 + m_y = ui_c̄_y * ui_e / c_speed^2 / 2 + + return [ui_e, m_x, m_y, ui_angular_σ] +end + +speed(x::Float64, y::Float64) = sqrt(x^2 + y^2) + +######### node--> Particle ########## + +""" +GetVariablesAtVertex(i_State::Vector{Float64}, x::Float64, y::Float64, Δn::Float64, Δφ_p::Float64 ) +i_state: [e, m_x, m_y] state vector at node +x, y: coordinates of the vertex + +""" +function GetVariablesAtVertex(i_State::Vector{Float64}, x::Float64, y::Float64) + e, m_x, m_y, sigma = i_State + m_amp = speed(m_x, m_y) + c_x = m_x * e / (2 * m_amp^2) + c_y = m_y * e / (2 * m_amp^2) + + return [log(e), c_x, c_y, x, y, sigma] +end + + +""" returns node state from shared array, given paritcle index """ +Get_u_FromShared(PI::AbstractParticleInstance, S::SharedArray,) = S[PI.position_ij[1], PI.position_ij[2], :] + +""" +GetGroupVelocity(i_State::Vector{Float64}) +i_state: [e, m_x, m_y] state vector at node +""" +function GetGroupVelocity(i_State) + e = i_State[:,:,1] + m_x = i_State[:,:,2] + m_y = i_State[:,:,3] + m_amp = speed.(m_x, m_y) + c_x = m_x .* e ./ (2 .* m_amp.^2) + c_y = m_y .* e ./ (2 .* m_amp.^2) + + return (c_x=c_x, c_y=c_y) +end + + +###### seed particles ##### + + +""" +InitParticleInstance(model::WaveGrowth1D, z_initials, pars, ij, boundary_flag, particle_on ; cbSets=nothing) +wrapper function to initalize a particle instance + inputs: + model is an initlized ODESytem + z_initials is the initial state of the ODESystem + pars are the parameters of the ODESystem + ij is the (i,j) tuple that of the initial position + boundary_flag is a boolean that indicates if the particle is on the boundary + particle_on is a boolean that indicates if the particle is on + chSet (optional) is the set of callbacks the ODE can have +""" +function InitParticleInstance(model, z_initials, ODE_settings, ij, boundary_flag, particle_on; cbSets=Nothing) + + # converty to ordered named tuple + ODE_parameters = NamedTuple{Tuple(Symbol.(keys(ODE_settings.Parameters)))}(values(ODE_settings.Parameters)) + + z_initials = initParticleDefaults(z_initials) + # create ODEProblem + problem = ODEProblem(model, z_initials, (0.0, ODE_settings.total_time), ODE_parameters) + # inialize problem + # works best with abstol = 1e-4,reltol=1e-3,maxiters=1e4, + integrator = init( + problem, + ODE_settings.solver, + saveat=ODE_settings.saving_step, + abstol=ODE_settings.abstol, + adaptive=ODE_settings.adaptive, + dt = ODE_settings.dt, + dtmin=ODE_settings.dtmin, + force_dtmin=ODE_settings.force_dtmin, + maxiters=ODE_settings.maxiters, + reltol=ODE_settings.reltol, + callback=ODE_settings.callbacks, + save_everystep=ODE_settings.save_everystep) + return StochasticParticleInstance2D(ij, (z_initials[4], z_initials[5]), integrator, boundary_flag, particle_on) +end + +""" +InitParticleInstance(model::ODESystem, z_initials, pars, ij, boundary_flag, particle_on ; cbSets=nothing) +wrapper function to initalize a particle instance + inputs: + model is an initlized ODESytem + z_initials is the initial state of the ODESystem + pars are the parameters of the ODESystem + ij is the (i,j) tuple that of the initial position + boundary_flag is a boolean that indicates if the particle is on the boundary + particle_on is a boolean that indicates if the particle is on + chSet (optional) is the set of callbacks the ODE can have +""" +function InitParticleInstance(model::ODESystem, z_initials, ODE_settings, ij, boundary_flag, particle_on; cbSets=Nothing) + + z_initials = initParticleDefaults(z_initials) + # create ODEProblem + problem = ODEProblem(model, z_initials, (0.0, ODE_settings.total_time), ODE_settings.Parameters) + # inialize problem + # works best with abstol = 1e-4,reltol=1e-3,maxiters=1e4, + integrator = init( + problem, + ODE_settings.solver, + saveat=ODE_settings.saving_step, + abstol=ODE_settings.abstol, + adaptive=ODE_settings.adaptive, + dt=ODE_settings.dt, + dtmin=ODE_settings.dtmin, + force_dtmin=ODE_settings.force_dtmin, + maxiters=ODE_settings.maxiters, + reltol=ODE_settings.reltol, + callback=ODE_settings.callbacks, + save_everystep=ODE_settings.save_everystep) + return StochasticParticleInstance2D(ij, (z_initials[4], z_initials[5]), integrator, boundary_flag, particle_on) +end + + + + +""" +InitParticleValues(defaults:: Union{Nothing,ParticleDefaults}, ij::Tuple, gridnote::OneDGridNotes, u, v, DT) + +Find initial conditions for particle. Used at the beginning of the experiment. + inputs: + defaults Dict with variables of the state vector, these will be replaced in this function + ij index of the grid point + gridnote grid to determine the position in the grid + winds NamedTuple (u,v) with interp. functions for wind values + DT time step of model, used to determine fetch laws +""" +function InitParticleValues( + defaults::PP, + xy::Tuple{Float64, Float64}, + uv::Tuple{Number, Number}, + DT) where {PP<:Union{Nothing,ParticleDefaults}} + + #i,j = ij + xx, yy = xy + + if defaults == nothing + + if speed(uv[1], uv[2]) > sqrt(2) + # defaults are not defined and there is wind + + # take in local wind velocities + WindSeaMin = FetchRelations.get_initial_windsea(uv[1], uv[2], DT) + # seed particle given fetch relations + lne = log(WindSeaMin["E"]) + c̄_x = WindSeaMin["cg_bar_x"] + c̄_y = WindSeaMin["cg_bar_y"] + angular_σ = π/8 + + particle_on = false + else + # defaults are not defined and there is no wind + u_min = FetchRelations.MinimalParticle(uv[1], uv[2], DT) + lne = u_min[1] + c̄_x = u_min[2] + c̄_y = u_min[3] + angular_σ = π/8 + + particle_on = false + end + + # initialize particle instance based on above devfined values + particle_defaults = ParticleDefaults(lne, c̄_x, c̄_y, xx, yy, angular_σ) + else + particle_defaults = defaults + particle_on = true + end + + #@show defaults + return particle_defaults, particle_on +end + +""" +InitParticleInstance(model::WaveGrowth1D, z_initials, pars, ij, boundary_flag, particle_on ; cbSets=nothing) +wrapper function to initalize a particle instance + inputs: + model is an initlized ODESytem + z_initials is the initial state of the ODESystem + pars are the parameters of the ODESystem + ij is the (i,j) tuple that of the particle position + xy is the (x,y) tuple that of the particle position + boundary_flag is a boolean that indicates if the particle is on the boundary + particle_on is a boolean that indicates if the particle is on + chSet (optional) is the set of callbacks the ODE can have +""" +function InitParticleInstance(model, z_initials, ODE_settings, ij, xy , boundary_flag, particle_on; cbSets=Nothing) + + ## to do's for add the Projection: + ## replace boundary_flag with mask that discrimiated by types [0,1,2,3] + ## add M as input to the function and add it to the ODE_parameters + ## modify ParticleInstance2D to deal with boundaries and mask + + # converty to ordered named tuple + ODE_parameters = NamedTuple{Tuple(Symbol.(keys(ODE_settings.Parameters)))}(values(ODE_settings.Parameters)) + ODE_parameters = (; ODE_parameters..., x = xy[1], y = xy[2]) + + z_initials = initParticleDefaults(z_initials) + # create ODEProblem + problem = ODEProblem(model, z_initials, (0.0, ODE_settings.total_time), ODE_parameters) + # inialize problem + # works best with abstol = 1e-4,reltol=1e-3,maxiters=1e4, + integrator = init( + problem, + ODE_settings.solver, + saveat=ODE_settings.saving_step, + abstol=ODE_settings.abstol, + adaptive=ODE_settings.adaptive, + dt = ODE_settings.dt, + dtmin=ODE_settings.dtmin, + force_dtmin=ODE_settings.force_dtmin, + maxiters=ODE_settings.maxiters, + reltol=ODE_settings.reltol, + callback=ODE_settings.callbacks, + save_everystep=ODE_settings.save_everystep) + + return StochasticParticleInstance2D(ij, (xy[1], xy[2]), integrator, boundary_flag, particle_on) +end + + +# --------------------------- up to here for now (remove MTK)--------------------------- # + +""" +ResetParticleValues(PI::AbstractParticleInstance, particle_defaults::Union{Nothing,ParticleDefaults,Vector{Float64}}, ODE_settings::ODESettings, gridnote::OneDGridNotes, winds, DT) + +Reset the state of a particle instance + inputs: + PI ParticleInstance + particle_defaults Dict with variables of the state vector, these will be replaced in this function + ODE_settings ODESettings type + gridnote grid to determine the position in the grid + winds NamedTuple (u,v) with interp. functions for wind values + DT time step of model, used to determine fetch laws +returns: + dict or vector +""" +function ResetParticleValues( + defaults::PP, + PI::AbstractParticleInstance, + wind_tuple, + DT, vector=true) where {PP<:Union{Nothing,ParticleDefaults,Vector{Float64}}} + + if defaults == nothing # this is boundary_defaults = "wind_sea" + #@info "init particles from fetch relations: $z_i" + #particle_defaults = Dict{Num,Float64}() + #xx, yy = PI.position_xy[1], PI.position_xy[2] + # take in local wind velocities + u_init, v_init = wind_tuple[1], wind_tuple[2] + #winds.u(xx, yy, 0), winds.v(xx, yy, 0) + + WindSeaMin = FetchRelations.get_initial_windsea(u_init, v_init, DT) + # seed particle given fetch relations + particle_defaults = ParticleDefaults(log(WindSeaMin["E"]), WindSeaMin["cg_bar_x"], WindSeaMin["cg_bar_y"], PI.position_xy[1], PI.position_xy[2], π/8) + elseif typeof(defaults) == Vector{Float64} # this is for the case of minimal wind sea + particle_defaults = defaults + particle_defaults[4] = PI.position_xy[1] + particle_defaults[5] = PI.position_xy[2] + else + particle_defaults = defaults + end + + #@show defaults + if vector + return initParticleDefaults(particle_defaults) + else + return particle_defaults + end + +end +""" +check_boundary_point(i, boundary, periodic_boundary) +checks where ever or not point is a boundary point. +returns Bool +""" +function check_boundary_point(i, boundary, periodic_boundary) + return periodic_boundary ? false : (i in boundary) +end + +function check_boundary_point(imesh, periodic_boundary) + if periodic_boundary + return imesh == 2 # only land points are treated as boundary points + else + return imesh >= 2 # land points and grid boundary points are treated as boundary points + end +end + + +# """ +# SeedParticle!(ParticleCollection ::Vector{Any}, State::SharedMatrix, i::Int64, +# particle_system::ODESystem, particle_defaults::Union{ParticleDefaults,Nothing}, ODE_settings, +# GridNotes, winds, DT:: Float64, Nx:: Int, boundary::Vector{Int}, periodic_boundary::Bool) + +# Seed Pickles to ParticleColletion and State +# """ +# function SeedParticle!( +# ParticleCollection::Vector{Any}, +# State::SharedArray, +# ij::Tuple{Int, Int}, + +# particle_system::SS, +# particle_defaults::PP, +# ODE_settings, #particle_waves_v3.ODESettings type + +# GridNotes, # ad type of grid note +# winds, # interp winds +# DT::Float64, + +# boundary::Vector{T}, +# periodic_boundary::Bool) where {T<:Union{Int,Any,Nothing,Int64},PP<:Union{ParticleDefaults,Nothing},SS<:Union{ODESystem,Any}} + +# xx, yy = GridNotes.x[ij[1]], GridNotes.y[ij[2]] + + +# #uv = winds.u(xx, yy, 0)::Union{Num,Float64}, winds.v(xx, yy, 0)::Union{Num,Float64} +# uv = winds.u(xx, yy, 0)::Float64, winds.v(xx, yy, 0)::Float64 + +# # define initial condition +# z_i, particle_on = InitParticleValues(particle_defaults, (xx,yy), uv, DT) +# # check if point is boundary point +# boundary_point = check_boundary_point(ij, boundary, periodic_boundary) +# #@info "boundary?", boundary_point + +# # if boundary_point +# # #@info "boundary point", boundary_point, particle_on +# # # if boundary point, then set particle to off +# # particle_on = false +# # end + +# # add initial state to State vector +# if particle_on +# init_z0_to_State!(State, ij, GetParticleEnergyMomentum(z_i)) +# end + +# # Push Inital condition to collection +# push!(ParticleCollection, +# InitParticleInstance( +# particle_system, +# z_i, +# ODE_settings, +# ij, +# boundary_point, +# particle_on)) +# nothing +# end + + +""" + function SeedParticle(State::StateTypeL1, ij:: (Int64, Int64) + particle_system::Any, particle_defaults::Union{ParticleDefaults,Nothing}, ODE_settings, + ij_mesh, ij_wind_tuple, DT:: Float64, boundary::Vector{Int}, periodic_boundary::Bool) + + returns ParicleInstance that can be pushed to ParticleColletion +""" +function SeedParticle( + State::StateTypeL1, + ij::II, # position tuple in grid + + particle_system::SS, + particle_defaults::PP, + ODE_settings, #particle_waves_v3.ODESettings type + + gridstats::AbstractGridStatistics, + ProjetionKernel::Function, + PropagationCorrection::Function, + + ij_mesh::NamedTuple, # local grid information + ij_wind::Tuple, # interp winds + + DT::Float64, + boundary::Vector{T}, periodic_boundary::Bool) where + {II<:Union{Tuple{Int,Int},CartesianIndex},T<:Union{Int,Any,Nothing,Int64},PP<:Union{ParticleDefaults,Nothing},SS<:Any} + + xy = (ij_mesh.x, ij_mesh.y) + # 1st check if particle is not in mask, Land points == 0 + if ij_mesh.mask == 0 # land point + # init dummy instance + return StochasticParticleInstance2D(ij, xy , nothing, false, false) + end + + # define initial condition + # particle initial condition is always (0,0) in relative coordinates not xy anymore + z_i, particle_on = InitParticleValues(particle_defaults, (0.0, 0.0) , ij_wind, DT) + + # check if point is boundary point <-- replace in the future with with mask: 0 = land, 1 = ocean, 2= land boundary, 3 = domain boundary + # boundary_point = check_boundary_point(ij, boundary, periodic_boundary) # old version that compares to list + boundary_point = check_boundary_point(ij_mesh.mask, periodic_boundary) + + # add initial state to State vector + if particle_on + init_z0_to_State!(State, ij, GetParticleEnergyMomentum(z_i)) + end + + # check if Propgation Correction is set in gridstats + # PropagationCorrection = gridstats.PropagationCorrection != nothing ? PropagationCorrection(ij_mesh, gridstats) : x -> 0.0 + + # set projection: + ODE_settings.Parameters = (; ODE_settings.Parameters..., M=ProjetionKernel(ij_mesh, gridstats), PC=PropagationCorrection(ij_mesh, gridstats)) + + return InitParticleInstance( + particle_system, + z_i, + ODE_settings, + ij, + xy, + boundary_point, + particle_on) + +end + +# """ +# SeedParticle!(ParticleCollection ::Vector{Any}, State::SharedMatrix, i::Int64, +# particle_system::ODESystem, particle_defaults::Union{ParticleDefaults,Nothing}, ODE_settings, +# GridNotes, winds, DT:: Float64, Nx:: Int, boundary::Vector{Int}, periodic_boundary::Bool) + +# Seed Pickles to ParticleColletion and State +# """ +# function SeedParticle!( +# ParticleCollection::Vector{Any}, +# State::SharedArray, +# ij::Tuple{Int, Int}, + +# particle_system::SS, +# particle_defaults::PP, +# ODE_settings, #particle_waves_v3.ODESettings type + +# GridNotes, # ad type of grid note +# winds, # interp winds +# DT::Float64, + +# boundary::Vector{T}, +# periodic_boundary::Bool) where {T<:Union{Int,Any,Nothing,Int64},PP<:Union{ParticleDefaults,Nothing},SS<:Union{ODESystem,Any}} + +# # Push Inital condition to collection +# push!(ParticleCollection, +# SeedParticle(State, ij, +# particle_system, particle_defaults, ODE_settings, #particle_waves_v3.ODESettings type +# GridNotes, winds, DT, +# boundary, periodic_boundary)) +# nothing +# end + + + +# end of module +end \ No newline at end of file diff --git a/src/Operators/initialize.jl b/src/Operators/initialize.jl index 4c48ff8..d54a2af 100644 --- a/src/Operators/initialize.jl +++ b/src/Operators/initialize.jl @@ -18,7 +18,7 @@ end """ sets particle state values to S. position is taking from particle """ -function set_u_to_shared!(S::StateTypeL1, PI::ParticleInstance2D) +function set_u_to_shared!(S::StateTypeL1, PI::AbstractParticleInstance) S[ PI.position_ij[1], PI.position_ij[2] , :] = PI.ODEIntegrator.u nothing end diff --git a/src/Operators/mapping_2D.jl b/src/Operators/mapping_2D.jl index a0452f1..396b47a 100644 --- a/src/Operators/mapping_2D.jl +++ b/src/Operators/mapping_2D.jl @@ -4,22 +4,27 @@ using SharedArrays using StaticArrays using DifferentialEquations using Printf +using Random, Distributions using ...ParticleMesh: TwoDGrid, TwoDGridNotes +using PiCLES.Grids.CartesianGrid: TwoDCartesianGridMesh, ProjetionKernel, TwoDCartesianGridStatistics import ...ParticleInCell as PIC using ...FetchRelations using ...custom_structures: ParticleInstance1D, ParticleInstance2D, MarkedParticleInstance -using ..core_2D: GetParticleEnergyMomentum, GetVariablesAtVertex, Get_u_FromShared, ResetParticleValues, ParticleDefaults +using ..core_2D_spread: GetParticleEnergyMomentum as StochasticGetParticleEnergyMomentum +using ..core_2D_spread: GetVariablesAtVertex as StochasticGetVariablesAtVertex +using ..core_2D_spread: Get_u_FromShared as StochasticGet_u_FromShared +using ..core_2D_spread: ResetParticleValues as StochasticResetParticleValues +using ..core_2D_spread: ParticleDefaults as StochasticParticleDefaults +using ..core_2D_spread: InitParticleInstance as StochasticInitParticleInstance +using ..core_2D: GetParticleEnergyMomentum, GetVariablesAtVertex, ParticleDefaults, InitParticleInstance, Get_u_FromShared, ResetParticleValues -using ...Architectures: AbstractParticleInstance, AbstractMarkedParticleInstance, AbstractODESettings, StateTypeL1 +using ...Architectures: AbstractParticleInstance, AbstractStochasticParticleInstance, AbstractMarkedParticleInstance, AbstractODESettings, StateTypeL1, Abstract2DModel, Abstract2DStochasticModel using ...Architectures: Grid2D, CartesianGrid, CartesianGridStatistics, CartesianGrid2D, CartesianGrid1D, AbstractGridStatistics, AbstractGrid, StandardRegular2D_old, MeshGrids, MeshGridStatistics - -# using ...ParticleMesh: TwoDGrid, TwoDGridNotes, TwoDGridMesh -# using ...Grids.CartesianGrid: TwoDCartesianGridMesh, TwoDCartesianGridStatistics ###### remeshing routines ############ @@ -40,6 +45,21 @@ S Shared array where particles are stored G (TwoDGrid) Grid that defines the nodepositions """ +function ParticleToNode!(PI::AbstractStochasticParticleInstance, particlesAtNode::Array{Array{Array{Any,1},1},1}, S::SharedArray, G::TwoDCartesianGridMesh, periodic_boundary::Bool) + + #u[4], u[5] are the x and y positions of the particle + #index_positions, weights = PIC.compute_weights_and_index(G, PI.ODEIntegrator.u[4], PI.ODEIntegrator.u[5]) + weights_and_index = PIC.compute_weights_and_index_mininal(G, PI.ODEIntegrator.u[4], PI.ODEIntegrator.u[5]) + + #ui[1:2] .= PI.position_xy + u_state = StochasticGetParticleEnergyMomentum(PI.ODEIntegrator.u) + #@show u_state + + #PIC.push_to_grid!(S, particlesAtNode, PI, u_state , index_positions, weights, G.stats.Nx.N, G.stats.Ny.N , periodic_boundary) + PIC.push_to_grid!(S, particlesAtNode, PI, u_state , weights_and_index, G.stats.Nx.N, G.stats.Ny.N , periodic_boundary) + nothing +end + function ParticleToNode!(PI::AbstractParticleInstance, S::StateTypeL1, G::TwoDGrid, periodic_boundary::Bool) #u[4], u[5] are the x and y positions of the particle @@ -51,8 +71,8 @@ function ParticleToNode!(PI::AbstractParticleInstance, S::StateTypeL1, G::TwoDGr u_state = GetParticleEnergyMomentum(PI.ODEIntegrator.u) #@show u_state - #PIC.push_to_grid!(S, u_state , index_positions, weights, G.Nx, G.Ny , periodic_boundary) - PIC.push_to_grid!(S, u_state , weights_and_index, G.Nx, G.Ny , periodic_boundary) + #PIC.push_to_grid!(S, u_state , index_positions, weights, G.stats.Nx.N, G.stats.Ny.N , periodic_boundary) + PIC.push_to_grid!(S, u_state , weights_and_index, G.stats.Nx.N, G.stats.Ny.N , periodic_boundary) nothing end @@ -67,7 +87,7 @@ function ParticleToNode!(PI::AbstractParticleInstance, S::StateTypeL1, G::MeshGr u_state = GetParticleEnergyMomentum(PI.ODEIntegrator.u) #@show u_state - #PIC.push_to_grid!(S, u_state , index_positions, weights, G.Nx, G.Ny , periodic_boundary) + #PIC.push_to_grid!(S, u_state , index_positions, weights, G.stats.Nx.N, G.stats.Ny.N , periodic_boundary) PIC.push_to_grid!(S, u_state, weights_and_index, G.stats.Nx, G.stats.Ny) nothing end @@ -82,6 +102,18 @@ function ParticleToNode!(PI::AbstractParticleInstance, S::SharedMatrix, u_state: nothing end +function set_u_and_t_stochas!(integrator, u_new::CC, t_new::Number) where CC <:Union{Vector{Float64},MVector} + # adding a small deviation due to directional spreading + d = Normal(0, u_new[6]^2) + delta_phi = rand(d,1) + c_x = u_new[2] * cos(delta_phi[1]) - u_new[3] * sin(delta_phi[1]) + c_y = u_new[2] * sin(delta_phi[1]) + u_new[3] * cos(delta_phi[1]) + u_new[2] = c_x + u_new[3] = c_y + integrator.u = u_new + integrator.t = t_new +end + function set_u_and_t!(integrator, u_new::CC, t_new::Number) where CC <:Union{Vector{Float64},MVector} integrator.u = u_new integrator.t = t_new @@ -96,9 +128,26 @@ function reset_PI_u!(PI::AbstractParticleInstance; ui::CC) where CC<:Union{Vecto end -function reset_PI_ut!(PI::AbstractParticleInstance; ui::CC, ti::Number) where CC <:Union{Vector{Float64},MVector} +function reset_PI_ut!(PI::AbstractStochasticParticleInstance; ui::CC, ti::Number, stochas::Bool) where CC <:Union{Vector{Float64},MVector} + # this method keeps the correct time for time varying forcing (~may 2023) + if stochas + set_u_and_t_stochas!(PI.ODEIntegrator, ui, ti) + else + @info "called the wrong version" + set_u_and_t!(PI.ODEIntegrator, ui, ti) + end + u_modified!(PI.ODEIntegrator, true) + auto_dt_reset!(PI.ODEIntegrator) +end + +function reset_PI_ut!(PI::AbstractParticleInstance; ui::CC, ti::Number, stochas::Bool) where CC <:Union{Vector{Float64},MVector} # this method keeps the correct time for time varying forcing (~may 2023) - set_u_and_t!(PI.ODEIntegrator, ui, ti) + if stochas + set_u_and_t!(PI.ODEIntegrator, ui, ti) + else + @info "called the wrong version" + set_u_and_t!(PI.ODEIntegrator, ui, ti) + end u_modified!(PI.ODEIntegrator, true) auto_dt_reset!(PI.ODEIntegrator) end @@ -113,8 +162,143 @@ end ######### Core routines for advancing and remeshing """ - advance!(PI::AbstractParticleInstance, S::SharedMatrix{Float64}, G::TwoDGrid, DT::Float64) + advance!(PI::AbstractStochasticParticleInstance, S::SharedMatrix{Float64}, G::TwoDGrid, DT::Float64) """ +function advance!(PI::AbstractStochasticParticleInstance, + particlesAtNode::Array{Array{Array{Any,1},1},1}, + S::StateTypeL1, + Failed::Vector{AbstractMarkedParticleInstance}, + Grid::Union{Grid2D,MeshGrids}, + winds::NamedTuple{(:u, :v)}, + DT::Float64, + log_energy_maximum::Float64, + wind_min_squared::Float64, + periodic_boundary::Bool, + default_particle::PP, + ) where {PP<:Union{StochasticParticleDefaults,Nothing}} + #@show PI.position_ij + + #if ~PI.boundary # if point is not a + t_start = copy(PI.ODEIntegrator.t) + add_saveat!(PI.ODEIntegrator, PI.ODEIntegrator.t ) + savevalues!(PI.ODEIntegrator) + + # set the position in particle state vector either to the node position or to the relative position in the CartesianGrid + if typeof(Grid) <: MeshGrids + xy = (0.0,0.0) + # @info "advance: CartesianGrid" + elseif typeof(Grid) <: StandardRegular2D_old + xy = PI.position_xy[1], PI.position_xy[2] + # @info "advance: StandardRegular2D_old" + else + @info "advance: no grid detected" + end + + + # advance particle + if PI.on #& ~PI.boundary # if Particle is on and not boundary + + try + # PI.ODEIntegrator.u[4] += DT*PI.ODEIntegrator.u[2] + # PI.ODEIntegrator.u[5] += DT*PI.ODEIntegrator.u[3] + before = sqrt(PI.ODEIntegrator.u[2]^2+PI.ODEIntegrator.u[2]^2) + @info PI + step!(PI.ODEIntegrator, DT, true) + after = sqrt(PI.ODEIntegrator.u[2]^2+PI.ODEIntegrator.u[2]^2) + # @info after/before + catch e + @printf "error on advancing ODE:\n" + print("- time after fail $(PI.ODEIntegrator.t)\n ") + print("- error message: $(e)\n") + print("- push to failed\n") + print("- state of particle: $(PI.ODEIntegrator.u)\n") + print("- winds are: $(winds.u( PI.ODEIntegrator.u[4], PI.ODEIntegrator.u[5], PI.ODEIntegrator.t))\n") + print("- winds are: $(winds.v( PI.ODEIntegrator.u[4], PI.ODEIntegrator.u[5], PI.ODEIntegrator.t))\n") + push!(Failed, + MarkedParticleInstance( + copy(PI), + copy(PI.ODEIntegrator.t), + copy(PI.ODEIntegrator.u), + PI.ODEIntegrator.sol.retcode + )) + return + + end + + elseif ~PI.on #& ~PI.boundary # particle is off, test if there was windsea + + t_end = t_start + DT + wind_end = convert(Tuple{Float64,Float64}, + (winds.u(PI.position_xy[1], PI.position_xy[2], t_end), + winds.v(PI.position_xy[1], PI.position_xy[2], t_end)))::Tuple{Float64,Float64} + + # test if winds where strong enough + if speed_square(wind_end[1], wind_end[2]) >= wind_min_squared + # winds are large eneough, reinit + ui = StochasticResetParticleValues(default_particle, xy, wind_end, DT) + reset_PI_u!(PI, ui =ui) + PI.on = true + end + + else #particle is on and boundary + + #@info "particle is on and boundary" + # particle stays off or is bounaday. do not advance + PI.on=false + return + end + + # # check if integration reached limits or is nan, or what ever. if so, reset + if sum(isnan.(PI.ODEIntegrator.u[1:3])) > 0 + @info "position or Energy is nan, reset" + @info PI.position_ij + @show PI + + t_end = t_start + DT + winds_start = convert( Tuple{Float64,Float64}, + (winds.u(PI.position_xy[1], PI.position_xy[2], t_end), + winds.v(PI.position_xy[1], PI.position_xy[2], t_end)))::Tuple{Float64,Float64} + @show winds_start + + ui = StochasticResetParticleValues(default_particle, xy, winds_start, DT) + @show PI.ODEIntegrator.u + reset_PI_u!(PI, ui=ui) + + elseif sum(isinf.(PI.ODEIntegrator.u[1:3])) > 0 + @info "position or Energy is inf" + @show PI + + winds_start = convert(Tuple{Float64,Float64}, + (winds.u(PI.position_xy[1], PI.position_xy[2], t_start), + winds.v(PI.position_xy[1], PI.position_xy[2], t_start)))::Tuple{Float64,Float64} + + ui = StochasticResetParticleValues(default_particle, xy, winds_start, DT) + reset_PI_u!(PI, ui=ui) + + elseif PI.ODEIntegrator.u[1] > log_energy_maximum + @info "e_max_log is reached" + #@show PI + + # winds_start = convert(Tuple{Float64,Float64}, + # (winds.u(PI.position_xy[1], PI.position_xy[2], t_start), + # winds.v(PI.position_xy[1], PI.position_xy[2], t_start)))::Tuple{Float64,Float64} + + ui = PI.ODEIntegrator.u + ui[1] = log_energy_maximum + # ui = StochasticResetParticleValues(default_particle, xy, winds_start, DT) + reset_PI_u!(PI, ui=ui) + + end + + #if PI.ODEIntegrator.u[1] > -13.0 #ODEs.log_energy_minimum # the minimum enerçy is distributed to 4 neighbouring particles + if PI.on + ParticleToNode!(PI, particlesAtNode, S, Grid, periodic_boundary) + # ParticleToNode!(PI, S, Grid, periodic_boundary) + end + + #return PI +end + function advance!(PI::AbstractParticleInstance, S::StateTypeL1, Failed::Vector{AbstractMarkedParticleInstance}, @@ -268,6 +452,140 @@ function remesh!(PI::ParticleInstance2D, S::StateTypeL1, #return PI end +""" + remesh!(PI::ParticleInstance2D, S::SharedMatrix{Float64, 3}) + Wrapper function that does everything necessary to remesh the particles. + - pushes the Node State to particle instance +""" +function remesh!(i::Int64, j::Int64, model::Abstract2DStochasticModel, DT::Float64) + winds = model.winds + ti = model.clock.time + grid = model.grid + + x = grid.xmin + grid.dx*(i-1) + y = grid.ymin + grid.dy*(j-1) + winds_i::Tuple{Float64,Float64} = winds.u(x, y, ti), winds.v(x, y, ti) + + NodeToParticle!(i, j, x, y, model, winds_i, DT) +end + +""" + NodeToParticle!(PI::AbstractParticleInstance, S::SharedMatrix) +Pushes node value to particle: +- If Node value is smaller than a minimal value, the particle is renintialized +- If Node value is okey, it is converted to state variable and pushed to particle. +- The particle position is set to the node positions +""" +function NodeToParticle!(i::Int64, j::Int64, x::Float64, y::Float64, model::Abstract2DStochasticModel, wind_tuple::Tuple{Float64,Float64}, DT::Float64) + S = model.State + minimal_particle = model.minimal_particle + minimal_state = model.minimal_state + wind_min_squared = model.ODEsettings.wind_min_squared + default_particle = model.ODEdefaults + e_min_log = model.ODEsettings.log_energy_minimum + current_is_boundary = i==0 || i==model.grid.stats.Nx || j==0 || j==model.grid.Ny + + # definition of a norm 2 + norm(vec) = sqrt(sum([vec[i]^2 for i in eachindex(vec)])) + one_norm(vec) = sum([vec[i] for i in eachindex(vec)]) + + + # load data from shared array + u_state = S[i, j, :] + + #last_t = ti + last_t = model.clock.time + # minimal_state[1] is the minmal Energy + # minimal_state[2] is the minmal momentum squared + if ~current_is_boundary & (u_state[1] >= minimal_state[1]) & (speed_square(u_state[2], u_state[3]) >= minimal_state[2]) # all interior nodes: convert note to particle values and push to ODEIntegrator + if model.angular_spreading_type == "geometrical" + nNewParticles = Int64((u_state[4] ÷ model.angular_spreading_thresh)*2+1) + if nNewParticles == 1 + nNewParticles = 3 + end + midParticleInd = Int64(nNewParticles ÷ 2 + 1) + ui = StochasticGetVariablesAtVertex(u_state, x, y) + energies = [exp(-(k-midParticleInd)^2/(2*π*midParticleInd)^2) for k in 1:nNewParticles] + energies = energies ./ sum(energies) .* exp(ui[1]) + + for k in 1:nNewParticles + delta_phi = (k-midParticleInd)/midParticleInd*u_state[4] + c_x = ui[2] * cos(delta_phi) - ui[3] * sin(delta_phi) + c_y = ui[2] * sin(delta_phi) + ui[3] * cos(delta_phi) + z_init = StochasticParticleDefaults(log(energies[k]), c_x, c_y, x, y, u_state[4]/nNewParticles) + push!(model.ParticleCollection,StochasticInitParticleInstance(model.ODEsystem, z_init, model.ODEsettings, (i,j), false, true)) + end + elseif model.angular_spreading_type == "nonparametric" + nNewParticles = model.n_particles_launch + nPartToDrawFrom = length(model.ParticlesAtNode[i][j]) + energiesNormed = [model.ParticlesAtNode[i][j][k][1]*exp(model.ParticlesAtNode[i][j][k][3].ODEIntegrator.u[1]) for k in 1:nPartToDrawFrom] + energiesNormed = energiesNormed / u_state[1] + a = Categorical(energiesNormed) + particlesDrawn = rand(a, nNewParticles) + new_energy_tot = sum([energiesNormed[k] for k in particlesDrawn]) + #@info "u_state = ", u_state[1] + #@info "new_energy_tot = ", new_energy_tot + #@info particlesDrawn + energy_factor = u_state[1] / new_energy_tot + #@info "new val = ", sum([energiesNormed[particlesDrawn[k]]*energy_factor for k in 1:nNewParticles]) + #@info model.ParticlesAtNode[i][j][particlesDrawn[1]][3].ODEIntegrator.u[4], model.ParticlesAtNode[i][j][particlesDrawn[1]][3].position_xy[1], model.ParticlesAtNode[i][j][particlesDrawn[1]][3].position_xy[1] + model.grid.dx + #@info model.ParticlesAtNode[i][j][particlesDrawn[1]][3].ODEIntegrator.u[5], model.ParticlesAtNode[i][j][particlesDrawn[1]][3].position_xy[2], model.ParticlesAtNode[i][j][particlesDrawn[1]][3].position_xy[2] + model.grid.dy + for k in 1:nNewParticles + + log_energy = log(energiesNormed[particlesDrawn[k]]*energy_factor) + c_x = model.ParticlesAtNode[i][j][particlesDrawn[k]][3].ODEIntegrator.u[2] + c_y = model.ParticlesAtNode[i][j][particlesDrawn[k]][3].ODEIntegrator.u[3] + spreading = model.ParticlesAtNode[i][j][particlesDrawn[k]][3].ODEIntegrator.u[6] + + z_init = StochasticParticleDefaults(log_energy, c_x, c_y, x, y, spreading) + push!(model.ParticleCollection,StochasticInitParticleInstance(model.ODEsystem, z_init, model.ODEsettings, (i,j), false, true)) + end + end + elseif ~PI.boundary & (speed_square(wind_tuple[1], wind_tuple[2]) >= wind_min_squared) #minimal windsea is not big enough but local winds are strong enough #(u_state[1] < exp(e_min_log)) | PI.boundary + # test if particle is below energy threshold, or + # if particle is at the boundary + + ui = StochasticResetParticleValues(default_particle, PI, wind_tuple, DT) # returns winds sea given DT and winds + reinit!(PI.ODEIntegrator, ui, erase_sol=false, reset_dt=true, reinit_cache=true)#, reinit_callbacks=true) + reset_PI_t!(PI, ti=last_t) + + PI.on = true + + elseif PI.boundary & (speed_square(wind_tuple[1], wind_tuple[2]) >= wind_min_squared) # at the boundary, reset particle if winds are strong enough + + ui = StochasticResetParticleValues(default_particle, PI, wind_tuple, DT) # returns winds sea given DT and winds + reinit!(PI.ODEIntegrator, ui, erase_sol=false, reset_dt=true, reinit_cache=true)#, reinit_callbacks=true) + reset_PI_t!(PI, ti=last_t) + + PI.on = true + + + else # particle is below energy threshold & on boundary + #PI.ODEIntegrator.u = StochasticResetParticleValues(minimal_particle, PI, wind_tuple, DT) + # if ~PI.boundary + # @info u_state + # end + PI.on = false + end + +end + +""" + remesh!(PI::ParticleInstance2D, S::SharedMatrix{Float64, 3}) + Wrapper function that does everything necessary to remesh the particles. + - pushes the Node State to particle instance +""" +function remesh!(i::Int64, j::Int64, model::Abstract2DStochasticModel, DT::Float64) + winds = model.winds + ti = model.clock.time + grid = model.grid + + x = grid.xmin + grid.dx*(i-1) + y = grid.ymin + grid.dy*(j-1) + winds_i::Tuple{Float64,Float64} = winds.u(x, y, ti), winds.v(x, y, ti) + + NodeToParticle!(i, j, x, y, model, winds_i, DT) +end """ NodeToParticle!(PI::AbstractParticleInstance, S::SharedMatrix) @@ -307,8 +625,8 @@ function NodeToParticle!(PI::AbstractParticleInstance, S::StateTypeL1, #@show "u_state", u_state ui = GetVariablesAtVertex(u_state, xy[1], xy[2]) - #@info exp(ui[1]), ui[2], ui[4]/1e3, ui[5]/1e3 - reset_PI_ut!(PI; ui=ui, ti=last_t) + #@info exp(ui[1]), ui[2], ui[4]/1e3, ui[5]/1e3 + reset_PI_ut!(PI; ui=ui, ti=last_t, stochas=true) PI.on = true # this method is more robust than the set_u! method (~february 2023) diff --git a/src/ParticleInCell.jl b/src/ParticleInCell.jl index cb87061..16bc854 100644 --- a/src/ParticleInCell.jl +++ b/src/ParticleInCell.jl @@ -2,6 +2,7 @@ module ParticleInCell using ..Architectures: StateTypeL1, AbstractBoundary using Statistics +using ..Grids.CartesianGrid: TwoDCartesianGridMesh, ProjetionKernel, TwoDCartesianGridStatistics using SharedArrays using StaticArrays @@ -10,6 +11,7 @@ using ..custom_structures: wni, N_Periodic, N_NonPeriodic, N_TripolarNorth # %% using ..ParticleMesh +using ...Architectures: AbstractParticleInstance # Tolerance for comparison of real numbers: set it here! # Set parameters @@ -99,16 +101,16 @@ end """ -compute_weights_and_index(g_pars::TwoDGrid, xp::Float64, yp:: Float64 ) +compute_weights_and_index(g_pars::TwoDCartesianGridMesh, xp::Float64, yp:: Float64 ) returns indexes and weights for in 2D for single x,y point """ -function compute_weights_and_index(g_pars::TwoDGrid, xp::Float64, yp:: Float64 ) +function compute_weights_and_index(g_pars::TwoDCartesianGridMesh, xp::Float64, yp:: Float64 ) """ 2d wrapper for 1d function """ - xp_normed = norm_distance(xp, g_pars.xmin, g_pars.dx) - yp_normed = norm_distance(yp, g_pars.ymin, g_pars.dy) + xp_normed = norm_distance(xp, g_pars.stats.xmin, g_pars.stats.dx) + yp_normed = norm_distance(yp, g_pars.stats.ymin, g_pars.stats.dy) xi, xw = get_absolute_i_and_w(xp_normed) yi, yw = get_absolute_i_and_w(yp_normed) @@ -121,16 +123,16 @@ end """ -compute_weights_and_index_mininal(g_pars::TwoDGrid, xp::Float64, yp:: Float64 ) +compute_weights_and_index_mininal(g_pars::TwoDCartesianGridMesh, xp::Float64, yp:: Float64 ) returns indexes and weights as FieldVector for in 2D for single x,y point """ -function compute_weights_and_index_mininal(g_pars::TwoDGrid, xp::Float64, yp::Float64) +function compute_weights_and_index_mininal(g_pars::TwoDCartesianGridMesh, xp::Float64, yp::Float64) """ 2d wrapper for 1d function """ - xp_normed = norm_distance(xp, g_pars.xmin, g_pars.dx) # multiples of grid spacing - yp_normed = norm_distance(yp, g_pars.ymin, g_pars.dy) # multiples of grid spacing + xp_normed = norm_distance(xp, g_pars.stats.xmin, g_pars.stats.dx) # multiples of grid spacing + yp_normed = norm_distance(yp, g_pars.stats.ymin, g_pars.stats.dy) # multiples of grid spacing xi, xw = get_absolute_i_and_w(xp_normed) yi, yw = get_absolute_i_and_w(yp_normed) @@ -174,10 +176,10 @@ end #indexes, weights = compute_weights_and_index(grid2d, 10.0, 9.9 ) """ -compute_weights_and_index(g_pars::TwoDGrid, xp::Float64, yp:: Float64 ) +compute_weights_and_index(g_pars::TwoDCartesianGridMesh, xp::Float64, yp:: Float64 ) returns indexes and weights for in 2D for vectors """ -function compute_weights_and_index(grid::TwoDGrid, xp::Vector{Float64}, yp:: Vector{Float64} ) +function compute_weights_and_index(grid::TwoDCartesianGridMesh, xp::Vector{Float64}, yp:: Vector{Float64} ) """ returns: weight list 2 N X 1 vector of tuples with indices and weights @@ -316,7 +318,7 @@ end ## most recent version 2D - not with Boundary types function push_to_grid!(grid::StateTypeL1, - charge::CC, + charge::CC, index_pos::II, weights::WW, Nx::Int, Ny::Int, @@ -337,6 +339,31 @@ function push_to_grid!(grid::StateTypeL1, end + +# ## most recent version 2D - not with Boundary types +function push_to_grid!(grid::SharedArray{Float64, 3}, + particlesAtNode::Array{Array{Array{Any,1},1},1}, + PI::AbstractParticleInstance, + charge::Vector{Float64}, + index_pos::Tuple{Int, Int}, + weights::Tuple{Float64, Float64}, + Nx::Int, Ny::Int, + periodic::Bool = true) + if periodic + grid[ wrap_index!(index_pos[1], Nx) , wrap_index!(index_pos[2], Ny), : ] += weights[1] * weights[2] * charge + else + if sum( test_domain(index_pos, Nx, Ny) ) != 2 + # position outside of domain + return + else + # position is inside the domain + push!(particlesAtNode[index_pos[1]][index_pos[2]], [weights[1] * weights[2], charge[4], PI]) + grid[ index_pos[1] , index_pos[2], 1:3 ] += weights[1] * weights[2] * charge[1:3] + end + end + +end + # Abstract Boundary Version function push_to_grid!(grid::StateTypeL1, charge::CC, @@ -488,6 +515,20 @@ function push_to_grid!(grid::StateTypeL1, end end +function push_to_grid!(grid::SharedArray{Float64, 3}, + particlesAtNode::Array{Array{Array{Any,1},1},1}, + PI::AbstractParticleInstance, + charge::Vector{Float64}, + index_pos::Vector{Tuple{Int, Int}}, + weights::Vector{Tuple{Float64, Float64}}, + Nx::Int, Ny::Int, + periodic::Bool=true) + #@info "this is version D" + for (i, w) in zip(index_pos, weights) + push_to_grid!(grid, particlesAtNode, PI, charge , i, w , Nx, Ny, periodic) + end +end + #push_to_grid!(charges_grid, 1.0 , index_positions[3], weights[3] , grid2d.Nx , grid2d.Ny ) ## allocation optimized: @@ -511,6 +552,19 @@ end """ wrapper over FieldVector weight&index (wni), """ +function push_to_grid!(grid::SharedArray{Float64, 3}, + particlesAtNode::Array{Array{Array{Any,1},1},1}, + PI::AbstractParticleInstance, + charge::CC, + wni::FieldVector, + Nx::Int, Ny::Int, + periodic::Bool=true) where CC <: Union{Vector{Float64}, SVector{3, Float64}} + #@info "this is version D" + for (i, w) in construct_loop(wni) + push_to_grid!(grid,particlesAtNode, PI, charge, i, w, Nx, Ny, periodic) + end +end + function push_to_grid!(grid::StateTypeL1, charge::CC, wni::FieldVector, @@ -556,6 +610,22 @@ function push_to_grid!(grid::SharedMatrix{Float64}, end end +# function push_to_grid!(grid::SharedMatrix{Float64}, +# particlesAtNode::Array{Array{Array{Any,1},1},1}, +# PI::AbstractParticleInstance, +# charge::CC, +# index_pos::Vector{Any}, +# weights::Vector{Any}, +# Nx::Int, +# periodic::Bool=true) where CC <: Union{Vector{Float64}, SVector{Float64}} +# #@info "this is version C" +# for (im, wm, c) in zip(index_pos, weights, charge) +# for (i, w) in zip(im, wm) +# push_to_grid!(grid, particlesAtNode, PI, c, i, w, Nx, periodic) +# end +# end +# end + # multiple index positions and 1 charge diff --git a/src/ParticleSystems/ParticleSystems.jl b/src/ParticleSystems/ParticleSystems.jl index 96d6521..10264b4 100644 --- a/src/ParticleSystems/ParticleSystems.jl +++ b/src/ParticleSystems/ParticleSystems.jl @@ -1,16 +1,18 @@ module ParticleSystems -export particle_waves_v3, particle_waves_v3beta, particle_waves_v4, particle_waves_v5 +export particle_waves_v3, particle_waves_v3beta, particle_waves_v4, particle_waves_v5, particle_waves_v6 # include("particle_waves_v3beta.jl") # include("particle_waves_v3.jl") # include("particle_waves_v4.jl") include("particle_waves_v5.jl") +include("particle_waves_v6.jl") # using .particle_waves_v3beta # using .particle_waves_v3 # using .particle_waves_v4 using .particle_waves_v5 +using .particle_waves_v6 end \ No newline at end of file diff --git a/src/ParticleSystems/particle_waves_v5.jl b/src/ParticleSystems/particle_waves_v5.jl index 4c3225f..2110003 100644 --- a/src/ParticleSystems/particle_waves_v5.jl +++ b/src/ParticleSystems/particle_waves_v5.jl @@ -4,7 +4,7 @@ using DifferentialEquations, IfElse using ...Architectures: AbstractODESettings, AbstractParticleSystem, IDConstantsInstance, ScgConstantsInstance -export particle_equations, ODESettings +export particle_equations, ODESettings, ODEParameters using LinearAlgebra using StaticArrays diff --git a/src/ParticleSystems/particle_waves_v6.jl b/src/ParticleSystems/particle_waves_v6.jl new file mode 100644 index 0000000..2871338 --- /dev/null +++ b/src/ParticleSystems/particle_waves_v6.jl @@ -0,0 +1,612 @@ +module particle_waves_v6 + +using DifferentialEquations, IfElse + +using ...Architectures: AbstractODESettings, AbstractParticleSystem + +export particle_equations, ODESettings +using LinearAlgebra +using StaticArrays + +using Parameters +using DocStringExtensions +#export t, x, y, e, c̄, φ_p, dist, Gₙ, u_10 + +# #using Plots +# #using PyPlot +# +# using PyCall +# #pygui(gui) #:tk, :gtk3, :gtk, :qt5, :qt4, :qt, or :wx +# using PyPlot + +# startupfile = joinpath(pwd(), "2022_particle_waves_startup.jl") +# isfile(startupfile) && include(startupfile) + +###### + + +""" +ODESettings +Structure to hold all information about the ODE system +# Fields +$(DocStringExtensions.FIELDS) +""" +@with_kw mutable struct ODESettings <: AbstractODESettings + "ODE parameters (Dict)" + Parameters::NamedTuple + "minimum allowed log energy on particle " + log_energy_minimum::Float64 + "maximum allowed log energy on particle " + log_energy_maximum::Float64 = log(17) + + "minimum needed squared wind velocity to seed particle" + wind_min_squared::Float64 = 4.0 + "solver method for ODE system" + #alternatives + #Rosenbrock23(), AutoVern7(Rodas4()) ,AutoTsit5(Rosenbrock23()) , Tsit5() + solver::Any = AutoTsit5(Rosenbrock23()) + "Internal saving timestep of the ODEs" + saving_step::Float64 + "remeshing time step, i.e. timestep of the model" + timestep::Float64 + + "Absolute allowed error" + abstol::Float64 = 1e-4 + "relative allowed error" + reltol::Float64 = 1e-3 + "max iteration for ODE solver (1e4)" + maxiters::Int = 1e4 + "Adaptive timestepping for ODE (true)" + adaptive::Bool = true + "timestep (if adaptive is false this is used), if adaptive is true this is the initial timestep" + dt::Float64 = 60 * 6 # seconds + "min timestep (if adaptive is true)" + dtmin::Float64 = 60 * 5 # seconds + "force min timestep (if adaptive is true)" + force_dtmin::Bool = false + + "Total time of the ODE integration, should not be needed, this is problematic .. " + total_time::Float64 + + "Callback function for ODE solver" + callbacks::Any = nothing + "save_everystep (false)" + save_everystep::Bool = false +end + + +# function init_vars() +# @variables t +# @variables x(t), y(t), c̄_x(t), c̄_y(t), lne(t), Δn(t), Δφ_p(t) +# @parameters r_g C_α C_φ g C_e +# return t, x, y, c̄_x, c̄_y, lne, Δn, Δφ_p, r_g, C_α, C_φ, g, C_e +# end + +# function init_vars_1D() +# @variables t +# @variables x(t), c̄_x(t), lne(t) +# @parameters r_g C_α g C_e +# return t, x, c̄_x, lne, r_g, C_α, g, C_e +# end + +# ------------------------------------------------------ +# Paramter functions +# ------------------------------------------------------ + +function magic_fractions(q::Float64=-1 / 4.0) + # returns universal exponent relations + p = (-1 - 10 * q) / 2 + n = 2 * q / (p + 4 * q) + [p, q, n] +end + +""" +get_I_D_constant(; c_D=2e-3, c_β=4e-2, c_e=1.3e-6, c_alpha= 11.8 , r_w = 2.35, q=-1/4) +this function returns a named tuple with the constants for the growth and dissipation +inputs: + c_D: drag coefficient + c_β: + c_e: + c_alpha: + r_w: + q: +returns: + c_D, c_β, c_e, c_alpha, r_w, C_e, γ, p, q, n +""" +function get_I_D_constant(; c_D=2e-3, c_β=4e-2, c_e=1.3e-6, c_alpha=11.8, r_w=2.35, q=-1 / 4) + p = (-1 - 10 * q) / 2 + n = 2 * q / (p + 4 * q) + + C_e = r_w * c_β * c_D + γ = (p - q) * c_alpha^(-4) * C_e^(-1) / 2 + return (c_D=c_D, c_β=c_β, c_e=c_e, c_alpha=c_alpha, r_w=r_w, C_e=C_e, γ=γ, p=p, q=q, n=n) +end + +""" +get_Scg_constants(C_alpha=1.41, C_varphi =1.81e-5) +this function returns a NamedTuple with constants for peak frequency shift. + C_alpha: 1.41 # constant for peak frequency shift (?) + C_varphi: 1.81e-5 # constant for peak frequency shift (?) +""" +function get_Scg_constants(; C_alpha=-1.41, C_varphi=1.81e-5) + return (C_alpha=C_alpha, C_varphi=C_varphi) +end + + + +# # test values +# u = (u=1, v=-1) +# cx = 1.0 +# cy = 1.0 + + + + +""" + αₚ(α::Number, φ::Number, φ_w::Number) + αₚ(α::Number, cφ_p::Number, sφ_p::Number, cφ_w::Number, sφ_w::Number) + αₚ(u::NamedTuple, cg::NamedTuple) + αₚ(u::NamedTuple, cx::Number, cy::Number ) + + returns angle between wave propagation direction and particle orientation. + cg, cx, and cy are the peak wave directions (!), not the mean wave direction! +""" +# αₚ(α::Number, φ::Number, φ_w::Number) = cos.(φ .- φ_w) .* α +#αₚ(α::Number, cφ_p::Number, sφ_p::Number, cφ_w::Number, sφ_w::Number) = (cφ_p .* cφ_w + sφ_p .* sφ_w) .* α +αₚ(u::NamedTuple, cg::NamedTuple) = αₚ(u.u, u.v, cg.cx, cg.cy) +αₚ(u::NamedTuple, cx::Number, cy::Number) = αₚ(u.u, u.v, cx, cy) +αₚ(u::Number, v::Number, cx::Number, cy::Number) = @. (u .* cx + v .* cy) ./ (2 .* max(speed(cx, cy), 1e-4)^2) + +#α_func(u_speed::Number, c_gp_speed::Number) = @. min( u_speed / (2.0 * c_gp_speed), 500) +function α_func(u_speed::Float64, c_gp_speed::Float64) + a = @. u_speed / (2.0 * c_gp_speed) + return IfElse.ifelse( a > 500, 500, a) + #return IfElse.ifelse(u_speed / (2.0 * c_gp_speed) > 500, 500, u_speed / (2.0 * c_gp_speed)) + # if u_speed / (2.0 * c_gp_speed) > 500 + # return 500.0 + # else + # return u_speed / (2.0 * c_gp_speed) + # end + # return min(u_speed / (2.0 * c_gp_speed), 500)::Float64 +end + +function α_func(u_speed, c_gp_speed) + a = @. u_speed / (2.0 * c_gp_speed) + return IfElse.ifelse(a > 500, 500, a) + #return IfElse.ifelse(u_speed / (2.0 * c_gp_speed) > 500, 500, u_speed / (2.0 * c_gp_speed)) + # if u_speed / (2.0 * c_gp_speed) > 500 + # return 500.0 + # else + # return u_speed / (2.0 * c_gp_speed) + # end + # return min(u_speed / (2.0 * c_gp_speed), 500)::Float64 +end + +#sin2_a_min_b(ca::Number, sa::Number, cb::Number, sb::Number) = 4 * sb * cb * (sa^2 -0.5) - 4 * sa * ca * (sb^2 -0.5) +# sign(cx) * +#sin2_a_min_b(ux::Number, uy::Number, cx::Number, cy::Number) = @. (2 / (speed(ux, uy) * speed(cx, cy))^2) * (ux * uy * (2 * cy^2 - speed(cx, cy)^2) - cx * cy * (2 * uy^2 - speed(ux, uy)^2)) +function sin2_a_min_b(ux::Number, uy::Number, cx::Number, cy::Number) + + IfElse.ifelse( + (speed(ux, uy) * speed(cx, cy)) == 0, + 0, + @. (2 / (speed(ux, uy) * speed(cx, cy))^2) * (ux * uy * (2 * cy^2 - speed(cx, cy)^2) - cx * cy * (2 * uy^2 - speed(ux, uy)^2)) + ) +end + +function sin2_a_min_b(u::NamedTuple, cx::Number, cy::Number) + sin2_a_min_b(u.u, u.v, cx, cy) +end +function sin2_a_min_b(u::NamedTuple, c::NamedTuple) + sin2_a_min_b(u.u, u.v, c.cx, c.cy) +end + +sin2_a_plus_b(ux::Number, uy::Number, cx::Number, cy::Number) = (2 / (speed(ux, uy) * speed(cx, cy))^2) * (cx * uy + cy * ux) * (cx * ux - cy * uy) +function sin2_a_plus_b(u::NamedTuple, cx::Number, cy::Number) + sin2_a_plus_b(u.u, u.v, cx, cy) +end +function sin2_a_plus_b(u::NamedTuple, c::NamedTuple) + sin2_a_plus_b(u.u, u.v, c.cx, c.cy) +end + + +#cos2_a_min_b(ca::Number, sa::Number, cb::Number, sb::Number) = (1 - 2 .* (sa .* cb - ca .* cb).^2 ) +e_T_func(γ::Float64, p::Float64, q::Float64, n::Float64; C_e::Number=2.16e-4, c_e::Float64=1.3e-6, c_α::Float64=11.8) = sqrt(c_e * c_α^(-p / q) / (γ * C_e)^(1 / n)) + + +H_β(α::Number, p::Float64; α_thresh::Float64=0.85) = @. 0.5 .* (1.0 + tanh.(p .* (α .- α_thresh))) # <--------------- tanh is slow +Δ_β(α::Number; α_thresh::Float64=0.85) = @. (1.0 .- 1.25 .* sech.(10.0 .* (α .- α_thresh)) .^ 2) + +""" +function c_g_conversions_vector(c̄::Number; g::Number=9.81, r_g::Number=0.9) +returns a vecotr with conversions between c̄, c_gp, kₚ, and ωₚ +""" +function c_g_conversions_vector(c̄::Number; g::Number=9.81, r_g::Number=0.9) # this is a slow function + c_gp = c_g_conversions(c̄, r_g=r_g) + kₚ = g / (4.0 * max(c_gp^2, 1e-2)) # < ---------- the power is slow + ωₚ = g / (2.0 * max(abs(c_gp), 0.1)) + #@SVector [c_gp, kₚ, ωₚ] + c_gp, kₚ, ωₚ +end + +function c_g_conversions(c̄::Float64; r_g::Number=0.9) # this is a slow function + c̄ / r_g +end + +function c_g_conversions(c̄; r_g::Number=0.9) # this is a slow function + c̄ / r_g +end + +speed(cx::Number, cy::Number) = @. sqrt(cx^2 + cy^2) + +function speed_and_angles(cx::Number, cy::Number) + #sqrt(cx.^2 + cy.^2), cx ./ sqrt(cx.^2 + cy.^2), cy ./sqrt(cx.^2 + cy.^2) + c = sqrt.(cx .^ 2 .+ cy .^ 2) + c, cx ./ c, cy ./ c +end + +function speed_and_angles(cx, cy) + #sqrt(cx.^2 + cy.^2), cx ./ sqrt(cx.^2 + cy.^2), cy ./sqrt(cx.^2 + cy.^2) + c = sqrt.(cx .^ 2 .+ cy .^ 2) + c, cx ./ c, cy ./ c +end + + +# ------------------------------- +# Forcing functions +# ------------------------------- + +# wind input +function Ĩ_func(alpha::Number, Hₚ::Number, C_e::Number) + # non-dimensional Wind energy input + # eq sec. 1.3 in the manual + C_e * Hₚ * alpha^2 +end + +# Dissipation +function D̃_func_e(e::Number, kₚ::Number, e_T::Number, n::Float64) + # non-dimensional Wind energy input + # eq sec. 1.3 in the manual + e .^ n .* (kₚ ./ e_T) .^ (2 * n) +end + +# disspation lne version +function D̃_func_lne(lne::Number, kₚ::Number, e_T::Number, n::Float64) + # non-dimensional Wind energy input + # eq sec. 1.3 in the manual + exp(n .* lne) .* (kₚ ./ e_T) .^ (2 * n) +end + + +# peak downshift +# C_α is negative in Kudravtec definition, here its a positive value +S_cg(lne::Number, Δₚ::Number, kₚ::Number, C_α::Number) = @. C_α * Δₚ * kₚ^4 * exp(2 * lne) + + +# Peak direction shift +function S_dir(u::Number, v::Number, cx::Number, cy::Number, C_φ::Number, Hₚ::Number) + return α_func(speed(u, v), speed(cx, cy))^2 * C_φ * Hₚ * sin2_a_min_b(u, v, cx, cy)#::Number + #alpha^2 * C_φ * Hₚ * sin2_a_min_b(cx, cy, u.u, u.v) + + #alpha^2 * C_φ * sin2_a_min_b(u, cx, cy) + #alpha^2 * C_φ * Hₚ * sin2_a_plus_b(u, cx, cy) +end + + +# function Gₙ( dphi_p::Number, dn::Number, dn_0::Number; dg_ratio::Float64 = 0.21 ) +# return (dphi_p / dn_0) * ( (dn /dn_0) / ( (dn/dn_0)^2 + (dg_ratio/2.0)^2 ) ) +# end + + +# tΔφ_w_func(c_x::Number, c_y::Number, u_x::Number, v_y::Number) = (c_x .* v_y - c_y .* u_x ) / ( c_x .* u_x + c_y .* v_y ) +# Δφ_p_RHS(tΔφ_w::Number, tΔφ_p::Number , T_inv::Number) = ( tΔφ_w - tΔφ_p ) * T_inv / ( 1 + tΔφ_w * tΔφ_p ) + + +""" +particle_equations(u , v ; γ::Number=0.88, q::Number=-1/4.0 ) +Particle wave equations in 2D +inputs: + u: zonal wind forcing field + v: meridional wind forcing field + γ: 0.88 + q: -1/4.0 + propagation: true + input: true + dissipation: true + peak_shift: true + direction: true + debug_output: false + +return an ODE system as function particle_system(dz, z, params, t) that provides 5 element state vector tendency of: + z = [lne, c̄_x, c̄_y, x, y] + params can be a named tuple with the parameters or a vector + params = [r_g, C_α, g, C_e] +""" +function particle_equations(u_wind, v_wind; γ::Number=0.88, q::Number=-1 / 4.0, + propagation=true, + input=true, + dissipation=true, + peak_shift=true, + direction=true, + debug_output=false, + static= false + ) + #t, x, y, c̄_x, c̄_y, lne, Δn, Δφ_p, r_g, C_α, C_φ, g, C_e = init_vars() + + p, q, n = magic_fractions(q) + e_T = e_T_func(γ, p, q, n)#, C_e=C_e) + + if static + + function particle_system_static(z, params, t) + + + # forcing fields + #lne, c̄_x, c̄_y, x, y = z + + r_g, C_α, g, C_e, C_φ = params.r_g, params.C_α, params.g, params.C_e, params.C_φ + #u = (u=u, v=v)::NamedTuple{(:u, :v),Tuple{Number,Number}} + u = u_wind(z[4], z[5], t)#::Number + v = v_wind(z[4], z[5], t)#::Number + + c̄ = speed(z[2], z[3]) + u_speed = speed(u, v) + + # peak parameters + c_gp_speed, kₚ, ωₚ = c_g_conversions_vector(abs(c̄), r_g=r_g) + c_gp_x = c_g_conversions(z[2], r_g=r_g) + c_gp_y = c_g_conversions(z[3], r_g=r_g) + + # direction equations + α = α_func(u_speed, c_gp_speed) + Hₚ = H_β(αₚ(u, v, c_gp_x, c_gp_y), p) + Δₚ = Δ_β(αₚ(u, v, c_gp_x, c_gp_y)) + + # Source terms + Ĩ = input ? Ĩ_func(α, Hₚ, C_e) : 0.0 + D̃ = dissipation ? D̃_func_lne(z[1], kₚ, e_T, n) : 0.0 + S_cg_tilde = peak_shift ? S_cg(z[1], Δₚ, kₚ, C_α) : 0.0 + S_dir_tilde = direction ? S_dir(u, v, c_gp_x, c_gp_y, C_φ, Hₚ) : 0.0 + + z = @SVector [ + # energy + +ωₚ .* r_g^2 .* S_cg_tilde + ωₚ .* (Ĩ - D̃), + + # peak group velocity vector + -z[2] .* ωₚ .* r_g^2 .* S_cg_tilde + z[3] .* S_dir_tilde, + -z[3] .* ωₚ .* r_g^2 .* S_cg_tilde - z[2] .* S_dir_tilde, + + # D(z[2]) ~ -z[2] .* ωₚ .* r_g^2 .* S_cg_tilde + (z[3] + 0.001) .* S_dir_tilde, #* (-1), + # D(z[3]) ~ -z[3] .* ωₚ .* r_g^2 .* S_cg_tilde - (z[2] + 0.001) .* S_dir_tilde, #* (1), + + # propagation + propagation ? z[2] : 0.0, + propagation ? z[3] : 0.0 + ] + + if debug_output + additional_output = @SVector [ + Ĩ, + -D̃, + r_g^2 * S_cg_tilde, + Hₚ, + S_dir_tilde, + Δₚ, + c_gp_y + ] + z = vcat(z, additional_output) + end + + return z + + end + return particle_system_static + + else + + function particle_system(dz, z, params, t)#::MVector{5, Number} + + # forcing fields + #u = (u=u(x, y, t), v=v(x, y, t))::NamedTuple{(:u, :v),Tuple{Number,Number}} + lne, c̄_x, c̄_y, x, y, angular_σ = z + + r_g, C_α, g, C_e, C_φ = params.r_g, params.C_α, params.g, params.C_e, params.C_φ + #u = (u=u, v=v)::NamedTuple{(:u, :v),Tuple{Number,Number}} + u = u_wind(z[4], y, t)#::Number + v = v_wind(x, y, t)#::Number + + c̄ = speed(c̄_x, c̄_y) + u_speed = speed(u, v) + + # peak parameters + c_gp_speed, kₚ, ωₚ = c_g_conversions_vector(abs(c̄), r_g=r_g) ## <--- this one is slow!! + c_gp_x = c_g_conversions(c̄_x, r_g=r_g) + c_gp_y = c_g_conversions(c̄_y, r_g=r_g) + + # direction equations + α = α_func(u_speed, c_gp_speed) + Hₚ = H_β(αₚ(u, v, c_gp_x, c_gp_y), p) + Δₚ = Δ_β(αₚ(u, v, c_gp_x, c_gp_y)) + + # Source terms + Ĩ = input ? Ĩ_func(α, Hₚ, C_e) : 0.0 + D̃ = dissipation ? D̃_func_lne(lne, kₚ, e_T, n) : 0.0 + S_cg_tilde = peak_shift ? S_cg(lne, Δₚ, kₚ, C_α) : 0.0 + S_dir_tilde = direction ? S_dir(u, v, c_gp_x, c_gp_y, C_φ, Hₚ) : 0.0 + + #particle_equations::Vector{Number} = [ + # energy + dz[1] = +ωₚ .* r_g .* S_cg_tilde + ωₚ .* (Ĩ - D̃) #- c̄ .* G_n, + + # peak group velocity vector + # @info "" + # @info "c̄_x = " * string(c̄_x) + # @info "ωₚ = " * string(ωₚ) + # @info "r_g = " * string(r_g) + # @info "S_cg_tilde = " * string(S_cg_tilde) + # @info "c̄_y = " * string(c̄_y) + # @info "S_dir_tilde = " * string(S_dir_tilde) + # total_x = -c̄_x .* ωₚ .* r_g .* S_cg_tilde + c̄_y .* S_dir_tilde + # total_y = -c̄_y .* ωₚ .* r_g .* S_cg_tilde - c̄_x .* S_dir_tilde + # total_tot = sqrt(total_x^2+total_y^2) + # @info "total_x = " * string(total_x) + # @info "total_y = " * string(total_y) + # @info "total_tot = " * string(total_tot) + # @info "" + dz[2] = -c̄_x .* ωₚ .* r_g .* S_cg_tilde + c̄_y .* S_dir_tilde #* (-1), + dz[3] = -c̄_y .* ωₚ .* r_g .* S_cg_tilde - c̄_x .* S_dir_tilde #* (1), + + # D(c̄_x) ~ -c̄_x .* ωₚ .* r_g^2 .* S_cg_tilde + (c̄_y + 0.001) .* S_dir_tilde, #* (-1), + # D(c̄_y) ~ -c̄_y .* ωₚ .* r_g^2 .* S_cg_tilde - (c̄_x + 0.001) .* S_dir_tilde, #* (1), + + # propagation + dz[4] = propagation ? c̄_x : 0.0 + dz[5] = propagation ? c̄_y : 0.0 + dz[6] = 0 + + + if debug_output + additional_output = [ + Ĩ, + -D̃, + r_g^2 * S_cg_tilde, + #alpha_p ~ αₚ(u, c_gp_x, c_gp_y), + Hₚ, + #alpha ~ α, + S_dir_tilde, + Δₚ, + c_gp_y] + append!(dz, additional_output) + end + + return dz + + end + + return particle_system + + end + + +end + +# ------------ 1D ------------ +""" +particle_equations(u ; γ::Number=0.88, q::Number=-1/4.0 ) +Particle wave equations in 2D +inputs: + u: forcing field + γ: 0.88 + q: -1/4.0 + propagation: true + input: true + dissipation: true + peak_shift: true + info: false + +returns an ODE system as function particle_system(dz, z, params, t) that provides 3 element state vector tendency of: + z = [lne, c̄_x, x] + params can be a named tuple with the parameters or a vector + params = [r_g, C_α, g, C_e] +""" +function particle_equations(u_wind; γ::Number=0.88, q::Number=-1 / 4.0, + propagation=true, + input=true, + dissipation=true, + peak_shift=true, info=false) + + #t, x, c̄_x, lne, r_g, C_α, g, C_e = init_vars_1D() + + # Define basic constants for wave equation (invariant throughout the simulation) + p, q, n = magic_fractions(q) + e_T = e_T_func(γ, p, q, n)#, C_e=C_e) + #D = Differential(t) + + ### ---------------- start function here + function partice_system(dz, z, params, t) #<: Vector{Number} + #unpack0 + lne, c̄_x, x = z + #lne, c̄_x, x = z.lne, z.c̄_x, z.x + #r_g, C_α, g, C_e = params + r_g, C_α, g, C_e = params.r_g, params.C_α, params.g, params.C_e + + # forcing fields, need to be global scope? + u = u_wind(x, t) + #u = (u=u2(x, y, t), v=v2(x, y, t)) + + # trig-values # we only use scalers, not vectors + c̄ = c̄_x + u_speed = abs(u) + + # peak parameters + c_gp_speed, kₚ, ωₚ = c_g_conversions_vector(abs(c̄), r_g=r_g) + + # direction equations + α = α_func(u_speed, c_gp_speed) + Hₚ = H_β(α, p) + Δₚ = Δ_β(α) + + # Source terms + Ĩ = input ? Ĩ_func(α, Hₚ, C_e) : 0.0 + D̃ = dissipation ? D̃_func_lne(lne, kₚ, e_T, n) : 0.0 + S_cg_tilde = peak_shift ? S_cg(lne, Δₚ, kₚ, C_α) : 0.0 + c_group_x = propagation ? c̄_x : 0.0 + # no directional changes in 1D! + if info + println("alpha = ", α) + println("Hₚ = ", Hₚ) + println("Δₚ = ", Δₚ) + println("I_tilde = ", Ĩ) + println("D_tilde = ", D̃) + println("S_cg_tilde = ", S_cg_tilde) + println("c_tilde_x = ", c_group_x) + end + + # energy + dz[1] = +ωₚ .* r_g^2 .* S_cg_tilde + ωₚ .* (Ĩ - D̃) #- c̄ .* G_n, + #e * ωₚ .* (Ĩ - D̃)- e^3 *ξ / c̄ , + + # peak group velocity vector + dz[2] = - c̄_x .* ωₚ .* r_g^2 .* S_cg_tilde + + # propagation + dz[3] = c_group_x + + return dz + + end + + return partice_system +end + + +""" +particle_rays(u ; γ::Number=0.88, q::Number=-1/4.0 ) +Particle wave equations in 2D + +returns 3 element state vector of the particle system as vector: + [lne, c̄_x, x], lne, c̄_x, are constants +""" +function particle_rays(info=false) + + # no directional changes in 1D! + if info + println("c_x = ", c̄_x) + end + #D = Differential(t) + function partice_system(dz, z, params, t) #<: Vector{Number} + lne, c̄_x, x = z + # energy + dz[1] = 0 + # peak group velocity vector + dz[2] = 0 + # propagation + dz[3] = c̄_x + end + + return partice_system +end + + +# end of module +end diff --git a/src/PiCLES.jl b/src/PiCLES.jl index b68ff9d..8069985 100644 --- a/src/PiCLES.jl +++ b/src/PiCLES.jl @@ -26,7 +26,7 @@ export Grids, # utils - Utils, Debugging, FetchRelations, ParticleTools, WindEmulator + Utils, Debugging, FetchRelations, ParticleTools, WindEmulator, visualization #externals @@ -43,14 +43,14 @@ using .ParticleSystems include("FetchRelations.jl") using .FetchRelations +include("Grids/Grids.jl") +using .Grids + include("ParticleMesh.jl") include("ParticleInCell.jl") using .ParticleMesh using .ParticleInCell -include("Grids/Grids.jl") -using .Grids - include("Operators/Operators.jl") include("Simulations/Simulations.jl") diff --git a/src/Simulations/run.jl b/src/Simulations/run.jl index 5776fcb..17c5289 100644 --- a/src/Simulations/run.jl +++ b/src/Simulations/run.jl @@ -1,10 +1,13 @@ +using ..Operators.core_2D_spread: SeedParticle as StochasticSeedParticle2D +using ..Operators.core_2D: SeedParticle as SeedParticle2D +using ..Operators.core_2D_spread: ParticleDefaults as ParticleDefaults2D using ..Operators.core_1D: ParticleDefaults as ParticleDefaults1D -using ..Operators.core_2D: ParticleDefaults as ParticleDefaults2D using ..Operators.core_1D: SeedParticle! as SeedParticle1D! -using ..Operators.core_2D: SeedParticle as SeedParticle2D +using ..Operators.core_2D_spread: SeedParticle! as SeedParticle2D! +# using ..Operators.core_2D: SeedParticle -using ..Architectures: Abstract2DModel, Abstract1DModel +using ..Architectures: Abstract2DModel, Abstract1DModel, Abstract2DStochasticModel using ..ParticleMesh: OneDGrid, OneDGridNotes, TwoDGrid, TwoDGridNotes #using WaveGrowthModels: init_particles! @@ -17,6 +20,11 @@ using Statistics using StructArrays +import Plots as plt + +using DataFrames, CSV +using Random, Distributions +using Dates #using ThreadsX @@ -28,12 +36,100 @@ function mean_of_state(model::Abstract1DModel) return mean(model.State[:, 1]) end +function plot_state_and_error_points(wave_simulation, gn) + plt.plot() + + X = (0:(gn.stats.Nx.N-1)).*gn.stats.dx .+ gn.stats.xmin + Y = (0:(gn.stats.Ny.N-1)).*gn.stats.dy .+ gn.stats.ymin + energy = get_tot_energy_domain(wave_simulation) + p1 = plt.heatmap(X, Y, transpose(wave_simulation.model.State[:, :, 1]), aspect_ratio=:equal, size=(1080, 1080))#,clim=(0,1)) + + plt.plot!(legend=:none, + title="total energy = "*string(round(energy,digits=3))*"; max = "*string(round(maximum(wave_simulation.model.State[:,:,1]), + digits=3))*"; pos = ("*string(argmax(wave_simulation.model.State[:,:,1])[1])*","* + string(argmax(wave_simulation.model.State[:,:,1])[2])*")", + ylabel="y position", + xlabel="x position", + xlims=(gn.stats.xmin, gn.stats.xmax), + ylims=(gn.stats.ymin, gn.stats.ymax)) |> display +end + +function write_particles_to_csv(wave_model) + sec=string(Int64(floor((wave_model.clock.time)/60))) + dec=string(Int64(floor(10*(wave_model.clock.time/60-floor((wave_model.clock.time)/60))))) + save_path = wave_model.plot_savepath + + nParticles = wave_model.n_particles_launch + # @info wave_model.clock.time/60 + + parts = wave_model.ParticleCollection[(end-nParticles+1):end] + + logE = zeros(nParticles) + cx = zeros(nParticles) + cy = zeros(nParticles) + x = zeros(nParticles) + y = zeros(nParticles) + for i in 1:nParticles + logE[i] = parts[i].ODEIntegrator[1] + cx[i] = parts[i].ODEIntegrator[2] + cy[i] = parts[i].ODEIntegrator[3] + x[i] = parts[i].ODEIntegrator[4] + y[i] = parts[i].ODEIntegrator[5] + end + + data = DataFrame(id=1:nParticles, logE = logE, cx = cx, cy = cy, x = x, y = y) + data2 = Tables.table(transpose(wave_model.State[:, :, 1])) + CSV.write(save_path*"/data/particles_"*sec*","*dec*".csv", data) + CSV.write(save_path*"/data/mesh_values_"*sec*","*dec*".csv", data2) + end + +function get_tot_energy_domain(wave_simulation) + return sum(wave_simulation.model.State[:,:,1]) +end + """ run!(sim::Simulation; store = false, pickup=false) main method to run the Simulation sim. Needs time_step! to be defined for the model, and push_state_to_storage! to be defined for the store. """ function run!(sim; store=false, pickup=false, cash_store=false, debug=false) + if sim.model isa Abstract2DStochasticModel + save_path = sim.model.plot_savepath + + if sim.model.save_particles + save_path = sim.model.plot_savepath + filename = save_path*"/data/simu_info.csv" + + Nx = sim.model.grid.stats.Nx.N + Ny = sim.model.grid.stats.Ny.N + xmin = sim.model.grid.stats.xmin + xmax = sim.model.grid.stats.xmax + ymin = sim.model.grid.stats.ymin + ymax = sim.model.grid.stats.ymax + lne_source = sim.model.ODEdefaults.lne + c_x_source = sim.model.ODEdefaults.c̄_x + c_y_source = sim.model.ODEdefaults.c̄_y + x_source = sim.model.ODEdefaults.x + y_source = sim.model.ODEdefaults.y + angular_spread_source = sim.model.ODEdefaults.angular_σ + Δt = sim.Δt + stop_time = sim.stop_time + + data = DataFrame(Nx=Nx, Ny=Ny, xmin=xmin, + xmax=xmax, ymin=ymin, ymax=ymax, + lne_source=lne_source, c_x_source=c_x_source, + c_y_source=c_y_source, x_source=x_source,y_source=y_source, + angular_spread_source=angular_spread_source, + Δt=Δt, stop_time=stop_time) + CSV.write(filename, data) + + filename2 = save_path*"/data/sigma.csv" + + covariance_init = sim.model.proba_covariance_init + data2 = DataFrame(covariance_init, :auto) + CSV.write(filename2, data2) + end + end start_time_step = time_ns() @@ -41,6 +137,12 @@ function run!(sim; store=false, pickup=false, cash_store=false, debug=false) initialize_simulation!(sim) end + if sim.model isa Abstract2DStochasticModel + if sim.model.save_particles && length(sim.model.ParticleCollection) > 0 + write_particles_to_csv(sim.model) + end + end + #sim.running = true sim.run_wall_time = 0.0 @@ -78,9 +180,17 @@ function run!(sim; store=false, pickup=false, cash_store=false, debug=false) sim.model.State .= 0.0 end # do time step - + + # l = length(sim.model.ParticleCollection) + # Uxs = [sim.model.ParticleCollection[k].ODEIntegrator.u[2] for k in 1:l] + # Uys = [sim.model.ParticleCollection[k].ODEIntegrator.u[3] for k in 1:l] + # maxUxs = maximum(Uxs) + # maxUys = maximum(Uys) + # @info "maximum(Uxs) = " * string(maxUxs) * " and maximum(Uys) = " * string(maxUys) + # @info "percentage of cell travelled : on x = " * string(maxUxs * sim.Δt / sim.model.grid.stats.dx) * " and on y = " * string(maxUys * sim.Δt / sim.model.grid.stats.dy) + time_step!(sim.model, sim.Δt, debug=debug) - + if debug & (length(sim.model.FailedCollection) > 0) @info "debug mode:" @info "found failed particles" @@ -111,7 +221,15 @@ function run!(sim; store=false, pickup=false, cash_store=false, debug=false) end sim.running = sim.stop_time >= sim.model.clock.time ? true : false - @info sim.model.clock + + if sim.model isa Abstract2DStochasticModel + if sim.model.plot_steps + plot_state_and_error_points(sim, sim.model.grid) + sec=string(Int64(floor((sim.model.clock.time)/60))) + dec=string(Int64(floor(10*(sim.model.clock.time/60-floor((sim.model.clock.time)/60))))) + plt.savefig(joinpath([save_path*"/plots/", "energy_plot_no_spread_"*sec*","*dec*".png"])) + end + end end end_time_step = time_ns() @@ -127,8 +245,8 @@ initialize_simulation!(sim::Simulation) initialize the simulation sim by calling init_particles! to initialize the model.ParticleCollection. -particle_initials::T=nothing was removed from arguments """ -function initialize_simulation!(sim::Simulation)# where {PP<:Union{ParticleDefaults,Nothing}} - # copy(ParticleDefaults(log(4e-8), 1e-2, 0.0))) +function initialize_simulation!(sim::Simulation)# where {PP<:Union{ParticleDefaults2D,Nothing}} + # copy(ParticleDefaults2D(log(4e-8), 1e-2, 0.0))) if sim.verbose @info "init particles..." @@ -151,7 +269,7 @@ reset_simulation!(sim::Simulation) reset the simulation sim by calling init_particles! to reinitialize the model.ParticleCollection, sets the model.clock.time, model.clock.iteration, and model.state to 0. - particle_initials::Dict{Num, Float64} was removed from arguments """ -function reset_simulation!(sim::Simulation)# where {PP<:Union{ParticleDefaults,Nothing}} +function reset_simulation!(sim::Simulation)# where {PP<:Union{ParticleDefaults2D,Nothing}} sim.running = false sim.run_wall_time = 0.0 @@ -196,6 +314,105 @@ initialize the model.ParticleCollection based on the model.grid and the defaults If defaults is nothing, then the model.ODEdev is used. usually the initilization uses wind constitions to seed the particles. """ +function init_particles!(model::Abstract2DStochasticModel; defaults::PP=nothing, verbose::Bool=false) where {PP<:Union{ParticleDefaults1D,ParticleDefaults2D,Array{Any,1},Nothing}} + #defaults = isnothing(defaults) ? model.ODEdev : defaults + if verbose + @info "seed PiCLES ... \n" + @info "defaults is $(defaults)" + if defaults isa Dict + @info "found particle initials, just replace position " + else + @info "no particle defaults found, use windsea to seed particles" + end + end + + # gridnotes = TwoDGridNotes(model.grid) + + # SeedParticle_i = SeedParticle_mapper(SeedParticle2D!, + # ParticleCollection, model.State, + # model.ODEsystem, nothing, model.ODEsettings, + # gridnotes, model.winds, model.ODEsettings.timestep, + # model.boundary, model.periodic_boundary) + + ParticleCollection = [] + model.ParticleCollection = ParticleCollection + + if defaults isa ParticleDefaults2D + i = Int64(floor((defaults.x - model.grid.stats.xmin) / model.grid.stats.dx)) + 1 + j = Int64(floor((defaults.y - model.grid.stats.ymin) / model.grid.stats.dy)) + 1 + # gridnotes = TwoDGridNotes(model.grid) + if model.angular_spreading_type == "nonparametric" + if sum(model.proba_covariance_init)==4e-50 + # if "nonparametric" is used, initialize a bigger number of particles at the original perturbation + n_part = model.n_particles_launch + ij = CartesianIndices((i,j)) + ij_mesh = model.grid.data[ij] + ij_wind = (model.winds.u(ij_mesh.x, ij_mesh.y, 0.0), + model.winds.v(ij_mesh.x, ij_mesh.y, 0.0) + ) + for _ in 1:n_part + defaults_temp = deepcopy(defaults) + delta_phi = rand() * defaults_temp.angular_σ - 0.5*defaults_temp.angular_σ + c_x = defaults_temp.c̄_x * cos(delta_phi) - defaults_temp.c̄_y * sin(delta_phi) + c_y = defaults_temp.c̄_x * sin(delta_phi) + defaults_temp.c̄_y * cos(delta_phi) + defaults_temp.lne += -log(n_part) + defaults_temp.c̄_x = c_x + defaults_temp.c̄_y = c_y + push!(ParticleCollection, StochasticSeedParticle2D(model.State, + (i,j), model.ODEsystem, defaults_temp, + model.ODEsettings, + model.grid.stats,model.grid.ProjetionKernel,model.grid.PropagationCorrection, + ij_mesh[i,j], ij_wind, model.ODEsettings.timestep, model.boundary, + model.periodic_boundary)) + end + else + n_part = model.n_particles_launch + for _ in 1:n_part + defaults_temp = deepcopy(defaults) + mu = [defaults_temp.c̄_x, defaults_temp.c̄_y, defaults_temp.x, defaults_temp.y] + d = MvNormal(mu, model.proba_covariance_init) + real = rand(d,1) + # delta_phi = real[1] + # delta_phi < 0 ? delta_phi = - delta_phi : delta_phi = delta_phi + # delta_phi < 0.01 ? delta_phi+=1 : delta_phi +=0 + # c_x = (real[2]+1)*(defaults_temp.c̄_x * cos(delta_phi) - defaults_temp.c̄_y * sin(delta_phi)) + # c_y = (real[2]+1)*(defaults_temp.c̄_x * sin(delta_phi) + defaults_temp.c̄_y * cos(delta_phi)) + c_x = real[1] + c_y = real[2] + defaults_temp.lne += -log(n_part) + defaults_temp.c̄_x = c_x + defaults_temp.c̄_y = c_y + defaults_temp.x = real[3] + defaults_temp.y = real[4] + push!(ParticleCollection, SeedParticle(model.State, + (i,j), model.ODEsystem, defaults_temp, + model.ODEsettings,gridnotes, model.winds, + model.ODEsettings.timestep, model.boundary, + model.periodic_boundary)) + end + end + else + push!(ParticleCollection, SeedParticle(model.State, + (i,j), model.ODEsystem, defaults, + model.ODEsettings,gridnotes, model.winds, + model.ODEsettings.timestep, model.boundary, + model.periodic_boundary)) + end + elseif defaults isa Array{Any,1} + for k in 1:length(defaults) + i = Int64(floor((defaults[k].x - model.grid.stats.xmin) / model.grid.stats.dx)) + 1 + j = Int64(floor((defaults[k].y - model.grid.stats.ymin) / model.grid.stats.dy)) + 1 + gridnotes = OneDGridNotes(model.grid) + push!(ParticleCollection, SeedParticle(model.State, + (i,j), model.ODEsystem, defaults[k], + model.ODEsettings,gridnotes, model.winds, + model.ODEsettings.timestep, model.boundary, + model.periodic_boundary)) + end + end + nothing +end + function init_particles!(model::Abstract2DModel; defaults::PP=nothing, verbose::Bool=false) where {PP<:Union{ParticleDefaults1D,ParticleDefaults2D,Nothing}} #defaults = isnothing(defaults) ? model.ODEdev : defaults if verbose @@ -248,8 +465,6 @@ end - - ### 1D version ### # """ # SeedParticle_mapper(f, p, s, b1, b2, b3, c1, c2, c3, c4, d1, d2 ) = x -> f( p, s, x, b1, b2, b3, c1, c2, c3, c4, d1, d2 ) @@ -286,7 +501,7 @@ function init_particles!(model::Abstract1DModel; defaults::PP=nothing, verbose:: gridnotes, model.winds, model.ODEsettings.timestep, model.boundary, model.periodic_boundary) - map(SeedParticle_i, range(1, length=model.grid.Nx)) + map(SeedParticle_i, range(1, length=model.grid.stats.Nx)) # print(defaults) diff --git a/src/custom_structures.jl b/src/custom_structures.jl index 9b0b63a..a78c47e 100644 --- a/src/custom_structures.jl +++ b/src/custom_structures.jl @@ -1,12 +1,12 @@ module custom_structures -export ParticleInstance1D, ParticleInstance2D, MarkedParticleInstance, AbstractParticleInstance, AbstractMarkedParticleInstance, wni +export ParticleInstance1D, ParticleInstance2D, MarkedParticleInstance, AbstractParticleInstance, AbstractMarkedParticleInstance, StochasticParticleInstance2D, wni using DifferentialEquations: OrdinaryDiffEq.ODEIntegrator using DocStringExtensions using StaticArrays -using ..Architectures: AbstractParticleInstance, AbstractMarkedParticleInstance +using ..Architectures: AbstractParticleInstance, AbstractMarkedParticleInstance, AbstractStochasticParticleInstance using ..Architectures: AbstractBoundary # ParticleInstance is the Stucture that carries each particle. @@ -18,6 +18,14 @@ mutable struct ParticleInstance2D <: AbstractParticleInstance on::Bool end +mutable struct StochasticParticleInstance2D <: AbstractStochasticParticleInstance + position_ij::Tuple{Int, Int} + position_xy::Tuple{Float64, Float64} + ODEIntegrator::Union{ODEIntegrator,Nothing} + boundary :: Bool + on::Bool +end + mutable struct ParticleInstance1D <: AbstractParticleInstance position_ij::Int position_xy::Float64 diff --git a/tests/T04_2D_box_rotation.jl b/tests/T04_2D_box_rotation.jl index 712444d..25f4896 100644 --- a/tests/T04_2D_box_rotation.jl +++ b/tests/T04_2D_box_rotation.jl @@ -23,6 +23,7 @@ import Oceananigans.Utils: prettytime using PiCLES.Architectures using GLMakie +using Revise using PiCLES.Operators.core_2D: GetGroupVelocity, speed using PiCLES.Plotting.movie: init_movie_2D_box_plot, init_movie_2D_box_plot_small diff --git a/tests/test_case_1.jl b/tests/test_case_1.jl new file mode 100644 index 0000000..152599a --- /dev/null +++ b/tests/test_case_1.jl @@ -0,0 +1,125 @@ +ENV["JULIA_INCREMENTAL_COMPILE"]=true +using Pkg +Pkg.activate(".") + +import Plots as plt +using Setfield, IfElse + +using PiCLES.ParticleSystems: particle_waves_v6 as PW + +import PiCLES: FetchRelations, ParticleTools +using PiCLES.Operators.core_2D_spread: ParticleDefaults, InitParticleInstance, GetGroupVelocity, PointSource +using PiCLES.Operators: TimeSteppers +using PiCLES.Simulations +using PiCLES.Operators.TimeSteppers: time_step!, movie_time_step! + +using PiCLES.ParticleMesh: TwoDGrid, TwoDGridNotes, TwoDGridMesh +using PiCLES.Grids.CartesianGrid: TwoDCartesianGridMesh, TwoDCartesianGridStatistics +using PiCLES.Models.GeometricalOpticsModels + +using Oceananigans.TimeSteppers: Clock, tick! +import Oceananigans: fields +using Oceananigans.Units +import Oceananigans.Utils: prettytime + +using PiCLES.Architectures + +using PiCLES.Operators.core_2D_spread: GetGroupVelocity, speed + +using Revise + +save_path = "plots/test_case_1/" +mkpath(save_path) +mkpath(save_path*"data/") +mkpath(save_path*"plots/") +touch(save_path*"data/simu_info.csv") +touch(save_path*"data/sigma.csv") + +# % Parameters +U10,V10 = 10.0, 10.0 +dt_ODE_save = 30minutes +DT = 2minutes +# version 3 +r_g0 = 0.85 + +# function to define constants +Const_ID = PW.get_I_D_constant() +@set Const_ID.γ = 0.88 +Const_Scg = PW.get_Scg_constants(C_alpha=- 1.41, C_varphi=1.81e-5) + +wind_angle = 3/8*pi +wind_magnitude = 20 + +u_func(x, y, t) = wind_magnitude*cos(wind_angle) +v_func(x, y, t) = wind_magnitude*sin(wind_angle) + +u(x, y, t) = u_func(x, y, t) +v(x, y, t) = v_func(x, y, t) +winds = (u=u, v=v) + +typeof(winds.u) +typeof(winds.u(1e3, 1e3, 11)) + +grid =TwoDCartesianGridMesh(22.5e4, 151, 22.5e4, 151; periodic_boundary=(false, true)) +# mesh = TwoDGridMesh(grid, skip=1); +# gn = TwoDGridNotes(grid); + +Revise.retry() + +# define variables based on particle equation + +particle_system = PW.particle_equations(u, v, γ=0.88, q=Const_ID.q, input=true, dissipation=false, peak_shift=true); + +# define V4 parameters absed on Const NamedTuple: +default_ODE_parameters = (r_g = r_g0, C_α = Const_Scg.C_alpha, + C_φ = Const_ID.c_β, C_e = Const_ID.C_e, g= 9.81 ); +#plt.scalefontsizes(1.75) + +# define setting and standard initial conditions +WindSeamin = FetchRelations.MinimalWindsea(U10, V10, DT ); +lne_local = log(WindSeamin["E"]) + +ODE_settings = PW.ODESettings( + Parameters=default_ODE_parameters, + # define mininum energy threshold + log_energy_minimum=lne_local, + #maximum energy threshold + log_energy_maximum=log(27), + saving_step=dt_ODE_save, + timestep=DT, + total_time=T=6days, + adaptive=true, + dt=1e-3, #60*10, + dtmin=1e-4, #60*5, + force_dtmin=true, + callbacks=nothing, + save_everystep=false) + + +default_particle = ParticleDefaults(5, 0.0, 5.0, 112500., 112500., 2*π) +Revise.retry() + +wave_model = GeometricalOpticsModels.GeometricalOptics(; grid=grid, + winds=winds, + ODEsys=particle_system, + ODEsets=ODE_settings, # ODE_settings + ODEinit_type=default_particle, + periodic_boundary=false, + boundary_type="same", + movie=true, + plot_steps=true, + save_particles=true, + plot_savepath="plots/test_case_1", + angular_spreading_type="nonparametric", + proba_covariance_init = [1e-50 0 0 0; + 0 1e-50 0 0; + 0 0 1e-50 0 + 0 0 0 1e-50], + n_particles_launch=150 + ) + +wave_simulation = Simulation(wave_model, Δt=0.75minutes, stop_time=75minutes) +initialize_simulation!(wave_simulation) + + +@time run!(wave_simulation, cash_store=false, debug=false) diff --git a/tests/test_case_original_method.jl b/tests/test_case_original_method.jl new file mode 100644 index 0000000..c45f9bd --- /dev/null +++ b/tests/test_case_original_method.jl @@ -0,0 +1,131 @@ +using Pkg +# This will be replaced by the module load in the future +Pkg.activate(".") # Activate the PiCLES package + +using PiCLES +using PiCLES.Operators.core_2D: ParticleDefaults +using PiCLES.Models.WaveGrowthModels2D: WaveGrowth2D +using PiCLES.Simulations +using PiCLES.Grids.CartesianGrid: TwoDCartesianGridMesh, ProjetionKernel, TwoDCartesianGridStatistics + +using PiCLES.ParticleSystems: particle_waves_v5 as PW +using Oceananigans.Units + +# just for simple plotting +import Plots as plt + +# Parameters +U10, V10 = 20.0, 10.0 +DT = 10minutes +r_g0 = 0.85 # ratio of c / c_g (phase velocity/ group velocity). + +# Define wind functions +function ind(x,a,b) + if x>= a && x