text:
using Parameters
using DiffEqBase
using ForwardDiff
using Sundials
using ModelingToolkit
abstract type AbstractReactor end
export AbstractReactor

struct Reactor{D,Q} <: AbstractReactor
    domain::D
    ode::ODEProblem
    recommendedsolver::Q
    forwardsensitivities::Bool
end

function Reactor(domain::T,y0::Array{W,1},tspan::Tuple,interfaces::Z=[];p::X=DiffEqBase.NullParameters(),forwardsensitivities=false,forwarddiff=false,modelingtoolkit=false) where {T<:AbstractDomain,W<:Real,Z<:AbstractArray,X}
    dydt(dy::X,y::T,p::V,t::Q) where {X,T,Q,V} = dydtreactor!(dy,y,t,domain,interfaces,p=p)
    jacy!(J::Q2,y::T,p::V,t::Q) where {Q2,T,Q,V} = jacobiany!(J,y,p,t,domain,interfaces,nothing)
    jacyforwarddiff!(J::Q2,y::T,p::V,t::Q) where {Q2,T,Q,V} = jacobianyforwarddiff!(J,y,p,t,domain,interfaces,nothing)
    jacp!(J::Q2,y::T,p::V,t::Q) where {Q2,T,Q,V} = jacobianp!(J,y,p,t,domain,interfaces,nothing)
    jacpforwarddiff!(J::Q2,y::T,p::V,t::Q) where {Q2,T,Q,V} = jacobianpforwarddiff!(J,y,p,t,domain,interfaces,nothing)

    if (forwardsensitivities || !forwarddiff) && domain isa Union{ConstantTPDomain,ConstantVDomain,ConstantPDomain,ParametrizedTPDomain,ParametrizedVDomain,ParametrizedPDomain,ConstantTVDomain,ParametrizedTConstantVDomain,ConstantTAPhiDomain}
        if !forwardsensitivities
            odefcn = ODEFunction(dydt;jac=jacy!,paramjac=jacp!)
        else
            odefcn = ODEFunction(dydt;paramjac=jacp!)
        end
    else
        odefcn = ODEFunction(dydt;jac=jacyforwarddiff!,paramjac=jacpforwarddiff!)
    end
    if forwardsensitivities
        ode = ODEForwardSensitivityProblem(odefcn,y0,tspan,p)
        recsolver = Sundials.CVODE_BDF(linear_solver=:GMRES)
    else
        ode = ODEProblem(odefcn,y0,tspan,p)
        recsolver  = Sundials.CVODE_BDF()
    end
    if modelingtoolkit
        sys = modelingtoolkitize(ode)
        jac = eval(ModelingToolkit.generate_jacobian(sys)[2])
        if (forwardsensitivities || !forwarddiff) && domain isa Union{ConstantTPDomain,ConstantVDomain,ConstantPDomain,ParametrizedTPDomain,ParametrizedVDomain,ParametrizedPDomain,ConstantTVDomain,ParametrizedTConstantVDomain,ConstantTAPhiDomain}
            odefcn = ODEFunction(dydt;jac=jac,paramjac=jacp!)
        else
            odefcn = ODEFunction(dydt;jac=jac,paramjac=jacpforwarddiff!)
        end
        if forwardsensitivities
            ode = ODEForwardSensitivityProblem(odefcn,y0,tspan,p)
        else
            ode = ODEProblem(odefcn,y0,tspan,p)
        end
    end
    return Reactor(domain,ode,recsolver,forwardsensitivities)
end
function Reactor(domains::T,y0s::W,tspan::W2,interfaces::Z=Tuple(),ps::X=DiffEqBase.NullParameters();forwardsensitivities=false,modelingtoolkit=false) where {T<:Tuple,W<:Tuple,Z,X,W2}
    #adjust indexing
    y0 = zeros(sum(length(y) for y in y0s))
    Nvars = 0
    for domain in domains
        Nvars += domain.indexes[end]
    end
    n = Nvars
    k = 1
    for (j,domain) in enumerate(domains)
        Nspcs = length(domain.phase.species)
        Ntherm = length(domain.indexes) - 2
        for i = 1:6, j = 1:size(domain.rxnarray)[2]
            if domain.rxnarray[i,j] != 0
                domain.rxnarray[i,j] += k-1
            end
        end
        for i in 1:length(domain.constantspeciesinds)
            domain.constantspeciesinds[i] += k-1
        end
        for (thermovar,ind) in domain.thermovariabledict
            domain.thermovariabledict[thermovar] += k-1
        end
        domain.indexes[1] = k
        k += Nspcs
        domain.indexes[2] = k-1
        y0[domain.indexes[1]:domain.indexes[2]] = y0s[j][1:Nspcs]
        for m = 3:2+Ntherm
            domain.indexes[m] = n
            y0[n] = y0s[j][Nspcs+m-2]
            n -= 1
        end
    end

    p = Array{Float64,1}()
    phases = []
    phaseinds = Array{Tuple,1}()
    for (i,domain) in enumerate(domains)
        if domain.phase in phases
            ind = findfirst(phases.==domain.phase)
            domain.parameterindexes[1] = phaseinds[ind][1]
            domain.parameterindexes[2] = phaseinds[ind][2]
            ds = [domains[ind] for ind in findall(x->x.phase==domain.phase,domains)]
            if any(.!isa.(ds,AbstractConstantKDomain)) && any(isa.(ds,AbstractConstantKDomain)) #handle different p formats for the same phase
                for d in ds
                    if !isa(d,AbstractConstantKDomain)
                        p[domain.parameterindexes[1]:domain.parameterindexes[2]] .= d.p
                        break
                    end
                end
                for d in ds
                    if isa(d,AbstractConstantKDomain)
                        if forwardsensitivities == true #error if running into sensitivity bug
                            error("forward sensitivities is not supported for domain combinations that share rate coefficients, but treat them differently ex: ConstantTPDomain and ConstantVDomain")
                        end
                        d.alternativepformat = true
                    end
                end
            end
        else
            domain.parameterindexes[1] = length(p)+1
            domain.parameterindexes[2] = length(p)+length(ps[i])
            push!(phaseinds,(length(p)+1,length(p)+length(ps[i])))
            push!(phases,domain.phase)
            p = vcat(p,ps[i])
        end
    end

    for (i,inter) in enumerate(interfaces)
        if isa(inter, AbstractReactiveInternalInterface)
            ind1 = findfirst(isequal(inter.domain1),domains)
            ind2 = findfirst(isequal(inter.domain2),domains)
            inter.domaininds[1] = ind1
            inter.domaininds[2] = ind2
            inter.parameterindexes[1] = length(p)+1
            inter.parameterindexes[2] = length(p)+length(ps[i+length(domains)])
            inter.rxnarray .= getinterfacereactioninds(inter.domain1,inter.domain2,inter.reactions)
            p = vcat(p,ps[i+length(domains)])
        end
    end

    dydt(dy::X,y::T,p::V,t::Q) where {X,T,Q,V} = dydtreactor!(dy,y,t,domains,interfaces,p=p)
    jacy!(J::Q2,y::T,p::V,t::Q) where {Q2,T,Q,V} = jacobianyforwarddiff!(J,y,p,t,domains,interfaces,nothing)
    jacp!(J::Q2,y::T,p::V,t::Q) where {Q2,T,Q,V} = jacobianpforwarddiff!(J,y,p,t,domains,interfaces,nothing)

    if forwardsensitivities
        odefcn = ODEFunction(dydt;paramjac=jacp!)
        ode = ODEForwardSensitivityProblem(odefcn,y0,tspan,p)
        recsolver = Sundials.CVODE_BDF(linear_solver=:GMRES)
        if modelingtoolkit
            sys = modelingtoolkitize(ode)
            jac = eval(ModelingToolkit.generate_jacobian(sys)[2])
            odefcn = ODEFunction(dydt;jac=jac,paramjac=jacp!)
            ode = ODEForwardSensitivityProblem(odefcn,y0,tspan,p)
        end
    else
        odefcn = ODEFunction(dydt;jac=jacy!,paramjac=jacp!)
        ode = ODEProblem(odefcn,y0,tspan,p)
        recsolver  = Sundials.CVODE_BDF()
        if modelingtoolkit
            sys = modelingtoolkitize(ode)
            jac = eval(ModelingToolkit.generate_jacobian(sys)[2])
            odefcn = ODEFunction(dydt;jac=jac,paramjac=jacp!)
            ode = ODEProblem(odefcn,y0,tspan,p)
        end
    end
    return Reactor(domains,ode,recsolver,forwardsensitivities),y0,p
end
export Reactor

@inline function getrate(rxn::T,cs::Array{W,1},kfs::Array{Q,1},krevs::Array{Q,1}) where {T<:AbstractReaction,Q,W<:Real}
    Nreact = length(rxn.reactantinds)
    Nprod = length(rxn.productinds)
    R = 0.0
    if Nreact == 1
        @fastmath @inbounds R += kfs[rxn.index]*cs[rxn.reactantinds[1]]
    elseif Nreact == 2
        @fastmath @inbounds R += kfs[rxn.index]*cs[rxn.reactantinds[1]]*cs[rxn.reactantinds[2]]
    elseif Nreact == 3
        @fastmath @inbounds R += kfs[rxn.index]*cs[rxn.reactantinds[1]]*cs[rxn.reactantinds[2]]*cs[rxn.reactantinds[3]]
    end

    if Nprod == 1
        @fastmath @inbounds R -= krevs[rxn.index]*cs[rxn.productinds[1]]
    elseif Nprod == 2
        @fastmath @inbounds R -= krevs[rxn.index]*cs[rxn.productinds[1]]*cs[rxn.productinds[2]]
    elseif Nprod == 3
        @fastmath @inbounds R -= krevs[rxn.index]*cs[rxn.productinds[1]]*cs[rxn.productinds[2]]*cs[rxn.productinds[3]]
    end

    return R
end
export getrate

@inline function getrates(rarray,cs,kfs,krevs)
    rts = zeros(length(kfs))
    for i = 1:length(rts)
        if @inbounds rarray[2,i] == 0
            @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]
        elseif @inbounds rarray[3,i] == 0
            @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]
        else
            @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]]
        end
        if @inbounds rarray[5,i] == 0
            @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]
        elseif @inbounds rarray[6,i] == 0
            @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]
        else
            @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]*cs[rarray[6,i]]
        end
        @fastmath R = fR - rR

        rts[i] = R
    end
    return rts
end
export getrates

@inline function addreactionratecontributions!(dydt::Q,rarray::Array{W2,2},cs::W,kfs::Z,krevs::Y) where {Q,Z,Y,T,W,W2}
    @inbounds @simd for i = 1:size(rarray)[2]
        if @inbounds rarray[2,i] == 0
            @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]
        elseif @inbounds rarray[3,i] == 0
            @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]
        else
            @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]]
        end
        if @inbounds rarray[5,i] == 0
            @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]
        elseif @inbounds rarray[6,i] == 0
            @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]
        else
            @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]*cs[rarray[6,i]]
        end
        @fastmath R = fR - rR
        @inbounds @fastmath dydt[rarray[1,i]] -= R
        if @inbounds rarray[2,i] != 0
            @inbounds @fastmath dydt[rarray[2,i]] -= R
            if @inbounds rarray[3,i] != 0
                @inbounds @fastmath dydt[rarray[3,i]] -= R
            end
        end
        @inbounds @fastmath dydt[rarray[4,i]] += R
        if @inbounds rarray[5,i] != 0
            @inbounds @fastmath dydt[rarray[5,i]] += R
            if @inbounds rarray[6,i] != 0
                @inbounds @fastmath dydt[rarray[6,i]] += R
            end
        end
    end
end

@inline function addreactionratecontributions!(dydt::Q,rarray::Array{W2,2},cs::W,kfs::Z,krevs::Y,V) where {Q,Z,Y,T,W,W2}
    @inbounds @simd for i = 1:size(rarray)[2]
        if @inbounds rarray[2,i] == 0
            @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]
        elseif @inbounds rarray[3,i] == 0
            @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]
        else
            @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]]
        end
        if @inbounds rarray[5,i] == 0
            @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]
        elseif @inbounds rarray[6,i] == 0
            @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]
        else
            @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]*cs[rarray[6,i]]
        end
        @fastmath R = (fR - rR)*V
        @inbounds @fastmath dydt[rarray[1,i]] -= R
        if @inbounds rarray[2,i] != 0
            @inbounds @fastmath dydt[rarray[2,i]] -= R
            if @inbounds rarray[3,i] != 0
                @inbounds @fastmath dydt[rarray[3,i]] -= R
            end
        end
        @inbounds @fastmath dydt[rarray[4,i]] += R
        if @inbounds rarray[5,i] != 0
            @inbounds @fastmath dydt[rarray[5,i]] += R
            if @inbounds rarray[6,i] != 0
                @inbounds @fastmath dydt[rarray[6,i]] += R
            end
        end
    end
