text:
#=
solves:
- Julia version:
- Author: qd4314
- Date: 2019-03-21
=#

function linFunction(phiVector::AbstractVector,KS::KineticSolver)

    # reshape vector input to tensor
    phi = reshape(phiVector,(KS.DecompNumberElements,KS.nx,KS.ny))

    OSigma = KS.OSigma
    M = KS.DecompM

    # compute sweeping rhs which is O*Sigma*phi
    OSigmaPhi = zeros(KS.nquadpoints,KS.nx,KS.ny)
    @einsum OSigmaPhi[q,i,j] = OSigma[q,moments] * phi[moments,i,j];

    # multiply with sigmaS
    @inbounds for j=KS.rangex,i=KS.rangey
        OSigmaPhi[:,i,j]=OSigmaPhi[:,i,j].*KS.SigmaS[i,j];
    end

    # perform sweeping
    Y= sweeping(KS,
				deepcopy(KS.psiC),
                OSigmaPhi,
                KS.ScatteringKernelConvolution);

    #Y= sweeping(KS,0.0*deepcopy(KS.psiC),OSigmaPhi,KS.ScatteringKernelConvolution);

    KS.psiC .= Y;

    # multiply with M
    MY = zeros(KS.DecompNumberElements,KS.nx,KS.ny)
    @einsum MY[m,i,j] = M[m,q]*Y[q,i,j]

    # substract from phi to get A*phi output
    Aphi = phi .- MY

    # reshape for vector output
    AphiVector = reshape(Aphi,prod((KS.DecompNumberElements,KS.nx,KS.ny)))

    return AphiVector;

end

function corysolve!(KS::KineticSolver)

    ### DEBUGGING
    # turn off scattering and absorption
    #KS.SigmaS = zeros(size(KS.SigmaS));
    #KS.SigmaT = zeros(size(KS.SigmaT));
    #KS.convolutionmagnitude = 0.0;
	close("all")
	io = open("mytestfile.txt", "a")

	writedlm(io,"start!!!!!!!!")
	close(io)
    # set convolution magnitude

    # set max number of sweeping steps
    KS.nSweepingSteps = 1000 # 1 is one-shot
    KS.epsSweeping = 1e-3#1e-7;

    # update SigmaT if convolution is active

    if KS.ictype == "linesource"
        println("Problem Type is ",KS.ictype)
        sigmaA = 0.0; # hard coded for linesource
    else
        sigmaA = KS.SigmaT - KS.SigmaS;
    end
    KS.SigmaT = KS.SigmaS+sigmaA.*ones(size(KS.SigmaT))+KS.convolutionmagnitude*ones(size(KS.SigmaT))

    # set estimate of contraction number
    L = KS.convolutionmagnitude / (1/KS.dt + minimum(KS.SigmaT)+sqrt(2)/KS.dx);
    KS.contractionEstimate = L/(1-L);

    # set moment order
    order = 0
    KS.DecompNumberElements = (order+1)*(order+1)

    # define total number of unknowns
    Nentries = KS.DecompNumberElements*KS.ny*KS.nx

    # decomposition order of kernel is Norder = (order+1)*(order+1)+1
    KS.DecompO,KS.DecompSigma, KS.DecompM = kerneldecomposition(KS,order) # Henyey Grenstein
    # precompute convolution matrix if needed
    if KS.convolutionmagnitude > 0
        KS.ScatteringKernelConvolution = computeScatteringKernelConvolution(KS);
    else
        KS.ScatteringKernelConvolution = zeros(KS.nquadpoints,KS.nquadpoints)
    end# precompute O*Sigma
    KS.OSigma = KS.DecompO*KS.DecompSigma;

    # write Decomposition matrices for physical scattering
    M = KS.DecompM
    Sigma = KS.DecompSigma
    O = KS.DecompO
    OSigma = KS.OSigma

    # Set initial condition
    psi = deepcopy(KS.psi0)
    psi0 = deepcopy(KS.psi0)
    KS.lastpsi = deepcopy(KS.psi0)
    KS.psiC = deepcopy(KS.psi0)

    phiNew = zeros(KS.DecompNumberElements,KS.ny,KS.nx)
    @einsum phiNew[m,i,j] = M[m,q]*psi[q,i,j]

    rho = computeMoments(KS,psi)
    println("Initial mass = $(sum(rho[:]))." )

    # define scattering kernel for source iteration
    S = deepcopy(KS.ScatteringKernelConvolution)
    for q = 1:size(S)[1]
        S[q,:] .*= KS.Q.weights[:]
    end

    KS.kernelTimesWeightsSparse = sparse(S);

    @showprogress 0.01 "Computing... "  for k=1:KS.nt

        # adapt dt if k*dt would be larger than tEnd
        dt = ifelse(k*KS.dt<KS.tEnd, KS.dt,KS.tEnd-k*KS.dt)

        # compute rhs
        rhsquads = sweeping(KS,deepcopy(KS.lastpsi),KS.lastpsi ./ KS.dt+KS.Source,KS.ScatteringKernelConvolution)

        rhsmoms = zeros(KS.DecompNumberElements,KS.ny,KS.nx)
        @einsum rhsmoms[m,i,j] = M[m,q]*rhsquads[q,i,j]

        ## reshape to vector
        rrhsmoms = reshape(rhsmoms,Nentries)

        # define linear Map for Krylov solver
        D = LinearMap((x)->linFunction(x,KS), Nentries,Nentries; ismutating=false) # muss in jeder Zeit iteration neu definiert werden, da KS sich aendert!

        # initial guess for Krylov solver
        #x0 = zeros(KS.DecompNumberElements,KS.nx,KS.ny)
        #@einsum x0[m,i,j] = M[m,q]*psi0[q,i,j]
        #x0 = reshape(x0,Nentries)

        # solve system with GMRES
        #println("############ starting krylov")
        println(sqrt(eps(real(eltype(rrhsmoms)))))
        phiNew = IterativeSolvers.gmres!(rrhsmoms,D,rrhsmoms,verbose=false)
		#println("done krylov")

        #phiNew = IterativeSolvers.gmres!(x0,D,rrhsmoms,verbose=false)
        phiNew = reshape(phiNew,(KS.DecompNumberElements,KS.nx,KS.ny))

        #visMoments(KS,phiNew,"OutputGMRES")

        # determine PsiNew
        OSigmaPhi = zeros(KS.nquadpoints,KS.nx,KS.ny)
        @einsum OSigmaPhi[q,i,j] = OSigma[q,moments] * phiNew[moments,i,j];

        # multiply with sigmaS
        @inbounds for j=KS.rangex,i=KS.rangey
            OSigmaPhi[:,i,j]=OSigmaPhi[:,i,j].*KS.SigmaS[i,j];
        end

        # perform sweeping to update psi
        psiNew = sweeping(KS,deepcopy(KS.lastpsi),KS.Source .+ KS.lastpsi / KS.dt+OSigmaPhi,KS.ScatteringKernelConvolution);
        KS.lastpsi .= psiNew
    end

    rho = computeMoments(KS,KS.lastpsi)
    println("Final mass = $(sum(rho[:]))." )

    vis(KS,KS.lastpsi)
	return rho
end

function solveWithSweeping!(KS::KineticSolver)

    # set max number of sweeping steps
    KS.nSweepingSteps = 1000
    KS.epsSweeping = 1e-3#1e-7;

    # update SigmaT if convolution is active
    if KS.ictype == "linesource"
        println("Problem Type is ",KS.ictype)
        sigmaA = 0.0; # hard coded for linesource
    else
        sigmaA = KS.SigmaT - KS.SigmaS;
    end
    KS.SigmaT = KS.SigmaS+sigmaA.*ones(size(KS.SigmaT))+KS.convolutionmagnitude*ones(size(KS.SigmaT))

    # set estimate of contraction number
    L = (KS.convolutionmagnitude+maximum(KS.SigmaS)) / (1/KS.dt + minimum(KS.SigmaT)+sqrt(2)/KS.dx);
    KS.contractionEstimate = L/(1-L);

    # compute artificial scattering kernel
    KS.ScatteringKernelConvolution = computeScatteringKernelConvolution(KS); # Convo

    # Set initial condition
    psi = zeros(KS.nquadpoints,KS.ny,KS.nx)
    psi = deepcopy(KS.psi0)
    KS.lastpsi = deepcopy(KS.psi0)
    rho = computeMoments(KS,psi)
    println("Initial mass = $(sum(rho[:]))." )

    nicewrite(string(KS.outputfolderprefix,"data/source.txt"),computeMoments(KS,KS.Source))

    # define scattering kernel for source iteration
    asKernel = deepcopy(KS.ScatteringKernelConvolution) # artificial inscattering
    physKernel =  ones(size(asKernel)) ./ 4.0 / pi; # add physical inscattering
    g = 0.01
    HG = computeScatteringKernelHenyeyGreenstein(KS,g)
    # decomposition order of kernel is Norder = (order+1)*(order+1)+1
    KS.DecompO,KS.DecompSigma, KS.DecompM = kerneldecomposition(KS,1) # Henyey Grenstein

    KS.kernelTimesWeightsSparse = sparse(asKernel);

    #println("HG = ",HG)
    #println("OSigmaM = ",KS.DecompO*KS.DecompSigma*KS.DecompM)
    #println("wQ = ",KS.Q.weights[:])
    for q = 1:size(physKernel)[1]
        asKernel[q,:] .*= KS.Q.weights[:];
        physKernel[q,:] .*= KS.Q.weights[:];
    end
    physKernel = KS.DecompO*KS.DecompSigma*KS.DecompM;

    @showprogress 0.01 "Computing... "  for k=1:KS.nt

        # adapt dt if k*dt would be larger than tEnd
        dt = ifelse(k*KS.dt<KS.tEnd, KS.dt,KS.tEnd-k*KS.dt)

        # first euler step
        psiNew = sweepingWithKernels(KS,deepcopy(KS.lastpsi),KS.Source .+ KS.lastpsi / KS.dt,asKernel,physKernel);
        KS.lastpsi .= psiNew
        #vis(KS,psi)
        #return;
    end
    vis(KS,KS.lastpsi) # do all post computational plotting
    rho = computeMoments(KS,KS.lastpsi)
    println("End mass = $(sum(rho[:]))." )
end # solve

function solveRSNWithSweeping!(KS::KineticSolver)

    # set max number of sweeping steps
    KS.nSweepingSteps = 1000
    KS.epsSweeping = 1e-3#1e-7;

    # update SigmaT if convolution is active
    if KS.ictype == "linesource"
        println("Problem Type is ",KS.ictype)
        sigmaA = 0.0; # hard coded for linesource
    else
        sigmaA = KS.SigmaT - KS.SigmaS;
    end
    KS.SigmaT = KS.SigmaS+sigmaA.*ones(size(KS.SigmaT))+KS.convolutionmagnitude*ones(size(KS.SigmaT))

    # set estimate of contraction number
    L = (KS.convolutionmagnitude+maximum(KS.SigmaS)) / (1/KS.dt + minimum(KS.SigmaT)+sqrt(2)/KS.dx);
    KS.contractionEstimate = L/(1-L);

    # compute artificial scattering kernel
    KS.ScatteringKernelConvolution = computeScatteringKernelConvolution(KS); # Convo

    # Set initial condition
    psi = zeros(KS.nquadpoints,KS.ny,KS.nx)
    psi = deepcopy(KS.psi0)
    KS.lastpsi = deepcopy(KS.psi0)
    rho = computeMoments(KS,psi)
    println("Initial mass = $(sum(rho[:]))." )

    nicewrite(string(KS.outputfolderprefix,"data/source.txt"),computeMoments(KS,KS.Source))

    # define scattering kernel for source iteration
    asKernel = deepcopy(KS.ScatteringKernelConvolution) # artificial inscattering
    physKernel =  ones(size(asKernel)) ./ 4.0 / pi; # add physical inscattering
    g = 0.01
    HG = computeScatteringKernelHenyeyGreenstein(KS,g)
    # decomposition order of kernel is Norder = (order+1)*(order+1)+1
    KS.DecompO,KS.DecompSigma, KS.DecompM = kerneldecomposition(KS,1) # Henyey Grenstein

    KS.kernelTimesWeightsSparse = sparse(asKernel);

    #println("HG = ",HG)
    #println("OSigmaM = ",KS.DecompO*KS.DecompSigma*KS.DecompM)
    #println("wQ = ",KS.Q.weights[:])
    for q = 1:size(physKernel)[1]
        asKernel[q,:] .*= KS.Q.weights[:];
        physKernel[q,:] .*= KS.Q.weights[:];
    end
    physKernel = KS.DecompO*KS.DecompSigma*KS.DecompM;

    @showprogress 0.01 "Computing... "  for k=1:KS.nt

        # adapt dt if k*dt would be larger than tEnd
        dt = ifelse(k*KS.dt<KS.tEnd, KS.dt,KS.tEnd-k*KS.dt)

        # first euler step
        psiNew = sweepingWithKernels(KS,deepcopy(KS.lastpsi),KS.Source .+ KS.lastpsi / KS.dt,asKernel,physKernel);
        KS.lastpsi .= psiNew

        # perform rotation and interpolation step, or convolution step, or other modifications
        methods!(KS,KS.lastpsi,psiNew,k)
        #vis(KS,psi)
        #return;
    end
    vis(KS,KS.lastpsi) # do all post computational plotting
    rho = computeMoments(KS,KS.lastpsi)
    println("End mass = $(sum(rho[:]))." )
