YingboMa

set-option -g history-limit 65536
set -g -a terminal-overrides ',*:Ss=\E[%p1%d q:Se=\E[2 q'
set -g mouse on

#bind -n C-k clear-history

# default shell
set-option -g default-shell /usr/bin/fish

# clipboard settings

YingboMa

using OrdinaryDiffEq, LabelledArrays

N = 40  # Number of heated units

Cu = N == 1 ? [2e7] : (ones(N) .+ range(0,1.348,length=N))*1e7 # "Heat capacity of heated units";
Cd = 2e6*N # "Heat capacity of distribution circuit";
Gh = 200 # "Thermal conductance of heating elements";
Gu = 150 # "Thermal conductance of heated units to the atmosphere";
Qmax = N*3000 # "Maximum power output of heat generation unit";
Teps = 0.5 # "Threshold of heated unit temperature controllers";

YingboMa

import Base.Broadcast: _broadcast_getindex, preprocess, preprocess_args, Broadcasted, broadcast_unalias, combine_axes, broadcast_shape, check_broadcast_axes, check_broadcast_shape
import Base: copyto!, tail, axes
struct DiffEqBC{T}
    x::T
end
@inline axes(b::DiffEqBC) = axes(b.x)
Base.@propagate_inbounds _broadcast_getindex(b::DiffEqBC, i) = _broadcast_getindex(b.x, i)
Base.@propagate_inbounds _broadcast_getindex(b::DiffEqBC{<:AbstractArray{<:Any,0}}, i) = b.x[]
Base.@propagate_inbounds _broadcast_getindex(b::DiffEqBC{<:AbstractVector}, i) = b.x[i[1]]
Base.@propagate_inbounds _broadcast_getindex(b::DiffEqBC{<:AbstractArray}, i) = b.x[i]

YingboMa

using OrdinaryDiffEq, DiffEqDevTools, Plots, ParameterizedFunctions, Sundials, ODEInterfaceDiffEq
using LinearAlgebra
BLAS.set_num_threads(1)
postfix = "PR"
setups = [
          Dict(:alg=>Kvaerno5()),
          Dict(:alg=>KenCarp4()),
          Dict(:alg=>radau5()),
          Dict(:alg=>Rodas4()),
          Dict(:alg=>Rodas5()),

YingboMa

using LinearAlgebra
lurec(A, blocksize=16) = lurec!(copy(A), Vector{LinearAlgebra.BlasInt}(undef, min(size(A)...)), blocksize)
function lurec!(A::AbstractMatrix{T}, ipiv, blocksize) where T
    info = Ref(zero(LinearAlgebra.BlasInt))
    m, n = size(A)
    mnmin = min(m, n)
    reckernel!(A, m, mnmin, ipiv, info, blocksize)
    LU{T, typeof(A)}(A, ipiv, info[])
end

YingboMa

using ReverseDiff, ForwardDiff
using Test

"""
    vecjac(f, x, v) -> u

``v'J(f(x))``
"""
function vecjac(f, x, v)
    tp = ReverseDiff.InstructionTape()

YingboMa

using OrdinaryDiffEq, Plots
plotly()
import GR: meshgrid

const N = 64
const xd, yd = linspace(0,1,N), linspace(0,1,N)
function brusselator_loop(du, u, p, t)
    @inbounds begin
        A, B, α = p
        # Interior

YingboMa

What did I do

My original goal that I proposed to my mentor Chris was solving boundary value problems (a.k.a. BVPs) which were determined from second order ordinary differential equations (a.k.a. ODEs). I started the BVP path, built a shooting method to solve BVPs from initial value problems (a.k.a. IVPs). Then, I built the beginning of a mono-implicit-Runge-Kutta (a.k.a. MIRK) method, but I don't have all of the bells and whistles. Both of these solvers are in the BoundaryValueDiffEq.jl repository. Instead of trying to jump directly to the end point, and talk about how to do every detail in MIRK, I went to explore how those details naturally arise in second order ODEs.

First, there the idea of symplecticity, because the Labatto MIRK (Lobatto IIIA-IIIB) tableaux are actually symplectic integrators. Symplecticity basically is another way to say energy conservation, so symplectic integrators are specialized for solving second order ODEs that are raise

YingboMa

using DifferentialEquations
using Optim

@public_abstract_type AbstractBVProblem{dType,bType,isinplace,F} <: DEProblem

type BVProblem{dType,bType,initType,F} <: AbstractBVProblem{dType,bType,F}
  f::F
  domin::dType
  bc::bType
  init::initType