end
export addreactionratecontributions!

@inline function addreactionratecontributionsforwardreverse!(dydt::Q,rarray::Array{W2,2},cs::W,kfs::Z,krevs::Y,V) where {Q,Z,Y,T,W,W2}
    frts = zeros(length(kfs))
    rrts = zeros(length(kfs))
    rts = zeros(length(kfs))
    @inbounds for i = 1:size(rarray)[2]
        if @inbounds rarray[2,i] == 0
            @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]
        elseif @inbounds rarray[3,i] == 0
            @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]
        else
            @inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]]
        end
        if @inbounds rarray[5,i] == 0
            @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]
        elseif @inbounds rarray[6,i] == 0
            @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]
        else
            @inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]*cs[rarray[6,i]]
        end
        @inbounds @fastmath frts[i] = fR*V
        @inbounds @fastmath rrts[i] = rR*V
        @inbounds @fastmath R = frts[i] - rrts[i]
        @inbounds rts[i] = R
        @inbounds @fastmath dydt[rarray[1,i]] -= R
        if @inbounds rarray[2,i] != 0
            @inbounds @fastmath dydt[rarray[2,i]] -= R
            if @inbounds rarray[3,i] != 0
                @inbounds @fastmath dydt[rarray[3,i]] -= R
            end
        end
        @inbounds @fastmath dydt[rarray[4,i]] += R
        if @inbounds rarray[5,i] != 0
            @inbounds @fastmath dydt[rarray[5,i]] += R
            if @inbounds rarray[6,i] != 0
                @inbounds @fastmath dydt[rarray[6,i]] += R
            end
        end
    end
    return rts,frts,rrts
end
export addreactionratecontributionsforwardreverse!

@inline function dydtreactor!(dydt::RC,y::U,t::Z,domain::Q,interfaces::B;p::RV=DiffEqBase.NullParameters(),sensitivity::Bool=true) where {RC,RV,B,Z,U,Q<:AbstractDomain}
    dydt .= 0.0
    ns,cs,T,P,V,C,N,mu,kfs,krevs,Hs,Us,Gs,diffs,Cvave,cpdivR,phi = calcthermo(domain,y,t,p)
    addreactionratecontributions!(dydt,domain.rxnarray,cs,kfs,krevs)
    dydt .*= V
    calcdomainderivatives!(domain,dydt,interfaces;t=t,T=T,P=P,Us=Us,Hs=Hs,V=V,C=C,ns=ns,N=N,Cvave=Cvave)
    return dydt
end
@inline function dydtreactor!(dydt::RC,y::U,t::Z,domains::Q,interfaces::B;p::RV=DiffEqBase.NullParameters(),sensitivity::Bool=true) where {RC,RV,B,Z,U,Q<:Tuple}
    cstot = similar(y)
    cstot .= 0.0
    dydt .= 0.0
    domain = domains[1]
    ns,cs,T,P,V,C,N,mu,kfs,krevs,Hs,Us,Gs,diffs,Cvave,cpdivR,phi = calcthermo(domain,y,t,p)
    vns = Array{Any,1}(undef,length(domains))
    vns[1] = ns
    vcs = Array{Any,1}(undef,length(domains))
    vcs[1] = cs
    cstot[domain.indexes[1]:domain.indexes[2]] = cs
    vT = Array{Any,1}(undef,length(domains))
    vT[1] = T
    vP = Array{Any,1}(undef,length(domains))
    vP[1] = P
    vV = Array{Any,1}(undef,length(domains))
    vV[1] = V
    vC = Array{Any,1}(undef,length(domains))
    vC[1] = C
    vN = Array{Any,1}(undef,length(domains))
    vN[1] = N
    vmu = Array{Any,1}(undef,length(domains))
    vmu[1] = mu
    vkfs = Array{Any,1}(undef,length(domains))
    vkfs[1] = kfs
    vkrevs = Array{Any,1}(undef,length(domains))
    vkrevs[1] = krevs
    vHs = Array{Any,1}(undef,length(domains))
    vHs[1] = Hs
    vUs = Array{Any,1}(undef,length(domains))
    vUs[1] = Us
    vGs = Array{Any,1}(undef,length(domains))
    vGs[1] = Gs
    vdiffs = Array{Any,1}(undef,length(domains))
    vdiffs[1] = diffs
    vCvave = Array{Any,1}(undef,length(domains))
    vCvave[1] = Cvave
    vcpdivR = Array{Any,1}(undef,length(domains))
    vcpdivR[1] = cpdivR
    vphi = Array{Any,1}(undef,length(domains))
    vphi[1] = phi
    addreactionratecontributions!(dydt,domain.rxnarray,cstot,kfs,krevs)
    @views dydt[domain.indexes[1]:domain.indexes[2]] .*= V
    for (i,domain) in enumerate(@views domains[2:end])
        k = i + 1
        vns[k],vcs[k],vT[k],vP[k],vV[k],vC[k],vN[k],vmu[k],vkfs[k],vkrevs[k],vHs[k],vUs[k],vGs[k],vdiffs[k],vCvave[k],vcpdivR[k],vphi[k] = calcthermo(domain,y,t,p)
        cstot[domain.indexes[1]:domain.indexes[2]] .= vcs[k]
        addreactionratecontributions!(dydt,domain.rxnarray,cstot,vkfs[k],vkrevs[k])
        @views dydt[domain.indexes[1]:domain.indexes[2]] .*= vV[k]
    end
    for (i,inter) in enumerate(interfaces)
        if isa(inter,AbstractReactiveInternalInterface)
            evaluate(inter,dydt,domains,vT[inter.domaininds[1]],vT[inter.domaininds[2]],vphi[inter.domaininds[1]],vphi[inter.domaininds[2]],vGs[inter.domaininds[1]],vGs[inter.domaininds[2]],cstot,p)
        end
    end
    for (i,domain) in enumerate(domains)
        calcdomainderivatives!(domain,dydt,interfaces;t=t,T=vT[i],P=vP[i],Us=vUs[i],Hs=vHs[i],V=vV[i],C=vC[i],ns=vns[i],N=vN[i],Cvave=vCvave[i])
    end
    return dydt
end
export dydtreactor!

function jacobianyforwarddiff!(J::Q,y::U,p::W,t::Z,domain::V,interfaces::Q3,colorvec::Q2=nothing) where {Q3,Q2,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,V<:AbstractDomain}
    f(dy::X,y::Array{T,1}) where {T<:Real,X} = dydtreactor!(dy,y,t,domain,interfaces;p=p,sensitivity=false)
    ForwardDiff.jacobian!(J,f,zeros(size(y)),y)
end
export jacobianyforwarddiff!

function jacobianyforwarddiff(y::U,p::W,t::Z,domain::V,interfaces::Q3,colorvec::Q2=nothing) where {Q3,Q2,U<:AbstractArray,W,Z<:Real,V<:AbstractDomain}
    J = zeros(length(y),length(y))
    jacobianyforwarddiff!(J,y,p,t,domain,interfaces,colorvec)
    return J
end
export jacobianyforwarddiff

function jacobianyforwarddiff!(J::Q,y::U,p::W,t::Z,domains::V,interfaces::Q3,colorvec::Q2=nothing) where {Q3,Q2,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,V<:Tuple}
    f(dy::X,y::Array{T,1}) where {T<:Real,X} = dydtreactor!(dy,y,t,domains,interfaces;p=p,sensitivity=false)
    ForwardDiff.jacobian!(J,f,zeros(size(y)),y)
end

function jacobianyforwarddiff(y::U,p::W,t::Z,domains::V,interfaces::Q3,colorvec::Q2=nothing) where {Q3,Q2,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,V<:Tuple}
    J = zeros(length(y),length(y))
    jacobianyforwarddiff!(J,y,p,t,domains,interfaces,colorvec)
    return J
end
# function jacobiany!(J::Q,y::U,p::W,t::Z,domain::V,interfaces::Q3,colorvec::Q2=nothing) where {Q3<:AbstractArray,Q2<:AbstractArray,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,V<:AbstractDomain}
#     f(y::Array{T,1}) where {T<:Real} = dydtreactor!(y,t,domain,interfaces;p=p,sensitivity=false)
#     forwarddiff_color_jacobian!(J,f,y,colorvec=colorvec)
# end
# function jacobiany!(J::Q,y::U,p::W,t::Z,domain::Q4,interfaces::Q3,colorvec::Q2=nothing) where {Q3<:AbstractArray,Q2,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,Q4<:Union{}}
#     ns,cs,T,P,V,C,N,mu,kfs,krevs,Hs,Us,Gs,diffs,Cvave = calcthermo(domain,y,t,p)
#     jacobiany!(y,t,domain,kfs,krevs,J;zero=true)
# end
# export jacobiany!

function jacobianpforwarddiff!(J::Q,y::U,p::W,t::Z,domain::V,interfaces::Q3,colorvec::Q2=nothing) where {Q3,Q2,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,V<:AbstractDomain}
    function f(dy::X,p::Array{T,1}) where {X,T<:Real}
        dydtreactor!(dy,y,t,domain,interfaces;p=p,sensitivity=false)
    end
    dy = zeros(length(y))
    ForwardDiff.jacobian!(J,f,dy,p)
end

function jacobianpforwarddiff!(J::Q,y::U,p::W,t::Z,domains::V,interfaces::Q3,colorvec::Q2=nothing) where {Q3,Q2,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,V<:Tuple}
    function f(dy::X,p::Array{T,1}) where {X,T<:Real}
        dydtreactor!(dy,y,t,domains,interfaces;p=p,sensitivity=false)
    end
    dy = zeros(length(y))
    ForwardDiff.jacobian!(J,f,dy,p)
end

export jacobianpforwarddiff!

function jacobianpforwarddiff(y::U,p::W,t::Z,domain::V,interfaces::Q3,colorvec::Q2=nothing) where {Q3,Q2,U<:AbstractArray,W,Z<:Real,V<:AbstractDomain}
    J = zeros(length(y),length(domain.phase.species)+length(domain.phase.reactions))
    jacobianpforwarddiff!(J,y,p,t,domain,interfaces,colorvec)
end

function jacobianpforwarddiff(y::U,p::W,t::Z,domains::V,interfaces::Q3,colorvec::Q2=nothing) where {Q3,Q2,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,V<:Tuple}
    J = zeros(length(y),length(domains.phase.species)+length(domains.phase.reactions))
    jacobianpforwarddiff!(J,y,p,t,domains,interfaces,colorvec)
end

export jacobianpforwarddiff

# function jacobianp!(J::Q,y::U,p::W,t::Z,domain::V,interfaces::Q3,colorvec::Q2=nothing) where {Q3<:AbstractArray,Q2<:AbstractArray,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,V<:AbstractDomain}
#     f(p::Array{T,1}) where {T<:Real} = dydtreactor!(y,domain.t[1],domain,interfaces;p=p,sensitivity=false)
#     forwarddiff_color_jacobian!(J,f,p,colorvec=colorvec)
# end
# function jacobianp!(J::Q,y::U,p::W,t::Z,domain::Q4,interfaces::Q3,colorvec::Q2=nothing) where {Q3<:AbstractArray,Q2,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,Q4<:Union{ConstantTPDomain,ConstantTVDomain}}
#     ns,cs,T,P,V,C,N,mu,kfs,krevs,Hs,Us,Gs,diffs,Cvave = calcthermo(domain,y,t,p)
#     dydt = zeros(length(y))
#     addreactionratecontributions!(dydt,domain.rxnarray,cs,kfs,krevs)
#     @views wV = dydt[domain.indexes[1]:domain.indexes[2]]
#     jacobianp!(domain;cs=cs,V=V,T=T,Us=Us,Cvave=Cvave,N=N,kfs=kfs,krevs=krevs,wV=wV,ratederiv=J)
# end

function jacobiany(y::U,p::W,t::Z,domain::V,interfaces::Q3,colorvec::Q2=nothing) where {Q3,Q2,U<:AbstractArray,W,Z<:Real,V<:AbstractDomain}
    J = zeros(length(y),length(y))
    jacobiany!(J,y,p,t,domain,interfaces,colorvec)
    return J
end
export jacobiany

function jacobianp(y::U,p::W,t::Z,domain::V,interfaces::Q3,colorvec::Q2=nothing) where {Q3,Q2,U<:AbstractArray,W,Z<:Real,V<:AbstractDomain}
    J = zeros(length(y),length(domain.phase.species)+length(domain.phase.reactions))
    jacobianp!(J,y,p,t,domain,interfaces,colorvec)
    return J
end
export jacobianp

