diff --git a/examples/exmpl_homogenous_box_mprocess.jl b/examples/exmpl_homogenous_box_mprocess.jl deleted file mode 100644 index 568cd09..0000000 --- a/examples/exmpl_homogenous_box_mprocess.jl +++ /dev/null @@ -1,315 +0,0 @@ -#JULIA_NUM_THREADS = 4 -#ENV["JULIA_NUM_THREADS"] = 4 - -using Base.Threads -nthreads() - -import Plots as plt - -# CPU initialization -using Distributed -addprocs(4) - - -@everywhere begin - using Setfield - using PiCLES.ParticleSystems: particle_waves_v5 as PW - - import PiCLES: FetchRelations, ParticleTools - using PiCLES.Operators.core_2D: ParticleDefaults, InitParticleInstance, GetGroupVelocity - using PiCLES.Simulations - using PiCLES.Operators.TimeSteppers: time_step!, movie_time_step! - - using PiCLES.ParticleMesh: TwoDGrid, TwoDGridNotes, TwoDGridMesh - using PiCLES.Models.WaveGrowthModels2D - - using Oceananigans.TimeSteppers: Clock, tick! - import Oceananigans: fields - using Oceananigans.Units - import Oceananigans.Utils: prettytime - - using PiCLES.Architectures - using GLMakie - using PiCLES.Plotting.movie: init_movie_2D_box_plot -end -# debugging: -#using ProfileView - -# %% -save_path = "plots/examples/homogenous_box/" -mkpath(save_path) - -@everywhere begin - # % Parameters - U10, V10 = 10.0, 10.0 - dt_ODE_save = 30minutes - DT = 30minutes - # 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) - - - u(x, y, t) = U10 #* sin(t / (6 * 60 * 60 * 2π)) * sin(x / 50e3) * sin(y / 50e3) - v(x, y, t) = V10 #* cos(t / (6 * 60 * 60 * 2π)) * sin(x / 50e3) * sin(y / 50e3) - - winds = (u=u, v=v) - - grid = TwoDGrid(100e3, 51, 100e3, 51) - #grid = TwoDGrid(20e3, 21, 20e3, 21) - mesh = TwoDGridMesh(grid, skip=1); - gn = TwoDGridNotes(grid); - -end - -Revise.retry() - -# define variables based on particle equation - -#ProfileView.@profview -#ProfileView.@profview -@everywhere begin - particle_system = PW.particle_equations(u, v, γ=0.88, q=Const_ID.q, input=true, dissipation=true); - #particle_equations = PW3.particle_equations_vec5(u, v, u, v, γ=0.88, q=Const_ID.q); - - # 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); - - # define setting and standard initial conditions - WindSeamin = FetchRelations.MinimalWindsea(U10, V10, DT); - #WindSeamin = FetchRelations.MinimalWindsea(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"] - - ODE_settings = PW.ODESettings( - Parameters=default_ODE_parameters, - # define mininum energy threshold - log_energy_minimum=lne_local,#log(FetchRelations.Eⱼ(0.1, DT)), - #maximum energy threshold - log_energy_maximum=log(27),#log(17), # correcsponds to Hs about 16 m - 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(lne_local, cg_u_local, cg_v_local, 0.0, 0.0) -end - -Revise.retry() -wave_model = WaveGrowthModels2D.WaveGrowth2D(; grid=grid, - winds=winds, - ODEsys=particle_system, - ODEsets=ODE_settings, # ODE_settings - ODEinit_type="wind_sea", # default_ODE_parameters - periodic_boundary=false, - boundary_type="same", - #minimal_particle=FetchRelations.MinimalParticle(U10, V10, DT), - movie=true) - -typeof(wave_model.State) - - -#global wave_model.ParticleCollection -### build Simulation -#wave_simulation = Simulation(wave_model, Δt=10minutes, stop_time=4hours)#1hours) -wave_simulation = Simulation(wave_model, Δt=10minutes, stop_time=6hour)#1hours) -initialize_simulation!(wave_simulation) - -# %% -using PiCLES.Operators: mapping_2D -using PiCLES.Operators.TimeSteppers: mean_of_state -import PiCLES.Operators.TimeSteppers: time_step! -# run single timestep - - -advance_wrapper(f, state, Fcol, grid, winds, dt, emax, windmin, boundary, defaults) = x -> f(x, state, Fcol, grid, winds, dt, emax, windmin, boundary, defaults) -remesh_wrapper(f, state, winds, time, sets, dt, minpar, minstate, defaults) = x -> f(x, state, winds, time, sets, dt, minpar, minstate, defaults) - -######## things to try: -# - ParticleCollection should be a global variable -# - try threads instead of processes for remeshing -# - try @async @sync for remeshing -- doesn't do much because task is still done by the same core. - - -""" -time_step!(model, Δt; callbacks=nothing) - -advances model by 1 time step: -1st) the model.ParticleCollection is advanced and then -2nd) the model.State is updated. -clock is ticked by Δt - -callbacks are not implimented yet - -""" -function time_step!(model::Abstract2DModel, Δt::Float64, workers::WorkerPool ; 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 - - @info "advance" - # @time @threads for a_particle in model.ParticleCollection - # #@info a_particle.position_ij - # mapping_2D.advance!(a_particle, model.State, FailedCollection, - # model.grid, model.winds, Δt, - # model.ODEsettings.log_energy_maximum, - # model.ODEsettings.wind_min_squared, - # model.periodic_boundary, - # model.ODEdefaults) - # end - - advance_wrap = advance_wrapper(mapping_2D.advance!, - model.State, FailedCollection, - model.grid, model.winds, Δt, - model.ODEsettings.log_energy_maximum, - model.ODEsettings.wind_min_squared, - model.periodic_boundary, - model.ODEdefaults) - - @time pmap(advance_wrap, workers, model.ParticleCollection) - - - print("mean energy after advance ", mean_of_state(model), "\n") - - if debug - @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 - - @info "re-mesh" - @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 - #model.ParticleCollection = fetch(model.ParticleCollection) - - # remesh_wrap = remesh_wrapper(mapping_2D.remesh!, - # model.State, - # model.winds, model.clock.time, - # model.ODEsettings, Δt, - # model.minimal_particle, - # model.minimal_state, - # model.ODEdefaults) - # pmap(remesh_wrap, workers, model.ParticleCollection); - - 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") - - tick!(model.clock, Δt) -end - -# external definitions -advance_workers = WorkerPool([2, 3, 4, 5]) -workers2 = WorkerPool([2, 3]) -workers2b = WorkerPool([1, 2]) -workers1 = WorkerPool([1]) - -# %% -for i in 1:3 - - wave_simulation.model.State .= 0.0 - - time_step!(wave_simulation.model, wave_simulation.Δt, workers1, debug=false) - #time_step!(wave_simulation.model, wave_simulation.Δt, debug=false) - istate = wave_simulation.model.State; - - p1 = plt.heatmap(gn.x / 1e3, gn.y / 1e3, istate[:, :, 1]) - display(p1) - sleep(0.5) -end -# %% -@fetchfrom 2 InteractiveUtils.varinfo() -using CUDA -@sync_gpu -# % - -@time time_step!(wave_simulation.model, wave_simulation.Δt, advance_workers, debug=false) - -@time time_step!(wave_simulation.model, wave_simulation.Δt, workers1, debug=false) - - -#wave_simulation.model.ParticleCollection -#fetch(wave_simulation.model.ParticleCollection) - -# code to make parallel map - - - -#init_state_store!(wave_simulation, save_path) -#wave_simulation.model.MovieState = wave_simulation.model.State - -@everywhere run!(wave_simulation, cash_store=true, debug=true) -#reset_simulation!(wave_simulation) -# run simulation -#ProfileView.@profview run!(wave_simulation, cash_store=true, debug=true) - - -istate = wave_simulation.store.store[end]; -p1 = plt.heatmap(gn.x / 1e3, gn.y / 1e3, istate[:, :, 1]) - - - - - -# %% Movie, just for fun ... -Revise.retry() -wave_model = WaveGrowthModels2D.WaveGrowth2D(; grid=grid, - winds=winds, - ODEsys=particle_system, - ODEsets=ODE_settings, # ODE_settings - ODEinit_type="wind_sea", # default_ODE_parameters - periodic_boundary=false, - boundary_type="same", - #minimal_particle=FetchRelations.MinimalParticle(U10, V10, DT), - movie=true) - -wave_simulation = Simulation(wave_model, Δt=5minutes, stop_time=4hours)#1hours) -initialize_simulation!(wave_simulation) - -#reset_simulation!(wave_simulation, particle_initials=copy(wave_model.ODEdefaults)) -fig, n = init_movie_2D_box_plot(wave_simulation, name_string="Rotating Stationary Winds") - -#wave_simulation.stop_time += 1hour -N = 10 -record(fig, save_path * "homogenous_box_nonper.gif", 1:N, framerate=10) do i - @info "Plotting frame $i of $N..." - #@time for _ = 1:10 - #run!(wave_simulation, store=false) - @info wave_simulation.model.clock - movie_time_step!(wave_simulation.model, wave_simulation.Δt) - #coupled_time_step!(ocean_simulation, ice_simulation) - #end - n[] = 1 -end diff --git a/src/Operators/core_1D_old.jl b/src/Operators/core_1D_old.jl deleted file mode 100644 index afb482f..0000000 --- a/src/Operators/core_1D_old.jl +++ /dev/null @@ -1,431 +0,0 @@ -module core_1D_old - -using SharedArrays -#using DifferentialEquations - -using DocStringExtensions - -using ParticleMesh - -using particle_waves_v3: init_vars_1D -t, x, c̄_x, lne, r_g, C_α, g, C_e = init_vars_1D() - -using ParticleInCell - -using DifferentialEquations -using ParticleMesh: OneDGrid, OneDGridNotes -# using ModelingToolkit: ODESystem - -using FetchRelations - -# initialize: -export init_z0_to_State - -# Callbacks -export wrap_pos!, periodic_BD_single_PI!, show_pos!, periodic_condition_x - -# Particle Node interaction -export GetParticleEnergyMomentum, GetVariablesAtVertex, Get_u_FromShared - -export ParticleInstance1D, MarkedParticleInstance, ParticleDefaults - -export InitParticleValues, check_boundary_point - -using custom_structures: ParticleInstance1D, MarkedParticleInstance - -# # ParticleInstance is the Stucture that carries each particle. -# mutable struct ParticleInstance -# position_ij::Int -# position_xy::Float64 -# ODEIntegrator::OrdinaryDiffEq.ODEIntegrator -# boundary::Bool -# end - -# # Debugging ParticleInstance -# mutable struct MarkedParticleInstance -# Particle::ParticleInstance -# time::Float64 -# state::Vector{Any} -# errorReturnCode -# end - -#Base.copy(s::ParticleInstance) = ParticleInstance(s.position_ij, s.position_xy, s.ODEIntegrator, s.boundary) - - -""" -ParticleDefaults -Structure holds default particles -# Fields -$(DocStringExtensions.FIELDS) -""" -struct ParticleDefaults - "log energy" - lne::Float64 - "horizontal velocity" - c̄_x::Float64 - "x Position" - x::Float64 -end - -""" -ParticleDefaults -Structure holds default particles -# Fields -$(DocStringExtensions.FIELDS) -""" -struct ParticleDefaults_old - "x Position" - x::Float64 - "horizontal velocity" - c̄_x::Float64 - "log energy" - lne::Float64 -end - - - -Base.copy(s::ParticleDefaults) = Dict(x => s.x, c̄_x => s.c̄_x, lne => s.lne) - - - -#### initialize ##### -""" sets node state values to S at tuple position ij to zi """ -function init_z0_to_State!(S::SharedMatrix, ij::Int, zi::Vector{Float64}) - S[ij[1], :] = zi - nothing -end - -""" sets particle state values to S. position is taking from particle """ -function set_u_to_shared!(S::SharedMatrix, PI::ParticleInstance1D) - S[PI.position_ij[1], PI.position_ij[2], :] = PI.ODEIntegrator.u - nothing -end - -""" takes node state values and pushes them to particle. position is taken from particle """ -function set_u_to_particle!(S::SharedMatrix, PI::ParticleInstance1D) - set_u!(PI.ODEIntegrator, Get_u_FromShared(S, PI)) - u_modified!(PI.ODEIntegrator, true) - nothing -end - -##### Callbacks ###### - -""" -wrap_pos!(xx::Float64, L::Float64) -makes periodic boundary conditions. If position xx exceeds 0 or L it is wrapped around -returns new position and if either or not the wrapping appends (bool) -""" -function wrap_pos!(xx::Float64, L::Float64) - wrap_flag = false - if xx > L - xx = xx - L - wrap_flag = true - elseif xx < 0 - xx = xx + L - wrap_flag = true - end - - return xx, wrap_flag -end - -""" -Checks if particle's PI position PI.ODEIntegrator.u[1,2] exceeeds the domain limits Lx, Ly - Lx, LY are global variables - If the PI's postions eceeeds the boundary they warpped around and the ODEIntegrator state is modified. -""" -function periodic_BD_single_PI!(PI::ParticleInstance1D, Lx::Float64) - - #@printf "PI periodic condition called" - ui = copy(PI.ODEIntegrator.u) - ui[3], wrap_pos_PI1 = wrap_pos!(ui[3], Lx) - #ui[2], wrap_pos_PI2 = wrap_pos!(ui[2], Ly) - - if wrap_pos_PI1 #|| wrap_pos_PI1 - #@show wrap_pos_PI1 , wrap_pos_PI2 - #@printf "wrap pos PI" - #@show PI.ODEIntegrator.u[1] - ui[1] - #@show PI.ODEIntegrator.u[2] - ui[2] - set_u!(PI.ODEIntegrator, ui) - u_modified!(PI.ODEIntegrator, true) - end - nothing #return PI -end - -# """ -# Checks if -# """ -# @everywhere function periodic_BD_single!(S::SharedMatrix) -# S.u[1], wrap_pos_PI = wrap_pos!(S.u[1], Lx) -# S.u[2], wrap_pos_PI = wrap_pos!(S.u[2], Ly) -# nothing -# end - -""" -return the mean position -""" -function show_pos!(PI) - @show "show current position" PI.position_ij, PI.position_xy / dx, PI.ODEIntegrator.u[3] / dx, PI.ODEIntegrator.u[2] -end - - -""" (debubugging function) returns the domains normalized position of the partivle """ -function show_pos!(integrator, G) - @show (integrator.ODEIntegrator.u[3] .- G.xmin) / grid1d.dx -end - -# the following function are another way to have wrapping boundary conditions. -function periodic_condition_x(u, t, integrator) - u[3] > 0 -end - -# function loop_x!(integrator) -# integrator.u[1] = integrator.u[1] - Lx -# @printf "wrap pos PI x" -# end -######### Particle --> Node ########## - -#### 1D code #### -""" -GetParticleEnergyMomentum(PI) - -""" -function GetParticleEnergyMomentum(PI) - #ui_x, ui_c̄_x, ui_lne = PI.ODEIntegrator.u - ui_lne, ui_c̄_x, ui_x = PI.ODEIntegrator.u - - ui_e = exp(ui_lne) - #c_speed = particle_waves_v3.speed(ui_c̄_x, ui_c̄_y) - m_x = ui_e / ui_c̄_x / 2 - #m_y = ui_c̄_y * ui_e / c_speed^2 / 2 - - return ui_e, m_x, 0 -end - -function GetParticleEnergyMomentum(PI::ParticleInstance1D) - return GetParticleEnergyMomentum(PI.ODEIntegrator.u) -end - -function GetParticleEnergyMomentum(zi::Dict) - return GetParticleEnergyMomentum([zi[lne], zi[c̄_x], zi[x]]) -end - -function GetParticleEnergyMomentum(zi::ParticleDefaults) - return GetParticleEnergyMomentum([zi.lne, zi.c̄_x, zi.x]) -end - - -""" -GetParticleEnergyMomentum(PI) - -""" -function GetParticleEnergyMomentum(z0::Vector{Float64}) - - ui_lne, ui_c̄_x, ui_x = z0 - ui_e = exp(ui_lne) - #c_speed = ui_c̄_x - m_x = ui_e / ui_c̄_x / 2 - #m_y = ui_c̄_y * ui_e / c_speed^2 / 2 - - return [ui_e, m_x, 0] -end - -# #### 2D code #### -# """ -# GetParticleEnergyMomentum(PI) - -# """ -# function GetParticleEnergyMomentum(PI) -# ui_x, ui_y, ui_c̄_x, ui_c̄_y, ui_lne, ui_Δn, ui_Δφ_p = PI.ODEIntegrator.u - -# ui_e = exp(ui_lne) -# c_speed = particle_waves_v3.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 -# end - -# function GetParticleEnergyMomentum(PI::ParticleInstance1D) -# 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[Δn], zi[Δφ_p]]) -# # end - -# function GetParticleEnergyMomentum(zi::Dict{Num, Float64}) -# return GetParticleEnergyMomentum([zi[lne], zi[c̄_x], zi[x]]) -# end - - -# """ -# GetParticleEnergyMomentum(PI) - -# """ -# function GetParticleEnergyMomentum(z0::Vector{Float64}) - -# ui_x, ui_c̄_x, ui_lne = z0 -# #ui_x, ui_y, ui_c̄_x, ui_c̄_y, ui_lne, ui_Δn, ui_Δφ_p = z0 -# ui_e = exp(ui_lne) -# c_speed = ui_c̄_x -# 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] -# end - - - - - - -######### node--> Particle ########## -""" -GetVariablesAtVertex(i_State::Vector{Float64}, x::Float64) -i_State: [e, m_x, m_y] state vector at node -x: coordinates of the vertex -""" -function GetVariablesAtVertex(i_State::Vector{Float64}, x::Float64) - e, m_x, m_y = i_State - #m_amp = particle_waves_v3.speed(m_x, m_y) - c_x = e / 2 / m_x - #c_y = m_y * e / (2 * m_amp^2) - - return [log(e), c_x, x] -end - - -# """ -# GetVariablesAtVertex(i_State::Vector{Float64}, x::Float64, y::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 = i_State -# m_amp = particle_waves_v3.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] -#end - - -""" returns node state from shared array, given paritcle index """ -Get_u_FromShared(PI::ParticleInstance1D, S::SharedMatrix,) = S[PI.position_ij[1], :] - - -###### seed particles ##### - -""" -InitParticleValues(defaults:: Dict{Num, Float64}, i::Int64, gridnote::OneDGridNotes, u, 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 - i index of the grid point - gridnote grid to determine the position in the grid - u interp. fucntion with wind values - DT time step of model, used to determine fetch laws -""" -function InitParticleValues( - defaults::Dict, - i::Int64, - gridnote::OneDGridNotes, - u, DT) - - # initalize state based on state vector - defaults[x] = gridnote.x[i] - # take in local wind velocities - u_init = u(defaults[x], 0) - - # seed particle given fetch relations - defaults[c̄_x] = FetchRelations.c_g_U_tau(abs(u_init), DT) - defaults[lne] = log(FetchRelations.Eⱼ(abs(u_init), DT)) - #@show defaults - return defaults -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 - - -""" -Init(model, z_initials, pars, ij ; 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 - chSet (optional) is the set of callbacks the ODE can have -""" -function InitParticleInstance(model, z_initials, ODE_settings, ij, boundary_flag; cbSets=Nothing) - - # 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, - maxiters=ODE_settings.maxiters, - reltol=ODE_settings.reltol, - callback=ODE_settings.callbacks, - save_everystep=ODE_settings.save_everystep) - return ParticleInstance1D(ij, z_initials[x], integrator, boundary_flag) -end - - - - -""" -SeedParticle!(ParticleCollection ::Vector{Any}, State::SharedMatrix, i::Int64, - particle_system::ODESystem, particle_defaults::Dict{Num, Float64}, ODE_defaults::Dict{Num, Float64}, - GridNotes, winds, DT:: Float64, Nx:: Int, boundary::Vector{Int}, periodic_boundary::Bool) - -Seed Pickles to ParticleColletion and State -""" -function SeedParticle!( - ParticleCollection::Vector{Any}, - State::SharedMatrix, - i::Int64, - particle_system::ODESystem, - particle_defaults, - ODE_settings, #particle_waves_v3.ODESettings type - GridNotes, # ad type of grid note - winds, # interp winds - DT::Float64, - Nx::Int, - boundary::Vector{Int}, - periodic_boundary::Bool) - - # define initial condition - z_i = InitParticleValues(particle_defaults, i, GridNotes, winds, DT) - # check if point is boundary point - boundary_point = check_boundary_point(i, boundary, periodic_boundary) - - # add initial state to State vector - init_z0_to_State!(State, i, GetParticleEnergyMomentum(z_i)) - - # Push Inital condition to collection - push!(ParticleCollection, - InitParticleInstance( - particle_system, - z_i, - ODE_settings, - i, - boundary_point)) - nothing -end - -# end of module -end diff --git a/src/old_structure/run_1D.jl b/src/old_structure/run_1D.jl deleted file mode 100644 index 8088f1c..0000000 --- a/src/old_structure/run_1D.jl +++ /dev/null @@ -1,528 +0,0 @@ - -#@show localARGS - -using Plots -using BenchmarkTools - -using Interpolations -using IfElse -using ModelingToolkit, DifferentialEquations, Statistics - - -# Particle Model modules -push!(LOAD_PATH, joinpath(pwd(), "code/")) -using ParticleMesh: OneDGrid, OneDGridNotes -using SharedArrays - -#include("./particle_waves_v1.jl") -using Revise - -import ParticleInCell -import FetchRelations -#includet("ParticleInCell.jl") -#includet("FetchRelations.jl") - -import particle_waves_v3 -using Printf - -using HDF5 -using JLD2 - - -""" -This load a parameter file, executes a 1D run, saves the data and run statistics. - -""" - -push!(LOAD_PATH, joinpath(pwd(), "code/")) -push!(LOAD_PATH, joinpath(pwd(), "code/Core")) -using core_1D_old: init_z0_to_State!, wrap_pos!, periodic_BD_single_PI!, show_pos!, periodic_condition_x -using core_1D_old: GetParticleEnergyMomentum, GetVariablesAtVertex, Get_u_FromShared -#using core_1D: InitParticleValues, check_boundary_point -using core_1D_old: InitParticleValues, SeedParticle! - - -using custom_structures: ParticleInstance1D, MarkedParticleInstance, AbstractParticleInstance - -using InputOutput: Argsettings, parse_args - -using ParticleTools: ParticleToDataframe - -includet("Core/mapping_1D.jl") -#using mapping_1D -using Debugging - - -Revise.retry() - -# Default values -save_path_base = "data/" -plot_path_base = "plots/static/" -parset = "1D_static/" - -T = 24 * 3 -Lx = 50 - -DT = 30 #60*30 # remeshing time stepping -dt_ODE_save = 10 # 3 min - -Nx = 40 - -### Boundary Conditions -periodic_boundary = false - -# parametric wind forcing -U10 = 20 - -# model parameters -r_g0 = 0.9 -c_beta = 4#e-2 # 4e-2 Growth rate constant # default value from Kudr. -gamma = γ = 0.88 # dissipation wind energy input - -@printf "Load passed arguments\n" -# arg_test = [ "--ID", "test1", -# "--Lx", "20", -# "--T", "24", -# "--DT", "30"] - -#localARGS = ["--ID", "Nx30_DT30_U15", "--Nx", "30", "--DT", "30", "--U10", "15", "--parset", "1D_static"]ß -#localARGS = ["--ID", "Nx50_DT2_U1", "--Nx", "50", "--c_beta", "2.000000", "--U10", "1", "--parset", "U10-c_beta", "--periodic"] -#localARGS = ["--ID", "Nx40_cbeta6.00_gamma1.32", "--c_beta", "6.00", "--Nx", "40", "--gamma", "1.32", "--parset", "gamma-c_beta"] -#localARGS = ["--ID", "rg0.85_cbeta6.00_gamma1.32", "--c_beta", "6.00", "--gamma", "1.32", "--rg", "0.85", "--parset", "gamma-c_beta-rg0.85"] -#localARGS = ["--ID", "rg1.25_gamma1.30_Nx40", "--Nx", "40", "--gamma", "1.32", "--rg", "1.25", "--parset", "gamma-rg0"] -#localARGS = ["--ID", "U25_DT60_Nx40", "--U10", "25", "--Nx", "40", "--DT", "60", "--parset", "U10-DT-periodic", "--periodic"] - -#localARGS = ["--ID", "U20_DT3_Nx40", "--U10", "20", "--Nx", "40", "--DT", "3", "--parset", "U10-DT"] -#arg_case = ["--ID", "U25_DT60_Nx40", "--U10", "25", "--Nx", "40", "--DT", "60", "--parset", "U10-DT-periodic", "--periodic"] -#localARGS = ["--ID", "U13_DT60_Nx40", "--U10", "13", "--Nx", "40", "--DT", "60", "--parset", "U10-DT"] -localARGS = ["--ID", "gamma1.32_rg1.20_Nx50", "--gamma", "1.32", "--r_g0", "1.25", "--Nx", "50", "--parset", "gamma-rg0_2"] -@show localARGS -passed_argument = parse_args(localARGS, Argsettings) - -# change here if more argument shuold be allowed to pass -#@unpack ID, parset, periodic, U10, Nx, DT = passed_argument #"U10-DT" -@unpack ID, parset, periodic, gamma, Nx, r_g0 = passed_argument #"gamma-rg_2" -#@unpack ID, parset, periodic, gamma, Nx, c_beta = passed_argument #"gamma-c_beta" - -periodic_boundary = periodic -c_β = c_beta * 1e-2 -C_e0 = (2.35 / r_g0) * 2e-3 * c_β -γ = Float64(gamma) -#r_g0 = rg -# -# if ~isnothing(ID) -# #ID = "Nx"*@sprintf("%i",Nx)*"_dt"*@sprintf("%i",DT)*"_U"*@sprintf("%i",U10) -# ID = "Nx"*@sprintf("%i",Nx)*"_cbeta"*@sprintf("%.1f",c_beta)*"_U"*@sprintf("%i",U10) -# end - -# rescale parameters for the right units. -T = T * 60 * 60 # seconds -Lx = Lx * 10e3 # km -DT = Float64(DT * 60) # seconds - - -# create ID and save name -#save_path = joinpath( "data/1D_static/", parsed_args["ID"] ) -save_path = save_path_base * parset * "/" * ID * "/" -plot_path = plot_path_base * parset * "/" - -# %% -@printf "Init Forcing Field\n" -# create wind test fucntion - -@register_symbolic u(x, t) -@register_symbolic u_x(x, t) - -#u_func(x, y, t) = y * 0 .+ 3 * exp(- ( ( x-25e3 + 20e3* t/T )./10e3).^2) #.+ 3 #.+ 3 * sin.(x *π/Ly/0.5 .+ y *π/Ly/0.5) -#u_func(x, t) = 3 * exp(- ( ( x-25e3 )./10e3).^2) .+ t *0 #.+ 3 #.+ 3 * sin.(x *π/Ly/0.5 -u_func(x, t) = x .* 0 + U10 + t * 0 #3 * exp(- ( ( x-25e3 )./10e3).^2) .+ t *0 #.+ 3 #.+ 3 * sin.(x *π/Ly/0.5 -#u_func(x, y, t) = y *0 .- x .*2/50e3 .+ 3 #.+ 3 * sin.(x *π/Ly/0.5 .+ y *π/Ly/0.5) -#u_func(x, y, t) = y *0 .- x .*0 .+ 3 #.+ 3 * sin.(x *π/Ly/0.5 .+ y *π/Ly/0.5) -#u_func(x, t) = x *0 .+ IfElse.ifelse.( x .< Lx*2.0/3.0, 3, 0.1) - -# u_func(x, y, t) = y *0 .+ IfElse.ifelse.( x .< Lx/2.0, -# 2 .+ 1, *sin.(x *π/Lx/2), -# 0 .* x + 0.1 -# ) - -# %% Define Grid, boundaries, and forcing field grid: - -### Define grid -grid1d = OneDGrid(1e3, Lx - 1e3, Nx) -grid1dnotes = OneDGridNotes(grid1d) -dx = grid1d.dx - -### storing matrix -Nstate = 3 -Nparticle = grid1d.Nx -State = SharedMatrix{Float64}(grid1d.Nx, Nstate) - -if ~periodic_boundary # if false, define boundary points here: - boundary = [1, Nx] -end - -# provide winds in the right format -xi = range(0, Lx, step=5e3) # note ': this is a row vector -ti = range(0, T, step=60 * 60) - -# % define wind forcing globally as interp functions -u_func_gridded = [u_func(xii, tii) for xii in xi, tii in ti] -u_grid = linear_interpolation((xi, ti), u_func_gridded, extrapolation_bc=Periodic()) -#u_grid = CubicSplineInterpolation( (xi, ti), u_func_gridded; bc=Line(OnGrid()), extrapolation_bc=Flat()) -#u_grid = interpolate((xi, ti), u_func_gridded, Gridded(Linear())) - - -u(x, t) = u_grid(x, t) -# x gradient only -u_x(x, t) = Interpolations.gradient(u_grid, x)[1] - -# %% Load Particle equations and derive ODE system - -particle_equations0 = particle_waves_v3.particle_equations(u, u_x, γ=γ) -@named particle_system = ODESystem(particle_equations0) - -# define variables based on particle equations -t, x, c̄_x, lne, r_g, C_α, g, C_e = particle_waves_v3.init_vars_1D() - -# % define storing stucture and populate inital conditions - -default_ODE_parameters = Dict( - r_g => (1 / r_g0), - C_α => -1.41, - g => 9.81, - C_e => C_e0, -) - -### limiters: -# define mininum energy threshold -e_min_log = log(FetchRelations.Eⱼ(0.1, DT)) -#e_min_log = log(FetchRelations.Eⱼ( 0.1 , 60.0 )) -e_max_log = log(17) # correcsponds to Hs about 16 m - -# Default values for particle -particle_defaults = Dict(x => 0.0, c̄_x => 1e-2, lne => e_min_log) - -# %%%% Define callbacks - -# terminate = PeriodicCallback(terminate_check!, 60*12 ) -show_mean = PeriodicCallback(show_pos!, DT / 2) -periodic = PeriodicCallback(periodic_BD_single_PI!, DT / 2) -cb = CallbackSet(periodic, show_mean)#,cb_terminate) - -# testing -#GetParticleEnergyMomentum(z0) -#ui_x, ui_c̄_x, ui_lne = PI.ODEIntegrator.u -#ie, imx, imy = GetParticleEnergyMomentum(PI.ODEIntegrator.u) - -# testing -#GetVariablesAtVertex( [ie, imx, imy], ui_x) - PI.ODEIntegrator.u - -# %% TEST 1 ODE solution -# -# z_i = copy(z0) -# z_i[x] = 20e4#grid1dnotes.x[20] -# -# ## seed particle given fetch relations -# u_init = u(z_i[x], 0) -# z_i[c̄_x] = FetchRelations.c_g_U_tau( u_init , DT ) -# z_i[lne] = log(FetchRelations.Eⱼ( u_init , DT )) -# -# #z_i[lne] = log(0.01) -# -# params_i = copy(params0) -# problem = ODEProblem(particle_system, z_i, (0.0, T) , params_i) -# problem -# -# sol = solve(problem, saveat=60*5) -# -# solver_method = AutoTsit5(Rosenbrock23()) -# #solver_method = Rosenbrock23() -# integrator = init(problem, solver_method , saveat =dt_ODE_save, adaptive =true, maxiters=1e3)#, reltol=1e-1) -# integrator -# -# step!(integrator, DT , true) -# -# # using PyPlot - -# %% - -# plt.figure() -# plt.subplot(3,1,1) -# plt.plot(sol[t], sol[c̄_x]) -# -# plt.subplot(3,1,2) -# plt.plot(sol[t], sol[x].- sol[x][1]) -# -# plt.subplot(3,1,3) -# plt.plot(sol[t], exp.(sol[lne]) ) -# -# display(plt.gcf()) -# %% seeding particles with initial states -@printf "seed particles ... \n" - -# """ -# InitParticleInstance(model, z_initials, pars, ij ; 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 -# chSet (optional) is the set of callbacks the ODE can have -# """ -# function InitParticleInstance(model, z_initials, pars, ij, boundary_flag ; cbSets=nothing) - -# # create ODEProblem -# problem = ODEProblem(model, z_initials, (0.0, T) , pars) - -# # inialize problem -# #solver_method = Rosenbrock23() -# #solver_method = AutoVern7(Rodas4()) -# solver_method = AutoTsit5(Rosenbrock23()) -# #solver_method =Tsit5() - -# # works best with abstol = 1e-4,reltol=1e-3,maxiters=1e4, -# integrator = init( -# problem, -# solver_method, -# saveat =dt_ODE_save, -# abstol = 1e-4, -# adaptive =true, -# maxiters=1e4, -# reltol=1e-3) -# #abstol=1e-4)#, -# # callbacks= cbSets, -# #save_everystep=false -# return ParticleInstance1D( ij , z_initials[x], integrator, boundary_flag ) -# end - - -ODE_settings = particle_waves_v3.ODESettings( - Parameters=particle_defaults, - # define mininum energy threshold - log_energy_minimum=e_min_log,#log(FetchRelations.Eⱼ(0.1, DT)), - #maximum energy threshold - log_energy_maximum=log(17), # correcsponds to Hs about 16 m - saving_step=dt_ODE_save, - timestep=DT, - total_time=T) - - -Revise.retry() -#seed particles -ParticleCollection = [] -for i in range(1, length=Nx) - SeedParticle!(ParticleCollection, State, i, - particle_system, default_ODE_parameters, ODE_settings, - grid1dnotes, u, DT, Nx, boundary, periodic_boundary) -end - -# test -PI = ParticleCollection[1] -PI.ODEIntegrator -#sol = solve(PI.ODEIntegrator) -PI.ODEIntegrator.sol.u -step!(PI.ODEIntegrator, DT, true) -PI.ODEIntegrator.sol.u - -#plot( PI.ODEIntegrator.sol.u[1], PI.ODEIntegrator.sol.u) -using ParticleTools -PID = ParticleTools.ParticleToDataframe(PI) - -# %% -# plit each row in PID and a figure -p1 = plot(PID[:, 2], exp.(PID[:, 4]), marker=2, title="e", xlabel="position", ylabel="e") #|> display -p2 = plot(PID[:, 2], PID[:, 3], marker=2, title="cg", ylabel="cg") #|> display -p3 = plot(PID[:, 1] / 60 / 60, PID[:, 4], marker=3, title="log e", xlabel="time", ylabel="postition") #|> display - -plot(p1, p2, p3, layout=(3, 1), legend=false, size=(600, 1200)) - - -# %% - -# # #seed particles -# ParticleCollection=[] -# local z_i -# local u_init -# for i in range(1,length = Nx) -# # there is SeedParticle! in the core module as alternative wrapper -# -# # define initial condition -# z_i = InitParticleValues(copy(z0), grid1dnotes, u, DT ) -# # check if point is boundary point -# boundary_point = check_boundary_point(i, boundary, periodic_boundary) -# -# -# #@show z_i, boundary_point -# -# # add initial state to State vector -# init_z0_to_State!(State, i, GetParticleEnergyMomentum(z_i) ) -# -# # Push Inital condition to collection -# push!( ParticleCollection, -# InitParticleInstance( -# particle_system, -# z_i , -# copy(params0), -# i , -# boundary_point)) -# -# end - - - -# %% single core -@printf "Start itteration ... \n" -# define state collector -State_collect = [] # storing array -State[:, :, :] .= 0 # reset current state # [ Energy, x momentum, not used] -push!(State_collect, copy(State)) # push initial state to storage - -# collecting failed particles -FailedCollection = Vector{MarkedParticleInstance}([]) - -# main time stepping: -t_range = range(0.0, DT * T / DT, step=DT) - -elapsed_time = @elapsed for t in t_range[2:end] - - State[:, :] .= 0 - #u_sum_m1 = mapping_1D.ShowTotalEnergyChange(ParticleCollection, u_sum_m1) - - for a_particle in ParticleCollection - #@show a_particle.position_ij - mapping_1D.advance!(a_particle, State, FailedCollection, grid1d, u, DT, e_min_log, e_max_log, periodic_boundary) - end - - #u_sum_m1 = mapping_1D.ShowTotalEnergyChange(ParticleCollection, u_sum_m1) - - #@show (State[:,:,1] .- grid1d.xmin)/grid1d.dx - #@printf "re-mesh" - for a_particle in ParticleCollection - mapping_1D.remesh!(a_particle, State, u, t, e_min_log, DT) - end - - #@printf "push" - push!(State_collect, copy(State)) - - #u_sum_m1 = mapping_1D.ShowTotalEnergyChange(ParticleCollection, u_sum_m1) - #@show (State[:,:,1] .- grid1d.xmin)/grid1d.dx -end - -@printf "... finished\n" - - -# %% plotting errors if applicaple - - -if ~isempty(FailedCollection) - @printf "plot failed particles\n" - using ParticleTools - mkpath(plot_path) - - ParticleTools.PlotFailedParticles(FailedCollection[1:min(4, size(FailedCollection)[1])], ID, DT, dx) - savefig(joinpath(plot_path, "failed_ov_" * ID * ".png")) - - ParticleTools.plot_cg_version1(FailedCollection) - savefig(joinpath(plot_path, "cg_failed_" * ID * ".png")) - - ParticleTools.plot_cg_version2(State_collect) - savefig(joinpath(plot_path, "cg_failed2_" * ID * ".png")) - -end - -# %% -PI = ParticleCollection[10] -ParticleInCell.compute_weights_and_index(grid1d, PI.ODEIntegrator.u[3]) - -# ParticleCollection[10].ODEIntegrator.u -# -# -energy = State_collect[end-1][:, 1] #GP.sel(state= 0 ).T # energy -m_x = State_collect[end-1][:, 2] #GP.sel(state= 1 ).T # energy -cg = energy ./ m_x ./ 2 # c_g -# @printf "\n max wave age %s" U10 / (2 * maximum( cg[.!isnan.(cg)] ) ) -# - - -# %% Checks - -@printf "\n size of collected states %s" size(State_collect) -@printf "\n size of each state %s" size(State_collect[1]) -@printf "\n length of timerange %s" size(t_range) - -# %% Saving Fields -@printf "save data \n" - -mkpath(save_path) -@printf "created model run folder %s \n" save_path - -rm(joinpath(save_path, "state.h5"), force=true) -file = h5open(joinpath(save_path, "state.h5"), "w") - -store_waves = create_group(file, "waves") -store_waves_data = create_dataset( - store_waves, "data", Float64, - (length(State_collect), (grid1dnotes.Nx), (size(State_collect[1])[2])) -)#, chunk=(2,)) - -for i in 1:length(State_collect) - store_waves_data[i, :, :] = State_collect[i] -end - -write_attribute(store_waves, "dims", ["time", "x", "state"]) -store_waves["time"] = Array(t_range) -store_waves["x"] = grid1dnotes.x -store_waves["var_names"] = ["e", "m_x", "m_y"] - - -store_winds = create_group(file, "wind") -store_winds_u = create_dataset(store_winds, "u", Float64, ((length(xi)), (length(ti))))#, chunk=(2,)) -store_winds_v = create_dataset(store_winds, "v", Float64, ((length(xi)), (length(ti))))#, chunk=(2,)) - -store_winds_u[:, :] = u_func_gridded -#store_winds_v[:,:,:] = v_func_gridded - -write_attribute(store_winds, "dims", ["x", "time"]) -store_winds["time"] = Array(ti) -store_winds["x"] = Array(xi) -#store_winds["y"] = Array(yi) - -close(file) - - -# %% Save Particles -# @printf "save particle \n" -# save_object(joinpath(save_path , "particles.jld2"), ParticleCollection) - -#PC2 = load_object(joinpath(save_path , "particles.jld2")) - -# %% -@printf "save parameters \n" -using JSON3 -params = Dict("ID" => ID, - "parset" => parset, - "elapsed_time" => elapsed_time, - "grid" => Dict( - "Nx" => Nx, - "dx" => dx, - "Nstate" => Nstate, - "T" => T, - "dt" => DT - ), - "ODE_sampling" => Dict( - "dt_ODE_save" => dt_ODE_save - ), - "forcing" => Dict( - "U10" => U10, - "max_alpha" => U10 / (2 * maximum([maximum(cg[.!isnan.(cg)]), 0.01])) - ) -) - -open(joinpath(save_path, "parameters.json"), "w") do io - JSON3.pretty(io, params) -end - -@printf "done.\n" - -# %% diff --git a/src/old_structure/run_1D_manual.jl b/src/old_structure/run_1D_manual.jl deleted file mode 100644 index 7220e2a..0000000 --- a/src/old_structure/run_1D_manual.jl +++ /dev/null @@ -1,594 +0,0 @@ - -using Plots -using BenchmarkTools - -# CPU initialization -using Distributed -#addprocs(2) -#advance_workers = WorkerPool([2, 3]) - -@everywhere using Interpolations -@everywhere using IfElse -@everywhere using ModelingToolkit, DifferentialEquations, Statistics - - -# Particle Model modules -@everywhere push!(LOAD_PATH, joinpath(pwd(), "code/") ) -@everywhere using ParticleMesh: OneDGrid, OneDGridNotes -@everywhere using SharedArrays - -#include("./particle_waves_v1.jl") -using Revise - -#@everywhere import ParticleInCell -includet("ParticleInCell.jl") -includet("FetchRelations.jl") - -@everywhere import particle_waves_v3 -@everywhere using Printf - -# modules for saveing: -using HDF5 -using JLD2 - -""" - -""" -push!(LOAD_PATH, joinpath(pwd(), "code/")) -push!(LOAD_PATH, joinpath(pwd(), "code/Core")) -using core_1D: init_z0_to_State!, wrap_pos!, periodic_BD_single_PI!, show_pos!, periodic_condition_x -using core_1D: GetParticleEnergyMomentum, GetVariablesAtVertex -using core_1D: Get_u_FromShared, ParticleDefaults -using custom_structures: AbstractParticleInstance, ParticleInstance1D, MarkedParticleInstance - -# %% - -@printf "Init Forcing Field\n" - -save_path_base= "data/1D_static/" - -@everywhere begin - - @register_symbolic u(x, t) - @register_symbolic u_x(x, t) - - # create fake wind data and its gradients - T = 60*60*24 *3 - Lx = 50e3 - xi = range(0,Lx,step=5e3) # note ': this is a row vector - ti = range(0,T ,step=60*60) - - DT = 60 * 20 #60*30 # remeshing time stepping - dt_ODE_save = 60 * 10 # 3 min - - U10 = 2 - periodic_boundary = true - - #u_func(x, y, t) = y * 0 .+ 3 * exp(- ( ( x-25e3 + 20e3* t/T )./10e3).^2) #.+ 3 #.+ 3 * sin.(x *π/Ly/0.5 .+ y *π/Ly/0.5) - - #u_func(x, t) = 3 * exp(- ( ( x-25e3 )./10e3).^2) .+ t *0 #.+ 3 #.+ 3 * sin.(x *π/Ly/0.5 - u_func(x, t) = x.*0+U10 + t *0 #3 * exp(- ( ( x-25e3 )./10e3).^2) .+ t *0 #.+ 3 #.+ 3 * sin.(x *π/Ly/0.5 - - #@everywhere u_func(x, y, t) = y *0 .- x .*2/50e3 .+ 3 #.+ 3 * sin.(x *π/Ly/0.5 .+ y *π/Ly/0.5) - #@everywhere u_func(x, y, t) = y *0 .- x .*0 .+ 3 #.+ 3 * sin.(x *π/Ly/0.5 .+ y *π/Ly/0.5) - - #u_func(x, t) = x *0 .+ IfElse.ifelse.( x .< Lx*2.0/3.0, 3, 0.1) - - # @everywhere u_func(x, y, t) = y *0 .+ IfElse.ifelse.( x .< Lx/2.0, - # 2 .+ 1, *sin.(x *π/Lx/2), - # 0 .* x + 0.1 - # ) - - # % define wind forcing globally as interp functions - nodes = (xi, ti) - - u_func_gridded = [ u_func(xii, tii) for xii in xi, tii in ti] - u_grid = LinearInterpolation( nodes , u_func_gridded ,extrapolation_bc=Periodic()) - #u_grid = CubicSplineInterpolation( nodes, u_func_gridded; bc=Line(OnGrid()), extrapolation_bc=Flat()) - #u_grid = interpolate(nodes, u_func_gridded, Gridded(Linear())) - u(x, t) = u_grid(x, t) - # only x gradient - u_x(x, t) = Interpolations.gradient(u_grid, x )[1] - - particle_equations0 = particle_waves_v3.particle_equations(u, u_x) - - @named particle_system0 = ODESystem(particle_equations0) - - # define variables based on particle equations - t, x, c̄_x, lne, r_g, C_α, g, C_e = particle_waves_v3.init_vars_1D() -end -# %% define storing stucture and populate inital conditions - -@everywhere Nx = 15 -@everywhere grid1d = OneDGrid(1e3, Lx-1e3, Nx) -@everywhere grid1dnotes = OneDGridNotes(grid1d) -@everywhere dx = grid1d.dx -@everywhere Nstate = 3 - - -Nparticle = grid1d.Nx -State = SharedMatrix{Float64}(grid1d.Nx, Nstate) -State.pids # check which cores can see state - -# parameters -params0 = Dict( - r_g => 0.9, - C_α => -1.41 , - g => 9.81 , - C_e => 2.16e-4 - ) - -@everywhere begin - - # define mininum energy threshold - e_0 = log(FetchRelations.Eⱼ( 0.5 , float(DT) )) - z0 = Dict( x => 0.0 , c̄_x=> 1e-2, lne => e_0) - -end -# %%%% Define callbacks - -# terminate = PeriodicCallback(terminate_check!, 60*12 ) -show_mean = PeriodicCallback(show_pos! , DT/2 ) -periodic = PeriodicCallback(periodic_BD_single_PI! , DT/2 ) -cbs = CallbackSet(periodic, show_mean )#,cb_terminate) - -# testing -#GetParticleEnergyMomentum(z0) -#ui_x, ui_c̄_x, ui_lne = PI.ODEIntegrator.u -#ie, imx, imy = GetParticleEnergyMomentum(PI.ODEIntegrator.u) - -# testing -#GetVariablesAtVertex( [ie, imx, imy], ui_x) - PI.ODEIntegrator.u - -# %% TEST 1 ODE solution -# -# z_i = copy(z0) -# z_i[x] = 20e4#grid1dnotes.x[20] -# -# ## seed particle given fetch relations -# u_init = u(z_i[x], 0) -# z_i[c̄_x] = FetchRelations.c_g_U_tau( u_init , float(DT) ) -# z_i[lne] = log(FetchRelations.Eⱼ( u_init , float(DT) )) -# -# #z_i[lne] = log(0.01) -# -# params_i = copy(params0) -# problem = ODEProblem(particle_system0, z_i, (0.0, T) , params_i) -# problem -# -# sol = solve(problem, saveat=60*5) -# -# solver_method = AutoTsit5(Rosenbrock23()) -# #solver_method = Rosenbrock23() -# integrator = init(problem, solver_method , saveat =dt_ODE_save, adaptive =true, maxiters=1e3)#, reltol=1e-1) -# integrator -# -# step!(integrator, DT , true) -# -# # using PyPlot - -# %% - -# plt.figure() -# plt.subplot(3,1,1) -# plt.plot(sol[t], sol[c̄_x]) -# -# plt.subplot(3,1,2) -# plt.plot(sol[t], sol[x].- sol[x][1]) -# -# plt.subplot(3,1,3) -# plt.plot(sol[t], exp.(sol[lne]) ) -# -# display(plt.gcf()) -# %% seeding particles with initial states -@printf "Init Particles\n" - -""" -InitParticleInstance(model, z_initials, pars, ij ; 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 - chSet (optional) is the set of callbacks the ODE can have -""" -function InitParticleInstance(model, z_initials, pars, ij, boundary_flag ; cbSets=nothing) - - # create ODEProblem - problem = ODEProblem(model, z_initials, (0.0, T) , pars) - # inialize problem - #solver_method = Rosenbrock23() - #solver_method = AutoVern7(Rodas4()) - solver_method = AutoTsit5(Rosenbrock23()) - #solver_method =Tsit5() - integrator = init(problem, solver_method , saveat =dt_ODE_save , abstol = 1e-3, adaptive =true, maxiters=1e3)#, reltol=1e-1) - #abstol=1e-4)#, reltol=1e-1) - - # callbacks= cbSets, - #save_everystep=false - return ParticleInstance1D( ij , z_initials[x], integrator, boundary_flag ) -end - -# define boundaries -boundary = [1 , Nx] - -ParticleCollection=[] -local z_i -local u_init -for i in range(1,length = Nx) - - # initalize state based on state vector - z_i = copy(z0) - z_i[x] = grid1dnotes.x[i] - - # take in local wind velocities - u_init = u(z_i[x], 0) - - # seed particle given fetch relations - z_i[c̄_x] = FetchRelations.c_g_U_tau( abs(u_init) , float(DT) ) - z_i[lne] = log(FetchRelations.Eⱼ( abs(u_init) , float(DT) )) - - #@show e_0 > z_i[lne] - #@show (z_i[x]-grid1d.xmin)/grid1d.dx, (z_i[y]-grid1d.ymin)/grid1d.dy - - # push z_i to shared array - init_z0_to_State!(State, i, GetParticleEnergyMomentum(z_i) ) - - # push z_i and params_i to particle system - # particle_system0 is an initilized ODESystem - - if ~periodic_boundary - boundary_points = (i in boundary) - else - boundary_points = false - end - - push!( ParticleCollection, - InitParticleInstance( - particle_system0, - z_i , - copy(params0) , - i , - boundary_points - ) - ) - #@show threadid() -end - - -###### Debugging functions ############ -@everywhere begin - - - """ - ResetParticle!(integrator) - (debubugging function) - Resets the integrator instance if the particle energy is nan or very high - resets the particles position to the domain center - resets the energy to a dfault value (e_0 is a global variable) - - """ - function ResetParticle!(integrator) - if isnan(integrator.ODEIntegrator.u[3]) || exp(integrator.ODEIntegrator.u[3]) >= 1e-3 - @show exp(integrator.ODEIntegrator.u[3]) - integrator.ODEIntegrator.u[1] = integrator.ODEIntegrator.u[1] - Lx/2 - #integrator.ODEIntegrator.u[2] = integrator.ODEIntegrator.u[2] - Ly/2 - #integrator.ODEIntegrator.u[4] = 1e-2 - integrator.ODEIntegrator.u[3] = e_0 - u_modified!(integrator.ODEIntegrator,true) - @show "rest particle" - end - nothing - end - - Lx_terminate_limit = 1 - - - """ - TerminateCheckSingle!(integrator) - (debubugging function) - - """ - function TerminateCheckSingle!(integrator) - if maximum(integrator.ODEIntegrator.u[3]) - Lx * Lx_terminate_limit >= 0 #|| maximum(exp.(integrator.u[3:N_state:end]) / e_0 ) >= 5 - terminate!(integrator.ODEIntegrator) - @show "terminate" - end - end -end - -###### remeshing routines ############ - - -""" - ParticleToNode!(PI::AbstractParticleInstance, S::SharedMatrix, G::TwoDGrid) -Pushes particle values to the neighboring nodes following the ParticleInCell rules. -1.) get weights and indexes of the neighboring notes, -2.) convert the particle state to nodestate -3.) push the calculated state to the shared arrays - -inputs: - -PI Particle instance -S Shared array where particles are stored -G (TwoDGrid) Grid that defines the nodepositions -""" -function ParticleToNode!(PI::AbstractParticleInstance, S::SharedMatrix, G::OneDGrid) - - index_positions, weights = ParticleInCell.compute_weights_and_index(G, PI.ODEIntegrator.u[1]) - #ui[1:2] .= PI.position_xy - #@show index_positions - u_state = GetParticleEnergyMomentum(PI.ODEIntegrator.u) - #@show u_state, index_positions, weights - ParticleInCell.push_to_grid!(S, u_state , index_positions, weights, G.Nx , periodic_boundary) - nothing -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!(PI::AbstractParticleInstance, S::SharedMatrix, ti::Number, e_0::Number) - u_state = Get_u_FromShared( PI, S) - - - last_t = PI.ODEIntegrator.t - # test if particle is below energy threshold, or - # if particle is at the boundary - if (u_state[1] < exp(e_0)) | PI.boundary # init new particle - - @show "re-init new particle" - # @show PI.position_ij, u_state - # z_i = copy(z0) - # z_i[x] = PI.position_xy[1] - # z_i[y] = PI.position_xy[2] - # ui = z_i - - u_init = u(PI.position_xy[1], ti) - # derive local wave speed and direction - cg_local = FetchRelations.c_g_U_tau( abs(u_init) , float(DT) ) - lne_local = log(FetchRelations.Eⱼ( abs(u_init) , float(DT) )) - - ui= [ PI.position_xy[1], cg_local, lne_local ] - reinit!(PI.ODEIntegrator, ui , erase_sol=false, reset_dt=true)#, reinit_callbacks=true) - #set_t!(PI.ODEIntegrator, last_t ) - u_modified!(PI.ODEIntegrator,true) - - else # load node value to particle - - #@show "get vertex variable" - #@show "u_state", u_state - ui= GetVariablesAtVertex( u_state, PI.position_xy[1] ) - #@show ui - set_u!(PI.ODEIntegrator, ui ) - #set_t!(PI.ODEIntegrator, last_t ) - u_modified!(PI.ODEIntegrator,true) - end - #@show ui - nothing - -end - -# %% -######### Core routines for advancing and remeshing -@everywhere begin - - """ - advance!(PI::AbstractParticleInstance, S::SharedMatrix{Float64}, G::OneDGrid, DT::Int) - """ - function advance!(PI::AbstractParticleInstance, S::SharedMatrix{Float64}, G::OneDGrid, DT::Int) - #@show PI.position_ij - - add_saveat!(PI.ODEIntegrator, PI.ODEIntegrator.t ) - savevalues!(PI.ODEIntegrator) - - step!(PI.ODEIntegrator, DT , true) - #@show PI.ODEIntegrator - - # use this for periodic boundary conditions - #periodic_BD_single_PI!(PI, Lx) #### !!!! not sure if the boundary condition has to be set there, it might have beeen set somewhere else as well .. - - # step!(PI.ODEIntegrator, DT/2 , true) - # PI = periodic_BD_single_PI!(PI ) - if isnan(PI.ODEIntegrator.u[3]) - @show "position is nan" - @show PI - PI.ODEIntegrator.u = [0,0,0] - end - - ParticleToNode!(PI, S, G) - return PI - end - - """ - remesh!(PI::AbstractParticleInstance, S::SharedMatrix{Float64, 3}) - Wrapper function that does everything necessary to remesh the particles. - - pushes the Node State to particle instance - """ - function remesh!(PI::AbstractParticleInstance, S::SharedMatrix{Float64}, ti::Number) - NodeToParticle!(PI, S, ti, e_0) - return PI - end - - - """ shows total energy of the all particles """ - function ShowTotalEnergyChange(ParticleCollection, u_sum_m1) - u_sum = zeros(Nstate) - for a_particle in ParticleCollection - u_sum += GetParticleEnergyMomentum(a_particle.ODEIntegrator.u) - end - @show u_sum_m1 - u_sum - return u_sum - end -end - - -# %% single core -@printf "Start itteration\n" -# define state collector -State_collect = [] # storing array -State[:,:,:] .= 0 # reset current state # [ Energy, x momentum, not used] -push!(State_collect, copy(State) ) # push initial state to storage - -# main time stepping: -t_range = range(0.0, DT*T/DT, step=DT) - - -for t in t_range[2:end] - #t = t_range[2] - @show t - - #@printf "advance" - State[:,:] .= 0 - #u_sum_m1 = ShowTotalEnergyChange(ParticleCollection, u_sum_m1) - - for a_particle in ParticleCollection - #@show a_particle.position_ij - #try - advance!(a_particle, State, grid1d, DT) - # catch err - # @show "err", a_particle.ODEIntegrator.u - # end - - # advance(a_particle, State) - # # check callback - #show_pos!( a_particle) - # ResetParticle!(a_particle) - end - - #u_sum_m1 = ShowTotalEnergyChange(ParticleCollection, u_sum_m1) - - #@show (State[:,:,1] .- grid1d.xmin)/grid1d.dx - #@printf "re-mesh" - for a_particle in ParticleCollection - #@show a_particle.position_ij - #@show "------------------------------------" - #show_pos!( a_particle)#, grid1d) - remesh!(a_particle, State, t) - #show_pos!( a_particle)#, grid1d) - #periodic_BD_single!(a_particle ) - #show_pos!( a_particle, grid1d) - end - - @printf "push" - push!(State_collect, copy(State) ) - - #u_sum_m1 = ShowTotalEnergyChange(ParticleCollection, u_sum_m1) - #@show (State[:,:,1] .- grid1d.xmin)/grid1d.dx -end - -# %% Parallel version - -# @everywhere carrier(f, y, g, d) = x -> f(x,y, g, d) -# @everywhere carrier(f, y, t) = x -> f(x,y, t) -# global ParticleCollection -# -# @time for t in t_range[2:end] -# -# @show t -# -# @show "advance" -# @time ParticleCollection = fetch(pmap(carrier( advance!, State, grid1d, DT) ,advance_workers, ParticleCollection)); -# -# @show "re-mesh" -# @time ParticleCollection = fetch(pmap(carrier( remesh!, State, t) ,advance_workers, ParticleCollection)); -# #@show mean(State) -# -# @show "push" -# push!(State_collect, State ) -# -# end - -# %% Checks - -@printf "\n size of collected states %s" size(State_collect) -@printf "\n size of each state %s" size(State_collect[1]) -@printf "\n length of timerange %s" size(t_range) - - -# %% Saving Fields - -@printf "save data \n" - -save_path = save_path_base*"dx"* @sprintf("%i",dx)*"_dt"*@sprintf("%i",DT)*"_U"*@sprintf("%i",U10)*"/" -mkpath(save_path) - -@printf "created model run folder %s" save_path - -rm(joinpath(save_path , "state.h5"), force=true) -file = h5open( joinpath(save_path , "state.h5") , "w" ) - -store_waves = create_group(file, "waves") -store_waves_data = create_dataset( - store_waves, "data", Float64, - (length(State_collect), (grid1dnotes.Nx) , (size(State_collect[1])[2]) ) - )#, chunk=(2,)) - -for i in 1:length(State_collect) - store_waves_data[i,:,:] = State_collect[i] -end - -#State_collect[1] - - -write_attribute(store_waves, "dims", ["time", "x", "state"]) -store_waves["time"] = Array(t_range) -store_waves["x"] = grid1dnotes.x -store_waves["var_names"] = ["e", "m_x", "m_y"] - - -store_winds = create_group( file, "wind") -store_winds_u = create_dataset(store_winds, "u", Float64, ( (length(xi)), (length(ti)) ) )#, chunk=(2,)) -store_winds_v = create_dataset(store_winds, "v", Float64, ( (length(xi)), (length(ti)) ) )#, chunk=(2,)) - -store_winds_u[:,:] = u_func_gridded -#store_winds_v[:,:,:] = v_func_gridded - - -write_attribute(store_winds, "dims", ["x", "time"]) -store_winds["time"] = Array(ti) -store_winds["x"] = Array(xi) -#store_winds["y"] = Array(yi) - -#attr = create_attribute(file, "dims", ["time", "x", "y", "state"]) -#HDF5.get_create_properties(state_store) - -close(file) - -# delete_object(parent, name) # for groups, datasets, and datatypes -# delete_attribute(parent, name) # for attributes - -# %% Save Particles -@printf "save particle \n" -save_object(joinpath(save_path , "particles.jld2"), ParticleCollection) -#PC2 = load_object(joinpath(save_path , "particles.jld2")) - -# %% -@printf "save parameters \n" -using JSON3 -params = Dict("grid" => Dict( - "Nx" => Nx, - "dx" => dx, - "Nstate" => Nstate, - "T" => T, - "dt" => DT - ), - "ODE_sampling" => Dict( - "dt_ODE_save" => dt_ODE_save - ), - "forcing" => Dict( - "U10" => U10 - ) - ) - -open( joinpath(save_path , "parameters.json") , "w") do io - JSON3.pretty(io, params) -end - -@printf "done.\n" diff --git a/src/old_structure/run_2d_manual.jl b/src/old_structure/run_2d_manual.jl deleted file mode 100644 index e56fb49..0000000 --- a/src/old_structure/run_2d_manual.jl +++ /dev/null @@ -1,786 +0,0 @@ - -using BenchmarkTools -#using Base.Threads - -# CPU initialization -using Distributed -addprocs(2) -advance_workers = WorkerPool([2, 3]) - - -@everywhere using Interpolations -@everywhere using IfElse -#@ redefining variables on all cores -@everywhere using ModelingToolkit, DifferentialEquations, Statistics - -@everywhere push!(LOAD_PATH, joinpath(pwd(), "MIC/") ) -@everywhere push!(LOAD_PATH, joinpath(pwd(), "code/") ) -@everywhere push!(LOAD_PATH, joinpath(pwd(), "analysis/") ) -#@everywhere import mesh -@everywhere using ParticleMesh -@everywhere using SharedArrays - -#include("./particle_waves_v1.jl") -@everywhere import PIC -@everywhere import particle_waves_v3 -@everywhere using Callbacks -@everywhere using Printf - - -@everywhere using ParticleMesh -using Revise -#using Plots -# %% - - -p_module = joinpath(pwd(), "code/particle_waves_v3.jl") -isfile(p_module) - -@everywhere begin - - @register_symbolic u(x, y, t) - @register_symbolic u_x(x, y, t) - @register_symbolic v(x, y, t) - @register_symbolic v_y(x, y, t) - - # %% create fake wind data and its gradients - T = 60*60*24 - Lx, Ly = 50e3, 50e3 - xi = range(0,Lx,step=3e3) # note ': this is a row vector - yi = range(0,Ly,step=3e3) - ti = range(0,T ,step=60*60) - - #u_func(x, y, t) = y * 0 .+ 3 * exp(- ( ( x-25e3 + 20e3* t/T )./10e3).^2) #.+ 3 #.+ 3 * sin.(x *π/Ly/0.5 .+ y *π/Ly/0.5) - u_func(x, y, t) = y * 0 .+ 3 * exp(- ( ( x-25e3 )./10e3).^2) #.+ 3 #.+ 3 * sin.(x *π/Ly/0.5 - #@everywhere u_func(x, y, t) = y *0 .- x .*2/50e3 .+ 3 #.+ 3 * sin.(x *π/Ly/0.5 .+ y *π/Ly/0.5) - #@everywhere u_func(x, y, t) = y *0 .- x .*0 .+ 3 #.+ 3 * sin.(x *π/Ly/0.5 .+ y *π/Ly/0.5) - - # @everywhere u_func(x, y, t) = y *0 .+ IfElse.ifelse.( x .< Lx/2.0, - # 2 .+ 1, *sin.(x *π/Lx/2), - # 0 .* x + 0.1 - # ) - - v_func(x, y, t) = y *0 .+ x .*0 .+ 2 * y *π/Ly/0.5 * cos.(y *π/Ly/0.5) - #v_func(x, y, t) = y * sin.(t * π ./T) * 2/Ly .+ x .*0 .+ 1* cos.( y *π/Ly/1 ) - #@everywhere v_func(x, y, t) = y *0 .+ x .*0 .+ 1* cos.(y *π/Ly/1 + t * 2 * pi /T ) - # @everywhere v_func(x, y, t) = x *0 .+ IfElse.ifelse.( x .< Lx/2.0, - # 2 * cos.(y *π/Ly/0.5) , - # 0 .* y + 0.1 - # ) - - - # % define wind forcing globally as interp functions - nodes = (xi, yi, ti) - - u_func_gridded = [ u_func(xii, yii, tii) for xii in xi, yii in yi, tii in ti] - u_grid = LinearInterpolation( nodes , u_func_gridded ,extrapolation_bc=Periodic()) - #u_grid = CubicSplineInterpolation( nodes, u_func_gridded; bc=Line(OnGrid()), extrapolation_bc=Flat()) - #u_grid = interpolate(nodes, u_func_gridded, Gridded(Linear())) - u(x, y, t) = u_grid(x, y, t) - - - v_func_gridded = [ v_func(xii, yii, tii) for xii in xi, yii in yi, tii in ti] - v_grid = LinearInterpolation( nodes,v_func_gridded ,extrapolation_bc=Periodic()) - #v_grid = CubicSplineInterpolation( nodes, v_func_gridded; bc=Line(OnGrid()), extrapolation_bc=Periodic()) - - #v_grid = interpolate(nodes, v_func_gridded, Gridded(Linear())) - #v_grid = extrapolate( interpolate(v_func_gridded, BSpline(Linear()) ) , nodes) - v(x, y, t) = v_grid(x, y, t)#v_grid(x, y) - -end -# % plotting focring field -#using GLMakie - - -# create mesh -xmesh = [ xii/1e3 for xii in xi, yii in yi] -ymesh = [ yii/1e3 for xii in xi, yii in yi] - -#arrows( xi/1e3, yi/1e3 , u_func_gridded[:,:,1], v_func_gridded[:,:,1] ,arrowsize = 10, lengthscale = 2 )# , color = "black" )# |> display -#arrows( xmesh, ymesh , u_func_gridded[:,:,1], v_func_gridded[:,:,1] )# , color = "black" )# |> display - - -# plot 1st and last wind field -#arrows( xmesh, ymesh , quiver = (u_func_gridded[:,:,1], v_func_gridded[:,:,1]) , color = "black" ) |> display - -############## -# quiver( xmesh, ymesh , quiver = (u_func_gridded[:,:,1], v_func_gridded[:,:,1]) , color = "black" ) |> display -# -# quiver!( xmesh, ymesh , quiver = (u_func_gridded[:,:,Int(floor(length(ti)))], v_func_gridded[:,:,Int(floor(length(ti)/2))]) , color = "red" ) |> display -############## - -#quiver!( xmesh, ymesh , quiver = (u_func_gridded[:,:,Int(floor(length(ti)/2))], v_func_gridded[:,:,Int(floor(length(ti)/2))]) , color = "red" ) |> display -# angles="xy", scale_units="width", scale=100 , headwidth=2, width= 0.005, headlength=3, color = "black", alpha= 0.5) - -# %% define all components to solve ODE -@everywhere begin - - # only x gradient - u_x(x, y, t) = Interpolations.gradient(u_grid, x, y, t )[1] - - # only y gradient - v_y(x, y, t) = Interpolations.gradient(v_grid, x,y, t )[2] - - particle_equations0 = particle_waves_v3.particle_equations(u, v, u_x, v_y) - - @named particle_system0 = ODESystem(particle_equations0) - - # define variables based on particle equations - t, x, y, c̄_x, c̄_y, lne, Δn, Δφ_p, r_g, C_α, C_φ, g, C_e = particle_waves_v3.init_vars() - - # ParticleInstance is the Stucture that carries each particle. - mutable struct ParticleInstance - position_ij :: Tuple - position_xy :: Tuple - ODEIntegrator::OrdinaryDiffEq.ODEIntegrator - - # dx::Float64 - # dy::Float64 - # function ParticleInstance(position_xy, ODE, dx::Float64, dy::Float64) - # new( ( Int(round( position_xy[1]/dx)), Int(round( position_xy[2]/dy)) ) , position_xy ,ODE) - # end - end - -end - - -# %% define storing stucture and populate inital conditions - - -@everywhere Nx, Ny = 30, 15 -@everywhere grid2d = ParticleMesh.TwoDGrid(1e3, Lx-1e3, Nx, 1e3, Ly-1e3 , Ny) -@everywhere grid2dnotes = ParticleMesh.TwoDGridNotes(grid2d) -@everywhere dx, dy = grid2d.dx , grid2d.dy -@everywhere Nstate = 3 - - -Nparticle = grid2d.Nx* grid2d.Ny -State = SharedArray{Float64}(grid2d.Nx, grid2d.Ny, Nstate) -State.pids # check which cores can see state - - -# function wrap_index(pos::Int, N::Int) -# if pos < 1 -# pos = pos + Int(N) -# #pos = N -# elseif pos > N -# pos = pos - Int(N) -# #pos = N -# end -# return pos -# end - -# function push_to_2d_grid_local!(grid::SharedArray{Float64, 3}, -# charge::Vector{Float64}, -# index_pos::Tuple{Int, Int}, -# weights::Tuple{Float64, Float64}, -# Nx::Int, Ny::Int ) -# -# grid[ wrap_index(index_pos[1], Nx) , wrap_index(index_pos[2], Ny), : ] += weights[1] * weights[2] * charge -# #grid[ index_pos[1] , index_pos[2] ] += weights[1] * weights[2] * charge -# -# end - - -#index_positions, weights = PIC.compute_weights_and_index(grid2d, 12.34e3, 13.4e3) -#PIC.push_to_grid!(State, [1.0, 1.0,1.0,1.0,1.0,1.0,1.0] , index_positions[1], weights[1], grid2d.Nx, grid2d.Ny ) - -# PIC.push_to_grid!(State, [1.0, 1.0,1.0,1.0,1.0,1.0,1.0] , index_positions, weights, grid2d.Nx, grid2d.Ny ) - -# %% - -# parameters -params0 = Dict( - r_g => 0.9, - C_α => -1.41 , - C_φ => 1.8e-5, - g => 9.81 , - C_e => 2.16e-4) - -@everywhere begin - -# define default initial state -e_0 = log(1e-4) -z0 = Dict( x => 0.0, y=> 0.0, c̄_x=> 1e-2, c̄_y=> -1e-9, lne=> e_0, Δn => dx, Δφ_p => 0.0) - - -# # copy State to shared array -# State[1, :,:] = [ z0[x] for i in range(0,length = Nparticle) ] -# -# State[1, :,1] = [ z0[x] for i in range(0,length = Nparticle) ] -# State[1, :,2] = [ 0 for i in y_range ] - - -# @everywhere function init_z0_to_State!(S::SharedArray, ij::Tuple{Int,Int}, zi::Vector{Float64} ) -# S[ ij[1], ij[2] , :] = [ zi[x],zi[y], zi[c̄_x],zi[c̄_y], zi[lne], zi[Δn],zi[Δφ_p] ]#collect(values(zi)) -# nothing -# end -""" sets node state values to S at tuple position ij to zi """ -function init_z0_to_State!(S::SharedArray, ij::Tuple{Int,Int}, zi::Vector{Float64} ) - S[ ij[1], ij[2] , :] = zi - nothing -end - -""" sets particle state values to S. position is taking from particle """ -function set_u_to_shared!(S::SharedArray, PI::ParticleInstance) - S[ PI.position_ij[1], PI.position_ij[2] , :] = PI.ODEIntegrator.u - nothing -end - -""" takes node state values and pushes them to particle. position is taken from particle """ -function set_u_to_particle!(S::SharedArray, PI::ParticleInstance) - set_u!(PI.ODEIntegrator, Get_u_FromShared(S, PI) ) - u_modified!(PI.ODEIntegrator,true) - nothing -end - -end -# %%%% Define callbacks - -@everywhere begin - -DT = 60*30 -dt_save = DT/3.0 - -""" -wrap_pos!(xx::Float64, L::Float64) -makes periodic boundary conditions. If position xx exceeds 0 or L it is wrapped around -returns new position and if either or not the wrapping appends (bool) -""" -function wrap_pos!(xx::Float64, L::Float64) - wrap_flag = false - if xx > L - xx = xx - L - wrap_flag = true - elseif xx < 0 - xx = xx + L - wrap_flag = true - end - - return xx, wrap_flag -end - -""" -Checks if particle's PI position PI.ODEIntegrator.u[1,2] exceeeds the domain limits Lx, Ly - Lx, LY are global variables - If the PI's postions eceeeds the boundary they warpped around and the ODEIntegrator state is modified. -""" -function periodic_BD_single_PI!(PI::ParticleInstance) - - #@printf "PI periodic condition called" - ui = copy(PI.ODEIntegrator.u) - ui[4], wrap_pos_PI1 = wrap_pos!(ui[4], Lx) - ui[5], wrap_pos_PI2 = wrap_pos!(ui[5], Ly) - - if wrap_pos_PI1 || wrap_pos_PI1 - #@show wrap_pos_PI1 , wrap_pos_PI2 - #@printf "wrap pos PI" - #@show PI.ODEIntegrator.u[4] - ui[4] - #@show PI.ODEIntegrator.u[5] - ui[5] - set_u!(PI.ODEIntegrator, ui ) - u_modified!(PI.ODEIntegrator,true) - end - nothing #return PI -end - - -# """ -# Checks if -# """ -# @everywhere function periodic_BD_single!(S::SharedArray) -# S.u[1], wrap_pos_PI = wrap_pos!(S.u[1], Lx) -# S.u[2], wrap_pos_PI = wrap_pos!(S.u[2], Ly) -# nothing -# end - -""" -return the mean position -""" -function show_pos!(PI) - @show "show current position" PI.position_xy, PI.ODEIntegrator.u[4], PI.ODEIntegrator.u[5] -end - -# the following function are another way to have wrapping boundary conditions. - -function periodic_condition_x(u,t,integrator) - u[4] > 0 -end - -function periodic_condition_y(u,t,integrator) - u[5] > Ly -end - -# function loop_x!(integrator) -# integrator.u[1] = integrator.u[1] - Lx -# @printf "wrap pos PI x" -# end -# -# function loop_y!(integrator) -# integrator.u[2] = integrator.u[2] - Ly -# @printf "wrap pos PI y" -# end - - -end # end of everywhere - -# terminate = PeriodicCallback(terminate_check!, 60*12 ) -show_mean = PeriodicCallback(show_pos! , DT/2 ) -# set_to_mean = PeriodicCallback(set_to_mean! , 60*60*6 ) -periodic = PeriodicCallback(periodic_BD_single_PI! , DT/2 ) -# -#periodic_x = DiscreteCallback(periodic_condition_x, loop_x!)# , interp_points=0) -#periodic_y = DiscreteCallback(periodic_condition_y, loop_y!)# , interp_points=0) - -cbs = CallbackSet(periodic, show_mean )#,cb_terminate) - - - - -######### Particle --> Node ########## -@everywhere begin - - """ - GetParticleEnergyMomentum(PI) - - """ - function GetParticleEnergyMomentum(PI) - ui_x, ui_y, ui_c̄_x, ui_c̄_y, ui_lne, ui_Δn, ui_Δφ_p = PI.ODEIntegrator.u - - ui_e = exp(ui_lne) - c_speed = particle_waves_v3.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 - end - - function GetParticleEnergyMomentum(PI::ParticleInstance) - 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[Δn], zi[Δφ_p]]) - end - - """ - GetParticleEnergyMomentum(PI) - - """ - function GetParticleEnergyMomentum(z0::Vector{Float64}) - - ui_x, ui_y, ui_c̄_x, ui_c̄_y, ui_lne, ui_Δn, ui_Δφ_p = z0 - ui_e = exp(ui_lne) - c_speed = particle_waves_v3.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] - end - - # testing - #GetParticleEnergyMomentum(z0) - #ui_x, ui_y, ui_c̄_x, ui_c̄_y, ui_lne, ui_Δn, ui_Δφ_p = PI.ODEIntegrator.u - #ie, imx, imy = GetParticleEnergyMomentum(a_particle.ODEIntegrator.u) - - ######### node--> Particle ########## - """ - GetVariablesAtVertex(i_State::Vector{Float64}, x::Float64, y::Float64, Δn::Float64, Δφ_p::Float64 ) - - """ - function GetVariablesAtVertex(i_State::Vector{Float64}, x::Float64, y::Float64, Δn::Float64, Δφ_p::Float64 ) - e, m_x, m_y = i_State - m_amp = particle_waves_v3.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, Δn, Δφ_p] - end -end -# testing -#GetVariablesAtVertex( [ie, imx, imy], 2e3, 4e3, dx, 0.0 ) - a_particle.ODEIntegrator.u - - -# %% seeding particles with initial states - -""" -InitParticleInstance(model, z_initials, pars, ij ; 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 - chSet (optional) is the set of callbacks the ODE can have -""" -function InitParticleInstance(model, z_initials, pars, ij ; cbSets=nothing) - - # create ODEProblem - problem = ODEProblem(model, z_initials, (0.0, T) , pars) - # inialize problem - #solver_method = Rosenbrock23() - #solver_method = AutoVern7(Rodas4()) - solver_method = AutoTsit5(Rosenbrock23()) - #solver_method =Tsit5() - integrator = init(problem, solver_method , saveat =dt_save , abstol = 1e-3, adaptive =true, maxiters=1e3)#, reltol=1e-1) - #abstol=1e-4)#, reltol=1e-1) - - # callbacks= cbSets, - #save_everystep=false - return ParticleInstance( ij , (z_initials[x], z_initials[y]), integrator) -end - - -ParticleCollection=[] -for i in range(1,length = Nx), j in range(1,length = Ny) - - - # initalize state based on state vector - z_i = copy(z0) - z_i[x] = grid2dnotes.x[i] - z_i[y] = grid2dnotes.y[j] - - # take in local wind velocities - u_init ,v_init = u(z_i[x], z_i[y], 0) , v(z_i[x], z_i[y], 0) - # derive local wave speed and direction - # ! replace by wind sea formula in the future - z_i[c̄_x] = 1e-1 * u_init / sqrt(u_init.^2 + v_init.^2) - z_i[c̄_y] = 1e-1 * v_init / sqrt(u_init.^2 + v_init.^2) - - #@show (z_i[x]-grid2d.xmin)/grid2d.dx, (z_i[y]-grid2d.ymin)/grid2d.dy - - # push z_i to shared array - init_z0_to_State!(State, (i,j), GetParticleEnergyMomentum(z_i) ) - - # copy initial parameters default - params_i = copy(params0) - - # push z_i and params_i to particle system - # particle_system0 is an initilized ODESystem - push!(ParticleCollection , InitParticleInstance(particle_system0, z_i , params_i, (i,j) )) - #@show threadid() -end - -###### Debugging functions ############ -@everywhere begin - - """ (debubugging function) returns the domains normalized position of the partivle """ - function show_pos!(integrator , G) - @show (integrator.ODEIntegrator.u[4] .- G.xmin)/grid2d.dx, (integrator.ODEIntegrator.u[5] .- G.ymin)/grid2d.dy - end - - - """ - ResetParticle!(integrator) - (debubugging function) - Resets the integrator instance if the particle energy is nan or very high - resets the particles position to the domain center - resets the energy to a dfault value (e_0 is a global variable) - - """ - function ResetParticle!(integrator) - if isnan(integrator.ODEIntegrator.u[1]) || exp(integrator.ODEIntegrator.u[1]) >= 1e-3 - @show exp(integrator.ODEIntegrator.u[1]) - integrator.ODEIntegrator.u[4] = integrator.ODEIntegrator.u[4] - Lx/2 - integrator.ODEIntegrator.u[5] = integrator.ODEIntegrator.u[5] - Ly/2 - #integrator.ODEIntegrator.u[4] = 1e-2 - integrator.ODEIntegrator.u[1] = e_0 - u_modified!(integrator.ODEIntegrator,true) - @show "rest particle" - end - nothing - end - - Lx_terminate_limit = 1 - - - """ - TerminateCheckSingle!(integrator) - (debubugging function) - - """ - function TerminateCheckSingle!(integrator) - if maximum(integrator.ODEIntegrator.u[4]) - Lx * Lx_terminate_limit >= 0 #|| maximum(exp.(integrator.u[3:N_state:end]) / e_0 ) >= 5 - terminate!(integrator.ODEIntegrator) - @show "terminate" - end - end - -end -#show_pos!(a_particle ) -#advance(a_particle, State) - -###### define remeshing routines ############ -@everywhere begin - - """ - ParticleToNode!(PI::ParticleInstance, S::SharedArray, G::TwoDGrid) - Pushes particle values to the neighboring nodes following the PIC rules. - 1.) get weights and indexes of the neighboring notes, - 2.) convert the particle state to nodestate - 3.) push the calculated state to the shared arrays - - inputs: - - PI Particle instance - S Shared array where particles are stored - G (TwoDGrid) Grid that defines the nodepositions - """ - function ParticleToNode!(PI::ParticleInstance, S::SharedArray, G::TwoDGrid) - - index_positions, weights = PIC.compute_weights_and_index(G, PI.ODEIntegrator.u[4], PI.ODEIntegrator.u[5]) - #ui[1:2] .= PI.position_xy - #@show index_positions - u_state = GetParticleEnergyMomentum(PI.ODEIntegrator.u) - #@show u_state - PIC.push_to_grid!(S, u_state , index_positions, weights, G.Nx, G.Ny ) - nothing - end - - """ returns node state from shared array, given paritcle index """ - Get_u_FromShared(PI::ParticleInstance, S::SharedArray, ) = S[ PI.position_ij[1], PI.position_ij[2] , :] - #u_state = Get_u_FromShared(a_particle, State ) - - """ - NodeToParticle!(PI::ParticleInstance, S::SharedArray) - 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!(PI::ParticleInstance, S::SharedArray, ti::Number) - u_state = Get_u_FromShared( PI, S) - - - last_t = PI.ODEIntegrator.t - if u_state[1] <= exp(e_0) # init new particle - - # @show "re-init new particle" - # @show PI.position_ij, u_state - # z_i = copy(z0) - # z_i[x] = PI.position_xy[1] - # z_i[y] = PI.position_xy[2] - # ui = z_i - - u_init ,v_init = u(PI.position_xy[1], PI.position_xy[2], ti) , v(PI.position_xy[1], PI.position_xy[2], ti) - # derive local wave speed and direction - # ! replace by wind sea formula in the future - u_local = 1e-1 * u_init / sqrt(u_init.^2 + v_init.^2) - v_local = 1e-1 * v_init / sqrt(u_init.^2 + v_init.^2) - - ui= [ PI.position_xy[1], PI.position_xy[2], u_local, v_local, z0[lne], z0[Δn], z0[Δφ_p] ] - reinit!(PI.ODEIntegrator, ui , erase_sol=false, reset_dt=true)#, reinit_callbacks=true) - #set_t!(PI.ODEIntegrator, last_t ) - u_modified!(PI.ODEIntegrator,true) - - else # load node value to particle - - ui= GetVariablesAtVertex( u_state, PI.position_xy[1], PI.position_xy[2], z0[Δn], z0[Δφ_p] ) - set_u!(PI.ODEIntegrator, ui ) - #set_t!(PI.ODEIntegrator, last_t ) - u_modified!(PI.ODEIntegrator,true) - end - #@show ui - nothing - end - - """ - advance!(PI::ParticleInstance, S::SharedArray) - """ - function advance!(PI::ParticleInstance, S::SharedArray{Float64, 3}, G::TwoDGrid, DT::Int) - #@show PI.position_ij - - add_saveat!(PI.ODEIntegrator, PI.ODEIntegrator.t ) - savevalues!(PI.ODEIntegrator) - - step!(PI.ODEIntegrator, DT , true) - periodic_BD_single_PI!(PI ) - - - # step!(PI.ODEIntegrator, DT/2 , true) - # PI = periodic_BD_single_PI!(PI ) - ParticleToNode!(PI, S, G) - return PI - end - - """ - remesh!(PI::ParticleInstance, S::SharedArray{Float64, 3}) - Wrapper function that does everything necessary to remesh the particles. - - pushes the Node State to particle instance - """ - function remesh!(PI::ParticleInstance, S::SharedArray{Float64, 3}, ti::Number) - NodeToParticle!(PI, S, ti) - return PI - end - - - """ shows total energy of the all particles """ - function ShowTotalEnergyChange(ParticleCollection, u_sum_m1) - u_sum = zeros(Nstate) - for a_particle in ParticleCollection - u_sum += GetParticleEnergyMomentum(a_particle.ODEIntegrator.u) - end - @show u_sum_m1 - u_sum - return u_sum - end -end - -State_collect = [] # storing array -State[:,:,:] .= 0 # reset current state -push!(State_collect, copy(State) ) # push initial state to storage - -u_sum_m1 = zeros(Nstate) - -T/DT -DT -# main time stepping: -t_range = range(0.0, DT*T/DT, step=DT) - - - -# %% single core -for t in t_range[2:end] - - @show t - - @show "advance" - State[:,:,:] .= 0 - #u_sum_m1 = ShowTotalEnergyChange(ParticleCollection, u_sum_m1) - - for a_particle in ParticleCollection - #@show a_particle.position_ij - #try - advance!(a_particle, State, grid2d, DT) - # catch err - # @show "err", a_particle.ODEIntegrator.u - # end - - # advance(a_particle, State) - # # check callback - # # show_pos!( a_particle) - # ResetParticle!(a_particle) - end - - #u_sum_m1 = ShowTotalEnergyChange(ParticleCollection, u_sum_m1) - - #@show (State[:,:,1] .- grid2d.xmin)/grid2d.dx - @show "re-mesh" - for a_particle in ParticleCollection - #@show a_particle.position_ij - #show_pos!( a_particle, grid2d) - remesh!(a_particle, State, t) - - #periodic_BD_single!(a_particle ) - #show_pos!( a_particle, grid2d) - end - - @show "push" - push!(State_collect, copy(State) ) - - #u_sum_m1 = ShowTotalEnergyChange(ParticleCollection, u_sum_m1) - #@show (State[:,:,1] .- grid2d.xmin)/grid2d.dx -end - - - -# %% Parallel version - -@everywhere carrier(f, y, g, d) = x -> f(x,y, g, d) -@everywhere carrier(f, y, t) = x -> f(x,y, t) -global ParticleCollection - -@time for t in t_range[2:end] - - @show t - - @show "advance" - @time ParticleCollection = fetch(pmap(carrier( advance!, State, grid2d, DT) ,advance_workers, ParticleCollection)); - - @show "re-mesh" - @time ParticleCollection = fetch(pmap(carrier( remesh!, State, t) ,advance_workers, ParticleCollection)); - #@show mean(State) - - @show "push" - push!(State_collect, State ) - -end - - -# %% - -heatmap( State_collect[2][:,:,1]' , levels = 20, color=:blues ) |> display - - -heatmap( State_collect[6][:,:,1]' , levels = 20, color=:blues ) |> display - -# %% -heatmap( State_collect[10][:,:,1]' , levels = 20, color=:blues ) |> display - - - - -# %% - -@printf "size of collected states %s" size(State_collect) -@printf "size of each state %s" size(State_collect[1]) -@printf "length of timerange %s" size(t_range) - - -# %% saving files - -@show "save data" -using HDF5 - -save_path = "data/first_try/" -mkpath(save_path) -rm(joinpath(save_path , "state.h5"), force=true) -file = h5open( joinpath(save_path , "state.h5") , "w" ) - -store_waves = create_group(file, "waves") -store_waves_data = create_dataset(store_waves, "data", Float64, (length(State_collect), (grid2dnotes.Nx), (grid2dnotes.Ny), (size(State_collect[1])[3]) ) )#, chunk=(2,)) - -for i in 1:length(State_collect) - store_waves_data[i,:,:,:] = State_collect[i] -end - -#State_collect[1] - - -write_attribute(store_waves, "dims", ["time", "x", "y", "state"]) -store_waves["time"] = Array(t_range) -store_waves["x"] = grid2dnotes.x -store_waves["y"] = grid2dnotes.y -store_waves["var_names"] = ["e", "m_x", "m_y"] - - -store_winds = create_group( file, "wind") -store_winds_u = create_dataset(store_winds, "u", Float64, ( (length(xi)), (length(yi)), (length(ti)) ) )#, chunk=(2,)) -store_winds_v = create_dataset(store_winds, "v", Float64, ( (length(xi)), (length(yi)), (length(ti)) ) )#, chunk=(2,)) - -store_winds_u[:,:,:] = u_func_gridded -store_winds_v[:,:,:] = v_func_gridded - - -write_attribute(store_winds, "dims", ["x", "y", "time"]) -store_winds["time"] = Array(ti) -store_winds["x"] = Array(xi) -store_winds["y"] = Array(yi) - - -#State_collect[1] -#attr = create_attribute(file, "dims", ["time", "x", "y", "state"]) -#HDF5.get_create_properties(state_store) - - -close(file) - -# delete_object(parent, name) # for groups, datasets, and datatypes -# delete_attribute(parent, name) # for attributes - - - -# %% save particles -@show "save particle" -using JLD2 - -save_object(joinpath(save_path , "particles.jld2"), ParticleCollection) -#PC2 = load_object(joinpath(save_path , "particles.jld2")) diff --git a/src/old_structure/run_collections.jl b/src/old_structure/run_collections.jl deleted file mode 100644 index e1f128c..0000000 --- a/src/old_structure/run_collections.jl +++ /dev/null @@ -1,188 +0,0 @@ -using Printf - - -# using Distributed -# addprocs(2) -# advance_workers = WorkerPool([2, 3]) - -push!(LOAD_PATH, joinpath(pwd(), "code/") ) -# %% - -function make_str(name::String, value::Number) - return name*@sprintf("%i",value) -end - - -function make_str(name::String, value::Number, digits::Number) - return name*@sprintf("%0.2f",round(value, digits= digits) ) -end - - -function concat_str(alist) - str_list ="" - for ai in alist - str_list *= ai * "_" - end - return str_list[1:end-1] -end - - -function concat_str_space(alist) - str_list ="" - for ai in alist - str_list *= ai * " " - end - return str_list[1:end-1] -end - -make_str - -# %% - -# parset_name = "U10-DT" -# #var1_list = [1, 3, 5, 8, 10, 13, 15, 20] -# var1_list = [20] -# var2_list = [3, 5, 8, 10, 15, 20, 30, 40, 60] -# var3_list = [50]#10, 20, 20] -# periodic = false - -# parset_name = "U10-DT-periodic" -# var1_list = [1, 3, 5, 8, 10, 13, 15, 18, 20, 25] # U10 -# var2_list = [3, 5, 8, 10, 15, 20, 30, 40, 60] # dt -# var3_list = [40]#10, 20, 20] # NZ -# periodic = true - -parset_name = "U10-DT" -var1_list = [1, 3, 5, 8, 10, 13, 15, 18, 20, 25] # U10 -var2_list = [3, 5, 8, 10, 15, 20, 30, 40, 45, 50, 60] # dt -var3_list = [40]#10, 20, 20] # NZ -periodic = false - -# parset_name = "Nx-DT" -# var1_list = [10] -# var2_list = [3, 5, 10, 15, 20, 30, 60] -# var3_list = [5, 8, 10, 20, 30, 40, 50, 80] -# periodic = false - - -# parset_name = "Nx-DT-periodic" -# var1_list = [10] -# var2_list = [3, 5, 10, 15, 20, 30, 60] -# var3_list = [5, 8, 10, 20, 30, 40, 50, 80] -# periodic = true - -# parset_name = "gamma-c_beta" -# var1_list = round.(range(0.88* 0.2 , 0.88* 1.5; length= 10), digits= 2) # gamma -# var2_list = [2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6] # c_beta_parameter -# var3_list = [40]#10, 20, 20] # Nx -# periodic = false - -# -# parset_name = "gamma-c_beta-rg0.85" -# var1_list = round.(range(0.88* 0.2 , 0.88* 1.5; length= 10), digits= 2) # gamma -# var2_list = [2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6] # c_beta_parameter -# var3_list = [0.85]#10, 20, 20] -# periodic = false -# -# -parset_name = "gamma-rg0_2" -var1_list = round.(range(0.88* 0.2 , 0.88* 1.5; length= 6), digits= 2) # gamma\ # gamma -#var1_list = round.(range(0.88* 0.2 , 0.88* 1.5; length= 10), digits= 2) # gamma\ # gamma -sort!(push!(var1_list, 0.88)) -var2_list = round.(range(0.80 , 1/0.80; length= 10), digits= 2) -#[2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6] # c_beta_parameter -var3_list = [50]#10, 20, 20] -periodic = false - -str_list =[] -arg_list = [] -for var1 in var1_list - for var2 in var2_list - for var3 in var3_list - push!(str_list , ) - - ID = concat_str([ make_str("U", var1), make_str("DT", var2),make_str("Nx", var3) ]) # "U10-DT" - #ID = concat_str([ make_str("rg", var2, 2), make_str("gamma", var1, 2),make_str("Nx", var3) ]) - - #ID = concat_str([ make_str("gamma", var1, 2), make_str("cbeta", var2, 1),make_str("Nx", var3) ]) # "gamma-c_beta" - ID = concat_str([ make_str("gamma", var1, 2), make_str("rg", var2, 1),make_str("Nx", var3) ]) # "gamma-rg_2" - - arg_case = [ "--ID", ID, - #"--U10", @sprintf("%i",var1), - #"--c_beta", @sprintf("%0.2f",var2), - "--gamma", @sprintf("%0.2f",var1), - "--r_g0", @sprintf("%0.2f",var2), - "--Nx", @sprintf("%i",var3), - #"--DT", @sprintf("%i",var2), - "--parset", parset_name] - if periodic - push!(arg_case, "--periodic") - end - push!(arg_list , arg_case) - @show arg_case - end - end -end - - - -#localARGS= arg_list[5] - -# -#aa = ["--ID", "Lx30_DT30_U10", "--Lx", 30, "--DT", 30, "--U", 10] -#push!(aa, "--periodic") -# ARGS = [ "--ID", "test1", -# "--Lx", "20", -# "--T", "24", -# "--DT", "30"] - -#include("./run_fake.jl")#, ARGS=["--foo", "123"]) -# %% -localARGS = arg_list[3] - -#mkdir(pwd()* "/data/" * parset * "/" ) -completed_runs = readdir( pwd()* "/data/" * parset_name * "/" ) - - -@show completed_runs -#@show arg_list -# include( joinpath(pwd(), "parameter_test/") * "/run_fake.jl") -# localARGS - -for args in arg_list - global localARGS = args - #push!(empty!(localARGS), args) - @show localARGS - - if args[2] in completed_runs - @printf "runner : case exist, obmitted " - else - @printf "runner : run code: " - #include( joinpath( pwd(), "parameter_test/") * "/run_fake.jl") - include( joinpath( pwd(), "code/") * "/run_1D.jl") - end - @printf "\n runner : next\n" - -end - -#run(pipeline( `julia -e 'print("hello")' ` , `grep el`)) -#localARGS -#run( `julia` ) - -# %% -# function do_run(args) -# localARGS = args -# #@show "before pass", localARGS -# include("./run_fake.jl") -# end - -#@everywhere carrier(argi) = x -> arg_list[x] - -# @everywhere carrier(f) = x -> f(x) -# -# -# carrier(do_run)(arg_list[1]) -# -# carrier( remesh!, State, t) -# -# ParticleCollection = fetch(pmap(carrier( remesh!, State, t) ,advance_workers, arg_list )); diff --git a/src/old_structure/run_collections_varying_forcing.jl b/src/old_structure/run_collections_varying_forcing.jl deleted file mode 100644 index 3c52b2f..0000000 --- a/src/old_structure/run_collections_varying_forcing.jl +++ /dev/null @@ -1,145 +0,0 @@ -using Printf - - -# using Distributed -# addprocs(2) -# advance_workers = WorkerPool([2, 3]) - -push!(LOAD_PATH, joinpath(pwd(), "code/") ) -# %% - -function make_str(name::String, value::Number) - return name*@sprintf("%i",value) -end - - -function make_str(name::String, value::Number, digits::Number) - return name*@sprintf("%0.2f",round(value, digits= digits) ) -end - - -function concat_str(alist) - str_list ="" - for ai in alist - str_list *= ai * "_" - end - return str_list[1:end-1] -end - - -function concat_str_space(alist) - str_list ="" - for ai in alist - str_list *= ai * " " - end - return str_list[1:end-1] -end - -make_str - -# %% - -base_parset="1D_gaussian" -parset_name = "Nx-DT" -var1_list = [20,50, 100, 150, 200, 300] -var2_list = [5, 10, 15, 20, 30, 40, 50, 60] -periodic = false - - -# parset_name = "Nx-DT-periodic" -# var1_list = [10] -# var2_list = [3, 5, 10, 15, 20, 30, 60] -# var3_list = [5, 8, 10, 20, 30, 40, 50, 80] -# periodic = true - - -str_list =[] -arg_list = [] -for var1 in var1_list - for var2 in var2_list - push!(str_list , ) - - ID = concat_str([ make_str("Nx", var1), make_str("DT", var2) ]) # "U10-DT" - #ID = concat_str([ make_str("rg", var2, 2), make_str("gamma", var1, 2),make_str("Nx", var3) ]) - - #ID = concat_str([ make_str("gamma", var1, 2), make_str("cbeta", var2, 1),make_str("Nx", var3) ]) # "gamma-c_beta" - - arg_case = [ "--ID", ID, - "--Nx", @sprintf("%i",var1), - #"--U10", @sprintf("%i",var1), - #"--c_beta", @sprintf("%0.2f",var2), - #"--gamma", @sprintf("%0.2f",var1), - #"--rg", @sprintf("%0.2f",var2), - "--DT", @sprintf("%i",var2), - "--parset", parset_name] - if periodic - push!(arg_case, "--periodic") - end - push!(arg_list , arg_case) - @show arg_case - end -end - - - -#localARGS= arg_list[5] - -# -#aa = ["--ID", "Lx30_DT30_U10", "--Lx", 30, "--DT", 30, "--U", 10] -#push!(aa, "--periodic") -# ARGS = [ "--ID", "test1", -# "--Lx", "20", -# "--T", "24", -# "--DT", "30"] - -#include("./run_fake.jl")#, ARGS=["--foo", "123"]) -mkdir(pwd()* "/data/"*base_parset*"/" * parset_name * "/") - -# %% -localARGS = arg_list[3] - -completed_runs = readdir( pwd()* "/data/"*base_parset*"/" * parset_name * "/" ) - - -@show completed_runs -#@show arg_list -# include( joinpath(pwd(), "parameter_test/") * "/run_fake.jl") -# localARGS - -for args in arg_list - localARGS = args - #push!(empty!(localARGS), args) - @show localARGS - - if args[2] in completed_runs - @printf "runner : case exist, obmitted " - else - @printf "runner : run code: " - #include( joinpath( pwd(), "parameter_test/") * "/run_fake.jl") - include( joinpath( pwd(), "code/") * "/run_1D_time_varing_forcing.jl") - end - @printf "\n runner : next\n" - -end - -#run(pipeline( `julia -e 'print("hello")' ` , `grep el`)) -#localARGS -#run( `julia` ) - -# %% -# function do_run(args) -# localARGS = args -# #@show "before pass", localARGS -# include("./run_fake.jl") -# end - -#@everywhere carrier(argi) = x -> arg_list[x] - -# @everywhere carrier(f) = x -> f(x) -# -# -# carrier(do_run)(arg_list[1]) -# -# carrier( remesh!, State, t) -# -# ParticleCollection = fetch(pmap(carrier( remesh!, State, t) ,advance_workers, arg_list ));