end # solve

function solve!(KS::KineticSolver)
    #close("all")
    psi = zeros(KS.nquadpoints,KS.ny,KS.nx)
    psiLastTimeStep = zeros(KS.nquadpoints,KS.ny,KS.nx)
    flux = zeros(KS.nquadpoints,KS.ny,KS.nx)
    psi = deepcopy(KS.psi0) # Set initial condition
    rho = computeMoments(KS,psi)
    println("Initial mass = $(sum(rho[:]))." )

    #nicewrite(string(KS.outputfolderprefix,"data/rho0.txt"),computeMoments(KS,psi))
    nicewrite(string(KS.outputfolderprefix,"data/source.txt"),computeMoments(KS,KS.Source))

    summarizeAsPlot(KS,psi); #close("all")

    # precompute convolution matrix if needed
    if KS.convolutionmagnitude > 0
        KS.ScatteringKernelConvolution = computeScatteringKernelConvolution(KS);
    else
        KS.ScatteringKernelConvolution = zeros(KS.nquadpoints,KS.nquadpoints)
	end

		# low rank computation?
    if KS.rank > 0
        dA = zeros(KS.nquadpoints,(KS.ny) * (KS.nx))
        psireshaped = reshape(psi,KS.nquadpoints,KS.ny*KS.nx)
        s0 = svd(psireshaped);
        U0 = s0.U[1:end, 1:KS.rank];
        S0 = diagm(0 => s0.S[1:KS.rank] );
        V0 = s0.V[1:end, 1:KS.rank];
        U1 = zeros(size(U0));    V1 = zeros(size(V0));        S1 = zeros(size(S0))
        K1 = zeros(size(U0*S0)); L1 = zeros(size(V0*S1')); Y0 = zeros(KS.nquadpoints,(KS.ny) * (KS.nx))
    end

    @showprogress 0.01 "Computing... "  for k=1:KS.nt

        # force periodicity (if wanted)
        forcePeriodicity!(KS,psi)

        # adapt dt if k*dt would be larger than tEnd
        dt = ifelse(k*KS.dt<KS.tEnd, KS.dt,KS.tEnd-k*KS.dt)

        # store old psi to average later on
        #@inbounds for q=1:length(psi);	psiLastTimeStep[q] = psi[q]	end
        psiLastTimeStep .= psi

        # first euler step
        solveFlux!(KS,psi,flux); Euler!(KS,psi,flux,dt)

        # second euler step
        solveFlux!(KS,psi,flux); Euler!(KS,psi,flux,dt)

        #average
        psi .= (psi .+ psiLastTimeStep) ./ 2.0

        # add source contribution
        psi .+= KS.Source .* KS.dt

        # Note: So far, low rank and rSN DO NOT WORK TOGETHER. Either lowrank or rSN.
        if KS.rank>0
            @inbounds for q=1:length(psi); dA[q] = psi[q]-psiLastTimeStep[q] end
            K1 = U0*S0 + dA*V0
            U1,S1 = qr(K1);            U1 = Matrix(U1)
            S1 = S1-U1'*dA*V0
            L1 = V0*S1'+dA'*U1;
            V1,S1 = qr(L1);            V1 = Matrix(V1)
            S1 = S1';
            Y1 = U1*S1*V1';
            for i=1:length(Y0),     Y0[i] = Y1[i];     end
            for i=1:length(U0),     U0[i] = U1[i];     end
            for i=1:length(V0),     V0[i] = V1[i];     end
            for i=1:length(S0),  S0[i] = S1[i];  end
            for i=1:length(psi), psi[i] = Y0[i]; end
        else
            # perform rotation and interpolation step, or convolution step, or other modifications
            methods!(KS,psi,psiLastTimeStep,k)
        end
    end
    vis(KS,psi) # do all post computational plotting
end # solve

max_stars_repo_name: camminady/SN

max_issues_repo_name: camminady/SN

max_forks_repo_name: camminady/SN


text:
using Test
using CSVFiles
using DataFrames
using Distributions

regenerate = false # do a large MC run, to regenerate information needed for std. errors
samplesize = 1000 # normal MC sample size (takes ~5 seconds)
confidence = 2.576 # 99% CI by default; use 1.96 to apply a 95% CI, but expect more spurious errors

# Monte Carlo distribution information
# Filled in from a run with regenerate = true
information = Dict(
    :td => Dict(:transform => x -> log(x), :mu => 20.814327104215817, :sigma => 0.7007588482770571),
    :tpc => Dict(:transform => x -> x, :mu => 3.23416123428041e7, :sigma => 6.374797685068551e7),
    :tac => Dict(:transform => x -> log(x), :mu => 14.944386813804241, :sigma => 0.46207147966982504),
    :te => Dict(:transform => x -> log(x), :mu => 20.84214660784536, :sigma => 0.6995498063474193),
    :c_co2concentration => Dict(:transform => x -> x, :mu => 651134.5999874412, :sigma => 46363.450549201516),
    :ft => Dict(:transform => x -> x, :mu => 6.216834969325794, :sigma => 0.4011489912917528),
    :rt_g => Dict(:transform => x -> x, :mu => 4.366704203451731, :sigma => 1.219601225681128),
    :sealevel => Dict(:transform => x -> x, :mu => 3.026419413461323, :sigma => 1.1459316083150264)
)

output_path = joinpath(@__DIR__, "../output")
compare = DataFrame(load(joinpath(@__DIR__, "validationdata/PAGE2020montecarloquantiles.csv")))

if regenerate
    println("Regenerating MC distribution information")

    # Perform a large MC run and extract statistics for std. errors
    do_monte_carlo_runs(100_000, "RCP4.5 & SSP2", output_path)
    df = DataFrame(load(joinpath(output_path, "mimipagemontecarlooutput.csv")))

    for ii in 1:nrow(compare)
        name = Symbol(compare[ii, :Variable_Name])
        if kurtosis(df[!, name]) > 2.9 # exponential distribution

            # in PAGE 2009 :tpc was negative across all quantiles so we took log(-x)
            # for the transform, in PAGE2020 they are not all negative so this
            # has been removed but may need to be revisited if we see high kurtosis
            # for :tpc since we cann't take log(negative value)

            # if name == :tpc # negative across all quantiles
                # print("    :$name => Dict(:transform => x -> log(-x), :mu => $(mean(log(-df[df[name] .< 0, name]))), :sigma => $(std(log(-df[df[name] .< 0, name]))))")
            # else
                # print("    :$name => Dict(:transform => x -> log(x), :mu => $(mean(log(df[df[name] .> 0, name]))), :sigma => $(std(log(df[df[name] .> 0, name]))))")
            # end

            print("    :$name => Dict(:transform => x -> log(x), :mu => $(mean(log.(df[df[!, name] .> 0, name]))), :sigma => $(std(log.(df[df[!, name] .> 0, name]))))")
        else
            print("    :$name => Dict(:transform => x -> x, :mu => $(mean(df[!, name])), :sigma => $(std(df[!, name])))")
        end
        if ii != nrow(compare)
            println(",")
        end
    end
else
    println("Performing MC sample")
    # Perform a small MC run
    do_monte_carlo_runs(samplesize, "RCP4.5 & SSP2", output_path)
    df = DataFrame(load(joinpath(output_path, "mimipagemontecarlooutput.csv")))
end

# Compare all known quantiles
for ii in 1:nrow(compare)
    name = Symbol(compare[ii, :Variable_Name])
    transform = information[name][:transform]
    distribution = Normal(information[name][:mu], information[name][:sigma])
    for qval in [.05, .10, .25, .50, .75, .90, .95]
        estimated = transform(quantile(collect(Missings.skipmissing(df[!, name])), qval)) # perform transform *after* quantile, so captures effect of all values
        stderr = sqrt(qval * (1 - qval) / (samplesize * pdf(distribution, estimated)^2))

        expected = transform(compare[ii, Symbol("perc_$(trunc(Int, qval * 100))")])

        # println("$name x $qval: $estimated ≈ $expected rtol=$(ceil(confidence * stderr; digits = -trunc(Int, log10(stderr))))")
        @test estimated ≈ expected rtol = ceil(confidence * stderr; digits=-trunc(Int, log10(stderr)))
    end
end

##
## Used to produce testing files on May 21, 2021 with master
##

# using DataFrames
# using CSV

# # get quantiles
# output_path = joinpath(@__DIR__, "../output")
# do_monte_carlo_runs(100_000, "RCP4.5 & SSP2", output_path)
# df_outputs = DataFrame(load(joinpath(output_path, "mimipagemontecarlooutput.csv")))

# # build and save quantiles dataframe
# Variable_Names = [:td, :tpc, :tac, :te, :c_co2concentration, :ft, :rt_g, :sealevel]
# Quantiles_Placeholder = Array{Any}(fill(missing, length(Variable_Names)))
# df_quantiles = DataFrame(Variable_Name = Variable_Names,
#                             perc_5 = Quantiles_Placeholder,
#                             perc_10 = Quantiles_Placeholder,
#                             perc_25 = Quantiles_Placeholder,
#                             perc_50 = Quantiles_Placeholder,
#                             perc_75 = Quantiles_Placeholder,
#                             perc_90 = Quantiles_Placeholder,
#                             perc_95 = Quantiles_Placeholder)
# quantiles = [5, 10, 25, 50, 75, 90, 95]

# for (i, variable) in enumerate(Variable_Names)
#     _name = Symbol(df_quantiles[i, :Variable_Name])
#     for quant in quantiles
#         df_quantiles[i, Symbol("perc_$(trunc(Int, quant))")] = quantile(collect(Missings.skipmissing(df_outputs[!, _name])), quant/100)
#     end
# end

# CSV.write(joinpath(@__DIR__, "validationdata/PAGE2020montecarloquantiles.csv"), df_quantiles)

# # save a smaller run with n = 1000
# do_monte_carlo_runs(1000, "RCP4.5 & SSP2", output_path)
# df = DataFrame(load(joinpath(output_path, "mimipagemontecarlooutput.csv")))
# CSV.write(joinpath(@__DIR__, "validationdata/mimipagemontecarlooutput.csv"), df)

ext: jl


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

"""
List of locations of North Magnetic Pole since 1590.

This is provided by World Magnetic Model
(https://www.ngdc.noaa.gov/geomag/data/poles/NP.xy).
"""
const POLELATLONG =
    try
        let
            local polelatlong, rows, floattype, temp
            polelatlong = readdlm(joinpath(dirname(@__FILE__),
                                           "..", "deps", "NP.xy"))
            rows = size(polelatlong, 1)
            floattype = typeof(polelatlong[1])
            temp = Dict{floattype, Tuple{floattype, floattype}}()
            for i = 1:rows
                merge!(temp, Dict(polelatlong[2rows + i]=>
                                  (polelatlong[rows + i],
                                   polelatlong[i])))
            end
            temp
        end
    catch
        warn("""Could not find file NP.xy, you will not be able to use
"geo2mag" and "mag2geo" functions.  Build again the package with
    Pkg.build("AstroLib")
then restart Julia session in order to fix this problem.""")
    end

const AU = 1.495978707e11 # Astronomical unit in meters
const J2000 = jdcnv(2000, 01, 01, 12) # 2451545
const JULIANYEAR = 365.25 # days in Julian year

# Constant used in ct2lst function, see Meeus, p.84.
const ct2lst_c  = (280.46061837, 360.98564736629, 0.000387933, 38710000.0)

# Used in "bprecess" and "jprecess".
const A_precess  = 1e-6*[-1.62557, -0.31919, -0.13843] # In radians
const A_dot_precess = 0.001*[1.244 , -1.579, -0.660] # In arc seconds per century

"""
List of observing sites.  The observatories have `Observatory` type.
"""
const observatories =
    Dict("mgio"=>Observatory("Mount Graham International Observatory",
                             (32,42,04.69), (-109,53,31.25), 3191.0, -7),
         "whitin"=>Observatory("Whitin Observatory, Wellesley College",
                               42.295, -71.305833, 32, -5),
         "fmo"=>Observatory("Fan Mountain Observatory",
                            (37,52,41), (-78,41,34), 556, -5),
         "lmo"=>Observatory("Leander McCormick Observatory",
                            (38,02,00), (-78,31,24), 264, -5),
         "holi"=>Observatory("Observatorium Hoher List (Universitaet Bonn) - Germany",
                             50.16276, 6.85, 541, 1),
         "ca"=>Observatory("Calar Alto Observatory",
                           (37,13,25), (-2,32,46.5), 2168, 1),
         "bgsuo"=>Observatory("Bowling Green State Univ Observatory",
                              (41,22,42), (-83,39,33), 225., -5),
         "irtf"=>Observatory("NASA Infrared Telescope Facility",
                             19.826218, -155.471999, 4168, -10),
         "rozhen"=>Observatory("National Astronomical Observatory Rozhen - Bulgaria",
                               (41,41,35), (24,44,38), 1759, 2),
         "bosque"=>Observatory("Estacion Astrofisica Bosque Alegre, Cordoba",
                               (-31,35,54), (-64,32,45), 1250, -3),
         "casleo"=>Observatory("Complejo Astronomico El Leoncito, San Juan",
                               (-31,47,57), (-69,18,00), 2552, -3),
         "saao"=>Observatory("South African Astronomical Observatory",
                             (-32,22,46), (20,48,38.5), 1798., 2),
         "lna"=>Observatory("Laboratorio Nacional de Astrofisica - Brazil",
                            (-22,32,04), -45.5825, 1864., -3),
         "oro"=>Observatory("Oak Ridge Observatory",
                            (42,30,18.94), (-71,33,29.32), 184., -5),
         "flwo"=>Observatory("Whipple Observatory",
                             (31,40,51.4), (-110,52,39), 2320., -7),
         "vbo"=>Observatory("Vainu Bappu Observatory",
                            12.57666, 78.8266, 725., 5.5),
         "lowell"=>Observatory("Lowell Observatory",
                               (35,05.8), (-111,32.1), 2198., -7),
         "apo"=>Observatory("Apache Point Observatory",
                            (32,46.8), (-105,49.2), 2798., -7),
         "loiano"=>Observatory("Bologna Astronomical Observatory, Loiano - Italy",
                               (44,15,33), (11,20,2), 785., 1),
         "ekar"=>Observatory("Mt. Ekar 182 cm. Telescope",
                             (45,50,54.92), (11,34,52.08), 1413.69, 1),
         "keck"=>Observatory("W. M. Keck Observatory",
                             (19,49.7), (-155,28.7), 4160., -10),
         "bao"=>Observatory("Beijing XingLong Observatory",
                            (40,23.6), (117,34.5), 950., 8),
         "bmo"=>Observatory("Black Moshannon Observatory",
                            (40,55.3), (-78,00.3), 738., -5),
         "nov"=>Observatory("National Observatory of Venezuela",
                            (8,47.4), (-70,52.0), 3610, -4),
         "mdm"=>Observatory("Michigan-Dartmouth-MIT Observatory",
                            (31,57.0), (-111,37.0), 1938.5, -7),
         "palomar"=>Observatory("The Hale Telescope",
                                (33,21,21.6), (-116,51,46.80), 1706., -8),
         # https://en.wikipedia.org/wiki/Tonantzintla_Observatory
         "tona"=>Observatory("Observatorio Astronomico Nacional, Tonantzintla",
                             (19,01,58), (-98,18,50), 2166., -6),
         "spm"=>Observatory("Observatorio Astronomico Nacional, San Pedro Martir",
                            (31,01,45), (-115,29,13), 2830., -7),
         "dao"=>Observatory("Dominion Astrophysical Observatory",
                            (48,31.3), (-123,25.0), 229., -8),
         "mtbigelow"=>Observatory("Catalina Observatory: 61 inch telescope",
                                  (32,25.0), (-110,43.9), 2510., -7),
         "lco"=>Observatory("Las Campanas Observatory",
                            (-29,0.2), (-70,42.1), 2282, -4),
         "mcdonald"=>Observatory("McDonald Observatory",
                                 30.6716667, -104.0216667, 2075, -6),
         "aao"=>Observatory("Anglo-Australian Observatory",
                            (-31,16,37.34), (149,3,57.91), 1164, 10),
         "sso"=>Observatory("Siding Spring Observatory",
                            (-31,16,24.10), (149,3,40.3), 1149, 10),
         "mso"=>Observatory("Mt. Stromlo Observatory",
                            (-35,19,14.34), (149,1,27.6), 767, 10),
         "lapalma"=>Observatory("Roque de los Muchachos, La Palma",
                                (28,45.5), (-17,52.8), 2327, 0),
         "cfht"=>Observatory("Canada-France-Hawaii Telescope",
                             (19,49.6), (-155,28.3), 4215., -10),
         "mmto"=>Observatory("MMT Observatory",
                             (31,41.3), (-110,53.1), 2600., -7),
         "lick"=>Observatory("Lick Observatory",
                             (37,20.6), (-121,38.2), 1290., -8),
         "eso"=>Observatory("European Southern Observatory",
                            (-29,15.4), (-70,43.8), 2347., -4),
         "ctio"=>Observatory("Cerro Tololo Interamerican Observatory",
                             -30.16527778, -70.815, 2215., -4),
         "kpno"=>Observatory("Kitt Peak National Observatory",
                             (31,57.8), (-111,36.0), 2120., -7))

"""
List of planets of the Solar System, from Mercury to Pluto.  The elements of the
list have `Planet` type.

Reference for most quantities is the Planetary Fact Sheet:
http://nssdc.gsfc.nasa.gov/planetary/factsheet/index.html
"""
const planets =
    Dict("mercury"=>Planet("mercury", 2439.7e3, 2439.7e3, 2439.7e3, 3.3011e23,
                           0.20563069, 57909050e3, 87.9691*86400),
         "venus"=>Planet("venus", 6051.8e3, 6051.8e3, 6051.8e3, 4.8675e24,
                         0.00677323, 108208000e3, 224.701*86400),
         "earth"=>Planet("earth", 6371e3, 6378.137e3,  6356.752e3, 5.97237e24,
                         0.0167086, 149598023e3, 365.256363004*86400),
         "mars"=>Planet("mars", 3389.5e3, 3396.2e3, 3376.2e3, 3.3011e23,
                        0.09341233, 57909050e3, 87.9691*86400),
         "jupiter"=>Planet("jupiter", 69911e3, 71492e3, 66854e3, 1898.19e24,
                           0.04839266, 5.20336301*AU, 4332.589*86400),
         "saturn"=>Planet("saturn", 58232e3, 60268e3, 54364e3, 568.34e23,
                          0.05415060, 9.53707032*AU, 10759.22*86400),
         "uranus"=>Planet("uranus", 25362e3, 25559e3, 24973e3, 86.813e24,
                          0.04716771, 19.19126393*AU, 30685.4*86400),
         "neptune"=>Planet("neptune", 24622e3, 24764e3, 24341e3, 102.413e24,
                           0.00858587, 30.06896348*AU, 60189*86400),
         "pluto"=>Planet("pluto", 1187e3, 1187e3, 1187e3, 0.01303e24,
                         0.24880766, 39.48168677*AU, 90560*86400))


text:
using BenchmarkTools
using ChainRulesCore
using FiniteDifferences
using LinearAlgebra
using Printf
using Random
using SparseArrays
using StaticArrays
using Test

Random.seed!(1)

@testset "FiniteDifferences" begin
    include("deprecated.jl")
    include("methods.jl")
    include("numerics.jl")
    include("to_vec.jl")
    include("grad.jl")
end


text:
using Base: has_nondefault_cmd_flags
import Base.+
import Base.!
import Base.count

function !(nucleotide::BioSymbols.DNA)::BioSymbols.RNA

    # ok: use watson-crick base pairing rules to return the opposite nucleotide
    if nucleotide == DNA_T
        return RNA_A
    elseif nucleotide == DNA_A
        return RNA_U
    elseif nucleotide == DNA_C
        return RNA_G
    elseif nucleotide == DNA_G
        return RNA_C
    else
        return RNA_N
    end
end

function +(buffer::Array{String,1}, content::String;
    prefix::Union{String,Nothing}=nothing,suffix::Union{String,Nothing}=nothing)

    # create a new content line -
    new_line = content

    # prefix -
    if (prefix !== nothing)
        new_line = prefix*new_line
    end

    # suffix -
    if (suffix !== nothing)
        new_line = new_line*suffix
    end

    # cache -
    push!(buffer,new_line)
end

function read_file_from_path(path_to_file::String)::Array{String,1}

    # initialize -
    buffer = String[]

    # Read in the file -
    open("$(path_to_file)", "r") do file
        for line in eachline(file)
            +(buffer,line)
        end
    end

    # return -
    return buffer
end

function +(buffer::Array{String,1}, content_array::Array{String,1})
    for line in content_array
        push!(buffer, line)
    end
end

function check(result::Some)::(Union{Nothing,T} where T <: Any)

    # ok, so check, do we have an error object?
    # Yes: log the error if we have a logger, then throw the error.
    # No: return the result.value

     # Error case -
    if (isa(something(result), Exception) == true)

        # get the error object -
        error_object = result.value

        # get the error message as a String -
        error_message = sprint(showerror, error_object, catch_backtrace())
        @error(error_message)

        # throw -
        throw(result.value)
    end

    # default -
    return result.value
end

function count(sequence::BioSequences.LongSequence, bioSymbol::BioSymbol)::Int64

    # initialize -
    number_of_biosymbols = 0

    # iterate -
    for test_symbol in sequence
        if (test_symbol == bioSymbol)
            number_of_biosymbols = number_of_biosymbols + 1
        end
    end

    # return -
    return number_of_biosymbols
end


text:
module TestBaseFilter

using Random
using Test
using AMLPipelineBase
using DataFrames: nrow

function test_basefilter()
  data = getiris()
  ohe  = OneHotEncoder()
  mptr = Imputer()
  @test fit_transform!(ohe,data) |> Matrix |> sum |> round == 2229.0
  @test fit_transform(ohe,data) |> Matrix |> sum |> round == 2229.0
  @test fit_transform!(mptr,data) |> Matrix |> x->x[:,1:4] |> sum |> round == 2079.0
  @test fit_transform(mptr,data) |> Matrix |> x->x[:,1:4] |> sum |> round == 2079.0
  Random.seed!(1)
  data.mss=rand([missing,(1:100)...],nrow(data))
  @test fit_transform!(mptr,data) |> Matrix |> x->x[:,[(1:4)...,6]]  |> sum |> round > 100.0
  @test fit_transform(mptr,data) |> Matrix |> x->x[:,[(1:4)...,6]]  |> sum |> round > 100.0
  wrp = Wrapper(Dict(:transformer => OneHotEncoder()))
  @test fit_transform!(wrp,data) |> Matrix |> sum |> round > 100.0
  Random.seed!(1)
  data.mss=rand([missing,(1:100)...],nrow(data))
  wrp = Wrapper(Dict(:transformer => Imputer()))
  @test fit_transform!(wrp,data) |> Matrix |> x->x[:,[(1:4)...,6]]  |> sum |> round > 100.0
  @test fit_transform(wrp,data) |> Matrix |> x->x[:,[(1:4)...,6]]  |> sum |> round > 100.0
end
@testset "BaseFilter" begin
  Random.seed!(123)
  test_basefilter()
end

end


text:
using Revise # reduce recompilation time
using Plots
# using Documenter
using LinearAlgebra
using SparseArrays
using BenchmarkTools
using UnPack

push!(LOAD_PATH, "./src")
using CommonUtils
using Basis1D
using Basis2DQuad
using UniformQuadMesh

using SetupDG

push!(LOAD_PATH, "./examples/EntropyStableEuler.jl/src")
include("../EntropyStableEuler.jl/src/logmean.jl")
#using EntropyStableEuler
#using EntropyStableEuler.Fluxes1D

function wavespeed_1D(rho,rhou,E)
    p = pfun_nd(rho,rhou,E)
    cvel = @. sqrt(γ*p/rho)
    return @. abs(rhou/rho) + cvel
end

function pfun_nd(rho, rhou, rhov, E)
    rhoUnorm2 = (rhou^2+rhov^2)/rho
    return @. (γ-1)*(E - .5*rhoUnorm2)
end

function pfun_nd(rho, rhou, E)
    rhoUnorm2 = (rhou^2)/rho
    return @. (γ-1)*(E - .5*rhoUnorm2)
end

function primitive_to_conservative_hardcode(rho,u,v,p)
    Unorm = u^2+v^2
    E = @. p/(γ-1) + .5*rho*Unorm
    return (rho,rho*u,rho*v,E)
end

function euler_fluxes_2D(rhoL,uL,vL,betaL,rhologL,betalogL,
    rhoR,uR,vR,betaR,rhologR,betalogR)

    rholog = logmean.(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean.(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = (@. .5*(rhoL+rhoR))
    uavg   = (@. .5*(uL+uR))
    vavg   = (@. .5*(vL+vR))

    unorm = (@. uL*uR + vL*vR)
    pa    = (@. rhoavg/(betaL+betaR))
    f4aux = (@. rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm)

    FxS1 = (@. rholog*uavg)
    FxS2 = (@. FxS1*uavg + pa)
    FxS3 = (@. FxS1*vavg)
    FxS4 = (@. f4aux*uavg)

    FyS1 = (@. rholog*vavg)
    FyS2 = (@. FxS3)
    FyS3 = (@. FyS1*vavg + pa)
    FyS4 = (@. f4aux*vavg)
    return (FxS1,FxS2,FxS3,FxS4),(FyS1,FyS2,FyS3,FyS4)
end

function euler_fluxes(rhoL,uL,vL,betaL,rhoR,uR,vR,betaR)
    rhologL,betalogL,rhologR,betalogR = map(x->log.(x),(rhoL,betaL,rhoR,betaR))
    return euler_fluxes_2D(rhoL,uL,vL,betaL,rhologL,betalogL,
                           rhoR,uR,vR,betaR,rhologR,betalogR)
end

function is_zero(a)
    return abs(a) < TOL
end

function Riemann_2D(x,y)
    if x >= 0 && y >= 0
        return .5313, 0.0, 0.0, 0.4
    elseif x < 0 && y >= 0
        return 1.0, .7276, 0.0, 1.0
    elseif x < 0 && y < 0
        return .8, 0.0, 0.0, 1.0
    else
        return 1.0, 0.0, .7276, 1.0
    end

    # if x >= 0 && y >= 0
    #     return .5313, 0.0, 0.0, 0.004#0.4
    # elseif x < 0 && y >= 0
    #     return 1.0, .7276, 0.0, 1.0
    # elseif x < 0 && y < 0
    #     return .8, 0.0, 0.0, 1.0
    # else
    #     return 1.0, 0.0, .7276, 1.0
    # end

    # if x >= 0 && y >= 0
    #     return 1.5, 0.0, 0.0, 1.5
    # elseif x < 0 && y >= 0
    #     return 0.5323, 1.206, 0.0, 0.3
    # elseif x < 0 && y < 0
    #     return .138, 1.206, 1.206, 0.029
    # else
    #     return 0.5323, 0.0, 1.206, 0.3
    # end
end

function vortex(x,y,t,γ=1.4)
    x0 = 0
    y0 = 0
    beta = 5
    r2 = @. (x-x0-t)^2 + (y-y0)^2

    u = @. 1 - beta*exp(1-r2)*(y-y0)/(2*pi)
    v = @. beta*exp(1-r2)*(x-x0-t)/(2*pi)
    rho = @. 1 - (1/(8*γ*pi^2))*(γ-1)/2*(beta*exp(1-r2))^2
    rho = @. rho^(1/(γ-1))
    p = @. rho^γ

    return (rho, u, v, p)
end

const γ = 1.4
const TOL = 1e-14
const Nc = 4 # number of components

"Approximation parameters"
# N = 3 # The order of approximation
# K1D = 20
N = 2
K1D = 32
# N = 3
# K1D = 64
T = 0.25
#T = 0.25

"Mesh related variables"
Kx = K1D
Ky = K1D
VX, VY, EToV = uniform_quad_mesh(Kx, Ky)

"Initialize reference element"
r,_ = gauss_lobatto_quad(0,0,N)   # Reference nodes
VDM = vandermonde_1D(N,r)         # modal to nodal
Dr = grad_vandermonde_1D(N,r)/VDM # nodal differentiation
V1 = vandermonde_1D(1,r)/vandermonde_1D(1,[-1;1]) # nodal linear interpolation

Nq = N
rq,wq = gauss_lobatto_quad(0,0,Nq)
M = diagm(wq)
Mlump = zeros(size(M))
Mlump_inv = zeros(size(M))
for i = 1:Nq+1
    Mlump[i,i] = sum(M[i,:])
    Mlump_inv[i,i] = 1.0/Mlump[i,i]
end

# 1D operators
M1D = Mlump
M1D_inv = Mlump_inv
Q1D = M*Dr
B1D = zeros(N+1,N+1)
B1D[1,1] = -1
B1D[end,end] = 1
L = Array(spdiagm(0=>-2*ones(N+1), 1=>ones(N), -1=>ones(N)))
L[1,1] = -1
L[end,end] = -1
psi = pinv(L)*-1/2*B1D*ones(N+1)
S0 = zeros(N+1,N+1)
for i = 1:N+1
    for j = 1:N+1
        if L[i,j] != 0
            S0[i,j] = psi[j] - psi[i]
        end
    end
end
Q01D = S0+1/2*B1D

# Drop zeros
Q1D = Matrix(droptol!(sparse(Q1D),TOL))
Q01D = Matrix(droptol!(sparse(Q01D),TOL))
B1D = Matrix(droptol!(sparse(B1D),TOL))
S01D = Matrix(droptol!(sparse(S0),TOL))
S1D = Matrix(droptol!(sparse((Q1D-Q1D')/2),TOL))

# Tensor product operators
M = droptol!(sparse(kron(M1D,M1D)),TOL)
Minv = spdiagm(0 => 1 ./ diag(M))
Qr = droptol!(sparse(kron(M1D,Q1D)),TOL)
Qs = droptol!(sparse(kron(Q1D,M1D)),TOL)
Qr0 = droptol!(sparse(kron(M1D,Q01D)),TOL)
Qs0 = droptol!(sparse(kron(Q01D,M1D)),TOL)
Sr = droptol!(sparse((Qr-Qr')/2),TOL)
Ss = droptol!(sparse((Qs-Qs')/2),TOL)
S0r = droptol!(sparse(kron(M1D,S01D)),TOL)
S0s = droptol!(sparse(kron(S01D,M1D)),TOL)
Br = droptol!(sparse(Qr+Qr'),TOL)
Bs = droptol!(sparse(Qs+Qs'),TOL)

# Make domain periodic, TODO: hardcoded
Jf = 1/K1D # Assume uniform quad mesh
K = K1D^2
Nfaces = 4
Nfp = Nfaces*N

# physical normal, assume standard uniform quad mesh
#=
nxJ = zeros(Nfp,K)
nyJ = zeros(Nfp,K)
for k = 1:K
    nxJ[1,k] = -Jf
    nxJ[2:N,k] .= zeros(N-1)
    nxJ[N+1,k] = Jf
    for i = 0:N-2
        nxJ[N+2+i*2,k] = -Jf
        nxJ[N+2+i*2+1,k] = Jf
    end
    nxJ[end-N,k] = -Jf
    nxJ[end-(N-1):end-1,k] .= zeros(N-1)
    nxJ[end,k] = Jf

    nyJ[1:N+1,k] .= -Jf*ones(N+1)
    nyJ[end-N:end,k] .= Jf*ones(N+1)
end
=#
nxJ = zeros((N+1)*(N+1),K)
nyJ = zeros((N+1)*(N+1),K)
for k = 1:K
    for i = 0:N
        nxJ[i*(N+1)+1,k] = -Jf
        nxJ[i*(N+1)+N+1,k] = Jf
    end
    nyJ[1:N+1,k] .= -Jf*ones(N+1)
    nyJ[end-N:end,k] .= Jf*ones(N+1)
end

# indices of nodal points on faces
face_idx = zeros(Nfp)
face_idx[1:N+1] .= 1:N+1
face_idx[N+2:2*N+2] .= N*(N+1)+1:(N+1)*(N+1)
for i = 1:N-1
    face_idx[2*N+3+2*(i-1)] = i*(N+1)+1
    face_idx[2*N+3+2*(i-1)+1] = i*(N+1)+N+1
end
sort!(face_idx)
face_idx = Int64.(face_idx)

# face connectivity via global node number
mapM = zeros(Int64,Nfp,K)
mapP = zeros(Int64,Nfp,K)
mapP_x = zeros(Int64,Nfp,K)
mapP_y = zeros(Int64,Nfp,K)
for kj = 1:K1D
    for ki = 1:K1D
        k = (kj-1)*K1D+ki
        mapM[:,k] .= (k-1)*(N+1)*(N+1).+face_idx
        mapP[1:N+1,k] = mod(k-1-K1D,K)*(N+1)*(N+1).+face_idx[end-N:end]
        mapP[end-N:end,k] = mod(k-1+K1D,K)*(N+1)*(N+1).+face_idx[1:N+1]
        mapP[N+2:2:end-N-2,k] = ((kj-1)*K1D+mod(ki-2,K1D))*(N+1)*(N+1).+face_idx[N+3:2:end-N-1] # left interface
        mapP[N+3:2:end-N-1,k] = ((kj-1)*K1D+mod(ki,K1D))*(N+1)*(N+1).+face_idx[N+2:2:end-N-2] # right interface
        mapP_x[[1;N+2:2:end-N],k] = ((kj-1)*K1D+mod(ki-2,K1D))*(N+1)*(N+1).+face_idx[[N+1:2:end-N-1;end]]
        mapP_x[[N+1:2:end-N-1;end],k] = ((kj-1)*K1D+mod(ki,K1D))*(N+1)*(N+1).+face_idx[[1;N+2:2:end-N]]
        mapP_y[1:N+1,k] = mod(k-1-K1D,K)*(N+1)*(N+1).+face_idx[end-N:end]
        mapP_y[end-N:end,k] = mod(k-1+K1D,K)*(N+1)*(N+1).+face_idx[1:N+1]
    end
end

rd = init_reference_quad(N,gauss_lobatto_quad(0,0,N))
md = init_mesh((VX,VY),EToV,rd)
@unpack x,y,rxJ,ryJ,sxJ,syJ = md

function translate_bump(x,y)
    rho = exp(1/10/(x^2+y^2+0.1))
    u = 1.0
    v = 0.0
    p = 1.0
    return rho,u,v,p
end
#U = @. translate_bump(x,y)
U = @. Riemann_2D(x,y)
#U = @. vortex(x,y,0.0)
rho,u,v,p = [(x->x[i]).(U) for i = 1:Nc]
U = @. primitive_to_conservative_hardcode(rho,u,v,p)
rho,rhou,rhov,E = [(x->x[i]).(U) for i = 1:Nc]
U = (rho,rhou,rhov,E)

function limiting_param(U_low, P_ij)
    l = 1.0
    # Limit density
    if U_low[1] + P_ij[1] < -TOL
        l = min(abs(U_low[1])/(abs(P_ij[1])+1e-14), 1.0)
    end

    # limiting internal energy (via quadratic function)
    # a = P_ij[1]*P_ij[3]-1.0/2.0*P_ij[2]^2
    # b = U_low[3]*P_ij[1]+U_low[1]*P_ij[3]-U_low[2]*P_ij[2]
    # c = U_low[3]*U_low[1]-1.0/2.0*U_low[2]^2
    a = P_ij[1]*P_ij[4]-1.0/2.0*(P_ij[2]^2+P_ij[3]^2)
    b = U_low[4]*P_ij[1]+U_low[1]*P_ij[4]-U_low[2]*P_ij[2]-U_low[3]*P_ij[3]
    c = U_low[4]*U_low[1]-1.0/2.0*(U_low[2]^2+U_low[3]^2)

    l_eps_ij = 1.0
    if b^2-4*a*c >= 0
        r1 = (-b+sqrt(b^2-4*a*c))/(2*a)
        r2 = (-b-sqrt(b^2-4*a*c))/(2*a)
        if r1 > TOL && r2 > TOL
            l_eps_ij = min(r1,r2)
        elseif r1 > TOL && r2 < -TOL
            l_eps_ij = r1
        elseif r2 > TOL && r1 < -TOL
            l_eps_ij = r2
        end
    end

    l = min(l,l_eps_ij)
    return l
end

function rhs_IDP(U,K1D,N,Nfaces,nxJ,nyJ,Minv,Sr,Ss,S0r,S0s,Br,Bs,face_idx,mapM,mapP_x,mapP_y,step)
    p = pfun_nd.(U[1],U[2],U[3],U[4])
    flux_x = zero.(U)
    flux_y = zero.(U)
    @. flux_x[1] = U[2]
    @. flux_x[2] = U[2]^2/U[1]+p
    @. flux_x[3] = U[2]*U[3]/U[1]
    @. flux_x[4] = U[4]*U[2]/U[1]+p*U[2]/U[1]
    @. flux_y[1] = U[3]
    @. flux_y[2] = U[2]*U[3]/U[1]
    @. flux_y[3] = U[3]^2/U[1]+p
    @. flux_y[4] = U[4]*U[3]/U[1]+p*U[3]/U[1]

    Ub = zero.(U)
    @. Ub[1] = U[1]
    @. Ub[2] = U[2]/U[1]
    @. Ub[3] = U[3]/U[1]
    @. Ub[4] = U[1]/(2*p)

    K = K1D^2
    Nfp = Nfaces*N
    Np = (N+1)*(N+1)
    J = (1/K1D)^2
    Jf = 1/K1D
    rxJ = 1/K1D
    syJ = 1/K1D
    sJ = 1/K1D

    F_low  = [zeros(Np,Np) for i = 1:Nc]
    F_high = [zeros(Np,Np) for i = 1:Nc]
    F_P    = [zeros(Nfp) for i = 1:Nc]
    rhsU   = [zeros(Np,K) for i = 1:Nc]
    U_low  = [zeros(Np) for i = 1:Nc]
    L = zeros(Np,Np)

    d_ii_arr = zeros(Np,K)
    # TODO: precompute
    for k = 1:K
        for i = 1:Np
            for j = 1:Np
                c_ij_norm = sqrt(S0r[i,j]^2+S0s[i,j]^2)
                if abs(c_ij_norm) >= TOL
                    n_ij = [S0r[i,j]; S0s[i,j]]./c_ij_norm
                    wavespd_i = wavespeed_1D(U[1][i,k],(n_ij[1]*U[2][i,k]+n_ij[2]*U[3][i,k]),U[4][i,k])
                    wavespd_j = wavespeed_1D(U[1][j,k],(n_ij[1]*U[2][j,k]+n_ij[2]*U[3][j,k]),U[4][j,k])
                    wavespd = max(wavespd_i,wavespd_j)
                    d_ij = wavespd*c_ij_norm
                    d_ii_arr[i,k] -= d_ij
                end
            end
        end

        S0r1 = sum(S0r,dims=2)
        S0s1 = sum(S0s,dims=2)
        for i = 1:Nfp
            fidM = mapM[i,k]
            S0rP = -S0r1[face_idx[i]]
            S0sP = -S0s1[face_idx[i]]

            # flux in x direction
            fidP = mapP_x[i,k]
            if fidP != 0
                wavespd_M = wavespeed_1D(U[1][fidM],U[2][fidM],U[4][fidM])
                wavespd_P = wavespeed_1D(U[1][fidP],U[2][fidP],U[4][fidP])
                wavespd = max(wavespd_M,wavespd_P)
                d_ij = wavespd*abs(S0rP)
                d_ii_arr[face_idx[i],k] -= d_ij
            end

            # flux in y direction
            fidP = mapP_y[i,k]
            if fidP != 0
                wavespd_M = wavespeed_1D(U[1][fidM],U[3][fidM],U[4][fidM])
                wavespd_P = wavespeed_1D(U[1][fidP],U[3][fidP],U[4][fidP])
                wavespd = max(wavespd_M,wavespd_P)
                d_ij = wavespd*abs(S0sP)
                d_ii_arr[face_idx[i],k] -= d_ij
            end
        end
    end

    dt = minimum(-J/2*M*(1 ./d_ii_arr))
    dt = dt

    for k = 1:K
        for c = 1:Nc
            for i = 1:Nfp
                F_P[c][i] = 0.0
            end
            for i = 1:Np
                U_low[c][i] = 0.0
            end
        end

        # Calculate interior low and high order flux
        for i = 1:Np
            for j = 1:Np
                c_ij_norm = sqrt(rxJ^2*S0r[i,j]^2+syJ^2*S0s[i,j]^2)
                if abs(c_ij_norm) >= TOL
                    n_ij = [rxJ*S0r[i,j]; syJ*S0s[i,j]]./c_ij_norm
                    wavespd_i = wavespeed_1D(U[1][i,k],(n_ij[1]*U[2][i,k]+n_ij[2]*U[3][i,k]),U[4][i,k])
                    wavespd_j = wavespeed_1D(U[1][j,k],(n_ij[1]*U[2][j,k]+n_ij[2]*U[3][j,k]),U[4][j,k])
                    wavespd = max(wavespd_i,wavespd_j)
                    d_ij = wavespd*c_ij_norm

                    for c = 1:Nc
                        F_low[c][i,j]  = (rxJ*S0r[i,j]*(flux_x[c][i,k]+flux_x[c][j,k])
                                        + syJ*S0s[i,j]*(flux_y[c][i,k]+flux_y[c][j,k])
                                        - d_ij*(U[c][j,k]-U[c][i,k]))
                    end
                end
            end
        end

        # High order flux
        for i = 1:Np-1
            for j = i+1:Np
                # TODO: can preallocate nonzero entries
                if Sr[i,j] != 0 || Ss[i,j] != 0
                    F1,F2 = euler_fluxes(Ub[1][i,k],Ub[2][i,k],Ub[3][i,k],Ub[4][i,k],Ub[1][j,k],Ub[2][j,k],Ub[3][j,k],Ub[4][j,k])
                    for c = 1:Nc
                        val = 2*rxJ*Sr[i,j]*F1[c]+2*syJ*Ss[i,j]*F2[c]
                        F_high[c][i,j] = val
                        F_high[c][j,i] = -val
                    end
                end
            end
        end

        # Calculate interface flux
        S0r1 = sum(S0r,dims=2)
        S0s1 = sum(S0s,dims=2)
        for i = 1:Nfp
            fidM = mapM[i,k]
            S0rP = -S0r1[face_idx[i]]
            S0sP = -S0s1[face_idx[i]]

            # flux in x direction
            fidP = mapP_x[i,k]
            if fidP != 0
                wavespd_M = wavespeed_1D(U[1][fidM],U[2][fidM],U[4][fidM])
                wavespd_P = wavespeed_1D(U[1][fidP],U[2][fidP],U[4][fidP])
                wavespd = max(wavespd_M,wavespd_P)
                d_ij = wavespd*abs(S0rP)
                for c = 1:Nc
                    F_P[c][i] += (Jf*S0rP*(flux_x[c][fidM]+flux_x[c][fidP])
                                - Jf*d_ij*(U[c][fidP]-U[c][fidM]))
                end
            end

            # flux in y direction
            fidP = mapP_y[i,k]
            if fidP != 0
                wavespd_M = wavespeed_1D(U[1][fidM],U[3][fidM],U[4][fidM])
                wavespd_P = wavespeed_1D(U[1][fidP],U[3][fidP],U[4][fidP])
                wavespd = max(wavespd_M,wavespd_P)
                d_ij = wavespd*abs(S0sP)
                for c = 1:Nc
                    F_P[c][i] += (Jf*S0sP*(flux_y[c][fidM]+flux_y[c][fidP])
                                - Jf*d_ij*(U[c][fidP]-U[c][fidM]))
                end
            end
        end

        # Calculate low order solution
        for i = 1:Np
            for c = 1:Nc
                U_low[c][i] = sum(F_low[c][i,:])
            end
        end
        # TODO: redundant
        for i = 1:Nfp
            for c = 1:Nc
                U_low[c][face_idx[i]] += F_P[c][i]
            end
        end
        for c = 1:Nc
            for i = 1:Np
                U_low[c][i] = U[c][i,k]-dt/J/M[i,i]*U_low[c][i]
            end
        end

        # Calculate limiting parameters
        P_ij = [0.0;0.0;0.0;0.0]
        for i = 1:Np
            m_i = J*M[i,i]
            lambda_j = 1/(Np-1)
            U_low_i = [U_low[c][i] for c = 1:Nc]
            for j = 1:Np
                if i != j
                    for c = 1:Nc
                        P_ij[c] = dt/(m_i*lambda_j)*(F_low[c][i,j]-F_high[c][i,j])
                    end
                    L[i,j] = limiting_param(U_low_i, P_ij)
                end
            end
        end

        # Symmetrize limiting parameters
        for i = 1:Np
            for j = 1:Np
                if i != j
                    l = min(L[i,j],L[j,i])
                    L[i,j] = l
                    L[j,i] = l
                end
            end
        end

        # L = ones(size(L))
        # if sum(L)-Np*Np+Np < 0
        #     @show sum(L)-Np*Np+Np
        # end

        for c = 1:Nc
            rhsU[c][:,k] = sum((L.-1).*F_low[c] - L.*F_high[c],dims=2)
            #rhsU[c][:,k] = sum((-1).*F_high[c],dims=2)
            #rhsU[c][:,k] = sum((-1).*F_low[c],dims=2)
            for i = 1:Nfp
                rhsU[c][face_idx[i],k] -= F_P[c][i]
            end
        end
    end

    for c = 1:Nc
        rhsU[c] = 1/J*Minv*rhsU[c]
    end

    return rhsU,dt
end

function rhs_IDP(U,K1D,N,Nfaces,nxJ,nyJ,Minv,Sr,Ss,S0r,S0s,Br,Bs,face_idx,mapM,mapP_x,mapP_y,step,dt)
    p = pfun_nd.(U[1],U[2],U[3],U[4])
    flux_x = zero.(U)
    flux_y = zero.(U)
    @. flux_x[1] = U[2]
    @. flux_x[2] = U[2]^2/U[1]+p
    @. flux_x[3] = U[2]*U[3]/U[1]
    @. flux_x[4] = U[4]*U[2]/U[1]+p*U[2]/U[1]
    @. flux_y[1] = U[3]
    @. flux_y[2] = U[2]*U[3]/U[1]
    @. flux_y[3] = U[3]^2/U[1]+p
    @. flux_y[4] = U[4]*U[3]/U[1]+p*U[3]/U[1]

    Ub = zero.(U)
    @. Ub[1] = U[1]
    @. Ub[2] = U[2]/U[1]
    @. Ub[3] = U[3]/U[1]
    @. Ub[4] = U[1]/(2*p)

    K = K1D^2
    Nfp = Nfaces*N
    Np = (N+1)*(N+1)
    J = (1/K1D)^2
    Jf = 1/K1D
    rxJ = 1/K1D
    syJ = 1/K1D
    sJ = 1/K1D

    F_low  = [zeros(Np,Np) for i = 1:Nc]
    F_high = [zeros(Np,Np) for i = 1:Nc]
    F_P    = [zeros(Nfp) for i = 1:Nc]
    rhsU   = [zeros(Np,K) for i = 1:Nc]
    U_low  = [zeros(Np) for i = 1:Nc]
    L = zeros(Np,Np)

    for k = 1:K
        for c = 1:Nc
            for i = 1:Nfp
                F_P[c][i] = 0.0
            end
            for i = 1:Np
                U_low[c][i] = 0.0
            end
        end

        # Calculate interior low and high order flux
        for i = 1:Np
            for j = 1:Np
                c_ij_norm = sqrt(rxJ^2*S0r[i,j]^2+syJ^2*S0s[i,j]^2)
                if abs(c_ij_norm) >= TOL
                    n_ij = [rxJ*S0r[i,j]; syJ*S0s[i,j]]./c_ij_norm
                    wavespd_i = wavespeed_1D(U[1][i,k],(n_ij[1]*U[2][i,k]+n_ij[2]*U[3][i,k]),U[4][i,k])
                    wavespd_j = wavespeed_1D(U[1][j,k],(n_ij[1]*U[2][j,k]+n_ij[2]*U[3][j,k]),U[4][j,k])
                    wavespd = max(wavespd_i,wavespd_j)
                    d_ij = wavespd*c_ij_norm

                    for c = 1:Nc
                        F_low[c][i,j]  = (rxJ*S0r[i,j]*(flux_x[c][i,k]+flux_x[c][j,k])
                                        + syJ*S0s[i,j]*(flux_y[c][i,k]+flux_y[c][j,k])
                                        - d_ij*(U[c][j,k]-U[c][i,k]))
                    end
                end
            end
        end

        # High order flux
        for i = 1:Np-1
            for j = i+1:Np
                # TODO: can preallocate nonzero entries
                if Sr[i,j] != 0 || Ss[i,j] != 0
                    F1,F2 = euler_fluxes(Ub[1][i,k],Ub[2][i,k],Ub[3][i,k],Ub[4][i,k],Ub[1][j,k],Ub[2][j,k],Ub[3][j,k],Ub[4][j,k])
                    for c = 1:Nc
                        val = 2*rxJ*Sr[i,j]*F1[c]+2*syJ*Ss[i,j]*F2[c]
                        F_high[c][i,j] = val
                        F_high[c][j,i] = -val
                    end
                end
            end
        end

        # Calculate interface flux
        S0r1 = sum(S0r,dims=2)
        S0s1 = sum(S0s,dims=2)
        for i = 1:Nfp
            fidM = mapM[i,k]
            S0rP = -S0r1[face_idx[i]]
            S0sP = -S0s1[face_idx[i]]

            # flux in x direction
            fidP = mapP_x[i,k]
            if fidP != 0
                wavespd_M = wavespeed_1D(U[1][fidM],U[2][fidM],U[4][fidM])
                wavespd_P = wavespeed_1D(U[1][fidP],U[2][fidP],U[4][fidP])
                wavespd = max(wavespd_M,wavespd_P)
                d_ij = wavespd*abs(S0rP)
                for c = 1:Nc
                    F_P[c][i] += (Jf*S0rP*(flux_x[c][fidM]+flux_x[c][fidP])
                                - Jf*d_ij*(U[c][fidP]-U[c][fidM]))
                end
            end

            # flux in y direction
            fidP = mapP_y[i,k]
            if fidP != 0
                wavespd_M = wavespeed_1D(U[1][fidM],U[3][fidM],U[4][fidM])
                wavespd_P = wavespeed_1D(U[1][fidP],U[3][fidP],U[4][fidP])
                wavespd = max(wavespd_M,wavespd_P)
                d_ij = wavespd*abs(S0sP)
                for c = 1:Nc
                    F_P[c][i] += (Jf*S0sP*(flux_y[c][fidM]+flux_y[c][fidP])
                                - Jf*d_ij*(U[c][fidP]-U[c][fidM]))
                end
            end
        end

        # Calculate low order solution
        for i = 1:Np
            for c = 1:Nc
                U_low[c][i] = sum(F_low[c][i,:])
            end
        end
        # TODO: redundant
        for i = 1:Nfp
            for c = 1:Nc
                U_low[c][face_idx[i]] += F_P[c][i]
            end
        end
        for c = 1:Nc
            for i = 1:Np
                U_low[c][i] = U[c][i,k]-dt/J/M[i,i]*U_low[c][i]
            end
        end

        # Calculate limiting parameters
        P_ij = [0.0;0.0;0.0;0.0]
        for i = 1:Np
            m_i = J*M[i,i]
            lambda_j = 1/(Np-1)
            U_low_i = [U_low[c][i] for c = 1:Nc]
            for j = 1:Np
                if i != j
                    for c = 1:Nc
                        P_ij[c] = dt/(m_i*lambda_j)*(F_low[c][i,j]-F_high[c][i,j])
                    end
                    L[i,j] = limiting_param(U_low_i, P_ij)
                end
            end
        end

        # Symmetrize limiting parameters
        for i = 1:Np
            for j = 1:Np
                if i != j
                    l = min(L[i,j],L[j,i])
                    L[i,j] = l
                    L[j,i] = l
                end
            end
        end

        # L = zeros(size(L))
        # if sum(L)-Np*Np+Np < 0
        #     @show sum(L)-Np*Np+Np
        # end

        for c = 1:Nc
            rhsU[c][:,k] = sum((L.-1).*F_low[c] - L.*F_high[c],dims=2)
            #rhsU[c][:,k] = sum((-1).*F_high[c],dims=2)
            #rhsU[c][:,k] = sum((-1).*F_low[c],dims=2)
            for i = 1:Nfp
                rhsU[c][face_idx[i],k] -= F_P[c][i]
            end
        end
    end

    for c = 1:Nc
        rhsU[c] = 1/J*Minv*rhsU[c]
    end

    return rhsU,dt
end

function rhs_IDPlow(U,K1D,N,Nfaces,nxJ,nyJ,Minv,Sr,Ss,S0r,S0s,Br,Bs,face_idx,mapM,mapP_x,mapP_y)
    p = pfun_nd.(U[1],U[2],U[3],U[4])
    flux_x = zero.(U)
    flux_y = zero.(U)
    @. flux_x[1] = U[2]
    @. flux_x[2] = U[2]^2/U[1]+p
    @. flux_x[3] = U[2]*U[3]/U[1]
    @. flux_x[4] = U[4]*U[2]/U[1]+p*U[2]/U[1]
    @. flux_y[1] = U[3]
    @. flux_y[2] = U[2]*U[3]/U[1]
    @. flux_y[3] = U[3]^2/U[1]+p
    @. flux_y[4] = U[4]*U[3]/U[1]+p*U[3]/U[1]

    K = K1D^2
    Nfp = Nfaces*N
    Np = (N+1)*(N+1)
    J = (1/K1D)^2
    Jf = 1/K1D
    rxJ = 1/K1D
    syJ = 1/K1D
    sJ = 1/K1D

    # Initialize sparse matrices to store numerical fluxes across interfaces
    I_arr = Int64[]
    J_arr = Int64[]
    for k = 1:K
        for n = 1:Nfp
            i = mapM[n,k]
            jx = mapP_x[n,k]
            jy = mapP_y[n,k]
            if jx != 0
                push!(I_arr,i)
                push!(J_arr,jx)
            end
            if jy != 0
                push!(I_arr,i)
                push!(J_arr,jy)
            end
        end
    end
    F_P_low = [sparse(I_arr,J_arr,zeros(size(I_arr))) for i = 1:Nc]
    # Low order and high order algebraic fluxes
    F_low   = [zeros(Np,Np) for i = 1:Nc]
    rhsU    = [zeros(Np,K) for i = 1:Nc]

    for k = 1:K
        # TODO: redundant iterations
        for i = 1:Np
            for j = 1:Np
                c_ij_norm = sqrt(rxJ^2*S0r[i,j]^2+syJ^2*S0s[i,j]^2)
                if abs(c_ij_norm) >= TOL
                    n_ij = [rxJ*S0r[i,j]; syJ*S0s[i,j]]./c_ij_norm
                    wavespd_i = wavespeed_1D(U[1][i,k],(n_ij[1]*U[2][i,k]+n_ij[2]*U[3][i,k]),U[4][i,k])
                    wavespd_j = wavespeed_1D(U[1][j,k],(n_ij[1]*U[2][j,k]+n_ij[2]*U[3][j,k]),U[4][j,k])
                    wavespd = max(wavespd_i,wavespd_j)
                    d_ij = wavespd*c_ij_norm

                    for c = 1:Nc
                        F_low[c][i,j] = (rxJ*S0r[i,j]*(flux_x[c][i,k]+flux_x[c][j,k])
                                       + syJ*S0s[i,j]*(flux_y[c][i,k]+flux_y[c][j,k])
                                       - d_ij*(U[c][j,k]-U[c][i,k]))
                    end
                end
            end
        end

        S0r1 = sum(S0r,dims=2)
        S0s1 = sum(S0s,dims=2)
        # Interface low order flux
        for i in 1:Nfp
            fidM = mapM[i,k]

            S0rP = -S0r1[face_idx[i]]
            S0sP = -S0s1[face_idx[i]]

            # numerical flux in x direction
            fidP = mapP_x[i,k]
            if fidP != 0
                c_ij_norm = abs(S0rP)
                if abs(c_ij_norm) >= TOL
                    n_ij_x = S0rP/c_ij_norm
                    wavespd_M = wavespeed_1D(U[1][fidM],n_ij_x*U[2][fidM],U[4][fidM])
                    wavespd_P = wavespeed_1D(U[1][fidP],n_ij_x*U[2][fidP],U[4][fidP])
                    wavespd = max(wavespd_M,wavespd_P)
                    d_ij = wavespd*c_ij_norm
                    for c = 1:Nc
                        F_P_low[c][fidM,fidP] += (Jf*S0rP*(flux_x[c][fidM]+flux_x[c][fidP])
                                                - Jf*d_ij*(U[c][fidP]-U[c][fidM]))
                    end
                end
            end

            # numerical flux in y direction
            fidP = mapP_y[i,k]
            if fidP != 0
                c_ij_norm = abs(S0sP)
                if abs(c_ij_norm) >= TOL
                    n_ij_y = S0sP/c_ij_norm
                    wavespd_M = wavespeed_1D(U[1][fidM],n_ij_y*U[3][fidM],U[4][fidM])
                    wavespd_P = wavespeed_1D(U[1][fidP],n_ij_y*U[3][fidP],U[4][fidP])
                    wavespd = max(wavespd_M,wavespd_P)
                    d_ij = wavespd*c_ij_norm
                    for c = 1:Nc
                        F_P_low[c][fidM,fidP] += (Jf*S0sP*(flux_y[c][fidM]+flux_y[c][fidP])
                                                - Jf*d_ij*(U[c][fidP]-U[c][fidM]))
                    end
                end
            end
        end

        for i = 1:Np
            for c = 1:Nc
                rhsU[c][i,k] += sum(F_low[c][i,:])
                rhsU[c][i,k] += sum(F_P_low[c][i+(k-1)*Np,:])
            end
        end

    end

    for c = 1:Nc
        rhsU[c] = -1/J*Minv*rhsU[c]
    end

    return rhsU
end

function rhs_ESDG(U,K1D,N,Nfaces,nxJ,nyJ,Minv,Sr,Ss,S0r,S0s,Br,Bs,face_idx,mapM,mapP_x,mapP_y)
    p = pfun_nd.(U[1],U[2],U[3],U[4])

    Ub = zero.(U)
    @. Ub[1] = U[1]
    @. Ub[2] = U[2]/U[1]
    @. Ub[3] = U[3]/U[1]
    @. Ub[4] = U[1]/(2*p)

    K = K1D^2
    Nfp = Nfaces*N
    Np = (N+1)*(N+1)
    J = (1/K1D)^2
    Jf = 1/K1D

    # Low order and high order algebraic fluxes
    # Initialize sparse matrices to store numerical fluxes across interfaces
    I_arr = Int64[]
    J_arr = Int64[]
    for k = 1:K
        for n = 1:Nfp
            i = mapM[n,k]
            jx = mapP_x[n,k]
            jy = mapP_y[n,k]
            if jx != 0
                push!(I_arr,i)
                push!(J_arr,jx)
            end
            if jy != 0
                push!(I_arr,i)
                push!(J_arr,jy)
            end
        end
    end
    F_P_high = [sparse(I_arr,J_arr,zeros(size(I_arr))) for i = 1:Nc]
    F_high   = [zeros(Np,Np) for i = 1:Nc]
    rhsU     = [zeros(Np,K) for i = 1:Nc]

    # TODO: assume mesh is uniform, so geometric factors constant
    #       and mesh is not rotated, i.e. ryJ = sxJ = 0.0
    rxJ = 1/K1D
    syJ = 1/K1D
    sJ = 1/K1D

    for k = 1:K
        for i = 1:Np-1
            for j = i+1:Np
                # TODO: can preallocate nonzero entries
                if Sr[i,j] != 0 || Ss[i,j] != 0
                    F1,F2 = euler_fluxes(Ub[1][i,k],Ub[2][i,k],Ub[3][i,k],Ub[4][i,k],Ub[1][j,k],Ub[2][j,k],Ub[3][j,k],Ub[4][j,k])
                    for c = 1:Nc
                        val = 2*rxJ*Sr[i,j]*F1[c]+2*syJ*Ss[i,j]*F2[c]
                        F_high[c][i,j] = val
                        F_high[c][j,i] = -val
                    end
                end
            end
        end

        # Interface high order flux
        for i in 1:Nfp
            fidM = mapM[i,k]
            nxJM = nxJ[fidM]
            nyJM = nyJ[fidM]
            rhou_nM_x = (U[2][fidM]*nxJM)/sJ
            rhou_nM_y = (U[3][fidM]*nyJM)/sJ
            wavespdM_x = wavespeed_1D(U[1][fidM],rhou_nM_x,U[4][fidM])
            wavespdM_y = wavespeed_1D(U[1][fidM],rhou_nM_y,U[4][fidM])

            # numerical flux in x direction
            fidP = mapP_x[i,k]
            if fidP != 0
                nxJP = nxJ[fidP]
                nyJP = nyJ[fidP]
                #rhou_nP = (U[2][fidP]*nxJP+U[3][fidP]*nyJP)/sJ
                rhou_nP = (U[2][fidP]*nxJP)/sJ
                wavespdP = wavespeed_1D(U[1][fidP],rhou_nP,U[4][fidP])
                wavespd = max(wavespdM_x,wavespdP)
                F1,_ = euler_fluxes(Ub[1][fidM],Ub[2][fidM],Ub[3][fidM],Ub[4][fidM],Ub[1][fidP],Ub[2][fidP],Ub[3][fidP],Ub[4][fidP])
                for c = 1:Nc
                    F_P_high[c][fidM,fidP] += Jf*Br[face_idx[i],face_idx[i]]*F1[c]-sJ*wavespd/2*(U[c][fidP]-U[c][fidM])#nxJM*Br[face_idx[i],face_idx[i]]*F1[c]#-wavespd/2*(U[c][fidP]-U[c][fidM])
                end
            end

            # numerical flux in y direction
            fidP = mapP_y[i,k]
            if fidP != 0
                nxJP = nxJ[fidP]
                nyJP = nyJ[fidP]
                #rhou_nP = (U[2][fidP]*nxJP+U[3][fidP]*nyJP)/sJ
                rhou_nP = (U[3][fidP]*nyJP)/sJ
                wavespdP = wavespeed_1D(U[1][fidP],rhou_nP,U[4][fidP])
                wavespd = max(wavespdM_y,wavespdP)
                _,F2 = euler_fluxes(Ub[1][fidM],Ub[2][fidM],Ub[3][fidM],Ub[4][fidM],Ub[1][fidP],Ub[2][fidP],Ub[3][fidP],Ub[4][fidP])
                for c = 1:Nc
                    F_P_high[c][fidM,fidP] += Jf*Bs[face_idx[i],face_idx[i]]*F2[c]-sJ*wavespd/2*(U[c][fidP]-U[c][fidM])#nyJM*Bs[face_idx[i],face_idx[i]]*F2[c]#-wavespd/2*(U[c][fidP]-U[c][fidM])
                end
            end
        end

        for i = 1:Np
            for c = 1:Nc
                rhsU[c][i,k] += sum(F_high[c][i,:])
                rhsU[c][i,k] += sum(F_P_high[c][i+(k-1)*Np,:])
            end
        end
    end

    for c = 1:Nc
        rhsU[c] = -1/J*Minv*rhsU[c]
    end

    return rhsU

end

# Time stepping
"Time integration"
t = 0.0
U = collect(U)
resU = [zeros(size(x)),zeros(size(x)),zeros(size(x)),zeros(size(x))]
resW = [zeros(size(x)),zeros(size(x)),zeros(size(x)),zeros(size(x))]
resZ = [zeros(size(x)),zeros(size(x)),zeros(size(x)),zeros(size(x))]

@unpack VDM = rd
rp,sp = equi_nodes_2D(20)
Vp = vandermonde_2D(N,rp,sp)/VDM
gr(aspect_ratio=:equal,legend=false,
   markerstrokewidth=0,markersize=2,axis=nothing)
const GIFINTERVAL = 100
#plot()

dt_hist = []
anim = Animation()
i = 1

while t < T
    #rhsU,dt = rhs_IDP(U,K1D,N,Nfaces,nxJ,nyJ,Minv,Sr,Ss,S0r,S0s,Br,Bs,face_idx,mapM,mapP_x,mapP_y,i)
    # @. U = U + dt*rhsU

    # SSPRK(3,3)
    dt = min(1e-4,T-t)
    rhsU,dt = rhs_IDP(U,K1D,N,Nfaces,nxJ,nyJ,Minv,Sr,Ss,S0r,S0s,Br,Bs,face_idx,mapM,mapP_x,mapP_y,i,dt)
    @. resW = U + dt*rhsU
    rhsU,_ = rhs_IDP(resW,K1D,N,Nfaces,nxJ,nyJ,Minv,Sr,Ss,S0r,S0s,Br,Bs,face_idx,mapM,mapP_x,mapP_y,i,dt)
    @. resZ = resW+dt*rhsU
    @. resW = 3/4*U+1/4*resZ
    rhsU,_ = rhs_IDP(resW,K1D,N,Nfaces,nxJ,nyJ,Minv,Sr,Ss,S0r,S0s,Br,Bs,face_idx,mapM,mapP_x,mapP_y,i,dt)
    @. resZ = resW+dt*rhsU
    @. U = 1/3*U+2/3*resZ

    push!(dt_hist,dt)
    global t = t + dt
    println("Current time $t with time step size $dt, and final time $T, at step $i")
    global i = i + 1
    # if mod(i,50) == 1
    #     scatter(Vp*x,Vp*y,Vp*U[1],zcolor=Vp*U[1],camera=(0,90))
    #     frame(anim)
    # end
end

xp = Vp*x
yp = Vp*y
rhop = Vp*U[1]
p = pfun_nd.(U[1],U[2],U[3],U[4])
pp = Vp*p
#plt = scatter(Vp*x,Vp*y,Vp*U[1],zcolor=Vp*U[1],camera=(0,90))
plt = scatter(xp,yp,rhop,zcolor=rhop,camera=(0,90))
display(plt)
#savefig(plt,"~/Desktop/tmp.png")
#gif(anim,"~/Desktop/tmp.gif",fps=15)

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
__precompile__()
module AutoAligns

import Base: position, print, println

export left, right, center, AutoAlign

abstract type Alignment end

struct Left <: Alignment end
const left = Left()
align_string(::Left, s::AbstractString, width::Int) = rpad(s, width)

struct Right <: Alignment end
const right = Right()
align_string(::Right, s::AbstractString, width::Int) = lpad(s, width)

struct Center <: Alignment end
const center = Center()
function align_string(::Center, s::AbstractString, width::Int)
    len = length(s)
    @assert len <= width
    lpad(rpad(s, len+floor(Int, (width-len)/2)), width)
end

get_alignment(a::Alignment, pos::Int) = a

function get_alignment(a::AbstractVector, pos::Int)
    if length(a) < pos
        a[end]
    else
        a[pos]
    end
end

function get_alignment(a::Dict, pos::Int)
    get(a, pos) do
        a[:default]
    end
end

type AutoAlign
    align
    table
    widths
    function AutoAlign(; align=left)
        table = Vector{Vector}()
        push!(table, Vector())
        new(align, table, Vector{Int}())
    end
end

"""
Get the position of the next item that would be printed.
"""
position(aa::AutoAlign) = length(aa.table[end])+1

_newline(aa::AutoAlign) =  (push!(aa.table, Vector()); nothing)

function _print(aa::AutoAlign, s::AbstractString, align=aa.align)
    pos = position(aa)
    if (length(aa.widths) < pos)
        push!(aa.widths, 0)
    end
    s = normalize_string(s, stripcc=true)
    aa.widths[pos] = max(aa.widths[pos], length(s))
    push!(aa.table[end], (get_alignment(align, pos),s))
    nothing
end

_print(aa::AutoAlign, x, align=aa.align) = _print(aa, string(x), align)

function print(aa::AutoAlign, xs...; align=aa.align)
    for x in xs
        _print(aa, x, align)
    end
end

function println(aa::AutoAlign, xs...; align=aa.align)
    print(aa, xs...; align=align)
    _newline(aa)
end

function print(io::IO, aa::AutoAlign)
    for (i,line) in enumerate(aa.table)
        if i != 1
            println(io)
        end
        for (as,w) in zip(line,aa.widths)
            print(io, align_string(as[1], as[2], w))
        end
    end
    nothing
end

end # module


text:
# Shallow ice approximation (SIA) implicit solver for Greenland (steady state)
# Runs on both CPU (threaded) and GPU (Nvidia Cuda)

const USE_GPU = false # switch here to use GPU
# enable plotting & saving by default
if !@isdefined do_visu; do_visu = true end
if !@isdefined do_save; do_save = true end

using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 2)
else
    @init_parallel_stencil(Threads, Float64, 2)
end
using JLD, Plots, Printf, LinearAlgebra

# CPU functions: finite difference stencil operation support functions
@views av(A)  = 0.25*(A[1:end-1,1:end-1].+A[2:end,1:end-1].+A[1:end-1,2:end].+A[2:end,2:end])
@views inn(A) = A[2:end-1,2:end-1]

# GPU function
@parallel_indices (iy) function bc_x!(A::Data.Array)
    A[1  , iy] = A[2    , iy]
    A[end, iy] = A[end-1, iy]
    return
end

@parallel_indices (ix) function bc_y!(A::Data.Array)
    A[ix, 1  ] = A[ix, 2    ]
    A[ix, end] = A[ix, end-1]
    return
end

@parallel function compute_Err1!(Err, H)
    @all(Err) = @all(H)
    return
end

@parallel function compute_Err2!(Err, H)
    @all(Err) = @all(Err) - @all(H)
    return
end

@parallel function compute_M_dS!(M, dSdx, dSdy, S, z_ELA, grad_b, b_max, dx, dy)
    @all(M)    = min(@all(grad_b)*(@all(S) - @all(z_ELA)), b_max)
    @all(dSdx) = @d_xa(S)/dx
    @all(dSdy) = @d_ya(S)/dy
    return
end

@parallel function compute_D!(D, H, dSdx, dSdy, a, npow)
    @all(D) = a*@av(H)^(npow+2) * sqrt(@av_ya(dSdx)*@av_ya(dSdx) + @av_xa(dSdy)*@av_xa(dSdy))^(npow-1)
    return
end

@parallel function compute_flux!(qHx, qHy, D, S, dx, dy)
    @all(qHx)  = -@av_ya(D)*@d_xi(S)/dx
    @all(qHy)  = -@av_xa(D)*@d_yi(S)/dy
    return
end

@parallel function compute_dHdt!(dtau, ResH, dHdt, D, qHx, qHy, M, dtausc, cfl, epsi, damp, dx, dy)
    @all(dtau) = dtausc*min(10.0, cfl/(epsi + @av(D)))
    @all(ResH) = -(@d_xa(qHx)/dx + @d_ya(qHy)/dy) + @inn(M)
    @all(dHdt) = @all(dHdt)*damp + @all(ResH)
    return
end

@parallel function compute_H!(H, dHdt, dtau)
    @inn(H) = max(0.0, @inn(H) + @all(dtau)*@all(dHdt))
    return
end

@parallel_indices (ix,iy) function compute_Mask_S!(H, S, B, Mask)
    if (ix<=size(H,1) && iy<=size(H,2)) if (Mask[ix,iy]==0) H[ix,iy] = 0.0 end end
    if (ix<=size(S,1) && iy<=size(S,2)) S[ix,iy] = B[ix,iy] + H[ix,iy] end
    return
end

@parallel function compute_Vel!(Vx, Vy, D, H, dSdx, dSdy, epsi)
    @all(Vx) = -@all(D)/(@av(H) + epsi)*@av_ya(dSdx)
    @all(Vy) = -@all(D)/(@av(H) + epsi)*@av_xa(dSdy)
    return
end

@views function iceflow(dx, dy, Zbed, Hice, Mask, grad_b, z_ELA, b_max)
    println("Initialising ice flow model ... ")
    # physics
    s2y      = 3600*24*365.25  # seconds to years
    rho_i    = 910.0           # ice density
    g        = 9.81            # gravity acceleration
    npow     = 3.0             # Glen's power law exponent
    a0       = 1.5e-24         # Glen's law enhancement term
    # numerics
    @assert (dx>0 && dy>0) "dx and dy need to be positive"
    nx, ny   = size(Zbed,1), size(Zbed,2) # numerical grid resolution
    @assert (nx, ny) == size(Zbed) == size(Hice) == size(Mask) "Sizes don't match"
    itMax    = 1e5             # number of iteration (max)
    nout     = 200             # error check frequency
    tolnl    = 1e-6            # nonlinear tolerance
    epsi     = 1e-4            # small number
    damp     = 0.85            # convergence accelerator (this is a tuning parameter, dependent on e.g. grid resolution)
    dtausc   = 1.0/3.0         # iterative dtau scaling (this is a tuning parameter, dependent on e.g. grid resolution)
    # derived physics
    a        = 2.0*a0/(npow+2)*(rho_i*g)^npow*s2y
    # derived numerics
    cfl      = max(dx^2,dy^2)/4.1
    # array initialization
    Err      = @zeros(nx  , ny  )
    dSdx     = @zeros(nx-1, ny  )
    dSdy     = @zeros(nx  , ny-1)
    D        = @zeros(nx-1, ny-1)
    qHx      = @zeros(nx-1, ny-2)
    qHy      = @zeros(nx-2, ny-1)
    dtau     = @zeros(nx-2, ny-2)
    ResH     = @zeros(nx-2, ny-2)
    dHdt     = @zeros(nx-2, ny-2)
    Vx       = @zeros(nx-1, ny-1)
    Vy       = @zeros(nx-1, ny-1)
    M        = @zeros(nx  , ny  )
    S        = @zeros(nx  , ny  )
    grad_b   = Data.Array(grad_b)
    z_ELA    = Data.Array(z_ELA)
    Mask     = Data.Array(Mask)
    # initial conditions
    B        = Data.Array(Zbed)
    H        = Data.Array(Hice)
    S       .= B .+ H
    println(" starting iteration loop:")
    # iteration loop
    iter = 1; err = 2*tolnl
    while err>tolnl && iter<itMax
        @parallel compute_Err1!(Err, H)
        @parallel compute_M_dS!(M, dSdx, dSdy, S, z_ELA, grad_b, b_max, dx, dy)
        @parallel compute_D!(D, H, dSdx, dSdy, a, npow)
        @parallel compute_flux!(qHx, qHy, D, S, dx, dy)
        @parallel compute_dHdt!(dtau, ResH, dHdt, D, qHx, qHy, M, dtausc, cfl, epsi, damp, dx, dy)
        @parallel compute_H!(H, dHdt, dtau)
        @parallel compute_Mask_S!(H, S, B, Mask)
        # error check
        if mod(iter, nout)==0
            @parallel compute_Err2!(Err, H)
            err = norm(Err)/length(Err)
            @printf(" iter = %d, error = %1.2e \n", iter, err)
            if isnan(err)
                error("""NaNs encountered.  Try a combination of:
                             decreasing `damp` and/or `dtausc`, more smoothing steps""")
            end
        end
        iter += 1
    end
    @parallel compute_Vel!(Vx, Vy, D, H, dSdx, dSdy, epsi)
    # return as GeoArrays
    return  as_geoarray(Array(H),  Zbed, name=:thickness),
            as_geoarray(Array(S),  Zbed, name=:surface),
            as_geoarray(Array(M),  Zbed, name=:smb),
            as_geoarray(Array(Vx), Zbed, name=:vel_x, staggerd=true),
            as_geoarray(Array(Vy), Zbed, name=:vel_y, staggerd=true)
end
# ------------------------------------------------------------------------------
include(joinpath(@__DIR__, "helpers.jl"))

# load the data
print("Loading the data ... ")
Zbed, Hice, Mask, dx, dy, xc, yc = load_data(; nx=96) # nx=96,160 are included in the repo
                                                      # other numbers will trigger a 2GB download
println("done.")

# apply some smoothing
print("Applying some smoothing ... ")
for is=1:2 # two smoothing steps
    smooth!(Zbed)
    smooth!(Hice)
end
println("done.")

# calculate mass balance coefficients for given spatial grid
grad_b, z_ELA, b_max = mass_balance_constants(xc, yc)

# run the SIA flow model
H, S, M, Vx, Vy = iceflow(dx, dy, Zbed, Hice, Mask, grad_b, z_ELA, b_max)

# visualization and save
nx, ny = size(H)
if do_visu
    !ispath("../output") && mkdir("../output")

    H_v = copy(H); H_v[Mask.==0].=NaN
    Hice_v = copy(Hice); Hice_v[Mask.==0].=NaN
    S_v = copy(S); S_v[Mask.==0].=NaN
    M_v = copy(M); M_v[Mask.==0].=NaN
    V_v = sqrt.(Vx.^2 .+ Vy.^2)

    # outputs
    fontsize  = 7
    opts = (aspect_ratio=1, yaxis=font(fontsize, "Courier"), xaxis=font(fontsize, "Courier"),
            ticks=nothing, framestyle=:box, titlefontsize=fontsize, titlefont="Courier", colorbar_title="",
            xlabel="", ylabel="", xlims=(dims(H_v)[1][1],dims(H_v)[1][end]), ylims=(dims(H_v)[2][end],dims(H_v)[2][1]) )
    p1 = heatmap(S_v; c=:davos, title="Surface elev. [m]", opts...)
    p2 = heatmap(H_v; c=:davos, title="Ice thickness [m]", opts...)
    p3 = heatmap(log10.(V_v); clims=(0.1, 2.0), title="log10(vel) [m/yr]", opts...)
    p4 = heatmap(M_v; c=:devon, title="Mass Bal. rate [m/yr]", opts...)
    p = plot(p1, p2, p3, p4, size=(400,400), dpi=200) #background_color=:transparent, foreground_color=:white
    ## uncomment if you want a pop-up plot pane showing:
    # display(p)
    savefig("../output/iceflow_out1_xpu_$(nx)x$(ny).png")

    # error
    H_diff = Hice_v.-H_v
    fontsize = 7
    p1 = heatmap(Hice_v; c=:davos, title="H data [m]", opts...)
    p2 = heatmap(H_v; c=:davos, title="H model [m]", opts...)
    clim = max(abs.(extrema(H_diff[.!isnan.(H_diff)]))...)
    p3 = heatmap(H_diff; title="H (data-model) [m]",  clims=(-clim,clim), seriescolor=:balance, opts...)
    p = plot(p1, p2, p3, layout=(1, 3), size=(500,160), dpi=200) #background_color=:transparent, foreground_color=:white
    ## uncomment if you want a pop-up plot pane showing:
    # display(p)
    savefig("../output/iceflow_out2_xpu_$(nx)x$(ny).png")
end

if do_save
    save("../output/iceflow_xpu_$(nx)x$(ny).jld", "Hice", Hice,
                                                  "Mask", Mask,
                                                  "H"   , H,
                                                  "S"   , S,
                                                  "M"   , M,
                                                  "Vx"  , Vx,
                                                  "Vy"  , Vy,
                                                  "xc", xc, "yc", yc)
end

println("... done.")


text:
# Append session to options.
# Given that options may be nothing, always use the return value as the new value for options.
function _join(options::Union{Nothing, BSON}, session::Session) :: BSON
    result = options == nothing ? Mongoc.BSON() : options
    err_ref = Ref{BSONError}()
    ok = mongoc_client_session_append(session.handle, result.handle, err_ref)
    if !ok
        throw(err_ref[])
    end
    return result
end

function set_casual_consistency!(session_options::SessionOptions, casual_consistency::Bool)
    mongoc_session_opts_set_causal_consistency(session_options.handle, casual_consistency)
    nothing
end

function get_casual_consistency(session_options::SessionOptions) :: Bool
    mongoc_session_opts_get_causal_consistency(session_options.handle)
end

get_session(database_session::DatabaseSession) = database_session.session
get_session(collection_session::CollectionSession) = get_session(collection_session.database_session)

function DatabaseSession(session::Session, db_name::String)
    DatabaseSession(Database(session.client, db_name), session)
end

function CollectionSession(database_session::DatabaseSession, collection_name::String)
    CollectionSession(database_session, Collection(database_session.database, collection_name))
end

Base.getindex(session::Session, db_name::String) = DatabaseSession(session, db_name)

function Base.getindex(database_session::DatabaseSession, collection_name::String)
    CollectionSession(database_session, collection_name)
end

#
# Overload
#
function find_collections(database::DatabaseSession; options::Union{Nothing, BSON}=nothing)
    options_with_session = _join(options, get_session(database))
    find_collections(database.database; options=options_with_session)
end

function get_collection_names(database::DatabaseSession; options::Union{Nothing, BSON}=nothing)
    options_with_session = _join(options, get_session(database))
    get_collection_names(database.database, options=options_with_session)
end

function insert_one(collection::CollectionSession, document::BSON;
                    options::Union{Nothing, BSON}=nothing)
    options_with_session = _join(options, get_session(collection))
    insert_one(collection.collection, document, options=options_with_session)
end

function delete_one(collection::CollectionSession, selector::BSON;
                    options::Union{Nothing, BSON}=nothing)
    options_with_session = _join(options, get_session(collection))
    delete_one(collection.collection, selector, options=options_with_session)
end

function delete_many(collection::CollectionSession, selector::BSON;
                     options::Union{Nothing, BSON}=nothing)
    options_with_session = _join(options, get_session(collection))
    delete_many(collection.collection, selector, options=options_with_session)
end

function update_one(collection::CollectionSession, selector::BSON, update::BSON;
                    options::Union{Nothing, BSON}=nothing)
    options_with_session = _join(options, get_session(collection))
    update_one(collection.collection, selector, update, options=options_with_session)
end

function update_many(collection::CollectionSession, selector::BSON, update::BSON;
                     options::Union{Nothing, BSON}=nothing)
    options_with_session = _join(options, get_session(collection))
    update_many(collection.collection, selector, update, options=options_with_session)
end

function insert_many(collection::CollectionSession, documents::Vector{BSON};
                     bulk_options::Union{Nothing, BSON}=nothing,
                     insert_options::Union{Nothing, BSON}=nothing)
    bulk_options_with_session = _join(bulk_options, get_session(collection))
    insert_many(collection.collection, documents,
                bulk_options=bulk_options_with_session, insert_options=insert_options)
end

function find(collection::CollectionSession, bson_filter::BSON=BSON();
              options::Union{Nothing, BSON}=nothing) :: Cursor
    options_with_session = _join(options, get_session(collection))
    find(collection.collection, bson_filter, options=options_with_session)
end

function count_documents(collection::CollectionSession, bson_filter::BSON=BSON();
                         options::Union{Nothing, BSON}=nothing)
    options_with_session = _join(options, get_session(collection))
    count_documents(collection.collection, bson_filter, options=options_with_session)
end

function find_one(collection::CollectionSession, bson_filter::BSON=BSON();
                  options::Union{Nothing, BSON}=nothing) :: Union{Nothing, BSON}
    options_with_session = _join(options, get_session(collection))
    find_one(collection.collection, bson_filter, options=options_with_session)
end

function aggregate(collection::CollectionSession, bson_pipeline::BSON;
                   flags::QueryFlags=QUERY_FLAG_NONE,
                   options::Union{Nothing, BSON}=nothing) :: Cursor
    options_with_session = _join(options, get_session(collection))
    aggregate(collection.collection, bson_pipeline, flags=flags, options=options_with_session)
end

function drop(collection::CollectionSession, opts::Union{Nothing, BSON}=nothing)
    drop(collection.collection, opts)
end

#
# Transaction
#

function start_transaction!(session::Session)
    err_ref = Ref{BSONError}()
    ok = mongoc_client_session_start_transaction(session.handle, C_NULL, err_ref)
    if !ok
        throw(err_ref[])
    end
    nothing
end

function abort_transaction!(session::Session)
    err_ref = Ref{BSONError}()
    ok = mongoc_client_session_abort_transaction(session.handle, err_ref)
    if !ok
        throw(err_ref[])
    end
    nothing
end

function commit_transaction!(session::Session) :: BSON
    reply = BSON()
    err_ref = Ref{BSONError}()
    ok = mongoc_client_session_commit_transaction(session.handle, reply.handle, err_ref)
    if !ok
        throw(err_ref[])
    end
    return reply
end

in_transaction(session::Session) :: Bool = mongoc_client_session_in_transaction(session.handle)

#
# High-level API
#

"""
    transaction(f::Function, client::Client; session_options::SessionOptions=SessionOptions())

Use *do-syntax* to execute a transaction.

Transaction will be commited automatically. If an error occurs, the transaction is aborted.

The `session` parameter should be treated the same way as a `Client`:
from a `session` you get a `database`,
and a `collection` that are bound to the session.

```julia
Mongoc.transaction(client) do session
    database = session["my_database"]
    collection = database["my_collection"]
    new_item = Mongoc.BSON()
    new_item["inserted"] = true
    push!(collection, new_item)
end
```
"""
function transaction(f::Function, client::Client; session_options::SessionOptions=SessionOptions())
    local result::Union{Nothing, BSON} = nothing
    local aborted::Bool = false
    session = Session(client, options=session_options)
    start_transaction!(session)

    try
        f(session)
    catch
        abort_transaction!(session)
        aborted = true
        rethrow()
    finally
        if !aborted
            result = commit_transaction!(session)
        end
    end

    return result
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg
import BinaryBuilderBase: nbits

julia_version = v"1.5.3"

name = "ConnectFourSolver"
version = v"0.2.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/findmyway/connect4.git", "d84c8a48b6c2d2cabe0a780d51758980f8084992"),
]

# Bash recipe for building across all platforms
script = raw"""
cd connect4*

install_license LICENSE
mkdir build
cd build
cmake \
    -DCMAKE_FIND_ROOT_PATH=${prefix} \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DJulia_PREFIX=${prefix} \
    ..
VERBOSE=ON cmake --build . --config Release --target install -- -j${nproc}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
include("../../L/libjulia/common.jl")
platforms = libjulia_platforms(julia_version)
platforms = expand_cxxstring_abis(platforms)
filter!(p -> nbits(p) != 32, platforms) # code use __int128 which is not available on 32bit targets

# The products that we will ensure are always built
products = [
    LibraryProduct("libconnect4jl", :libconnect4jl),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency(PackageSpec(name="libjulia_jll", version=julia_version)),
    Dependency("libcxxwrap_julia_jll")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
    preferred_gcc_version=v"8",
    julia_compat = "$(julia_version.major).$(julia_version.minor)")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
abstract type PSS <: GeneratorComponent end

@def pss_ports begin
      state_input = [:ω]
      inner_input = [τe_var, VR_gen_var]
  end

"""
Parameters of a PSS that returns a fixed voltage to add to
the reference for the AVR

# Conmutable structor
```julia
PSSFixed(Vs)
```

# Arguments
* `Vs`::Float64 : Fixed voltage stabilization signal

"""
mutable struct PSSFixed <: PSS
      Vs::Float64
      n_states::Int64
      states::Vector{Symbol}
      ports::Ports

      function PSSFixed(Vs::Float64)
            n_states = 0
            states = Vector{Symbol}()
            @pss_ports
            new(Vs, n_states, states, Ports(state_input, inner_input))
      end
end

"""
Parameters of a PSS that returns a proportional droop voltage to add to
the reference for the AVR

# Conmutable structor
```julia
PSSSimple(K_ω, K_p)
```

# Arguments
* `K_ω`::Float64 : Proportional gain for frequency
* `K_p`::Float64 : Proportional gain for active power

"""
mutable struct PSSSimple <: PSS
      K_ω::Float64
      K_p::Float64
      n_states::Int64
      states::Vector{Symbol}
      ports::Ports

      function PSSSimple(K_ω::Float64, K_p::Float64)
            n_states = 0
            states = Vector{Symbol}()
            @pss_ports
            new(K_ω, K_p, n_states, states, Ports(state_input, inner_input))
      end
end