@inline function _spreadreactantpartials!(jac::S,deriv::Float64,rxnarray::Array{Int64,2},rxnind::Int64,ind::Int64) where {S<:AbstractArray}
    @inbounds jac[rxnarray[4,rxnind],ind] += deriv
    if @inbounds rxnarray[5,rxnind] !== 0
        @inbounds jac[rxnarray[5,rxnind],ind] += deriv
        if @inbounds rxnarray[6,rxnind] !== 0
            @inbounds jac[rxnarray[6,rxnind],ind] += deriv
        end
    end
end
@inline function _spreadproductpartials!(jac::S,deriv::Float64,rxnarray::Array{Int64,2},rxnind::Int64,ind::Int64) where {S<:AbstractArray}
    @inbounds jac[rxnarray[1,rxnind],ind] += deriv
    if @inbounds rxnarray[2,rxnind] !== 0
        @inbounds jac[rxnarray[2,rxnind],ind] += deriv
        if @inbounds rxnarray[3,rxnind] !== 0
            @inbounds jac[rxnarray[3,rxnind],ind] += deriv
        end
    end
end

@inline function _jacobianynswrtns!(jac::S,rxnarray::Array{Int64,2},rxnind::Int64,cs::Array{Float64,1},kf::Float64,krev::Float64) where {S<:AbstractArray}
    k=kf
    if rxnarray[2,rxnind] == 0
        deriv = k
        @inbounds jac[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= deriv
        @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[1,rxnind])
    elseif rxnarray[3,rxnind] == 0
        if rxnarray[1,rxnind] == rxnarray[2,rxnind]
            @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[1,rxnind]]
            @inbounds jac[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= 2.0*deriv
            @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[1,rxnind])
        else
            @inbounds @fastmath deriv = k*cs[rxnarray[2,rxnind]]
            @inbounds jac[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= deriv
            @inbounds jac[rxnarray[2,rxnind],rxnarray[1,rxnind]] -= deriv
            @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[1,rxnind])
            @inbounds @fastmath deriv = k*cs[rxnarray[1,rxnind]]
            @inbounds jac[rxnarray[1,rxnind],rxnarray[2,rxnind]] -= deriv
            @inbounds jac[rxnarray[2,rxnind],rxnarray[2,rxnind]] -= deriv
            @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[2,rxnind])
        end
    else
        if rxnarray[1,rxnind]==rxnarray[2,rxnind] && rxnarray[1,rxnind]==rxnarray[3,rxnind]
            @inbounds @fastmath deriv = 3.0*k*cs[rxnarray[1,rxnind]]*cs[rxnarray[1,rxnind]]
            @inbounds jac[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= 3.0*deriv
            @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[1,rxnind])
        elseif rxnarray[1,rxnind]==rxnarray[2,rxnind]
            @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[1,rxnind]]*cs[rxnarray[3,rxnind]]
            @inbounds jac[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= 2.0*deriv
            @inbounds jac[rxnarray[3,rxnind],rxnarray[1,rxnind]] -= deriv
            @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[1,rxnind])
            @inbounds @fastmath deriv = k*cs[rxnarray[1,rxnind]]*cs[rxnarray[1,rxnind]]
            @inbounds jac[rxnarray[1,rxnind],rxnarray[3,rxnind]] -= 2.0*deriv
            @inbounds jac[rxnarray[3,rxnind],rxnarray[3,rxnind]] -= deriv
            @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[3,rxnind])
        elseif rxnarray[2,rxnind]==rxnarray[3,rxnind]
            @inbounds @fastmath deriv = k*cs[rxnarray[2,rxnind]]*cs[rxnarray[2,rxnind]]
            @inbounds jac[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= deriv
            @inbounds jac[rxnarray[2,rxnind],rxnarray[1,rxnind]] -= 2.0*deriv
            @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[1,rxnind])
            @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]
            @inbounds jac[rxnarray[1,rxnind],rxnarray[2,rxnind]] -= deriv
            @inbounds jac[rxnarray[2,rxnind],rxnarray[2,rxnind]] -= 2.0*deriv
            @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[2,rxnind])
        elseif rxnarray[1,rxnind]==rxnarray[3,rxnind]
            @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]
            @inbounds jac[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= 2.0*deriv
            @inbounds jac[rxnarray[2,rxnind],rxnarray[1,rxnind]] -= deriv
            @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[1,rxnind])
            @inbounds @fastmath deriv = k*cs[rxnarray[1,rxnind]]*cs[rxnarray[1,rxnind]]
            @inbounds jac[rxnarray[1,rxnind],rxnarray[2,rxnind]] -= 2.0*deriv
            @inbounds jac[rxnarray[2,rxnind],rxnarray[2,rxnind]] -= deriv
            @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[2,rxnind])
        else
            @inbounds @fastmath deriv = k*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]]
            @inbounds jac[rxnarray[1,rxnind],rxnarray[1,rxnind]] -= deriv
            @inbounds jac[rxnarray[2,rxnind],rxnarray[1,rxnind]] -= deriv
            @inbounds jac[rxnarray[3,rxnind],rxnarray[1,rxnind]] -= deriv
            @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[1,rxnind])
            @inbounds @fastmath deriv = k*cs[rxnarray[1,rxnind]]*cs[rxnarray[3,rxnind]]
            @inbounds jac[rxnarray[1,rxnind],rxnarray[2,rxnind]] -= deriv
            @inbounds jac[rxnarray[2,rxnind],rxnarray[2,rxnind]] -= deriv
            @inbounds jac[rxnarray[3,rxnind],rxnarray[2,rxnind]] -= deriv
            @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[2,rxnind])
            @inbounds @fastmath deriv = k*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]
            @inbounds jac[rxnarray[1,rxnind],rxnarray[3,rxnind]] -= deriv
            @inbounds jac[rxnarray[2,rxnind],rxnarray[3,rxnind]] -= deriv
            @inbounds jac[rxnarray[3,rxnind],rxnarray[3,rxnind]] -= deriv
            @inbounds _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,rxnarray[3,rxnind])
        end
    end
    k=krev
    if rxnarray[5,rxnind] == 0
        deriv = k
        @inbounds jac[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= deriv
        @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[4,rxnind])
    elseif rxnarray[6,rxnind] == 0
        if rxnarray[4,rxnind] == rxnarray[5,rxnind]
            @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[4,rxnind]]
            @inbounds jac[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= 2.0*deriv
            @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[4,rxnind])
        else
            @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]
            @inbounds jac[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= deriv
            @inbounds jac[rxnarray[5,rxnind],rxnarray[4,rxnind]] -= deriv
            @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[4,rxnind])
            @inbounds @fastmath deriv = k*cs[rxnarray[4,rxnind]]
            @inbounds jac[rxnarray[4,rxnind],rxnarray[5,rxnind]] -= deriv
            @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= deriv
            @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind])
        end
    else
        if rxnarray[4,rxnind]==rxnarray[5,rxnind] && rxnarray[4,rxnind]==rxnarray[6,rxnind]
            @inbounds @fastmath deriv = 3.0*k*cs[rxnarray[4,rxnind]]*cs[rxnarray[4,rxnind]]
            @inbounds jac[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= 3.0*deriv
            @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[4,rxnind])
        elseif rxnarray[4,rxnind]==rxnarray[5,rxnind]
            @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[4,rxnind]]*cs[rxnarray[6,rxnind]]
            @inbounds jac[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= 2.0*deriv
            @inbounds jac[rxnarray[6,rxnind],rxnarray[4,rxnind]] -= deriv
            @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[4,rxnind])
            @inbounds @fastmath deriv = k*cs[rxnarray[4,rxnind]]*cs[rxnarray[4,rxnind]]
            @inbounds jac[rxnarray[4,rxnind],rxnarray[6,rxnind]] -= 2.0*deriv
            @inbounds jac[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= deriv
            @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[6,rxnind])
        elseif rxnarray[5,rxnind]==rxnarray[6,rxnind]
            @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[5,rxnind]]
            @inbounds jac[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= deriv
            @inbounds jac[rxnarray[5,rxnind],rxnarray[4,rxnind]] -= 2.0*deriv
            @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[4,rxnind])
            @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]]
            @inbounds jac[rxnarray[4,rxnind],rxnarray[5,rxnind]] -= deriv
            @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= 2.0*deriv
            @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind])
        elseif rxnarray[4,rxnind]==rxnarray[6,rxnind]
            @inbounds @fastmath deriv = 2.0*k*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]]
            @inbounds jac[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= 2.0*deriv
            @inbounds jac[rxnarray[5,rxnind],rxnarray[4,rxnind]] -= deriv
            @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[4,rxnind])
            @inbounds @fastmath deriv = k*cs[rxnarray[4,rxnind]]*cs[rxnarray[4,rxnind]]
            @inbounds jac[rxnarray[4,rxnind],rxnarray[5,rxnind]] -= 2.0*deriv
            @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= deriv
            @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind])
        else
            @inbounds @fastmath deriv = k*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]
            @inbounds jac[rxnarray[4,rxnind],rxnarray[4,rxnind]] -= deriv
            @inbounds jac[rxnarray[5,rxnind],rxnarray[4,rxnind]] -= deriv
            @inbounds jac[rxnarray[6,rxnind],rxnarray[4,rxnind]] -= deriv
            @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[4,rxnind])
            @inbounds @fastmath deriv = k*cs[rxnarray[4,rxnind]]*cs[rxnarray[6,rxnind]]
            @inbounds jac[rxnarray[4,rxnind],rxnarray[5,rxnind]] -= deriv
            @inbounds jac[rxnarray[5,rxnind],rxnarray[5,rxnind]] -= deriv
            @inbounds jac[rxnarray[6,rxnind],rxnarray[5,rxnind]] -= deriv
            @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[5,rxnind])
            @inbounds @fastmath deriv = k*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]]
            @inbounds jac[rxnarray[4,rxnind],rxnarray[6,rxnind]] -= deriv
            @inbounds jac[rxnarray[5,rxnind],rxnarray[6,rxnind]] -= deriv
            @inbounds jac[rxnarray[6,rxnind],rxnarray[6,rxnind]] -= deriv
            @inbounds _spreadproductpartials!(jac,deriv,rxnarray,rxnind,rxnarray[6,rxnind])
        end
    end
end

@inline function _jacobianynswrtV!(jac::S,Vind::Int64,rxnarray::Array{Int64,2},rxnind::Int64,cs::Array{Float64,1},kf::Float64,krev::Float64) where {S<:AbstractArray}
    k=kf
    if rxnarray[2,rxnind]==0
        nothing
    elseif rxnarray[3,rxnind]== 0
        @inbounds @fastmath deriv = -k*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]
        @inbounds jac[rxnarray[1,rxnind],Vind] -= deriv
        @inbounds jac[rxnarray[2,rxnind],Vind] -= deriv
        _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,Vind)
    else
        @inbounds @fastmath deriv = -2.0*k*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]]
        @inbounds jac[rxnarray[1,rxnind],Vind] -= deriv
        @inbounds jac[rxnarray[2,rxnind],Vind] -= deriv
        @inbounds jac[rxnarray[3,rxnind],Vind] -= deriv
        _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,Vind)
    end
    k=krev
    if rxnarray[5,rxnind]==0
        nothing
    elseif rxnarray[6,rxnind] == 0
        @inbounds @fastmath deriv = -k*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]]
        @inbounds jac[rxnarray[4,rxnind],Vind] -= deriv
        @inbounds jac[rxnarray[5,rxnind],Vind] -= deriv
        _spreadproductpartials!(jac,deriv,rxnarray,rxnind,Vind)
    else
        @inbounds @fastmath deriv = -2.0*k*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]
        @inbounds jac[rxnarray[4,rxnind],Vind] -= deriv
        @inbounds jac[rxnarray[5,rxnind],Vind] -= deriv
        @inbounds jac[rxnarray[6,rxnind],Vind] -= deriv
        _spreadproductpartials!(jac,deriv,rxnarray,rxnind,Vind)
    end
end

