first commit

.gitignore

@@ -1 +1,3 @@
 /Manifest.toml
+log/
+videos/

Project.toml

@@ -1,7 +1,22 @@
 name = "RDE"
 uuid = "1d57ab14-0807-4ee5-93aa-665821a12dfd"
 authors = ["Kristian Holme"]
-version = "1.0.0-DEV"
+version = "0.0.1-DEV"
+
+[deps]
+CommonRLInterface = "d842c3ba-07a1-494f-bbec-f5741b0a3e98"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+Observables = "510215fc-4207-5dde-b226-833fc4488ee2"
+POMDPTools = "7588e00f-9cae-40de-98dc-e0c70c48cdd7"
+POMDPs = "a93abf59-7444-517b-a68a-c42f96afdd7d"
+PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
+ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
 
 [compat]
 julia = "1.6.7"

README.md

@@ -1,3 +1,45 @@
 # RDE
 
 [![Build Status](https://github.com/KristianHolme/RDE.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/KristianHolme/RDE.jl/actions/workflows/CI.yml?query=branch%3Amain)
+
+This package provides a solver for the rotating detonation engine (RDE) model equations presented in SOURCES:
+
+$u_{t}+ uu_{x} = (1-\lambda)\omega(u)q_0 + \nu_1 u_{xx} + \epsilon \xi (u, u_0)$
+
+$\lambda_t = (1-\lambda)\omega(u) - \beta (u, u_p, s)\lambda + \nu_{2}\lambda_{xx}$.
+
+The solver uses a pseudospectral approach.
+
+This package also provides an interface to the PDE solver using [CommonRLInterface.jl](https://github.com/JuliaReinforcementLearning/CommonRLInterface.jl).
+
+
+# References
+```bibtex
+@article{PhysRevE.101.013106,
+  title = {Mode-locked rotating detonation waves: Experiments and a model equation},
+  author = {Koch, James and Kurosaka, Mitsuru and Knowlen, Carl and Kutz, J. Nathan},
+  journal = {Phys. Rev. E},
+  volume = {101},
+  issue = {1},
+  pages = {013106},
+  numpages = {11},
+  year = {2020},
+  month = {Jan},
+  publisher = {American Physical Society},
+  doi = {10.1103/PhysRevE.101.013106},
+  url = {https://link.aps.org/doi/10.1103/PhysRevE.101.013106}
+}
+
+@article{Koch_2021,
+   title={Multiscale physics of rotating detonation waves: Autosolitons and modulational instabilities},
+   volume={104},
+   ISSN={2470-0053},
+   url={http://dx.doi.org/10.1103/PhysRevE.104.024210},
+   DOI={10.1103/physreve.104.024210},
+   number={2},
+   journal={Physical Review E},
+   publisher={American Physical Society (APS)},
+   author={Koch, James and Kurosaka, Mitsuru and Knowlen, Carl and Kutz, J. Nathan},
+   year={2021},
+   month=aug }
+```

examples/example_DRL.jl

@@ -0,0 +1,28 @@
+using POMDPs
+using POMDPTools
+using Crux
+using Flux
+
+mdp = convert(MDP, RDEEnv())
+
+as = POMDPs.actions(mdp)
+S = ContinuousSpace((mdp.env.prob.params.N*2,))
+layer_size = 32
+A() = DiscreteNetwork(Chain(Dense(Crux.dim(S)..., layer_size, relu),
+                            Dense(layer_size, layer_size, relu),
+                            Dense(layer_size, length(as))), as)
+V() = ContinuousNetwork(Chain(Dense(Crux.dim(S)..., layer_size, relu),
+                              Dense(layer_size, layer_size, relu),
+                              Dense(layer_size, 1)))
+
+
+
+𝒮_ppo = PPO(π=ActorCritic(A(), V()), S=S, N=5000, ΔN=200, max_steps=10000)
+@time π_ppo = POMDPs.solve(𝒮_ppo, mdp)
+
+
+
+p = plot_learning(𝒮_ppo)
+
+PolRunData = run_policy(π_ppo, mdp.env;tmax = 26.0)
+animate_policy_data(PolRunData, mdp.env; fname="ppo")

examples/example_animation.jl

@@ -0,0 +1,10 @@
+
+π_const = ConstantRDEPolicy(RDEEnv());
+
+constPolRunData = run_policy(π_const, RDEEnv();tmax=26.0);
+
+animate_policy_data(constPolRunData, RDEEnv();fname="constant_policy")
+
+
+
+

examples/standard_solver.jl

@@ -0,0 +1,9 @@
+# Create an instance of RDEProblem with custom parameters if desired
+params = RDEParam(;N=128, tmax=26.0)
+rde_prob = RDEProblem(params);
+
+# Solve the PDE
+@time out = solve_pde!(rde_prob; progress=true);
+
+# Plot solution
+plot_solution(rde_prob);

src/RDE.jl

@@ -1,5 +1,43 @@
 module RDE
 
 # Write your package code here.
+    using DifferentialEquations
+    using FFTW
+    using LinearAlgebra
+    using LoopVectorization
+    using Makie
+    using Observables
+    using POMDPs
+    using POMDPTools
+    using CommonRLInterface
+    using Interpolations
+    using ProgressMeter
+
+    using PrecompileTools
+
+
+    export RDEParam, RDEProblem, RDEEnv, solve_pde!
+    export ConstantRDEPolicy, run_policy, PolicyRunData, plot_policy, animate_policy, animate_RDE
+
+
+
+    include("solver.jl")
+    include("RLenv.jl")
+    include("plotting.jl")
+    include("animations.jl")
+    include("interactive_control.jl")
+    include("utils.jl")
+
+
+    @compile_workload begin
+        try
+            #simulate tiny case for a short time
+            prob = RDEProblem(RDEParam(;N=32, tmax = 0.01));
+            solve_pde!(prob, progress=true);
+        catch e
+            @warn "Precompilation failure: $e"
+        end
+    end
+
 
 end

src/RLenv.jl

@@ -0,0 +1,189 @@
+using CommonRLInterface
+using Interpolations
+# using IntervalSets
+# using DomainSets
+using POMDPs
+
+mutable struct RDEEnv <: AbstractEnv
+    prob::RDEProblem                  # RDE problem
+    state::Vector{Float64}  
+    observation::Vector{Float64}          
+    dt::Float64                       # time step
+    t::Float64                        # Current time
+    done::Bool                        # Termination flag
+    reward::Float64
+    smax::Float64
+    u_pmax::Float64
+    observation_samples::Int64
+    max_state_val::Float64
+    reward_func::Function
+    action_num::Int64
+
+    function RDEEnv(;
+        dt = 0.1,
+        smax = 5.0,
+        u_pmax = 3.0,
+        observation_samples = 10,
+        params = RDEParam(),
+        reward_func = RDE_reward_energy_balance!,
+        action_num = 10,
+        kwargs...)
+        prob = RDEProblem(params;kwargs...)
+        initial_state = vcat(prob.u0, prob.λ0)
+        init_observation = Vector{Float64}(undef, observation_samples*2)
+        return new(prob, initial_state, init_observation, dt, 0.0, false, 0.0, 
+                   smax, u_pmax, observation_samples, 100.0, reward_func, action_num)
+    end
+end
+
+function interpolate_state(env::RDEEnv)
+    N = env.prob.params.N
+    n = env.observation_samples
+    L = env.prob.params.L
+    dx = n/L
+    xs_sample = LinRange(0.0, L, n+1)[1:end-1]
+    u = env.state[1:N]
+    λ = env.state[N+1:end]
+
+    itp_u = LinearInterpolation(env.prob.x, u)
+    itp_λ = LinearInterpolation(env.prob.x, λ)
+
+    env.observation[1:n] = itp_u(xs_sample)
+    env.observation[n+1:end] = itp_λ(xs_sample)
+    return env.observation
+end
+
+# CommonRLInterface.reward(env::RDEEnv) = env.reward
+CommonRLInterface.state(env::RDEEnv) = env.state
+CommonRLInterface.terminated(env::RDEEnv) = env.done
+CommonRLInterface.observe(env::RDEEnv) = interpolate_state(env)
+
+function CommonRLInterface.actions(env::RDEEnv)
+    single_actions = LinRange(0, 1, env.action_num)
+    collect(Iterators.product(single_actions, single_actions))[:]
+end
+
+#TODO test that this works
+function CommonRLInterface.clone(env::RDEEnv)
+    env2 = deepcopy(env)
+    @debug "env is copied!"
+    return env2
+end
+
+function CommonRLInterface.setstate!(env::RDEEnv, s)
+    env.state = s
+end
+
+# function CommonRLInterface.state_space(env::RDEEnv)
+    # (0.0 .. env.max_state_val)^env.observation_samples × (0.0 .. 1.1)^env.observation_samples 
+# end
+
+function CommonRLInterface.reset!(env::RDEEnv)
+    # if env.t > 0 && env.t < 40
+    #     error("resetting early?")
+    # end
+    env.t = 0
+    set_init_state!(env.prob)
+    env.state = vcat(env.prob.u0, env.prob.λ0)
+    env.reward = 0.0
+    env.done = false
+    nothing
+end
+
+function CommonRLInterface.act!(env::RDEEnv, action)
+    t_span = (env.t, env.t + env.dt)
+
+    env.prob.params.s = action[1]*env.smax #actions between 0 and 1
+    env.prob.params.u_p = action[2]*env.u_pmax
+
+    prob_ode = ODEProblem(RDE_RHS!, env.state, t_span, env.prob)
+
+    sol = DifferentialEquations.solve(prob_ode)
+    env.prob.sol = sol
+    env.t = sol.t[end]
+    env.state = sol.u[end]
+
+    if env.t ≥ env.prob.params.tmax || sol.retcode != :Success || maximum(abs.(env.state)) > 100
+        env.done = true
+    end
+    env.reward_func(env)
+    return env.reward
+end
+
+function RDE_reward_max!(env::RDEEnv) #just to have something
+    prob = env.prob
+    N = prob.params.N
+    u = env.state[1:N]
+    env.reward = maximum(u)
+    nothing
+end
+
+function RDE_reward_energy_balance!(env::RDEEnv)
+    env.reward = -1*energy_balance(env.state, env.prob.params)
+    nothing
+end
+
+function run_policy(π::P, env::RDEEnv; sparse_skip=5, tmax=26.0, T=Float64) where P <: Policy
+    reset!(env)
+    env.prob.params.tmax = tmax
+    dt = env.dt
+    N = ceil(tmax/dt)+1 |> Int
+
+    ts = Vector{T}(undef,N)
+    ss = Vector{T}(undef,N)
+    u_ps = Vector{T}(undef,N)
+    energy_bal = Vector{T}(undef,N)
+    sparse_N = Int(ceil(N/sparse_skip))
+    sparse_logged = 0
+    sparse_states = Vector{Vector{T}}(undef, sparse_N)
+    sparse_ts = Vector{T}(undef, sparse_N)
+
+    function log!(step)
+        ts[step] = env.t
+        ss[step] = env.prob.params.s
+        u_ps[step] = env.prob.params.u_p
+        energy_bal[step] = energy_balance(state(env), env.prob.params)
+        if (step-1)%sparse_skip == 0
+            sparse_step = (step-1) ÷ sparse_skip + 1
+            sparse_states[sparse_step] = state(env)
+            sparse_ts[sparse_step] = env.t
+            sparse_logged += 1
+        elseif step == N
+            sparse_logged += 1
+            sparse_states[end] = state(env)
+            sparse_ts[end] = env.t
+        end
+    end
+
+    for step = 1:N
+        log!(step)
+        action = POMDPs.action(π, state(env))[1]
+        act!(env, action)
+    end
+
+    if sparse_logged != sparse_N
+        @warn "sparse logged $(sparse_logged) not equal to sparse N $(sparse_N)"
+        pop!(sparse_ts)
+        pop!(sparse_states)
+    end
+
+    return PolicyRunData(ts, ss, u_ps, energy_bal, sparse_ts, sparse_states)
+end
+
+struct PolicyRunData
+    ts
+    ss
+    u_ps
+    energy_bal
+    sparse_ts
+    sparse_states
+end
+
+struct ConstantRDEPolicy <: Policy
+    env::RDEEnv
+    ConstantRDEPolicy(env::RDEEnv = RDEEnv()) = new(env)
+end
+
+function POMDPs.action(π::ConstantRDEPolicy, s)
+    return [(π.env.prob.params.s/π.env.smax, π.env.prob.params.u_p/π.env.u_pmax)]
+end