"""
This function calculates the ns partials in jacobiany involving k derivatives. dkdx is either dkdni and dkdV. x is either ni or V.
"""
@inline function _jacobianykderiv!(jac::S,xind::Int64,dkfdx::Float64,dkrevdx::Float64,rxnarray::Array{Int64,2},rxnind::Int64,cs::Array{Float64,1},V::Float64) where {S<:AbstractArray}
    dkdx = dkfdx
    if rxnarray[2,rxnind] == 0
        @inbounds @fastmath deriv = dkdx*cs[rxnarray[1,rxnind]]*V
        @inbounds jac[rxnarray[1,rxnind],xind] -= deriv
        _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,xind)
    elseif rxnarray[3,rxnind] == 0
        @inbounds @fastmath deriv = dkdx*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*V
        @inbounds jac[rxnarray[1,rxnind],xind] -= deriv
        @inbounds jac[rxnarray[2,rxnind],xind] -= deriv
        _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,xind)
    else
        @inbounds @fastmath deriv = dkdx*cs[rxnarray[1,rxnind]]*cs[rxnarray[2,rxnind]]*cs[rxnarray[3,rxnind]]*V
        @inbounds jac[rxnarray[1,rxnind],xind] -= deriv
        @inbounds jac[rxnarray[2,rxnind],xind] -= deriv
        @inbounds jac[rxnarray[3,rxnind],xind] -= deriv
        _spreadreactantpartials!(jac,deriv,rxnarray,rxnind,xind)
    end
    dkdx = dkrevdx
    if rxnarray[5,rxnind] == 0
        @inbounds @fastmath deriv = dkdx*cs[rxnarray[4,rxnind]]*V
        @inbounds jac[rxnarray[4,rxnind],xind] -= deriv
        _spreadproductpartials!(jac,deriv,rxnarray,rxnind,xind)
    elseif rxnarray[6,rxnind] == 0
        @inbounds @fastmath deriv = dkdx*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]]*V
        @inbounds jac[rxnarray[4,rxnind],xind] -= deriv
        @inbounds jac[rxnarray[5,rxnind],xind] -= deriv
        _spreadproductpartials!(jac,deriv,rxnarray,rxnind,xind)
    else
        @inbounds @fastmath deriv = dkdx*cs[rxnarray[4,rxnind]]*cs[rxnarray[5,rxnind]]*cs[rxnarray[6,rxnind]]*V
        @inbounds jac[rxnarray[4,rxnind],xind] -= deriv
        @inbounds jac[rxnarray[5,rxnind],xind] -= deriv
        @inbounds jac[rxnarray[6,rxnind],xind] -= deriv
        _spreadproductpartials!(jac,deriv,rxnarray,rxnind,xind)
    end
end

function _dydttherm(dy::X,x::T,y::Q,p::W,t::Z,domain::D,interfaces::Q3,ind::T1) where {X,Q3<:AbstractArray,T<:Real,Q<:AbstractArray,Z<:Real,D<:AbstractDomain,T1<:Integer,W}
    v = [ i != ind ? convert(typeof(x),z) : x for (i,z) in enumerate(y)]
    return dydtreactor!(dy,v,t,domain,interfaces;p=p,sensitivity=false)
end

function jacobianytherm!(jac::Q,y::U,p::W,t::Z,domain::D,interfaces::Q3,ind::I,x::F,colorvec::Q2=nothing) where {Q3<:AbstractArray,Q2,Q<:AbstractArray,U<:AbstractArray,W,Z<:Real,D<:AbstractDomain,I<:Int64,F<:Float64}
    f(dy::X,x::Y) where {Y<:Real,X} = _dydttherm(dy,x,y,p,t,domain,interfaces,ind)
    jac[:,ind] = ForwardDiff.derivative(f,zeros(size(y)),x)
end

# function jacobianp!(d::W;cs::Q,V::Y,T::Y2,Us::Z3,Cvave::Z4,N::Z5,kfs::Z,krevs::X,wV::Q2,ratederiv::Q3) where {Q3,W<:Union{ConstantTPDomain,ConstantTVDomain},Z4<:Real,Z5<:Real,Z3<:AbstractArray,Q2<:AbstractArray,Q<:AbstractArray,Y2<:Real,Y<:Real,Z<:AbstractArray,X<:AbstractArray}
#     Nspcs = length(cs)
#     rxns = d.phase.reactions
#     Nrxns = length(rxns)
#     RTinv = 1.0/(R*T)
#     ratederiv .= 0.0
#
#     for (j,rxn) in enumerate(rxns)
#         Nreact = length(rxn.reactantinds)
#         Nprod = length(rxn.productinds)
#
#         if Nreact == 1
#             rind1 = rxn.reactantinds[1]
#             fderiv = cs[rind1]
#         elseif Nreact == 2
#             rind1,rind2 = rxn.reactantinds
#             fderiv = cs[rind1]*cs[rind2]
#         else
#             rind1,rind2,rind3 = rxn.reactantinds
#             fderiv = cs[rind1]*cs[rind2]*cs[rind3]
#         end
#
#         if Nprod == 1
#             pind1 = rxn.productinds[1]
#             rderiv = krevs[j]/kfs[j]*cs[pind1]
#         elseif Nprod == 2
#             pind1,pind2 = rxn.productinds
#             rderiv = krevs[j]/kfs[j]*cs[pind1]*cs[pind2]
#         else
#             pind1,pind2,pind3 = rxn.productinds
#             rderiv = krevs[j]/kfs[j]*cs[pind1]*cs[pind2]*cs[pind3]
#         end
#
#         flux = fderiv-rderiv
#         gderiv = rderiv*kfs[j]*RTinv
#
#         deriv = zeros(Nspcs)
#
#         deriv[rind1] += gderiv
#         if Nreact > 1
#             deriv[rind2] += gderiv
#             if Nreact > 2
#                 deriv[rind3] == gderiv
#             end
#         end
#
#         deriv[pind1] -= gderiv
#         if Nprod > 1
#             deriv[pind2] -= gderiv
#             if Nprod > 2
#                 deriv[pind3] -= gderiv
#             end
#         end
#
#         ratederiv[rind1,j] -= flux
#         ratederiv[rind1,Nrxns+1:Nrxns+Nspcs] .-= deriv
#         if Nreact > 1
#             ratederiv[rind2,j] -= flux
#             ratederiv[rind2,Nrxns+1:Nrxns+Nspcs] .-= deriv
#             if Nreact > 2
#                 ratederiv[rind3,j] -= flux
#                 ratederiv[rind3,Nrxns+1:Nrxns+Nspcs] .-= deriv
#             end
#         end
#
#         ratederiv[pind1,j] += flux
#         ratederiv[pind1,Nrxns+1:Nrxns+Nspcs] .+= deriv
#         if Nprod > 1
#             ratederiv[pind2,j] += flux
#             ratederiv[pind2,Nrxns+1:Nrxns+Nspcs] .+= deriv
#             if Nprod > 2
#                 ratederiv[pind3,j] += flux
#                 ratederiv[pind3,Nrxns+1:Nrxns+Nspcs] .+= deriv
#             end
#         end
#     end
#     return V*ratederiv
# end
#
# function jacobianp!(d::W; cs::Q,V::Y,T::Y2,Us::Z3,Cvave::Y3,N::Y2,kfs::Z,krevs::X,wV::Q2,ratederiv::Q3) where {W<:Union{ConstantVDomain,ParametrizedVDomain},Q3,Z3<:AbstractArray,Q<:AbstractArray,Q2<:AbstractArray,Y3<:Real,Y2<:Real,Y<:Real,Z<:AbstractArray,X<:AbstractArray}
#     Nspcs = length(cs)
#     rxns = d.phase.reactions
#     Nrxns = length(rxns)
#     RTinv = 1.0/(R*T)
#     ratederiv .= 0.0
#
#     for (j,rxn) in enumerate(rxns)
#         Nreact = length(rxn.reactantinds)
#         Nprod = length(rxn.productinds)
#
#         if Nreact == 1
#             rind1 = rxn.reactantinds[1]
#             fderiv = cs[rind1]
#         elseif Nreact == 2
#             rind1,rind2 = rxn.reactantinds
#             fderiv = cs[rind1]*cs[rind2]
#         else
#             rind1,rind2,rind3 = rxn.reactantinds
#             fderiv = cs[rind1]*cs[rind2]*cs[rind3]
#         end
#
#         if Nprod == 1
#             pind1 = rxn.productinds[1]
#             rderiv = krevs[j]/kfs[j]*cs[pind1]
#         elseif Nprod == 2
#             pind1,pind2 = rxn.productinds
#             rderiv = krevs[j]/kfs[j]*cs[pind1]*cs[pind2]
#         else
#             pind1,pind2,pind3 = rxn.productinds
#             rderiv = krevs[j]/kfs[j]*cs[pind1]*cs[pind2]*cs[pind3]
#         end
#
#         flux = fderiv-rderiv
#         gderiv = rderiv*kfs[j]*RTinv
#
#         deriv = zeros(Nspcs)
#
#         deriv[rind1] += gderiv
#         if Nreact > 1
#             deriv[rind2] += gderiv
#             if Nreact > 2
#                 deriv[rind3] == gderiv
#             end
#         end
#
#         deriv[pind1] -= gderiv
#         if Nprod > 1
#             deriv[pind2] -= gderiv
#             if Nprod > 2
#                 deriv[pind3] -= gderiv
#             end
#         end
#
#         ratederiv[rind1,j] -= flux
#         ratederiv[rind1,Nrxns+1:Nrxns+Nspcs] .-= deriv
#         if Nreact > 1
#             ratederiv[rind2,j] -= flux
#             ratederiv[rind2,Nrxns+1:Nrxns+Nspcs] .-= deriv
#             if Nreact > 2
#                 ratederiv[rind3,j] -= flux
#                 ratederiv[rind3,Nrxns+1:Nrxns+Nspcs] .-= deriv
#             end
#         end
#
#         ratederiv[pind1,j] += flux
#         ratederiv[pind1,Nrxns+1:Nrxns+Nspcs] .+= deriv
#         if Nprod > 1
#             ratederiv[pind2,j] += flux
#             ratederiv[pind2,Nrxns+1:Nrxns+Nspcs] .+= deriv
#             if Nprod > 2
#                 ratederiv[pind3,j] += flux
#                 ratederiv[pind3,Nrxns+1:Nrxns+Nspcs] .+= deriv
#             end
#         end
#     end
#     ratederiv *= V
#     #Temperature stuff
#     @views ratederiv[end,:] += (Us'*ratederiv[1:end-1,:])[1,:]
#     ratederiv[end,1:Nspcs] += wV
#     ratederiv[end,:] /= (N*Cvave)
#     return ratederiv
# end

# @inline function spreadpartials!(jac::S,deriv::T,inds::V,ind::Q,N::Q) where {S<:AbstractArray, T<:Real, V<:AbstractArray, Q<:Integer}
#     if N == 1
#         jac[inds[1],ind] += deriv
#     elseif N == 2
#         jac[inds[1],ind] += deriv
#         jac[inds[2],ind] += deriv
#     elseif N == 3
#         jac[inds[1],ind] += deriv
#         jac[inds[2],ind] += deriv
#         jac[inds[3],ind] += deriv
#     end
# end
#
# @inline function spreadpartials!(jac::S,deriv::T,inds::V,ind::Q,N::Q) where {S<:AbstractArray, T<:Real, V<:AbstractArray, Q<:Integer}
#     if N == 1
#         jac[inds[1],ind] += deriv
#     elseif N == 2
#         jac[inds[1],ind] += deriv
#         jac[inds[2],ind] += deriv
#     elseif N == 3
#         jac[inds[1],ind] += deriv
#         jac[inds[2],ind] += deriv
#         jac[inds[3],ind] += deriv
#     end
# end
#
# function jacobiany!(y::Array{T,1},t::T,domain::ConstantTPDomain,kfs::Array{T,1},krevs::Array{T,1},jac::P;zero::Bool=true) where {P<:AbstractArray,T<:Real,J<:Integer}
#     if zero
#         jac .= 0
#     end
#     N = sum(y)
#     V = N*R*domain.T/domain.P
#     cs = y./V
#     C = N/V
#     rxnarray = domain.rxnarray
#     Nrxns = size(rxnarray)[2]
#     Nspcs = length(y)
#     for i in 1:Nrxns
#         kf = kfs[i]
#         krev = krevs[i]
#         if rxnarray[2,i] == 0
#             jac[rxnarray[1,i],rxnarray[1,i]] -= kf
#             if rxnarray[5,i] == 0
#                 jac[rxnarray[4,i],rxnarray[1,i]] += kf
#             elseif rxnarray[6,i] == 0
#                 jac[rxnarray[4,i],rxnarray[1,i]] += kf
#                 jac[rxnarray[5,i],rxnarray[1,i]] += kf
#             else
#                 jac[rxnarray[4,i],rxnarray[1,i]] += kf
#                 jac[rxnarray[5,i],rxnarray[1,i]] += kf
#                 jac[rxnarray[6,i],rxnarray[1,i]] += kf
#             end
#         elseif rxnarray[3,i] == 0
#             corr = -kf*cs[rxnarray[1,i]]*cs[rxnarray[2,i]]/C #correction for the partial of the volume term
#             if rxnarray[1,i] == rxnarray[2,i]
#                 deriv = 2*kf*cs[rxnarray[1,i]]
#                 jac[rxnarray[1,i],rxnarray[1,i]] -= 2.0*deriv
#                 for j in 1:Nspcs
#                     jac[rxnarray[1,i],j] -= 2.0*corr
#                 end
#                 jac[rxnarray[4,i],rxnarray[1,i]] += deriv
#                 for j in 1:Nspcs
#                     jac[rxnarray[4,i],j] += corr
#                 end
#                 if rxnarray[5,i] != 0
#                     jac[rxnarray[5,i],rxnarray[1,i]] += deriv
#                     for j = 1:Nspcs
#                         jac[rxnarray[5,i],j] += corr
#                     end
#                     if rxnarray[6,i] != 0
#                         jac[rxnarray[6,i],rxnarray[1,i]] += deriv
#                         for j = 1:Nspcs
#                             jac[rxnarray[6,i],j] += corr
#                         end
#                     end
#                 end
#             else
#                 #derivative with respect to reactant 1
#                 deriv = kf*cs[rxnarray[2,i]]
#                 jac[rxnarray[1,i],rxnarray[1,i]] -= deriv
#                 jac[rxnarray[2,i],rxnarray[1,i]] -= deriv
#
#                 jac[rxnarray[4,i],rxnarray[1,i]] += deriv
#                 if rxnarray[5,i] != 0
#                     jac[rxnarray[5,i],rxnarray[1,i]] += deriv
#                     if rxnarray[6,i] != 0
#                         jac[rxnarray[6,i],rxnarray[1,i]] += deriv
#                     end
#                 end
#
#                 #derivative with respect to reactant 2
#                 deriv = kf*cs[rxnarray[1,i]]
#                 jac[rxnarray[1,i],rxnarray[2,i]] -= deriv
#                 jac[rxnarray[2,i],rxnarray[2,i]] -= deriv
#                 for j = 1:Nspcs
#                     jac[rxnarray[1,i],j] -= corr
#                     jac[rxnarray[2,i],j] -= corr
#                 end
#                 jac[rxnarray[4,i],rxnarray[2,j]] += deriv
#                 if rxnarray[5,i] != 0
#                     jac[rxnarray[5,i],rxnarray[2,i]] += deriv
#                     for j = 1:Nspcs
#                         jac[rxnarray[5,i],j] += corr
#                     end
#                     if rxnarray[6,i] != 0
#                         jac[rxnarray[6,i],rxnarray[2,i]] += deriv
#                         for j = 1:Nspcs
#                             jac[rxnarray[6,i],j] += corr
#                         end
#                     end
#                 end
#             end
#         else
#             corr = -2.0*kf*cs[rxnarray[1,i]]*cs[rxnarray[2,i]]*cs[rxnarray[3,i]]/C
#             if (rxnarray[1,i] == rxnarray[2,i] && rxnarray[1,i] == rxnarray[3,i])
#                 deriv = 3.0*kf*cs[rxnarray[1,i]]*cs[rxnarray[1,i]]
#                 jac[rxnarray[1,i],rxnarray[1,i]] -= 3.0*deriv
#                 for j = 1:Nspcs
#                     jac[rxnarray[1,i],j] -= 3.0*corr
#                 end
#                 jac[rxnarray[4,i],rxnarray[1,i]] += deriv
#                 for j = 1:Nspcs
#                     jac[rxnarray[4,i],j] += corr
#                 end
#                 if rxnarray[5,i] != 0
#                     jac[rxnarray[5,i],rxnarray[1,i]] += deriv
#                     for j = 1:Nspcs
#                         jac[rxnarray[5,i],j] += corr
#                     end
#                     if rxnarray[6,i] != 0
#                         jac[rxnarray[6,i],rxnarray[1,i]] += deriv
#                         for j = 1:Nspcs
#                             jac[rxnarray[6,i],j] += corr
#                         end
#                     end
#                 end
#             elseif rxnarray[1,i] == rxnarray[2,i]
#                 #derivative with respect to reactant 1
#                 deriv = 2.0*kf*cs[rxnarray[1,i]]*cs[rxnarray[3,i]]
#                 jac[rxnarray[1,i],rxnarray[1,i]] -= 2.0*deriv
#                 jac[rxnarray[3,i],rxnarray[1,i]] -= deriv
#
#                 jac[rxnarray[4,i],rxnarray[1,i]] += deriv
#                 if rxnarray[5,i] != 0
#                     jac[rxnarray[5,i],rxnarray[1,i]] += deriv
#                     if rxnarray[6,i] != 0
#                         jac[rxnarray[6,i],rxnarray[1,i]] += deriv
#                     end
#                 end
#
#                 #derivative with respect to reactant 3
#                 deriv = kf*cs[rxnarray]
#
#             ind1,ind2,ind3 = rxn.reactantinds
#             corr = -2.0*kf*cs[ind1]*cs[ind2]*cs[ind3]/C
#             deriv = kf*cs[ind1]*cs[ind2]
#             jac[ind1,ind3] -= deriv
#             jac[ind2,ind3] -= deriv
#             jac[ind3,ind3] -= deriv
#             spreadpartials!(jac,deriv,rxn.productinds,ind3,Nprod)
#             deriv = kf*cs[ind1]*cs[ind3]
#             jac[ind1,ind2] -= deriv
#             jac[ind2,ind2] -= deriv
#             jac[ind3,ind2] -= deriv
#             spreadpartials!(jac,deriv,rxn.productinds,ind2,Nprod)
#             deriv = kf*cs[ind3]*cs[ind2]
#             jac[ind1,ind1] -= deriv
#             jac[ind2,ind1] -= deriv
#             jac[ind3,ind1] -= deriv
#             spreadpartials!(jac,deriv,rxn.productinds,ind1,Nprod)
#             for i in 1:length(domain.phase.species)
#                 jac[ind1,i] -= corr
#                 jac[ind2,i] -= corr
#                 jac[ind3,i] -= corr
#                 spreadpartials!(jac,corr,rxn.productinds,i,Nprod)
#             end
#         end
#         #reverse direction
#         if Nprod == 1
#             ind1 = rxn.productinds[1]
#             jac[ind1,ind1] -= krev
#             if Nprod == 1
#                 jac[rxn.reactantinds[1],ind1] += krev
#             elseif Nprod == 2
#                 jac[rxn.reactantinds[1],ind1] += krev
#                 jac[rxn.reactantinds[2],ind1] += krev
#             elseif Nprod == 3
#                 jac[rxn.reactantinds[1],ind1] += krev
#                 jac[rxn.reactantinds[2],ind1] += krev
#                 jac[rxn.reactantinds[3],ind1] += krev
#             end
#         elseif Nprod == 2
#             ind1,ind2 = rxn.productinds
#             corr = -krev*cs[ind1]*cs[ind2]/C #correction for the partial of the volume term
#             deriv = krev*cs[ind1]
#             jac[ind1,ind2] -= deriv
#             jac[ind2,ind2] -= deriv
#             spreadpartials!(jac,deriv,rxn.reactantinds,ind2,Nreact)
#             deriv = kf*cs[ind2]
#             jac[ind1,ind1] -= deriv
#             jac[ind2,ind1] -= deriv
#             spreadpartials!(jac,deriv,rxn.reactantinds,ind1,Nreact)
#             for i in 1:length(domain.phase.species)
#                 jac[ind1,i] -= corr
#                 jac[ind2,i] -= corr
#                 spreadpartials!(jac,corr,rxn.reactantinds,i,Nreact)
#             end
#         elseif Nprod == 3
#             ind1,ind2,ind3 = rxn.productinds
#             corr = -2.0*krev*cs[ind1]*cs[ind2]*cs[ind3]/C
#             deriv = krev*cs[ind1]*cs[ind2]
#             jac[ind1,ind3] -= deriv
#             jac[ind2,ind3] -= deriv
#             jac[ind3,ind3] -= deriv
#             spreadpartials!(jac,deriv,rxn.reactantinds,ind3,Nreact)
#             deriv = kf*cs[ind1]*cs[ind3]
#             jac[ind1,ind2] -= deriv
#             jac[ind2,ind2] -= deriv
#             jac[ind3,ind2] -= deriv
#             spreadpartials!(jac,deriv,rxn.reactantinds,ind2,Nreact)
#             deriv = kf*cs[ind3]*cs[ind2]
#             jac[ind1,ind1] -= deriv
#             jac[ind2,ind1] -= deriv
#             jac[ind3,ind1] -= deriv
#             spreadpartials!(jac,deriv,rxn.reactantinds,ind1,Nreact)
#             for i in 1:length(domain.phase.species)
#                 jac[ind1,i] -= corr
#                 jac[ind2,i] -= corr
#                 jac[ind3,i] -= corr
#                 spreadpartials!(jac,corr,rxn.reactantinds,i,Nreact)
#             end
#         end
#     end
#     for ind in domain.constantspeciesinds
#         jac[ind,:] .= 0
#     end
#     return jac
# end
#
# function jacobiany!(y::Array{T,1},t::T,domain::ConstantTPDomain,kfs::Array{T,1},krevs::Array{T,1},jac::P;zero::Bool=true) where {P<:AbstractArray,T<:Real,J<:Integer}
#     if zero
#         jac .= 0
#     end
#     N = sum(y)
#     V = N*R*domain.T/domain.P
#     cs = y./V
#     C = N/V
#     for (i,rxn) in enumerate(domain.phase.reactions)
#         Nreact = length(rxn.reactantinds)
#         Nprod = length(rxn.productinds)
#         kf = kfs[i]
#         krev = krevs[i]
#         if Nreact == 1
#             ind1 = rxn.reactantinds[1]
#             jac[ind1,ind1] -= kf
#             if Nprod == 1
#                 jac[rxn.productinds[1],ind1] += kf
#             elseif Nprod == 2
#                 jac[rxn.productinds[1],ind1] += kf
#                 jac[rxn.productinds[2],ind1] += kf
#             elseif Nprod == 3
#                 jac[rxn.productinds[1],ind1] += kf
#                 jac[rxn.productinds[2],ind1] += kf
#                 jac[rxn.productinds[3],ind1] += kf
#             end
#         elseif Nreact == 2
#             ind1,ind2 = rxn.reactantinds
#             corr = -kf*cs[ind1]*cs[ind2]/C #correction for the partial of the volume term
#             deriv = kf*cs[ind1]
#             jac[ind1,ind2] -= deriv
#             jac[ind2,ind2] -= deriv
#             spreadpartials!(jac,deriv,rxn.productinds,ind2,Nprod)
#             deriv = kf*cs[ind2]
#             jac[ind1,ind1] -= deriv
#             jac[ind2,ind1] -= deriv
#             spreadpartials!(jac,deriv,rxn.productinds,ind1,Nprod)
#             for i in 1:length(domain.phase.species)
#                 jac[ind1,i] -= corr
#                 jac[ind2,i] -= corr
#                 spreadpartials!(jac,corr,rxn.productinds,i,Nprod)
#             end
#         elseif Nreact == 3
#             ind1,ind2,ind3 = rxn.reactantinds
#             corr = -2.0*kf*cs[ind1]*cs[ind2]*cs[ind3]/C
#             deriv = kf*cs[ind1]*cs[ind2]
#             jac[ind1,ind3] -= deriv
#             jac[ind2,ind3] -= deriv
#             jac[ind3,ind3] -= deriv
#             spreadpartials!(jac,deriv,rxn.productinds,ind3,Nprod)
#             deriv = kf*cs[ind1]*cs[ind3]
#             jac[ind1,ind2] -= deriv
#             jac[ind2,ind2] -= deriv
#             jac[ind3,ind2] -= deriv
#             spreadpartials!(jac,deriv,rxn.productinds,ind2,Nprod)
#             deriv = kf*cs[ind3]*cs[ind2]
#             jac[ind1,ind1] -= deriv
#             jac[ind2,ind1] -= deriv
#             jac[ind3,ind1] -= deriv
#             spreadpartials!(jac,deriv,rxn.productinds,ind1,Nprod)
#             for i in 1:length(domain.phase.species)
#                 jac[ind1,i] -= corr
#                 jac[ind2,i] -= corr
#                 jac[ind3,i] -= corr
#                 spreadpartials!(jac,corr,rxn.productinds,i,Nprod)
#             end
#         end
#         #reverse direction
#         if Nprod == 1
#             ind1 = rxn.productinds[1]
#             jac[ind1,ind1] -= krev
#             if Nprod == 1
#                 jac[rxn.reactantinds[1],ind1] += krev
#             elseif Nprod == 2
#                 jac[rxn.reactantinds[1],ind1] += krev
#                 jac[rxn.reactantinds[2],ind1] += krev
#             elseif Nprod == 3
#                 jac[rxn.reactantinds[1],ind1] += krev
#                 jac[rxn.reactantinds[2],ind1] += krev
#                 jac[rxn.reactantinds[3],ind1] += krev
#             end
#         elseif Nprod == 2
#             ind1,ind2 = rxn.productinds
#             corr = -krev*cs[ind1]*cs[ind2]/C #correction for the partial of the volume term
#             deriv = krev*cs[ind1]
#             jac[ind1,ind2] -= deriv
#             jac[ind2,ind2] -= deriv
#             spreadpartials!(jac,deriv,rxn.reactantinds,ind2,Nreact)
#             deriv = kf*cs[ind2]
#             jac[ind1,ind1] -= deriv
#             jac[ind2,ind1] -= deriv
#             spreadpartials!(jac,deriv,rxn.reactantinds,ind1,Nreact)
#             for i in 1:length(domain.phase.species)
#                 jac[ind1,i] -= corr
#                 jac[ind2,i] -= corr
#                 spreadpartials!(jac,corr,rxn.reactantinds,i,Nreact)
#             end
#         elseif Nprod == 3
#             ind1,ind2,ind3 = rxn.productinds
#             corr = -2.0*krev*cs[ind1]*cs[ind2]*cs[ind3]/C
#             deriv = krev*cs[ind1]*cs[ind2]
#             jac[ind1,ind3] -= deriv
#             jac[ind2,ind3] -= deriv
#             jac[ind3,ind3] -= deriv
#             spreadpartials!(jac,deriv,rxn.reactantinds,ind3,Nreact)
#             deriv = kf*cs[ind1]*cs[ind3]
#             jac[ind1,ind2] -= deriv
#             jac[ind2,ind2] -= deriv
#             jac[ind3,ind2] -= deriv
#             spreadpartials!(jac,deriv,rxn.reactantinds,ind2,Nreact)
#             deriv = kf*cs[ind3]*cs[ind2]
#             jac[ind1,ind1] -= deriv
#             jac[ind2,ind1] -= deriv
#             jac[ind3,ind1] -= deriv
#             spreadpartials!(jac,deriv,rxn.reactantinds,ind1,Nreact)
#             for i in 1:length(domain.phase.species)
#                 jac[ind1,i] -= corr
#                 jac[ind2,i] -= corr
#                 jac[ind3,i] -= corr
#                 spreadpartials!(jac,corr,rxn.reactantinds,i,Nreact)
#             end
#         end
#     end
#     for ind in domain.constantspeciesinds
#         jac[ind,:] .= 0
#     end
#     return jac
# end
# function jacobiany!(y::Array{T,1},t::T,domain::ConstantTVDomain,kfs::Array{T,1},krevs::Array{T,1},jac::P;zero::Bool=true) where {P<:AbstractArray,T<:Real,J<:Integer}
#     if zero
#         jac .= 0
#     end
#     cs = y./domain.V
#     for (i,rxn) in enumerate(domain.phase.reactions)
#         Nreact = length(rxn.reactantinds)
#         Nprod = length(rxn.productinds)
#         kf = kfs[i]
#         krev = krevs[i]
#         if Nreact == 1
#             ind1 = rxn.reactantinds[1]
#             jac[ind1,ind1] -= kf
#             if Nprod == 1
#                 jac[rxn.productinds[1],ind1] += kf
#             elseif Nprod == 2
#                 jac[rxn.productinds[1],ind1] += kf
#                 jac[rxn.productinds[2],ind1] += kf
#             elseif Nprod == 3
#                 jac[rxn.productinds[1],ind1] += kf
#                 jac[rxn.productinds[2],ind1] += kf
#                 jac[rxn.productinds[3],ind1] += kf
#             end
#         elseif Nreact == 2
#             ind1,ind2 = rxn.reactantinds
#             deriv = kf*cs[ind1]
#             jac[ind1,ind2] -= deriv
#             jac[ind2,ind2] -= deriv
#             spreadpartials!(jac,deriv,rxn.productinds,ind2,Nprod)
#             deriv = kf*cs[ind2]
#             jac[ind1,ind1] -= deriv
#             jac[ind2,ind1] -= deriv
#             spreadpartials!(jac,deriv,rxn.productinds,ind1,Nprod)
#         elseif Nreact == 3
#             ind1,ind2,ind3 = rxn.reactantinds
#             deriv = kf*state.cs[ind1]*cs[ind2]
#             jac[ind1,ind3] -= deriv
#             jac[ind2,ind3] -= deriv
#             jac[ind3,ind3] -= deriv
#             spreadpartials!(jac,deriv,rxn.productinds,ind3,Nprod)
#             deriv = kf*cs[ind1]*cs[ind3]
#             jac[ind1,ind2] -= deriv
#             jac[ind2,ind2] -= deriv
#             jac[ind3,ind2] -= deriv
#             spreadpartials!(jac,deriv,rxn.productinds,ind2,Nprod)
#             deriv = kf*cs[ind3]*cs[ind2]
#             jac[ind1,ind1] -= deriv
#             jac[ind2,ind1] -= deriv
#             jac[ind3,ind1] -= deriv
#             spreadpartials!(jac,deriv,rxn.productinds,ind1,Nprod)
#         end
#         #reverse direction
#         if Nprod == 1
#             ind1 = rxn.productinds[1]
#             jac[ind1,ind1] -= krev
#             if Nprod == 1
#                 jac[rxn.reactantinds[1],ind1] += krev
#             elseif Nprod == 2
#                 jac[rxn.reactantinds[1],ind1] += krev
#                 jac[rxn.reactantinds[2],ind1] += krev
#             elseif Nprod == 3
#                 jac[rxn.reactantinds[1],ind1] += krev
#                 jac[rxn.reactantinds[2],ind1] += krev
#                 jac[rxn.reactantinds[3],ind1] += krev
#             end
#         elseif Nprod == 2
#             ind1,ind2 = rxn.productinds
#             deriv = krev*cs[ind1]
#             jac[ind1,ind2] -= deriv
#             jac[ind2,ind2] -= deriv
#             spreadpartials!(jac,deriv,rxn.reactantinds,ind2,Nreact)
#             deriv = kf*cs[ind2]
#             jac[ind1,ind1] -= deriv
#             jac[ind2,ind1] -= deriv
#             spreadpartials!(jac,deriv,rxn.reactantinds,ind1,Nreact)
#         elseif Nprod == 3
#             ind1,ind2,ind3 = rxn.productinds
#             deriv = krev*cs[ind1]*cs[ind2]
#             jac[ind1,ind3] -= deriv
#             jac[ind2,ind3] -= deriv
#             jac[ind3,ind3] -= deriv
#             spreadpartials!(jac,deriv,rxn.reactantinds,ind3,Nreact)
#             deriv = kf*cs[ind1]*cs[ind3]
#             jac[ind1,ind2] -= deriv
#             jac[ind2,ind2] -= deriv
#             jac[ind3,ind2] -= deriv
#             spreadpartials!(jac,deriv,rxn.reactantinds,ind2,Nreact)
#             deriv = kf*cs[ind3]*cs[ind2]
#             jac[ind1,ind1] -= deriv
#             jac[ind2,ind1] -= deriv
#             jac[ind3,ind1] -= deriv
#             spreadpartials!(jac,deriv,rxn.reactantinds,ind1,Nreact)
#         end
#     end
#     for ind in domain.constantspeciesinds
#         jac[ind,:] .= 0
#     end
#     return jac
# end


text:
using CUDAnative: device!
device!(0)
using BenchmarkTools, OMEinsum, CuArrays
CuArrays.allowscalar(false)

function bfunc_star(N::Int)
    ca = CuArrays.curandn(Float32, N, N)
    res = CuArrays.@sync ein"ji,kl,li->jkl"(ca,ca,ca)
    #res ≈ ein"ji,kl,li->jkl"(Matrix(ca),Matrix(ca), Matrix(ca))
end

function bfunc_t3(N::Int)
    ca = CuArrays.curandn(Float32, N, N)
    ct = CuArrays.curandn(Float32, N, N, N)
    CuArrays.@sync ein"ji,kli,li->jkl"(ca,ct,ca)
end

function bfunc_psum(N::Int)
    cq = CuArrays.curandn(Float32, N, N, N, N)
    CuArrays.@sync ein"iikl->"(cq)
end

display(@benchmark bfunc_star(300) seconds = 1)
display(@benchmark bfunc_t3(300) seconds = 1)
display(@benchmark bfunc_psum(100) seconds = 1)

ext: jl


text:
convert(::Type{JString}, str::AbstractString) = JString(str)
convert(::Type{JObject}, str::AbstractString) = convert(JObject, JString(str))

#Cast java object from S to T . Needed for polymorphic calls
function convert(::Type{JavaObject{T}}, obj::JavaObject{S}) where {T,S}
    if isConvertible(T, S)   #Safe static cast
        ptr = ccall(jnifunc.NewLocalRef, Ptr{Nothing}, (Ptr{JNIEnv}, Ptr{Nothing}), penv, obj.ptr)
        ptr === C_NULL && geterror()
        return JavaObject{T}(ptr)
    end
    isnull(obj) && throw(ArgumentError("Cannot convert NULL"))
    realClass = ccall(jnifunc.GetObjectClass, Ptr{Nothing}, (Ptr{JNIEnv}, Ptr{Nothing} ), penv, obj.ptr)
    if isConvertible(T, realClass)  #dynamic cast
        ptr = ccall(jnifunc.NewLocalRef, Ptr{Nothing}, (Ptr{JNIEnv}, Ptr{Nothing}), penv, obj.ptr)
        ptr === C_NULL && geterror()
        return JavaObject{T}(ptr)
    end
    throw(JavaCallError("Cannot cast java object from $S to $T"))
end

#Is java type convertible from S to T.
isConvertible(T, S) = (ccall(jnifunc.IsAssignableFrom, jboolean,
                             (Ptr{JNIEnv}, Ptr{Nothing}, Ptr{Nothing}), penv, metaclass(S),
                             metaclass(T)) == JNI_TRUE)
isConvertible(T, S::Ptr{Nothing}) = (ccall(jnifunc.IsAssignableFrom, jboolean,
                                           (Ptr{JNIEnv}, Ptr{Nothing}, Ptr{Nothing}), penv, S,
                                           metaclass(T)) == JNI_TRUE)

unsafe_convert(::Type{Ptr{Nothing}}, cls::JavaMetaClass) = cls.ptr

# Get the JNI/C type for a particular Java type
function real_jtype(rettype)
    if rettype <: JavaObject || rettype <: Array || rettype <: JavaMetaClass
        jnitype = Ptr{Nothing}
    else
        jnitype = rettype
    end
    return jnitype
end

function convert_args(argtypes::Tuple, args...)
    convertedArgs = Array{Int64}(undef, length(args))
    savedArgs = Array{Any}(undef, length(args))
    for i in 1:length(args)
        r = convert_arg(argtypes[i], args[i])
        savedArgs[i] = r[1]
        convertedArgs[i] = jvalue(r[2])
    end
    return savedArgs, convertedArgs
end

function convert_arg(argtype::Type{JString}, arg)
    x = convert(JString, arg)
    return x, x.ptr
end

function convert_arg(argtype::Type, arg)
    x = convert(argtype, arg)
    return x,x
end
function convert_arg(argtype::Type{T}, arg) where T<:JavaObject
    x = convert(T, arg)::T
    return x, x.ptr
end

for (x, y, z) in [ (:jboolean, :(jnifunc.NewBooleanArray), :(jnifunc.SetBooleanArrayRegion)),
                  (:jchar, :(jnifunc.NewCharArray), :(jnifunc.SetCharArrayRegion)),
                  (:jbyte, :(jnifunc.NewByteArray), :(jnifunc.SetByteArrayRegion)),
                  (:jshort, :(jnifunc.NewShortArray), :(jnifunc.SetShortArrayRegion)),
                  (:jint, :(jnifunc.NewIntArray), :(jnifunc.SetIntArrayRegion)),
                  (:jlong, :(jnifunc.NewLongArray), :(jnifunc.SetLongArrayRegion)),
                  (:jfloat, :(jnifunc.NewFloatArray), :(jnifunc.SetFloatArrayRegion)),
                  (:jdouble, :(jnifunc.NewDoubleArray), :(jnifunc.SetDoubleArrayRegion)) ]
    m = quote
        function convert_arg(argtype::Type{Array{$x,1}}, arg)
            carg = convert(argtype, arg)
            sz=length(carg)
            arrayptr = ccall($y, Ptr{Nothing}, (Ptr{JNIEnv}, jint), penv, sz)
            ccall($z, Nothing, (Ptr{JNIEnv}, Ptr{Nothing}, jint, jint, Ptr{$x}), penv, arrayptr, 0, sz,
                  carg)
            return carg, arrayptr
        end
    end
    eval( m)
end

function convert_arg(argtype::Type{Array{T,1}}, arg) where T<:JavaObject
    carg = convert(argtype, arg)
    sz = length(carg)
    init = carg[1]
    arrayptr = ccall(jnifunc.NewObjectArray, Ptr{Nothing},
                     (Ptr{JNIEnv}, jint, Ptr{Nothing}, Ptr{Nothing}), penv, sz, metaclass(T), init.ptr)
    for i=2:sz
        ccall(jnifunc.SetObjectArrayElement, Nothing, (Ptr{JNIEnv}, Ptr{Nothing}, jint, Ptr{Nothing}),
              penv, arrayptr, i-1, carg[i].ptr)
    end
    return carg, arrayptr
end

convert_result(rettype::Type{T}, result) where {T<:JString} = unsafe_string(JString(result))
convert_result(rettype::Type{T}, result) where {T<:JavaObject} = T(result)
convert_result(rettype, result) = result

for (x, y, z) in [(:jboolean, :(jnifunc.GetBooleanArrayElements), :(jnifunc.ReleaseBooleanArrayElements)),
                  (:jchar, :(jnifunc.GetCharArrayElements), :(jnifunc.ReleaseCharArrayElements)),
                  (:jbyte, :(jnifunc.GetByteArrayElements), :(jnifunc.ReleaseByteArrayElements)),
                  (:jshort, :(jnifunc.GetShortArrayElements), :(jnifunc.ReleaseShortArrayElements)),
                  (:jint, :(jnifunc.GetIntArrayElements), :(jnifunc.ReleaseIntArrayElements)),
                  (:jlong, :(jnifunc.GetLongArrayElements), :(jnifunc.ReleaseLongArrayElements)),
                  (:jfloat, :(jnifunc.GetFloatArrayElements), :(jnifunc.ReleaseFloatArrayElements)),
                  (:jdouble, :(jnifunc.GetDoubleArrayElements), :(jnifunc.ReleaseDoubleArrayElements)) ]
    m = quote
        function convert_result(rettype::Type{Array{$(x),1}}, result)
            sz = ccall(jnifunc.GetArrayLength, jint, (Ptr{JNIEnv}, Ptr{Nothing}), penv, result)
            arr = ccall($(y), Ptr{$(x)}, (Ptr{JNIEnv}, Ptr{Nothing}, Ptr{jboolean} ), penv, result,
                        C_NULL)
            jl_arr::Array = unsafe_wrap(Array, arr, Int(sz))
            jl_arr = deepcopy(jl_arr)
            ccall($(z), Nothing, (Ptr{JNIEnv},Ptr{Nothing}, Ptr{$(x)}, jint), penv, result, arr, 0)
            return jl_arr
        end
    end
    eval(m)
end

function convert_result(rettype::Type{Array{JavaObject{T},1}}, result) where T
    sz = ccall(jnifunc.GetArrayLength, jint, (Ptr{JNIEnv}, Ptr{Nothing}), penv, result)

    ret = Array{JavaObject{T}}(undef, sz)

    for i=1:sz
        a = ccall(jnifunc.GetObjectArrayElement, Ptr{Nothing}, (Ptr{JNIEnv},Ptr{Nothing}, jint), penv,
                  result, i-1)
        ret[i] = JavaObject{T}(a)
    end
    return ret
end

# covers return types like Vector{Vector{T}}
function convert_result(rettype::Type{Array{T,1}}, result) where T
    sz = ccall(jnifunc.GetArrayLength, jint, (Ptr{JNIEnv}, Ptr{Nothing}), penv, result)

    ret = Array{T}(undef, sz)

    for i=1:sz
        a=ccall(jnifunc.GetObjectArrayElement, Ptr{Nothing}, (Ptr{JNIEnv},Ptr{Nothing}, jint), penv, result, i-1)
        ret[i] = convert_result(T, a)
    end
    return ret
end

function convert_result(rettype::Type{Array{JavaObject{T},2}}, result) where T
    sz = ccall(jnifunc.GetArrayLength, jint, (Ptr{JNIEnv}, Ptr{Nothing}), penv, result)
    if sz == 0
        return Array{T}(undef, 0,0)
    end
    a_1 = ccall(jnifunc.GetObjectArrayElement, Ptr{Nothing}, (Ptr{JNIEnv},Ptr{Nothing}, jint), penv,
                result, 0)
    sz_1 = ccall(jnifunc.GetArrayLength, jint, (Ptr{JNIEnv}, Ptr{Nothing}), penv, a_1)
    ret = Array{JavaObject{T}}(undef, sz, sz_1)
    for i=1:sz
        a = ccall(jnifunc.GetObjectArrayElement, Ptr{Nothing}, (Ptr{JNIEnv},Ptr{Nothing}, jint), penv,
                  result, i-1)
        # check that size of the current subarray is the same as for the first one
        sz_a = ccall(jnifunc.GetArrayLength, jint, (Ptr{JNIEnv}, Ptr{Nothing}), penv, a)
        @assert(sz_a == sz_1, "Size of $(i)th subrarray is $sz_a, but size of the 1st subarray was $sz_1")
        for j=1:sz_1
            x = ccall(jnifunc.GetObjectArrayElement, Ptr{Nothing}, (Ptr{JNIEnv},Ptr{Nothing}, jint),
                      penv, a, j-1)
            ret[i, j] = JavaObject{T}(x)
        end
    end
    return ret
end

# matrices of primitive types and other arrays
function convert_result(rettype::Type{Array{T,2}}, result) where T
    sz = ccall(jnifunc.GetArrayLength, jint, (Ptr{JNIEnv}, Ptr{Nothing}), penv, result)
    if sz == 0
        return Array{T}(undef, 0,0)
    end
    a_1 = ccall(jnifunc.GetObjectArrayElement, Ptr{Nothing}, (Ptr{JNIEnv},Ptr{Nothing}, jint), penv,
                result, 0)
    sz_1 = ccall(jnifunc.GetArrayLength, jint, (Ptr{JNIEnv}, Ptr{Nothing}), penv, a_1)
    ret = Array{T}(undef, sz, sz_1)
    for i=1:sz
        a = ccall(jnifunc.GetObjectArrayElement, Ptr{Nothing}, (Ptr{JNIEnv},Ptr{Nothing}, jint), penv,
                  result, i-1)
        # check that size of the current subarray is the same as for the first one
        sz_a = ccall(jnifunc.GetArrayLength, jint, (Ptr{JNIEnv}, Ptr{Nothing}), penv, a)
        @assert(sz_a == sz_1, "Size of $(i)th subrarray is $sz_a, but size of the 1st subarray was $sz_1")
        ret[i, :] = convert_result(Vector{T}, a)
    end
    return ret
end

convert(::Type{jlong}, obj::JavaObject{Symbol("java.lang.Long")}) = jcall(obj, "longValue", jlong, ())
convert(::Type{jint}, obj::JavaObject{Symbol("java.lang.Integer")}) = jcall(obj, "intValue", jint, ())
convert(::Type{jdouble}, obj::JavaObject{Symbol("java.lang.Double")}) = jcall(obj, "doubleValue", jdouble, ())
convert(::Type{jfloat}, obj::JavaObject{Symbol("java.lang.Float")}) = jcall(obj, "floatValue", jfloat, ())
convert(::Type{jboolean}, obj::JavaObject{Symbol("java.lang.Boolean")}) = jcall(obj, "booleanValue", jboolean, ())

#The second term in this addition is due to the fact that Java converts all times to local time
function convert(::Type{DateTime}, x::@jimport(java.util.Date))
    if isnull(x)
        Dates.DateTime(1970,1,1,0,0,0)
    else
        Dates.unix2datetime(jcall(x, "getTime", jlong, ())/1000) +
            Second(round(div(Dates.value(now() - now(Dates.UTC)),1000)/900)*(900))
    end
end

function convert(::Type{DateTime}, x::JavaObject)
    isnull(x) && return Dates.DateTime(1970,1,1,0,0,0)
    JDate = @jimport(java.util.Date)
    if isConvertible(JDate, x)
        return convert(DateTime, convert(JDate, x))
    elseif isConvertible(@jimport(java.util.Calendar), x)
        return convert(DateTime, jcall(x, "getTime", JDate, ()))
    end
end

function convert(::Type{@jimport(java.util.Properties)}, x::Dict)
    Properties = @jimport(java.util.Properties)
    p = Properties(())
    for (n,v) in x
        jcall(p, "setProperty", @jimport(java.lang.Object), (JString, JString), n, v)
    end
    return p
end

function convert(::Type{@jimport(java.util.HashMap)}, K::Type{JavaObject{X}}, V::Type{JavaObject{Y}},
                 x::Dict) where {X,Y}
    Hashmap = @jimport(java.util.HashMap)
    p = Hashmap(())
    for (n,v) in x
        jcall(p, "put", @jimport(java.lang.Object), (JObject, JObject), n, v)
    end
    return p
end

function convert(::Type{@jimport(java.util.Map)}, K::Type{JavaObject{X}}, V::Type{JavaObject{Y}},
                 x::Dict) where {X,Y}
    convert(@jimport(java.util.Map), convert(@jimport(java.util.HashMap), K, V, x))
end

function convert(::Type{@jimport(java.util.ArrayList)}, x::Vector, V::Type{JavaObject{X}}=JObject) where X
    ArrayList = @jimport(java.util.ArrayList)
    a = ArrayList(())
    for v in x
        jcall(a, "add", jboolean, (JObject,), convert(V, v))
    end
    return a
end

function convert(::Type{@jimport(java.util.List)}, x::Vector, V::Type{JavaObject{X}}=JObject) where X
    convert(@jimport(java.util.ArrayList), x, V)
end

# Convert a reference to a java.lang.String into a Julia string. Copies the underlying byte buffer
function unsafe_string(jstr::JString)  #jstr must be a jstring obtained via a JNI call
    if isnull(jstr); return ""; end #Return empty string to keep type stability. But this is questionable
    pIsCopy = Array{jboolean}(undef, 1)
    buf::Ptr{UInt8} = ccall(jnifunc.GetStringUTFChars, Ptr{UInt8},
                            (Ptr{JNIEnv}, Ptr{Nothing}, Ptr{jboolean}), penv, jstr.ptr, pIsCopy)
    s = unsafe_string(buf)
    ccall(jnifunc.ReleaseStringUTFChars, Nothing, (Ptr{JNIEnv}, Ptr{Nothing}, Ptr{UInt8}), penv,
          jstr.ptr, buf)
    return s
end

for (x, y, z) in [(:jboolean, :(jnifunc.GetBooleanArrayElements), :(jnifunc.ReleaseBooleanArrayElements)),
                  (:jchar, :(jnifunc.GetCharArrayElements), :(jnifunc.ReleaseCharArrayElements)),
                  (:jbyte, :(jnifunc.GetByteArrayElements), :(jnifunc.ReleaseByteArrayElements)),
                  (:jshort, :(jnifunc.GetShortArrayElements), :(jnifunc.ReleaseShortArrayElements)),
                  (:jint, :(jnifunc.GetIntArrayElements), :(jnifunc.ReleaseIntArrayElements)),
                  (:jlong, :(jnifunc.GetLongArrayElements), :(jnifunc.ReleaseLongArrayElements)),
                  (:jfloat, :(jnifunc.GetFloatArrayElements), :(jnifunc.ReleaseFloatArrayElements)),
                  (:jdouble, :(jnifunc.GetDoubleArrayElements), :(jnifunc.ReleaseDoubleArrayElements)) ]
    m = quote
        function convert(::Type{Array{$(x),1}}, obj::JObject)
            sz = ccall(jnifunc.GetArrayLength, jint, (Ptr{JNIEnv}, Ptr{Nothing}), penv, obj.ptr)
            arr = ccall($(y), Ptr{$(x)}, (Ptr{JNIEnv}, Ptr{Nothing}, Ptr{jboolean}), penv, obj.ptr,
                        C_NULL)
            jl_arr::Array = unsafe_wrap(Array, arr, Int(sz))
            jl_arr = deepcopy(jl_arr)
            ccall($(z), Nothing, (Ptr{JNIEnv},Ptr{Nothing}, Ptr{$(x)}, jint), penv, obj.ptr, arr, 0)
            return jl_arr
        end
    end
    eval(m)
end

function convert(::Type{Array{T, 1}}, obj::JObject) where T
    sz = ccall(jnifunc.GetArrayLength, jint,
               (Ptr{JNIEnv}, Ptr{Nothing}), penv, obj.ptr)
    ret = Array{T}(undef, sz)
    for i=1:sz
        ptr = ccall(jnifunc.GetObjectArrayElement, Ptr{Nothing},
                  (Ptr{JNIEnv}, Ptr{Nothing}, jint), penv, obj.ptr, i-1)
        ret[i] = convert(T, JObject(ptr))
    end
    return ret
end

##Iterator
iterator(obj::JavaObject) = jcall(obj, "iterator", @jimport(java.util.Iterator), ())

"""
Given a `JavaObject{T}` narrows down `T` to a real class of the underlying object.
For example, `JavaObject{:java.lang.Object}` pointing to `java.lang.String`
will be narrowed down to `JavaObject{:java.lang.String}`
"""
function narrow(obj::JavaObject)
    c = jcall(obj,"getClass", @jimport(java.lang.Class), ())
    t = jcall(c, "getName", JString, ())
    return convert(JavaObject{Symbol(t)}, obj)
end

has_next(itr::JavaObject) = (jcall(itr, "hasNext", jboolean, ()) == JNI_TRUE)

function Base.iterate(itr::JavaObject, state=nothing)
    if has_next(itr)
        o = jcall(itr, "next", @jimport(java.lang.Object), ())
        return (narrow(o), state)
    else
        return nothing
    end
end


text:
using Measurements

d(a_in, a_out) = (1 - a_in*a_in) / (1 - a_in*a_out) * a_out - a_in
d̃(a_in, a_out) = a_out - a_in

# Au(111)
# a_in, a_out...

# a = [ 0.0   -0.31
#      -0.22  -0.35
#       0.22  -0.29 ]

# σ = [ 0.05  0.04
#       0.05  0.03
#       0.05  0.03 ]

# Polycrystalline
# a_in, a_out...

a = [
     0.020  -0.380  -0.355  -0.572
     0.063  -0.351  -0.315  -0.545
    -0.041  -0.401  -0.388  -0.608
]

σ = [
    0.058  0.065  0.055  0.059
    0.044  0.083  0.057  0.070
    0.063  0.055  0.059  0.058
]

ā = a .± σ

D = d.(ā[:, 1], ā[:, 2:end])
D̃ = d̃.(ā[:, 1], ā[:, 2:end])

Δ = D[1, :] .- (D[2, :] .+ D[3, :]) ./ 2
Δ̃ = D̃[1, :] .- (D̃[2, :] .+ D̃[3, :]) ./ 2


text:
#!/usr/bin/env julia --project
using ClimateMachine
using ArgParse
using UnPack

s = ArgParseSettings()
@add_arg_table! s begin
    "--number-of-tracers"
    help = "Number of dummy tracers"
    metavar = "<number>"
    arg_type = Int
    default = 0
end

parsed_args = ClimateMachine.init(parse_clargs = true, custom_clargs = s)
const number_of_tracers = parsed_args["number-of-tracers"]

using ClimateMachine.Atmos
using ClimateMachine.Orientations
using ClimateMachine.ConfigTypes
using ClimateMachine.Diagnostics
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.StdDiagnostics
using ClimateMachine.SystemSolvers: ManyColumnLU
using ClimateMachine.TurbulenceClosures
using ClimateMachine.Mesh.Filters
using ClimateMachine.Mesh.Grids
using ClimateMachine.TemperatureProfiles
using ClimateMachine.Thermodynamics:
    air_pressure, air_density, air_temperature, total_energy, internal_energy
using ClimateMachine.VariableTemplates

using ClimateMachine.BalanceLaws
import ClimateMachine.BalanceLaws: source, prognostic_vars

using LinearAlgebra
using StaticArrays
using Test

using CLIMAParameters
using CLIMAParameters.Planet:
    MSLP, R_d, day, cp_d, cv_d, grav, Omega, planet_radius
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

function init_heldsuarez!(problem, bl, state, aux, localgeo, t)
    FT = eltype(state)

    param_set = parameter_set(bl)
    # parameters
    _a::FT = planet_radius(param_set)

    z_t::FT = 15e3
    λ_c::FT = π / 9
    φ_c::FT = 2 * π / 9
    d_0::FT = _a / 6
    V_p::FT = 10

    # grid
    φ = latitude(bl.orientation, aux)
    λ = longitude(bl.orientation, aux)
    z = altitude(bl.orientation, param_set, aux)

    # deterministic velocity perturbation
    F_z::FT = 1 - 3 * (z / z_t)^2 + 2 * (z / z_t)^3
    if z > z_t
        F_z = FT(0)
    end
    d::FT = _a * acos(sin(φ) * sin(φ_c) + cos(φ) * cos(φ_c) * cos(λ - λ_c))
    c3::FT = cos(π * d / 2 / d_0)^3
    s1::FT = sin(π * d / 2 / d_0)
    if 0 < d < d_0 && d != FT(_a * π)
        u′::FT =
            -16 * V_p / 3 / sqrt(3) *
            F_z *
            c3 *
            s1 *
            (-sin(φ_c) * cos(φ) + cos(φ_c) * sin(φ) * cos(λ - λ_c)) /
            sin(d / _a)
        v′::FT =
            16 * V_p / 3 / sqrt(3) * F_z * c3 * s1 * cos(φ_c) * sin(λ - λ_c) /
            sin(d / _a)
    else
        u′ = FT(0)
        v′ = FT(0)
    end
    w′::FT = 0
    u_sphere = SVector{3, FT}(u′, v′, w′)
    u_cart = sphr_to_cart_vec(bl.orientation, u_sphere, aux)

    ## potential & kinetic energy
    e_kin::FT = 0.5 * u_cart' * u_cart

    ## Assign state variables
    state.ρ = aux.ref_state.ρ
    state.ρu = state.ρ * u_cart
    state.energy.ρe = aux.ref_state.ρe + state.ρ * e_kin
    if number_of_tracers > 0
        state.tracers.ρχ = @SVector [FT(ii) for ii in 1:number_of_tracers]
    end

    nothing
end

"""
    HeldSuarezForcing <: TendencyDef{Source}

Defines a forcing that parametrises radiative and frictional effects using
Newtonian relaxation and Rayleigh friction, following Held and Suarez (1994)
"""
struct HeldSuarezForcing <: TendencyDef{Source} end

prognostic_vars(::HeldSuarezForcing) = (Momentum(), Energy())

function held_suarez_forcing_coefficients(bl, args)
    @unpack state, aux = args
    @unpack ts = args.precomputed
    FT = eltype(state)

    # Parameters
    T_ref = FT(255)
    param_set = parameter_set(bl)
    _R_d = FT(R_d(param_set))
    _day = FT(day(param_set))
    _grav = FT(grav(param_set))
    _cp_d = FT(cp_d(param_set))
    _p0 = FT(MSLP(param_set))

    # Held-Suarez parameters
    k_a = FT(1 / (40 * _day))
    k_f = FT(1 / _day)
    k_s = FT(1 / (4 * _day))
    ΔT_y = FT(60)
    Δθ_z = FT(10)
    T_equator = FT(315)
    T_min = FT(200)
    σ_b = FT(7 / 10)

    # Held-Suarez forcing
    φ = latitude(bl, aux)
    p = air_pressure(ts)

    #TODO: replace _p0 with dynamic surface pressure in Δσ calculations to account
    #for topography, but leave unchanged for calculations of σ involved in T_equil
    σ = p / _p0
    exner_p = σ^(_R_d / _cp_d)
    Δσ = (σ - σ_b) / (1 - σ_b)
    height_factor = max(0, Δσ)
    T_equil = (T_equator - ΔT_y * sin(φ)^2 - Δθ_z * log(σ) * cos(φ)^2) * exner_p
    T_equil = max(T_min, T_equil)
    k_T = k_a + (k_s - k_a) * height_factor * cos(φ)^4
    k_v = k_f * height_factor
    return (k_v = k_v, k_T = k_T, T_equil = T_equil)
end

function source(::Energy, s::HeldSuarezForcing, m, args)
    @unpack state = args
    @unpack ts = args.precomputed
    nt = held_suarez_forcing_coefficients(m, args)
    FT = eltype(state)
    param_set = parameter_set(m)
    _cv_d = FT(cv_d(param_set))
    @unpack k_T, T_equil = nt
    T = air_temperature(ts)
    return -k_T * state.ρ * _cv_d * (T - T_equil)
end

function source(::Momentum, s::HeldSuarezForcing, m, args)
    nt = held_suarez_forcing_coefficients(m, args)
    return -nt.k_v * projection_tangential(m, args.aux, args.state.ρu)
end

function config_heldsuarez(FT, poly_order, resolution)
    # Set up a reference state for linearization of equations
    temp_profile_ref =
        DecayingTemperatureProfile{FT}(param_set, FT(290), FT(220), FT(8e3))
    ref_state = HydrostaticState(temp_profile_ref)

    # Set up the atmosphere model
    exp_name = "HeldSuarez"
    domain_height::FT = 30e3 # distance between surface and top of atmosphere (m)

    if number_of_tracers > 0
        δ_χ = @SVector [FT(ii) for ii in 1:number_of_tracers]
        tracers = NTracers{number_of_tracers, FT}(δ_χ)
    else
        tracers = NoTracers()
    end

    physics = AtmosPhysics{FT}(
        param_set;
        ref_state = ref_state,
        turbulence = ConstantKinematicViscosity(FT(0)),
        hyperdiffusion = DryBiharmonic(FT(8 * 3600)),
        moisture = DryModel(),
        tracers = tracers,
    )

    model = AtmosModel{FT}(
        AtmosGCMConfigType,
        physics;
        init_state_prognostic = init_heldsuarez!,
        source = (Gravity(), Coriolis(), HeldSuarezForcing()),
    )

    config = ClimateMachine.AtmosGCMConfiguration(
        exp_name,
        poly_order,
        resolution,
        domain_height,
        param_set,
        init_heldsuarez!;
        model = model,
    )

    return config
end

function main()
    # Driver configuration parameters
    FT = Float64                             # floating type precision
    poly_order = 5                           # discontinuous Galerkin polynomial order
    n_horz = 8                              # horizontal element number
    n_vert = 4                               # vertical element number
    n_days::FT = 1
    timestart::FT = 0                        # start time (s)
    timeend::FT = n_days * day(param_set)    # end time (s)

    # Set up driver configuration
    driver_config = config_heldsuarez(FT, poly_order, (n_horz, n_vert))

    # Set up experiment
    ode_solver_type = ClimateMachine.IMEXSolverType(
        implicit_model = AtmosAcousticGravityLinearModel,
        implicit_solver = ManyColumnLU,
        solver_method = ARK2GiraldoKellyConstantinescu,
        split_explicit_implicit = true,
        discrete_splitting = false,
    )

    CFL = FT(0.1) # target acoustic CFL number

    # time step is computed such that the horizontal acoustic Courant number is CFL
    solver_config = ClimateMachine.SolverConfiguration(
        timestart,
        timeend,
        driver_config,
        Courant_number = CFL,
        ode_solver_type = ode_solver_type,
        CFL_direction = HorizontalDirection(),
        diffdir = HorizontalDirection(),
    )

    # Set up diagnostics
    dgn_ssecs = cld(timeend, 2) + 30
    dgn_interval = "$(dgn_ssecs)ssecs"
    dgn_config = config_diagnostics(FT, driver_config, dgn_interval)

    # Set up user-defined callbacks
    filterorder = 20
    filter = ExponentialFilter(solver_config.dg.grid, 0, filterorder)
    cbfilter = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(
            solver_config.Q,
            AtmosFilterPerturbations(driver_config.bl),
            solver_config.dg.grid,
            filter,
            state_auxiliary = solver_config.dg.state_auxiliary,
        )
        nothing
    end

    # Run the model
    result = ClimateMachine.invoke!(
        solver_config;
        diagnostics_config = dgn_config,
        user_callbacks = (cbfilter,),
        check_euclidean_distance = true,
    )
end

function config_diagnostics(FT, driver_config, interval)
    _planet_radius = FT(planet_radius(param_set))

    info = driver_config.config_info
    boundaries = [
        FT(-90.0) FT(-180.0) _planet_radius
        FT(90.0) FT(180.0) FT(_planet_radius + info.domain_height)
    ]
    resolution = (FT(1), FT(1), FT(1000)) # in (deg, deg, m)
    interpol = ClimateMachine.InterpolationConfiguration(
        driver_config,
        boundaries,
        resolution,
    )

    dgngrp = StdDiagnostics.AtmosGCMDefault(
        interval,
        driver_config.name,
        interpol = interpol,
    )

    return ClimateMachine.DiagnosticsConfiguration([dgngrp])
end

main()


text:
module GL

const lib = "libGL"

################################################################################
#   Types
################################################################################

include("type.jl")

################################################################################
#   Shared Enums
################################################################################

# Program and Shader GetInfoLog()
const DELETE_STATUS   = 0x8B80
const INFO_LOG_LENGTH = 0x8B84

# Data Types
const BYTE           = 0x1400
const UNSIGNED_BYTE  = 0x1401
const SHORT          = 0x1402
const UNSIGNED_SHORT = 0x1403
const INT            = 0x1404
const UNSIGNED_INT   = 0x1405
const FLOAT          = 0x1406
const DOUBLE         = 0x140A

const GLtype = [
	Int8    => BYTE,
	Uint8   => UNSIGNED_BYTE,
	Int16   => SHORT,
	Uint16  => UNSIGNED_SHORT,
	Int32   => INT,
	Uint32  => UNSIGNED_INT,
	Float32 => FLOAT,
	Float64 => DOUBLE
]

################################################################################
#   Includes
################################################################################

include("attribute.jl")
include("buffer.jl")
include("draw.jl")
include("program.jl")
include("shader.jl")
include("texture.jl")
include("uniform.jl")
include("utils.jl")
include("vertexarray.jl")

end # module